Surface refl

src/core/geodesy/geodetic.cpp

@@ -283,12 +283,9 @@ Vector3 approx_geometrical_tangent_point(Vector3 ecef,
 
     const Numeric a2 = refellipsoid[0] * refellipsoid[0];
     const Numeric b2 = refellipsoid[1] * refellipsoid[1];
-    Vector3 yunit, zunit;
 
-    cross3(zunit, decef, ecef);
-    zunit /= norm2(zunit);
-    cross3(yunit, zunit, decef);
-    yunit /= norm2(yunit);
+    const Vector3 zunit = normalized(cross(decef, ecef));
+    const Vector3 yunit = normalized(cross(zunit, decef));
 
     const Numeric yr = dot(ecef, yunit);
     const Numeric xr = dot(ecef, decef);

src/core/matpack/matpack_mdspan_cdata_t.h

@@ -275,6 +275,17 @@ struct [[nodiscard]] cdata_t {
   }
 };
 
+template <any_cdata out, exact_md<value_type<out>, rank<out>()> in>
+constexpr out to(const in& x) {
+  out o{};
+
+  assert(x.shape() == o.shape());
+
+  stdr::copy(x | by_elem, o.data.begin());
+
+  return o;
+}
+
 template <any_cdata T>
 constexpr T operator+(const T& x, const T& y) {
   T z  = x;

src/core/matpack/matpack_mdspan_common_select.h

@@ -15,20 +15,23 @@
 #include "matpack_mdspan_common_types.h"
 
 namespace matpack {
-namespace stdx = std::experimental;
-
-using Joker                 = stdx::full_extent_t;
-constexpr inline auto joker = stdx::full_extent;
+namespace stdx               = std::experimental;
+using Joker                  = stdx::full_extent_t;
+constexpr inline Joker joker = stdx::full_extent;
 
 struct [[nodiscard]] StridedRange {
-  Index offset;
-  Index nelem;
-  Index stride;
-
-  template <integral Offset = Index,
-            integral Nelem  = Index,
-            integral Stride = Index>
-  constexpr StridedRange(Offset i0 = 0, Nelem n = 0, Stride d = 1)
+  Index offset{0};
+  Index nelem{0};
+  Index stride{1};
+
+  constexpr StridedRange()                                   = default;
+  constexpr StridedRange(const StridedRange&)                = default;
+  constexpr StridedRange(StridedRange&&) noexcept            = default;
+  constexpr StridedRange& operator=(const StridedRange&)     = default;
+  constexpr StridedRange& operator=(StridedRange&&) noexcept = default;
+
+  template <integral Offset, integral Nelem, integral Stride>
+  constexpr StridedRange(Offset i0, Nelem n, Stride d)
       : offset(static_cast<Index>(i0)),
         nelem(static_cast<Index>(n)),
         stride(static_cast<Index>(d)) {}
@@ -41,11 +44,17 @@ struct [[nodiscard]] StridedRange {
 };
 
 struct [[nodiscard]] Range {
-  Index offset;
-  Index nelem;
+  Index offset{0};
+  Index nelem{0};
+
+  constexpr Range()                            = default;
+  constexpr Range(const Range&)                = default;
+  constexpr Range(Range&&) noexcept            = default;
+  constexpr Range& operator=(const Range&)     = default;
+  constexpr Range& operator=(Range&&) noexcept = default;
 
-  template <integral Offset = Index, integral Nelem = Index>
-  constexpr Range(Offset i0 = 0, Nelem n = 0)
+  template <integral Offset, integral Nelem>
+  constexpr Range(Offset i0, Nelem n)
       : offset(static_cast<Index>(i0)), nelem(static_cast<Index>(n)) {}
 
   [[nodiscard]] constexpr stdx::
@@ -55,6 +64,25 @@ struct [[nodiscard]] Range {
   }
 };
 
+struct PartialRange {
+  Index n;
+
+  template <integral Offset>
+  constexpr PartialRange(Offset i) : n(static_cast<Index>(i)) {}
+
+  friend constexpr Range operator|(const PartialRange& pr, Index n) {
+    return Range{pr.n, n};
+  }
+
+  friend constexpr Range operator|(Index n, const PartialRange& pr) {
+    return Range{n, pr.n};
+  }
+};
+
+namespace literals {
+constexpr PartialRange operator""_i0(unsigned long long int i) { return i; }
+}  // namespace literals
+
 template <typename T>
 concept exhaustive_access_operator =
     integral<T> or std::is_same_v<std::remove_cvref_t<T>, Joker> or
@@ -331,3 +359,5 @@ struct std::formatter<Joker> {
     return tags.format(ctx, "joker"sv);
   }
 };
+
+using namespace matpack::literals;

src/core/matpack/matpack_mdspan_helpers_reduce.h

@@ -231,7 +231,17 @@ constexpr T dot(const Self& self, const Other& other, T init = {}) {
  */
 template <any_md Self>
 constexpr auto hypot(const Self& self) {
-  return std::sqrt(dot(self, self));
+  if constexpr (any_cdata<Self> and ranked<Self, 1>) {
+    if constexpr (Self::ndata == 3) {
+      return std::hypot(self[0], self[1], self[2]);
+    } else if constexpr (Self::ndata == 2) {
+      return std::hypot(self[0], self[1]);
+    } else {
+      return std::sqrt(dot(self, self));
+    }
+  } else {
+    return std::sqrt(dot(self, self));
+  }
 }
 
 template <any_md Self, any_md Other>

src/core/matpack/matpack_mdspan_helpers_vector.h

@@ -6,22 +6,13 @@ namespace matpack {
 Vector uniform_grid(Numeric x0, Index N, Numeric dx);
 ComplexVector uniform_grid(Complex x0, Index N, Complex dx);
 
-template <ranked_md<1> VEC1, ranked_md<1> VEC2>
-constexpr Vector3 cross(const VEC1 &a, const VEC2 &b) noexcept {
+constexpr Vector3 cross(const Vector3 &a, const Vector3 &b) noexcept {
   assert(a.size() == 3 and b.size() == 3);
   return {a[1] * b[2] - a[2] * b[1],
           a[2] * b[0] - a[0] * b[2],
           a[0] * b[1] - a[1] * b[0]};
 }
 
-template <mut_ranked_md<1> VEC1, ranked_md<1> VEC2, ranked_md<1> VEC3>
-constexpr void cross3(VEC1 &&x, const VEC2 &a, const VEC3 &b) noexcept {
-  assert(x.size() == 3 and a.size() == 3 and b.size() == 3);
-  x[0] = a[1] * b[2] - a[2] * b[1];
-  x[1] = a[2] * b[0] - a[0] * b[2];
-  x[2] = a[0] * b[1] - a[1] * b[0];
-}
-
 void gaussian_grid(Vector &y,
                    const Vector &x,
                    const Numeric &x0,

src/core/montecarlo/mc_antenna.cc

@@ -175,16 +175,14 @@ std::pair<Vector2, Matrix33> MCAntenna::draw_los(
         R_los[joker, 1]    = R_ant2enu[1, joker];
         sampled_rte_los[1] = bore_sight_los[1];
       } else {
-        const Vector uhat{0.0, 0.0, 1.0};
-        Numeric magh;
+        const Vector3 uhat{0.0, 0.0, 1.0};
         sampled_rte_los[1] = atan2d(R_los[0, 2], R_los[1, 2]);
-        cross3(R_los[joker, 1], R_los[joker, 2], uhat);
-        magh             = hypot(R_los[joker, 1]);
-        R_los[joker, 1] /= magh;
+        R_los[joker, 1] = normalized(cross(to<Vector3>(R_los[joker, 2]), uhat));
       }
 
       // Vertical polarization basis
-      cross3(R_los[joker, 0], R_los[joker, 1], R_los[joker, 2]);
+      R_los[joker, 0] =
+          cross(to<Vector3>(R_los[joker, 1]), to<Vector3>(R_los[joker, 2]));
 
       break;
 

src/core/rtepack/rtepack_surface.cc

@@ -29,6 +29,160 @@ muelmat fresnel_reflectance(Complex Rv, Complex Rh) {
   return out;
 }
 
+namespace {
+muelmat stokes_rotation(Numeric cos2psi, Numeric sin2psi) {
+  muelmat L{};
+  L[0, 0] = 1.0;
+  L[1, 1] = cos2psi;
+  L[1, 2] = sin2psi;
+  L[2, 1] = -sin2psi;
+  L[2, 2] = cos2psi;
+  L[3, 3] = 1.0;
+  return L;
+}
+
+muelmat stokes_rotation_refl(Numeric cos2psi, Numeric sin2psi) {
+  muelmat L{};
+  L[0, 0] = 1.0;
+  L[1, 1] = cos2psi;
+  L[1, 2] = -sin2psi;
+  L[2, 1] = sin2psi;
+  L[2, 2] = -cos2psi;
+  L[3, 3] = -1.0;
+  return L;
+}
+
+/** Compute polarization basis vectors (v, h) for a propagation direction k.
+ *
+ * Uses the local vertical z = (0, 0, 1) as reference.
+ * h = k x z / |k x z|,  v = h x k
+ */
+static std::pair<Vector3, Vector3> pol_basis(const Vector3& k) {
+  const Vector3 z{0.0, 0.0, 1.0};
+  Vector3 h = cross(k, z);
+  const Numeric norm_h = hypot(h);
+
+  if (norm_h < 1e-12) {
+    h = Vector3{1.0, 0.0, 0.0};
+  } else {
+    h[0] /= norm_h;
+    h[1] /= norm_h;
+    h[2] /= norm_h;
+  }
+
+  const Vector3 v = cross(h, k);
+  return {v, h};
+}
+}  // namespace
+
+/** Fresnel reflection Mueller matrix for specular reflection off a
+ *  surface with arbitrary tilt.
+ *
+ *  @param Rv        Complex Fresnel coefficient for vertical polarization
+ *  @param Rh        Complex Fresnel coefficient for horizontal polarization
+ *  @param k_inc     Unit propagation vector of incident beam (toward surface)
+ *  @param n_surface Outward unit surface normal
+ *  @return          4x4 Mueller matrix in the incident beam's polarization frame
+ *
+ *  For a flat horizontal surface (n_surface = (0,0,1)) this reduces to
+ *  fresnel_reflectance(Rv, Rh) with U and V sign flips.
+ */
+muelmat fresnel_reflectance_specular(Complex Rv,
+                                     Complex Rh,
+                                     const Vector3& k_inc,
+                                     const Vector3& n_surface) {
+  const muelmat M_fresnel = fresnel_reflectance(Rv, Rh);
+
+  Vector3 m = cross(k_inc, n_surface);
+  const Numeric norm_m = hypot(m);
+
+  if (norm_m < 1e-12) {
+    muelmat F = muelmat::id();
+    F[2, 2] = -1.0;
+    F[3, 3] = -1.0;
+    return F * M_fresnel;
+  }
+
+  m[0] /= norm_m;
+  m[1] /= norm_m;
+  m[2] /= norm_m;
+
+  const auto [v_i, h_i] = pol_basis(k_inc);
+
+  const Numeric cos_psi1 = dot(h_i, m);
+  const Numeric sin_psi1 = dot(v_i, m);
+  const Numeric cos2psi1 = 2.0 * cos_psi1 * cos_psi1 - 1.0;
+  const Numeric sin2psi1 = 2.0 * sin_psi1 * cos_psi1;
+
+  const muelmat L1 = stokes_rotation(cos2psi1, sin2psi1);
+  const muelmat L2 = stokes_rotation_refl(cos2psi1, -sin2psi1);
+
+  return L2 * M_fresnel * L1;
+}
+
+/** Fresnel reflection Mueller matrix for non-specular (e.g. Lambertian,
+ *  BRDF, or tilted-facet) reflection where incident and outgoing directions
+ *  are independent.
+ *
+ *  @param Rv        Complex Fresnel coefficient for vertical polarization
+ *  @param Rh        Complex Fresnel coefficient for horizontal polarization
+ *  @param k_inc     Unit propagation vector of incident beam (toward surface)
+ *  @param k_out     Unit propagation vector of outgoing beam (toward observer)
+ *  @param n_surface Outward unit surface normal
+ *  @return          4x4 Mueller matrix mapping incident Stokes to outgoing Stokes
+ */
+muelmat fresnel_reflectance_nonspecular(Complex Rv,
+                                        Complex Rh,
+                                        const Vector3& k_inc,
+                                        const Vector3& k_out,
+                                        const Vector3& n_surface) {
+  // Fresnel matrix in the plane-of-incidence frame
+  const muelmat M_fresnel = fresnel_reflectance(Rv, Rh);
+
+  // Plane-of-incidence normal: m = k_inc x n / |k_inc x n|
+  Vector3 m = cross(k_inc, n_surface);
+  const Numeric norm_m = hypot(m);
+
+  if (norm_m < 1e-12) {
+    // Normal incidence: plane of incidence undefined
+    muelmat F = muelmat::id();
+    F[2, 2] = -1.0;
+    F[3, 3] = -1.0;
+    return F * M_fresnel;
+  }
+
+  m[0] /= norm_m;
+  m[1] /= norm_m;
+  m[2] /= norm_m;
+
+  // Incident beam polarization basis and rotation angle psi_1
+  const auto [v_i, h_i] = pol_basis(k_inc);
+  const Numeric cos_psi1 = dot(h_i, m);
+  const Numeric sin_psi1 = dot(v_i, m);
+  const Numeric cos2psi1 = 2.0 * cos_psi1 * cos_psi1 - 1.0;
+  const Numeric sin2psi1 = 2.0 * sin_psi1 * cos_psi1;
+
+  const muelmat L1 = stokes_rotation(cos2psi1, sin2psi1);
+
+  // Outgoing beam polarization basis and rotation angle psi_2
+  const auto [v_r, h_r] = pol_basis(k_out);
+  const Numeric cos_psi2 = dot(m, h_r);
+  const Numeric sin_psi2 = dot(m, v_r);
+  const Numeric cos2psi2 = 2.0 * cos_psi2 * cos_psi2 - 1.0;
+  const Numeric sin2psi2 = 2.0 * sin_psi2 * cos_psi2;
+
+  // Reflected-frame rotation: includes handedness change
+  const muelmat L2 = stokes_rotation_refl(cos2psi2, sin2psi2);
+
+  return L2 * M_fresnel * L1;
+}
+
+Vector3 specular_reflected_direction(const Vector3& k_inc,
+                                     const Vector3& n_surface) {
+  const Vector3 out = k_inc - 2.0 * dot(k_inc, n_surface) * n_surface;
+  return normalized(out);
+}
+
 stokvec flat_scalar_reflection(stokvec I, const Numeric R, const stokvec B) {
   I      = I * R;
   I.V() *= -1.0;  //! NOTE: The elementwise multiplication is [R, R, R, -R]
@@ -44,14 +198,14 @@ stokvec dflat_scalar_reflection_dr(stokvec I, const stokvec B) {
 
 stokvec reflection(stokvec I, const muelmat R, const stokvec B) {
   I      = R * I;
-  I.V() *= -1.0;  //! FIXME: Is this correct?
+  I.V() *= -1.0;
   I     += (1.0 - R) * B;
   return I;
 }
 
 stokvec dreflection(stokvec I, const muelmat dR, const stokvec B) {
   I      = dR * I;
-  I.V() *= -1.0;  //! FIXME: Is this correct?
+  I.V() *= -1.0;
   I     -= dR * B;
   return I;
 }