diff --git a/src/core/geodesy/geodetic.cpp b/src/core/geodesy/geodetic.cpp index daa23de985..f73c533948 100644 --- a/src/core/geodesy/geodetic.cpp +++ b/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); diff --git a/src/core/matpack/matpack_mdspan_cdata_t.h b/src/core/matpack/matpack_mdspan_cdata_t.h index b940942f7b..22445dae03 100644 --- a/src/core/matpack/matpack_mdspan_cdata_t.h +++ b/src/core/matpack/matpack_mdspan_cdata_t.h @@ -275,6 +275,17 @@ struct [[nodiscard]] cdata_t { } }; +template , rank()> 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 constexpr T operator+(const T& x, const T& y) { T z = x; diff --git a/src/core/matpack/matpack_mdspan_common_select.h b/src/core/matpack/matpack_mdspan_common_select.h index 8b8615732a..5ad83b033c 100644 --- a/src/core/matpack/matpack_mdspan_common_select.h +++ b/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 - 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 + constexpr StridedRange(Offset i0, Nelem n, Stride d) : offset(static_cast(i0)), nelem(static_cast(n)), stride(static_cast(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 - constexpr Range(Offset i0 = 0, Nelem n = 0) + template + constexpr Range(Offset i0, Nelem n) : offset(static_cast(i0)), nelem(static_cast(n)) {} [[nodiscard]] constexpr stdx:: @@ -55,6 +64,25 @@ struct [[nodiscard]] Range { } }; +struct PartialRange { + Index n; + + template + constexpr PartialRange(Offset i) : n(static_cast(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 concept exhaustive_access_operator = integral or std::is_same_v, Joker> or @@ -331,3 +359,5 @@ struct std::formatter { return tags.format(ctx, "joker"sv); } }; + +using namespace matpack::literals; diff --git a/src/core/matpack/matpack_mdspan_helpers_reduce.h b/src/core/matpack/matpack_mdspan_helpers_reduce.h index 39d5affeca..2791332bd0 100644 --- a/src/core/matpack/matpack_mdspan_helpers_reduce.h +++ b/src/core/matpack/matpack_mdspan_helpers_reduce.h @@ -231,7 +231,17 @@ constexpr T dot(const Self& self, const Other& other, T init = {}) { */ template constexpr auto hypot(const Self& self) { - return std::sqrt(dot(self, self)); + if constexpr (any_cdata and ranked) { + 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 diff --git a/src/core/matpack/matpack_mdspan_helpers_vector.h b/src/core/matpack/matpack_mdspan_helpers_vector.h index e360d9e149..11bb65b7d3 100644 --- a/src/core/matpack/matpack_mdspan_helpers_vector.h +++ b/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 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 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, diff --git a/src/core/montecarlo/mc_antenna.cc b/src/core/montecarlo/mc_antenna.cc index baa0ad35f0..989635f697 100644 --- a/src/core/montecarlo/mc_antenna.cc +++ b/src/core/montecarlo/mc_antenna.cc @@ -175,16 +175,14 @@ std::pair 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(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(R_los[joker, 1]), to(R_los[joker, 2])); break; diff --git a/src/core/rtepack/rtepack_surface.cc b/src/core/rtepack/rtepack_surface.cc index 76b362a2a1..a0c1b85e45 100644 --- a/src/core/rtepack/rtepack_surface.cc +++ b/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 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; } diff --git a/src/core/rtepack/rtepack_surface.h b/src/core/rtepack/rtepack_surface.h index 95847b3d13..c08efc67be 100644 --- a/src/core/rtepack/rtepack_surface.h +++ b/src/core/rtepack/rtepack_surface.h @@ -29,11 +29,49 @@ stokvec dflat_scalar_reflection_dr(stokvec I, const stokvec B); /** Fresnel style reflectance matrix * - * @param[out] Rv Reflection AMPLITUDE coefficient for vertical polarisation. - * @param[out] Rh Reflection AMPLITUDE coefficient for vertical polarisation. + * @param Rv Reflection AMPLITUDE coefficient for vertical polarisation. + * @param Rh Reflection AMPLITUDE coefficient for horizontal polarisation. */ muelmat fresnel_reflectance(Complex Rv, Complex Rh); +/** Fresnel reflection Mueller matrix for specular reflection off a + * surface with arbitrary tilt. + * + * UNTESTED + * + * @param Rv Reflection AMPLITUDE coefficient for vertical polarisation. + * @param Rh Reflection AMPLITUDE coefficient for horizontal polarisation. + * @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); + +/** Fresnel reflection Mueller matrix for non-specular (e.g. Lambertian, + * BRDF, or tilted-facet) reflection where incident and outgoing directions + * are independent. + * + * UNTESTED + * + * @param Rv Reflection AMPLITUDE coefficient for vertical polarisation. + * @param Rh Reflection AMPLITUDE coefficient for horizontal polarisation. + * @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); + /** Reflect and emit a Stokes vector * * The reflection is scalar, i.e. the same for all Stokes components. @@ -49,4 +87,12 @@ muelmat fresnel_reflectance(Complex Rv, Complex Rh); */ stokvec reflection(stokvec I, const muelmat R, const stokvec B); stokvec dreflection(stokvec I, const muelmat dR, const stokvec B); + +/** Return the outgoing propagation direction for specular reflection. + * + * Both `k_inc` and `n_surface` are assumed to be unit vectors with + * `k_inc` pointing toward the surface and `n_surface` pointing outwards. + */ +Vector3 specular_reflected_direction(const Vector3& k_inc, + const Vector3& n_surface); } // namespace rtepack diff --git a/src/python_interface/py_rtepack.cpp b/src/python_interface/py_rtepack.cpp index 70425ed211..9d69d9d04c 100644 --- a/src/python_interface/py_rtepack.cpp +++ b/src/python_interface/py_rtepack.cpp @@ -15,7 +15,6 @@ #include "hpy_arts.h" #include "hpy_numpy.h" #include "hpy_vector.h" -#include "rtepack_transmission.h" namespace Python { template @@ -608,7 +607,37 @@ Returns ------- Propmat The logarithm of the Mueller matrix as a Propagation matrix. -)"); +)") + .def( + "specular_reflected_direction", + &rtepack::specular_reflected_direction, + "k_inc"_a, + "n_surface"_a, + R"(Return the direction of the specularly reflected propagation vector.)") + .def("fresnel_reflectance", + &rtepack::fresnel_reflectance, + "Rv"_a, + "Rh"_a, + R"(Return the Fresnel Mueller matrix for the complex amplitude +coefficients `Rv` and `Rh`.)") + .def("fresnel_reflectance_specular", + &rtepack::fresnel_reflectance_specular, + "Rv"_a, + "Rh"_a, + "k_inc"_a, + "n_surface"_a, + R"(Return the Fresnel Mueller matrix for specular reflection. `k_inc` +is the incident propagation vector toward the surface and `n_surface` is the +outward surface normal.)") + .def("fresnel_reflectance_nonspecular", + &rtepack::fresnel_reflectance_nonspecular, + "Rv"_a, + "Rh"_a, + "k_inc"_a, + "k_out"_a, + "n_surface"_a, + R"(Return the Fresnel Mueller matrix for non-specular reflection with +independent incident and outgoing directions.)"); } catch (std::exception &e) { throw std::runtime_error( std::format("DEV ERROR:\nCannot initialize rtepack\n{}", e.what())); diff --git a/tests/python/math/test_surface_helpers.py b/tests/python/math/test_surface_helpers.py new file mode 100644 index 0000000000..62b8af648f --- /dev/null +++ b/tests/python/math/test_surface_helpers.py @@ -0,0 +1,107 @@ +import numpy as np + +from pyarts3 import arts + + +def _expected_fresnel_matrix(Rv: complex, Rh: complex) -> np.ndarray: + rv = np.abs(Rv) ** 2 + rh = np.abs(Rh) ** 2 + rmean = 0.5 * (rv + rh) + rdiff = 0.5 * (rv - rh) + a = Rh * np.conjugate(Rv) + b = Rv * np.conjugate(Rh) + c = 0.5 * np.real(a + b) + d = 0.5 * np.imag(a - b) + + expected = np.zeros((4, 4), dtype=np.float64) + expected[0, 0] = rmean + expected[1, 0] = rdiff + expected[0, 1] = rdiff + expected[1, 1] = rmean + expected[2, 2] = c + expected[2, 3] = d + expected[3, 2] = -d + expected[3, 3] = c + return expected + + +def _to_array(value) -> np.ndarray: + return np.asarray(value) + + +results = {} + +n_surface = [0.0, 0.0, 1.0] + +Rv = 0.3 + 0.2j +Rh = -0.1 + 0.7j + +actual = _to_array(arts.rtepack.fresnel_reflectance(Rv, Rh)) +expected = _expected_fresnel_matrix(Rv, Rh) + +print("Fresnel reflectance matrix (expected):\n", expected) +print("Fresnel reflectance matrix (actual):\n", actual) + +assert np.allclose(actual, expected, rtol=1e-12, atol=1e-15) +results["fresnel_matrix"] = (actual, expected) + + +straight_in = [0.0, 0.0, -1.0] +straight_out = _to_array( + arts.rtepack.specular_reflected_direction(straight_in, n_surface) +) +print("Specular reflection (normal incidence) computed:", straight_out) +print("Specular reflection (normal incidence) expected: [0.0, 0.0, 1.0]") +assert np.allclose(straight_out, [0.0, 0.0, 1.0]) +results["specular_normal"] = straight_out + +angle = np.pi / 6 +oblique_in = [np.sin(angle), 0.0, -np.cos(angle)] +oblique_out = _to_array( + arts.rtepack.specular_reflected_direction(oblique_in, n_surface) +) +expected_oblique = [np.sin(angle), 0.0, np.cos(angle)] +print("Specular reflection (oblique incidence) computed:", oblique_out) +print("Specular reflection (oblique incidence) expected:", expected_oblique) +assert np.allclose(oblique_out, expected_oblique, rtol=1e-12) +results["specular_oblique"] = oblique_out + +Rv = 0.4 + 0.3j +Rh = 0.1 + 0.5j +k_inc = [0.0, 0.0, -1.0] + +specular = _to_array( + arts.rtepack.fresnel_reflectance_specular(Rv, Rh, k_inc, n_surface) +) +fresnel = _expected_fresnel_matrix(Rv, Rh) + +print("Specular Fresnel rows (computed):\n", specular) +print("Specular Fresnel rows (expected with UV flip):\n", fresnel) +assert np.allclose(specular[0], fresnel[0], rtol=1e-12, atol=1e-15) +assert np.allclose(specular[1], fresnel[1], rtol=1e-12, atol=1e-15) +assert np.allclose(specular[2], -fresnel[2], rtol=1e-12, atol=1e-15) +assert np.allclose(specular[3], -fresnel[3], rtol=1e-12, atol=1e-15) +results["specular_uv_flip"] = specular + +Rv = 0.2 + 0.4j +Rh = -0.3 + 0.1j +angle = np.pi / 4 +k_inc = [np.sin(angle), 0.0, -np.cos(angle)] + + +k_out = arts.rtepack.specular_reflected_direction(k_inc, n_surface) + +specular = _to_array( + arts.rtepack.fresnel_reflectance_specular(Rv, Rh, k_inc, n_surface) +) +nonspecular = _to_array( + arts.rtepack.fresnel_reflectance_nonspecular( + Rv, Rh, k_inc, k_out, n_surface + ) +) + +print("Specular Fresnel Mueller (computed):\n", specular) +print("Nonspecular Fresnel Mueller (computed):\n", nonspecular) +assert np.allclose(specular, nonspecular, rtol=1e-12, atol=1e-15) +results["nonspecular_agreement"] = specular +