From 7ec49b1cc53c4dcfbf69587b469651b8ff880ff4 Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Wed, 17 Dec 2025 19:06:27 +0900 Subject: [PATCH 1/5] Add inverse to matrix exponent for Mueller matrix --- src/core/rtepack/rtepack_mueller_matrix.h | 28 ++- src/core/rtepack/rtepack_transmission.cc | 282 ++++++++++++++++++---- src/core/rtepack/rtepack_transmission.h | 2 + src/python_interface/py_rtepack.cpp | 43 ++++ 4 files changed, 298 insertions(+), 57 deletions(-) diff --git a/src/core/rtepack/rtepack_mueller_matrix.h b/src/core/rtepack/rtepack_mueller_matrix.h index 62b15450ee..86f0280d81 100644 --- a/src/core/rtepack/rtepack_mueller_matrix.h +++ b/src/core/rtepack/rtepack_mueller_matrix.h @@ -167,16 +167,30 @@ constexpr muelmat avg(muelmat a, const muelmat &b) { return a; } +constexpr Numeric det(const muelmat &A) { + const auto [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p] = A; + + return a * f * k * p - a * f * l * o - a * g * j * p + a * g * l * n + + a * h * j * o - a * h * k * n - b * e * k * p + b * e * l * o + + b * g * i * p - b * g * l * m - b * h * i * o + b * h * k * m + + c * e * j * p - c * e * l * n - c * f * i * p + c * f * l * m + + c * h * i * n - c * h * j * m - d * e * j * o + d * e * k * n + + d * f * i * o - d * f * k * m - d * g * i * n + d * g * j * m; +} + +constexpr Numeric tr(const muelmat &A) { + return A[0, 0] + A[1, 1] + A[2, 2] + A[3, 3]; +} + +constexpr Numeric midtr(const muelmat &A) { + return std::midpoint(std::midpoint(A[0, 0], A[1, 1]), + std::midpoint(A[2, 2], A[3, 3])); +} + constexpr muelmat inv(const muelmat &A) { const auto [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p] = A; - const Numeric div = - 1.0 / (a * f * k * p - a * f * l * o - a * g * j * p + a * g * l * n + - a * h * j * o - a * h * k * n - b * e * k * p + b * e * l * o + - b * g * i * p - b * g * l * m - b * h * i * o + b * h * k * m + - c * e * j * p - c * e * l * n - c * f * i * p + c * f * l * m + - c * h * i * n - c * h * j * m - d * e * j * o + d * e * k * n + - d * f * i * o - d * f * k * m - d * g * i * n + d * g * j * m); + const Numeric div = 1.0 / det(A); return div * muelmat{(f * k * p - f * l * o - g * j * p + g * l * n + h * j * o - h * k * n), diff --git a/src/core/rtepack/rtepack_transmission.cc b/src/core/rtepack/rtepack_transmission.cc index 87fd57eb18..6ac8c3acf0 100644 --- a/src/core/rtepack/rtepack_transmission.cc +++ b/src/core/rtepack/rtepack_transmission.cc @@ -4,6 +4,7 @@ #include #include +#include #include "rtepack_mueller_matrix.h" #include "rtepack_propagation_matrix.h" @@ -31,7 +32,7 @@ tran::tran(const propmat &k1, const propmat &k2, const Numeric r) v2 = v * v; w2 = w * w; - /* Solve: + /** Solve: 0 = L^4 + B L^2 + C B = U^2+V^2+W^2-B^2-C^2-D^2 C = - (DU - CV + BW)^2 @@ -40,7 +41,7 @@ tran::tran(const propmat &k1, const propmat &k2, const Numeric r) C = -Math::pow2(d * u - c * v + b * w); S = std::sqrt(B * B - 4 * C); - /* + /** We define: x2 and y2 are the squares of x and y x and y are real and positive @@ -95,58 +96,48 @@ tran::tran(const propmat &k1, const propmat &k2, const Numeric r) C0 = either_zero ? 1.0 : (cy * x2 + cx * y2) * inv_x2y2; C1 = either_zero ? 1.0 : (sy * x2 * iy + sx * y2 * ix) * inv_x2y2; C2 = both_zero ? 0.5 : (cx - cy) * inv_x2y2; - C3 = both_zero ? 1.0 / 6.0 - : (x_zero ? 1.0 - sy * iy - : y_zero ? sx * ix - 1.0 - : sx * ix - sy * iy) * - inv_x2y2; + C3 = both_zero + ? 1.0 / 6.0 + : ((x_zero ? 1.0 : sx * ix) - (y_zero ? 1.0 : sy * iy)) * inv_x2y2; } muelmat tran::operator()() const noexcept { + if (not polarized) return exp_a; + // Do the calculation exp(a) * (C0 * I + C1 * K + C2 * K^2 + C3 * K^3) - return not polarized - ? muelmat{exp_a} - : muelmat{ - exp_a * (C0 + C2 * (b2 + c2 + d2)), - -exp_a * (-C1 * b + C2 * (c * u + d * v) + - C3 * (u * (b * u - d * w) - b * (b2 + c2 + d2) + - v * (b * v + c * w))), - exp_a * (C1 * c + C2 * (b * u - d * w) + - C3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) - - w * (b * v + c * w))), - exp_a * (C1 * d + C2 * (b * v + c * w) + - C3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) + - w * (b * u - d * w))), - exp_a * (C1 * b + C2 * (c * u + d * v) + - C3 * (c * (b * c - v * w) - b * (-b2 + u2 + v2) + - d * (b * d + u * w))), - exp_a * (C0 + C2 * (b2 - u2 - v2)), - exp_a * (C1 * u + C2 * (b * c - v * w) + - C3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) - - w * (b * d + u * w))), - exp_a * (C1 * v + C2 * (b * d + u * w) + - C3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) + - w * (b * c - v * w))), - exp_a * (C1 * c + C2 * (d * w - b * u) + - C3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) + - d * (c * d - u * v))), - -exp_a * (C1 * u + C2 * (v * w - b * c) + - C3 * (b * (b * u - d * w) - u * (-c2 + u2 + w2) + - v * (c * d - u * v))), - exp_a * (C0 + C2 * (c2 - u2 - w2)), - exp_a * (C1 * w + C2 * (c * d - u * v) + - C3 * (v * (b * c - v * w) - d * (b * u - d * w) - - w * (-c2 + u2 + w2))), - exp_a * (C1 * d - C2 * (b * v + c * w) + - C3 * (b * (b * d + u * w) + c * (c * d - u * v) - - d * (-d2 + v2 + w2))), - -exp_a * (C1 * v - C2 * (b * d + u * w) + - C3 * (b * (b * v + c * w) + u * (c * d - u * v) - - v * (-d2 + v2 + w2))), - -exp_a * (C1 * w + C2 * (u * v - c * d) + - C3 * (c * (b * v + c * w) - u * (b * d + u * w) - - w * (-d2 + v2 + w2))), - exp_a * (C0 + C2 * (d2 - v2 - w2))}; + + const Numeric C2b = C2 * (c * u + d * v); + const Numeric C2c = C2 * (b * u - d * w); + const Numeric C2d = C2 * (b * v + c * w); + const Numeric C2u = C2 * (b * c - v * w); + const Numeric C2v = C2 * (b * d + u * w); + const Numeric C2w = C2 * (c * d - u * v); + + const Numeric C3b = + C3 * (u * (b * u - d * w) - b * (b2 + c2 + d2) + v * (b * v + c * w)); + const Numeric C3c = + C3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) - w * (b * v + c * w)); + const Numeric C3d = + C3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) + w * (b * u - d * w)); + const Numeric C3u = + C3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) - w * (b * d + u * w)); + const Numeric C3v = + C3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) + w * (b * c - v * w)); + const Numeric C3w = + C3 * (v * (b * c - v * w) - d * (b * u - d * w) - w * (-c2 + u2 + w2)); + + const Numeric M00 = C0 + C2 * (b2 + c2 + d2); + const Numeric M11 = C0 + C2 * (b2 - u2 - v2); + const Numeric M22 = C0 + C2 * (c2 - u2 - w2); + const Numeric M33 = C0 + C2 * (d2 - v2 - w2); + + // clang-format off + return exp_a * muelmat{ + M00, C1 * b - C2b - C3b, C1 * c + C2c + C3c, C1 * d + C2d + C3d, + C1 * b + C2b - C3b, M11, C1 * u + C2u + C3u, C1 * v + C2v + C3v, + C1 * c - C2c + C3c, - C1 * u + C2u - C3u, M22, C1 * w + C2w + C3w, + C1 * d - C2d + C3d, - C1 * v + C2v - C3v, - C1 * w + C2w - C3w, M33}; + // clang-format on } muelmat tran::deriv(const muelmat &t, @@ -499,4 +490,195 @@ void two_level_exp(std::vector &T, dr[1, i - 1]); } } + +propmat logK(const muelmat &m) { + if (not m.is_polarized()) return std::log(midtr(m)); + + /** + The code tries to retrieve K from exp(K) = M input as muelmat m, + where K is a propmat. + + K is unknown but, as a reminder, it looks like: + K = [ a b c d + b a u v + c -u a w + d -v -w a ], + + where a, b, c, d, u, v, w are real numbers. + + Since det(e^A) = e^tr(A), we can find a from the determinant of m, + using tr(K) = 4 a. + */ + + const Numeric a = 0.25 * std::log(det(m)); + + /** + We use this knowing + m = exp(K) + = exp(a) * (C0 I + C1 K + C2 K^2 + C3 K^3) + + where C0, C1, C2, C3 are the Cayley Hamilton coefficients derived in tran::tran. + + This means that we can extract + T = m * exp(-a) + = C0 I + C1 K + C2 K^2 + C3 K^3 + using the same logic and notation. + */ + + const muelmat T = m * std::exp(-a); + + /** + + These coefficients C0, C1, C2, C3 are found from the eigenvalues of K, which + have the form: lambda = exp(x), exp(-x), exp(iy), exp(-iy), where x and y are real. + + Since the trace of a matrix is the sum of its eigenvalues, and the eigenvalues + of a squared matrix are the squares of the eigenvalues of the original matrix, + we can form these two relations from our known T matrix: + S1 = tr(T) + = e^x + e^-x + e^iy + e^-iy + = 2 cosh(x) + 2 cos(y) + and + S2 = tr(T^2) + = (e^x + e^-x)^2 + (e^iy + e^-iy)^2 + = 2 cosh(2x) + 2 cos(2y) + = 2 (2 cosh^2(x) - 1) + 2 (2 cos^2(y) - 1) + + Letting u = cosh(x), v = cos(y), we can write: + S1 = 2 u + 2 v -> + u + v = S1 / 2 + and + S2 = 2 (2 u^2 - 1) + 2 (2 v^2 - 1) + = 4 u^2 + 4 v^2 - 4 -> + u^2 + v^2 = (S2 + 4) / 4 + + Since + (u - v)^2 = 2 (u^2 + v^2) - (u + v)^2 + we can solve for u and v: + (u - v)^2 = 2 * (S2 + 4) / 4 - (S1 / 2)^2 + = (S2 + 8 - S1^2) / 4 + and write + D = (u - v)^2 + = (S2 + 8 - S1^2) / 4 + */ + + const Numeric S1 = tr(T); + const Numeric D = 0.25 * (2 * tr(T * T) + 8 - S1 * S1); + + /** + We can use D to find u and v: + u = (u + v + u - v) / 2 + = (S1 / 2 + sqrt(D)) / 2 + and + v = (u + v - (u - v)) / 2 + = (S1 / 2 - sqrt(D)) / 2 + Since u - v = sqrt(D) + */ + + const Numeric sqrtD = D > 0 ? std::sqrt(D) : 0.0; + const Numeric u = std::max((S1 * 0.5 + sqrtD) * 0.5, 1.0); + const Numeric v = std::clamp((S1 * 0.5 - sqrtD) * 0.5, -1.0, 1.0); + + /** + From the definition of u and v we can now define x and y: + x = acosh(u) + y = acos(v) + */ + + const Numeric x = std::acosh(u); + const Numeric y = std::acos(v); + + // From here, we simply compute the coefficients C0, C1, C2, C3 as in tran::tran + const Numeric x2 = x * x; + const Numeric y2 = y * y; + const bool x_zero = x < 1e-4; + const bool y_zero = y < 1e-4; + const bool both_zero = x_zero && y_zero; + const bool either_zero = x_zero || y_zero; + const Numeric cy = v; + const Numeric sy = std::sin(y); + const Numeric cx = u; + const Numeric sx = std::sinh(x); + const Numeric ix = x_zero ? 0.0 : 1.0 / x; + const Numeric iy = y_zero ? 0.0 : 1.0 / y; + const Numeric inv_x2y2 = both_zero ? 1.0 : 1.0 / (x2 + y2); + const Numeric C0 = either_zero ? 1.0 : (cy * x2 + cx * y2) * inv_x2y2; + const Numeric C1 = + either_zero ? 1.0 : (sy * x2 * iy + sx * y2 * ix) * inv_x2y2; + // Skipping C2 and C3 + + /** + Now we are just one step away from reconstructing log(T) and thus K from exp(K). + + We have: + T = C0 I + C1 K + C2 K^2 + C3 K^3 + + We write: + T02 = C0 I + C2 K^2 -> + K^2 = (T02 - C0 I) / C2 + and + T13 = C1 K + C3 K^3 -> + T13 = K (C1 I + C3 K^2) + + So that: + K = T13 * inv(C1 I + C3 K^2) + = T13 * inv(C1 I + C3 (T02 - C0 I) / C2) + + OK. We just have to be able to construct T02 and T13 from T. + This requires some separation of terms in T. + + For T02, + [ T[0,0] (T[0,1]-T[1,0])/2 (T[0,2]-T[2,0])/2 (T[0,3]-T[3,0])/2 + -(T[1,0]-T[0,1])/2 T[1,1] (T[1,2]+T[2,1])/2 (T[1,3]+T[3,1])/2 + -(T[2,0]-T[0,2])/2 (T[2,1]+T[1,2])/2 T[2,2] (T[2,3]+T[3,2])/2 + -(T[3,0]-T[0,3])/2 (T[3,1]+T[1,3])/2 (T[3,2]+T[2,3])/2 T[3,3] ] + For T13, + [ 0 (T[0,1]+T[1,0])/2 (T[0,2]+T[2,0])/2 (T[0,3]+T[3,0])/2 + (T[1,0]+T[0,1])/2 0 (T[1,2]-T[2,1])/2 (T[1,3]-T[3,1])/2 + (T[2,0]+T[0,2])/2 -(T[2,1]-T[1,2])/2 0 (T[2,3]-T[3,2])/2 + (T[3,0]+T[0,3])/2 -(T[3,1]-T[1,3])/2 -(T[3,2]-T[2,3])/2 0 ] + + (To make the above symmetries and anti-symmetries clear, the return value of + tran::operator() has been formatted to demonstrate them.) + */ + + const Numeric factor = + both_zero ? 1.0 / 3.0 + : ((x_zero ? 1.0 : sx * ix) - (y_zero ? 1.0 : sy * iy)) / sqrtD; + const Numeric z01 = factor * std::midpoint(T[0, 1], -T[1, 0]); + const Numeric z02 = factor * std::midpoint(T[0, 2], -T[2, 0]); + const Numeric z03 = factor * std::midpoint(T[0, 3], -T[3, 0]); + const Numeric z12 = factor * std::midpoint(T[1, 2], T[2, 1]); + const Numeric z13 = factor * std::midpoint(T[1, 3], T[3, 1]); + const Numeric z23 = factor * std::midpoint(T[2, 3], T[3, 2]); + const Numeric z00 = factor * (T[0, 0] - C0) + C1; + const Numeric z11 = factor * (T[1, 1] - C0) + C1; + const Numeric z22 = factor * (T[2, 2] - C0) + C1; + const Numeric z33 = factor * (T[3, 3] - C0) + C1; + const Numeric b01 = std::midpoint(T[0, 1], T[1, 0]); + const Numeric b02 = std::midpoint(T[0, 2], T[2, 0]); + const Numeric b03 = std::midpoint(T[0, 3], T[3, 0]); + const Numeric b12 = std::midpoint(T[1, 2], -T[2, 1]); + const Numeric b13 = std::midpoint(T[1, 3], -T[3, 1]); + const Numeric b23 = std::midpoint(T[2, 3], -T[3, 2]); + + // clang-format off + const muelmat Kp = muelmat{ 0.0, b01, b02, b03, + b01, 0.0, b12, b13, + b02, -b12, 0.0, b23, + b03, -b13, -b23, 0.0} * inv(muelmat + { z00, z01, z02, z03, + -z01, z11, z12, z13, + -z02, z12, z22, z23, + -z03, z13, z23, z33}); + // clang-format on + + return {a, + std::midpoint(Kp[0, 1], Kp[1, 0]), + std::midpoint(Kp[0, 2], Kp[2, 0]), + std::midpoint(Kp[0, 3], Kp[3, 0]), + std::midpoint(Kp[1, 2], -Kp[2, 1]), + std::midpoint(Kp[1, 3], -Kp[3, 1]), + std::midpoint(Kp[2, 3], -Kp[3, 2])}; +} } // namespace rtepack diff --git a/src/core/rtepack/rtepack_transmission.h b/src/core/rtepack/rtepack_transmission.h index 8fec15f2c2..fe8dcc07e9 100644 --- a/src/core/rtepack/rtepack_transmission.h +++ b/src/core/rtepack/rtepack_transmission.h @@ -41,6 +41,8 @@ void two_level_exp(muelmat &t, muelmat exp(propmat k, Numeric r = 1.0); +propmat logK(const muelmat& m); + void two_level_exp(muelmat_vector_view t, const propmat_vector_const_view &k1, const propmat_vector_const_view &k2, diff --git a/src/python_interface/py_rtepack.cpp b/src/python_interface/py_rtepack.cpp index 4b386809d6..3d44de7e12 100644 --- a/src/python_interface/py_rtepack.cpp +++ b/src/python_interface/py_rtepack.cpp @@ -501,6 +501,49 @@ void py_rtepack(py::module_ &m) try { "dJs"_a, "I0"_a, "Returns the two-level radiative transfer of the input matrices"); + + auto rp = m.def_submodule("rtepack"); + auto tr = py::class_(rp, "tran"); + generic_interface(tr); + tr.def(py::init(), "k1"_a, "k2"_a, "r"_a) + .def("__call__", &rtepack::tran::operator(), "Returns the Mueller matrix") + .def("deriv", + &rtepack::tran::deriv, + "t"_a, + "k1"_a, + "k2"_a, + "dk"_a, + "r"_a, + "dr"_a, + "Returns the derivative of the Mueller matrix"); + rp.def("logK", + &rtepack::logK, + "m"_a, + R"(Returns the logarithm of the Mueller matrix as a Propagation matrix + +This only works for a Mueller matrix that has been generated from a Propagation matrix using + +.. math:: + \mathrm{M} = e^{-\mathrm{K} r}, + +where :math:`\mathrm{M}` is the Mueller matrix, :math:`\mathrm{K}` is the Propagation matrix, and :math:`r` is the distance. + +The return value of this method is then not the original Propagation matrix, but rather the logarithm of the Mueller matrix +as if it were generated from a Propagation matrix using the same equation. That is + +.. math:: + -\mathrm{K}r = \log(\mathrm{M} + +Parameters +---------- +m : Muelmat + The Mueller matrix to compute the logarithm of. + +Returns +------- +Propmat + The logarithm of the Mueller matrix as a Propagation matrix. +)"); } catch (std::exception &e) { throw std::runtime_error( std::format("DEV ERROR:\nCannot initialize rtepack\n{}", e.what())); From aa296dc965b5c00770481a7beabd9af50c6179bb Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Wed, 17 Dec 2025 19:43:33 +0900 Subject: [PATCH 2/5] Add test and fix limit --- src/core/rtepack/rtepack_transmission.cc | 10 ++- tests/python/math/logK.py | 81 ++++++++++++++++++++++++ 2 files changed, 85 insertions(+), 6 deletions(-) create mode 100644 tests/python/math/logK.py diff --git a/src/core/rtepack/rtepack_transmission.cc b/src/core/rtepack/rtepack_transmission.cc index 6ac8c3acf0..befe3d19ac 100644 --- a/src/core/rtepack/rtepack_transmission.cc +++ b/src/core/rtepack/rtepack_transmission.cc @@ -528,7 +528,6 @@ propmat logK(const muelmat &m) { const muelmat T = m * std::exp(-a); /** - These coefficients C0, C1, C2, C3 are found from the eigenvalues of K, which have the form: lambda = exp(x), exp(-x), exp(iy), exp(-iy), where x and y are real. @@ -634,16 +633,15 @@ propmat logK(const muelmat &m) { -(T[3,0]-T[0,3])/2 (T[3,1]+T[1,3])/2 (T[3,2]+T[2,3])/2 T[3,3] ] For T13, [ 0 (T[0,1]+T[1,0])/2 (T[0,2]+T[2,0])/2 (T[0,3]+T[3,0])/2 - (T[1,0]+T[0,1])/2 0 (T[1,2]-T[2,1])/2 (T[1,3]-T[3,1])/2 + (T[1,0]+T[0,1])/2 0 (T[1,2]-T[2,1])/2 (T[1,3]-T[3,1])/2 (T[2,0]+T[0,2])/2 -(T[2,1]-T[1,2])/2 0 (T[2,3]-T[3,2])/2 (T[3,0]+T[0,3])/2 -(T[3,1]-T[1,3])/2 -(T[3,2]-T[2,3])/2 0 ] (To make the above symmetries and anti-symmetries clear, the return value of tran::operator() has been formatted to demonstrate them.) */ - const Numeric factor = - both_zero ? 1.0 / 3.0 + (both_zero or sqrtD <= 0) ? 1.0 / 3.0 : ((x_zero ? 1.0 : sx * ix) - (y_zero ? 1.0 : sy * iy)) / sqrtD; const Numeric z01 = factor * std::midpoint(T[0, 1], -T[1, 0]); const Numeric z02 = factor * std::midpoint(T[0, 2], -T[2, 0]); @@ -666,8 +664,8 @@ propmat logK(const muelmat &m) { const muelmat Kp = muelmat{ 0.0, b01, b02, b03, b01, 0.0, b12, b13, b02, -b12, 0.0, b23, - b03, -b13, -b23, 0.0} * inv(muelmat - { z00, z01, z02, z03, + b03, -b13, -b23, 0.0} * + inv(muelmat{ z00, z01, z02, z03, -z01, z11, z12, z13, -z02, z12, z22, z23, -z03, z13, z23, z33}); diff --git a/tests/python/math/logK.py b/tests/python/math/logK.py new file mode 100644 index 0000000000..13d5fbbce2 --- /dev/null +++ b/tests/python/math/logK.py @@ -0,0 +1,81 @@ +import pyarts3 as pyarts +import numpy as np + +ln = pyarts.arts.rtepack.logK +tra = pyarts.arts.rtepack.tran + +made_it = False + +# Since random noise, we try multiple times to avoid occasional failures +for i in range(10): + try: + arr = np.array([1, 0, 0, 0, 0, 0, 0]) * 1e-2 + + K = pyarts.arts.Propmat(arr) + expK = tra(K, K, 1.0)() + logK = -ln(expK) + print(f"Testing logK with arr=\n{arr}") + print(logK) + assert np.allclose(K, logK, atol=1e-5) + + arr[1] = np.random.random() * 1e-3 + + K = pyarts.arts.Propmat(arr) + expK = tra(K, K, 1.0)() + logK = -ln(expK) + print(f"Testing logK with arr=\n{arr}") + print(logK) + assert np.allclose(K, logK, atol=1e-5) + + arr[2] = np.random.random() * 1e-3 + + K = pyarts.arts.Propmat(arr) + expK = tra(K, K, 1.0)() + logK = -ln(expK) + print(f"Testing logK with arr=\n{arr}") + print(logK) + assert np.allclose(K, logK, atol=1e-5) + + arr[3] = np.random.random() * 1e-3 + + K = pyarts.arts.Propmat(arr) + expK = tra(K, K, 1.0)() + logK = -ln(expK) + print(f"Testing logK with arr=\n{arr}") + print(logK) + assert np.allclose(K, logK, atol=1e-5) + + arr[4] = np.random.random() * 1e-3 + + K = pyarts.arts.Propmat(arr) + expK = tra(K, K, 1.0)() + logK = -ln(expK) + print(f"Testing logK with arr=\n{arr}") + print(logK) + assert np.allclose(K, logK, atol=1e-5) + + arr[5] = np.random.random() * 1e-3 + + K = pyarts.arts.Propmat(arr) + expK = tra(K, K, 1.0)() + logK = -ln(expK) + print(f"Testing logK with arr=\n{arr}") + print(logK) + assert np.allclose(K, logK, atol=1e-5) + + arr[6] = np.random.random() * 1e-3 + + K = pyarts.arts.Propmat(arr) + expK = tra(K, K, 1.0)() + logK = -ln(expK) + print(f"Testing logK with arr=\n{arr}") + print(logK) + assert np.allclose(K, logK, atol=1e-5) + made_it = True + print("logK tests passed at attempt", i+1) + break + except AssertionError: + print("logK test failed at attempt", i+1) + pass + +assert made_it, "logK tests failed after 10 attempts" From aa267e5d7fb6fff41694746eeb0189db07cd37f8 Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Wed, 17 Dec 2025 20:24:28 +0900 Subject: [PATCH 3/5] Simplify some of the core maths --- src/core/rtepack/rtepack_mueller_matrix.h | 62 ++++++++--------------- 1 file changed, 22 insertions(+), 40 deletions(-) diff --git a/src/core/rtepack/rtepack_mueller_matrix.h b/src/core/rtepack/rtepack_mueller_matrix.h index 86f0280d81..79663c6713 100644 --- a/src/core/rtepack/rtepack_mueller_matrix.h +++ b/src/core/rtepack/rtepack_mueller_matrix.h @@ -170,12 +170,10 @@ constexpr muelmat avg(muelmat a, const muelmat &b) { constexpr Numeric det(const muelmat &A) { const auto [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p] = A; - return a * f * k * p - a * f * l * o - a * g * j * p + a * g * l * n + - a * h * j * o - a * h * k * n - b * e * k * p + b * e * l * o + - b * g * i * p - b * g * l * m - b * h * i * o + b * h * k * m + - c * e * j * p - c * e * l * n - c * f * i * p + c * f * l * m + - c * h * i * n - c * h * j * m - d * e * j * o + d * e * k * n + - d * f * i * o - d * f * k * m - d * g * i * n + d * g * j * m; + return a * (f * (k * p - l * o) + g * (l * n - j * p) + h * (j * o - k * n)) + + b * (e * (l * o - k * p) + g * (i * p - l * m) + h * (k * m - i * o)) + + c * (e * (j * p - l * n) + f * (l * m - i * p) + h * (i * n - j * m)) + + d * (e * (k * n - j * o) + f * (i * o - k * m) + g * (j * m - i * n)); } constexpr Numeric tr(const muelmat &A) { @@ -190,40 +188,24 @@ constexpr Numeric midtr(const muelmat &A) { constexpr muelmat inv(const muelmat &A) { const auto [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p] = A; - const Numeric div = 1.0 / det(A); - - return div * muelmat{(f * k * p - f * l * o - g * j * p + g * l * n + - h * j * o - h * k * n), - (-b * k * p + b * l * o + c * j * p - c * l * n - - d * j * o + d * k * n), - (b * g * p - b * h * o - c * f * p + c * h * n + - d * f * o - d * g * n), - (-b * g * l + b * h * k + c * f * l - c * h * j - - d * f * k + d * g * j), - (-e * k * p + e * l * o + g * i * p - g * l * m - - h * i * o + h * k * m), - (a * k * p - a * l * o - c * i * p + c * l * m + - d * i * o - d * k * m), - (-a * g * p + a * h * o + c * e * p - c * h * m - - d * e * o + d * g * m), - (a * g * l - a * h * k - c * e * l + c * h * i + - d * e * k - d * g * i), - (e * j * p - e * l * n - f * i * p + f * l * m + - h * i * n - h * j * m), - (-a * j * p + a * l * n + b * i * p - b * l * m - - d * i * n + d * j * m), - (a * f * p - a * h * n - b * e * p + b * h * m + - d * e * n - d * f * m), - (-a * f * l + a * h * j + b * e * l - b * h * i - - d * e * j + d * f * i), - (-e * j * o + e * k * n + f * i * o - f * k * m - - g * i * n + g * j * m), - (a * j * o - a * k * n - b * i * o + b * k * m + - c * i * n - c * j * m), - (-a * f * o + a * g * n + b * e * o - b * g * m - - c * e * n + c * f * m), - (a * f * k - a * g * j - b * e * k + b * g * i + - c * e * j - c * f * i)}; + return muelmat{ + f * (k * p - l * o) + g * (l * n - j * p) + h * (j * o - k * n), + b * (l * o - k * p) + c * (j * p - l * n) + d * (k * n - j * o), + b * (g * p - h * o) + c * (h * n - f * p) + d * (f * o - g * n), + b * (h * k - g * l) + c * (f * l - h * j) + d * (g * j - f * k), + e * (l * o - k * p) + g * (i * p - l * m) + h * (k * m - i * o), + a * (k * p - l * o) + c * (l * m - i * p) + d * (i * o - k * m), + a * (h * o - g * p) + c * (e * p - h * m) + d * (g * m - e * o), + a * (g * l - h * k) + c * (h * i - e * l) + d * (e * k - g * i), + e * (j * p - l * n) + f * (l * m - i * p) + h * (i * n - j * m), + a * (l * n - j * p) + b * (i * p - l * m) + d * (j * m - i * n), + a * (f * p - h * n) + b * (h * m - e * p) + d * (e * n - f * m), + a * (h * j - f * l) + b * (e * l - h * i) + d * (f * i - e * j), + e * (k * n - j * o) + f * (i * o - k * m) + g * (j * m - i * n), + a * (j * o - k * n) + b * (k * m - i * o) + c * (i * n - j * m), + a * (g * n - f * o) + b * (e * o - g * m) + c * (f * m - e * n), + a * (f * k - g * j) + b * (g * i - e * k) + c * (e * j - f * i)} / + det(A); } using muelmat_vector = matpack::data_t; From 4eb2aa11f01f6205cb79998f711dfe1744c7abf3 Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Wed, 17 Dec 2025 20:29:08 +0900 Subject: [PATCH 4/5] Fix paranthesis --- src/python_interface/py_rtepack.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/python_interface/py_rtepack.cpp b/src/python_interface/py_rtepack.cpp index 3d44de7e12..33ed5c13a4 100644 --- a/src/python_interface/py_rtepack.cpp +++ b/src/python_interface/py_rtepack.cpp @@ -532,7 +532,7 @@ The return value of this method is then not the original Propagation matrix, but as if it were generated from a Propagation matrix using the same equation. That is .. math:: - -\mathrm{K}r = \log(\mathrm{M} + -\mathrm{K}r = \log(\mathrm{M}) Parameters ---------- From ff23bcd192e3d4a9a787b839965005cdeee6ce31 Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Wed, 17 Dec 2025 20:29:56 +0900 Subject: [PATCH 5/5] Add doc-string --- src/python_interface/py_rtepack.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/python_interface/py_rtepack.cpp b/src/python_interface/py_rtepack.cpp index 33ed5c13a4..eff0b4aedb 100644 --- a/src/python_interface/py_rtepack.cpp +++ b/src/python_interface/py_rtepack.cpp @@ -505,6 +505,7 @@ void py_rtepack(py::module_ &m) try { auto rp = m.def_submodule("rtepack"); auto tr = py::class_(rp, "tran"); generic_interface(tr); + tr.doc() = "Class for computing the transmission Mueller matrix and its derivative"; tr.def(py::init(), "k1"_a, "k2"_a, "r"_a) .def("__call__", &rtepack::tran::operator(), "Returns the Mueller matrix") .def("deriv",