Implement and explain linear in tau RT approximation as WIP (#1075)

This is still a bit of a WIP. But I have implemented linear in tau
approximation in fully polarized form.

I also took the time to clean up some documentation. I specifically
fixed some weirdities in the description of the Jacobian calculations.
(It's now more 1-to-1 with how the code works.)

I also describes the new implementation, but this part specifically is
not tested. The forwards is, but only for a few cases.

I need this merged because I need to be able to think of abstractions
around the RTE. Currently, this adds variables, which breaks workspace.
I could introduce a class to wrap a lot of the RTE, like we did for the
atmosphere, but it's still not beautiful to me how that would work.

Author: riclarsson

doc/arts/concept.radiative_transfer.rst

@@ -127,6 +127,25 @@ struct muelmat final : mat44 {
                     a30 * b03 + a31 * b13 + a32 * b23 + a33 * b33};
   }
 
+  constexpr muelmat operator-() const {
+    return muelmat{-data[0],
+                   -data[1],
+                   -data[2],
+                   -data[3],
+                   -data[4],
+                   -data[5],
+                   -data[6],
+                   -data[7],
+                   -data[8],
+                   -data[9],
+                   -data[10],
+                   -data[11],
+                   -data[12],
+                   -data[13],
+                   -data[14],
+                   -data[15]};
+  }
+
   [[nodiscard]] constexpr bool is_polarized() const noexcept {
     return data[1] != 0.0 or data[2] != 0.0 or data[3] != 0.0 or
            data[4] != 0.0 or data[6] != 0.0 or data[7] != 0.0 or

src/core/rtepack/rtepack_mueller_matrix.h

@@ -1,5 +1,6 @@
 #pragma once
 
+#include "configtypes.h"
 #include "rtepack_concepts.h"
 #include "rtepack_mueller_matrix.h"
 #include "rtepack_propagation_matrix.h"
@@ -20,33 +21,40 @@ constexpr stokvec absvec(const propmat &k) {
 //! Treat a list of propagation matrices as a list of Stokes (absorption) vectors
 stokvec_vector absvec(const propmat_vector_const_view &k);
 
-//! Get the inverse of a propagation matrix as a Mueller matrix
-constexpr muelmat inv(const propmat &k) {
+constexpr muelmat adj(const propmat &k) {
   const auto &[a, b, c, d, u, v, w] = k.data;
 
-  const Numeric div =
-      1.0 / (a * a * (a * a - b * b - c * c - d * d + u * u + v * v + w * w) -
-             b * w * (b * w - 2 * c * v + 2 * d * u) +
-             c * v * (-c * v + 2 * d * u) - d * d * u * u);
-
-  return {a * (a * a + u * u + v * v + w * w) * div,
-          -(a * (a * b + c * u + d * v) + w * (b * w - c * v + d * u)) * div,
-          (a * (-a * c + b * u - d * w) + v * (b * w - c * v + d * u)) * div,
-          (a * (-a * d + b * v + c * w) - u * (b * w - c * v + d * u)) * div,
-          (a * (-a * b + c * u + d * v) - w * (b * w - c * v + d * u)) * div,
-          a * (a * a - c * c - d * d + w * w) * div,
-          (d * (b * w - c * v + d * u) - a * (a * u - b * c + v * w)) * div,
-          (a * (-a * v + b * d + u * w) - c * (b * w - c * v + d * u)) * div,
-          (v * (b * w - c * v + d * u) - a * (a * c + b * u - d * w)) * div,
-          (a * (a * u + b * c - v * w) - d * (b * w - c * v + d * u)) * div,
-          a * (a * a - b * b - d * d + v * v) * div,
-          (b * (b * w - c * v + d * u) - a * (a * w - c * d + u * v)) * div,
-          -(a * (a * d + b * v + c * w) + u * (b * w - c * v + d * u)) * div,
-          (a * (a * v + b * d + u * w) + c * (b * w - c * v + d * u)) * div,
-          -(a * (-a * w - c * d + u * v) + b * (b * w - c * v + d * u)) * div,
-          a * (a * a - b * b - c * c + u * u) * div};
+  const Numeric a2 = a * a;
+  const Numeric b2 = b * b;
+  const Numeric c2 = c * c;
+  const Numeric d2 = d * d;
+  const Numeric u2 = u * u;
+  const Numeric v2 = v * v;
+  const Numeric w2 = w * w;
+
+  const Numeric C = b * w - c * v + d * u;
+
+  return muelmat{a * (a2 + u2 + v2 + w2),
+                 -a * (a * b + c * u + d * v) - w * C,
+                 a * (-a * c + b * u - d * w) + v * C,
+                 a * (-a * d + b * v + c * w) - u * C,
+                 a * (-a * b + c * u + d * v) - w * C,
+                 a * (a2 - c2 - d2 + w2),
+                 d * C - a * (a * u - b * c + v * w),
+                 a * (-a * v + b * d + u * w) - c * C,
+                 v * C - a * (a * c + b * u - d * w),
+                 a * (a * u + b * c - v * w) - d * C,
+                 a * (a2 - b2 - d2 + v2),
+                 b * C - a * (a * w - c * d + u * v),
+                 -a * (a * d + b * v + c * w) - u * C,
+                 a * (a * v + b * d + u * w) + c * C,
+                 -a * (-a * w - c * d + u * v) - b * C,
+                 a * (a2 - b2 - c2 + u2)};
 }
 
+//! Get the inverse of a propagation matrix as a Mueller matrix
+constexpr muelmat inv(const propmat &k) { return adj(k) / det(k); }
+
 //! muelmat matrix multiplied by a stokvec vector
 constexpr stokvec operator*(const muelmat &a, const stokvec &b) {
   const auto

src/core/rtepack/rtepack_multitype.h

@@ -74,6 +74,22 @@ constexpr propmat operator/(propmat a, const propmat &b) {
   return a;
 }
 
+constexpr Numeric det(const propmat &k) {
+  const auto &[a, b, c, d, u, v, w] = k.data;
+
+  const Numeric a2 = a * a;
+  const Numeric b2 = b * b;
+  const Numeric c2 = c * c;
+  const Numeric d2 = d * d;
+  const Numeric u2 = u * u;
+  const Numeric v2 = v * v;
+  const Numeric w2 = w * w;
+
+  const Numeric C = b * w - c * v + d * u;
+
+  return a2 * (a2 - b2 - c2 - d2 + u2 + v2 + w2) - C * C;
+}
+
 //! Take the average of two propmat matrixes
 constexpr propmat avg(propmat a, const propmat &b) {
   a += b;

src/core/rtepack/rtepack_propagation_matrix.h

@@ -76,22 +76,77 @@ void two_level_linear_emission_step_by_step_full(
   for (Size iv = 0; iv < nv; iv++) {
     stokvec &Iv = I[iv];
     for (Size i = N - 2; i < N; i--) {
+      const muelmat &T = Ts[i + 1][iv];
       const stokvec Jv = avg(Js[i][iv], Js[i + 1][iv]);
 
       Iv -= Jv;
 
       for (Index iq = 0; iq < nq; iq++) {
         dI[i][iq, iv] +=
-            Pi[i][iv] *
-            (dTs[i][0, iq, iv] * Iv +
-             0.5 * (dJs[i][iq, iv] - Ts[i + 1][iv] * dJs[i][iq, iv]));
+            Pi[i][iv] * (dTs[i][0, iq, iv] * Iv +
+                         0.5 * (dJs[i][iq, iv] - T * dJs[i][iq, iv]));
+        dI[i + 1][iq, iv] +=
+            Pi[i][iv] * (dTs[i + 1][1, iq, iv] * Iv +
+                         0.5 * (dJs[i + 1][iq, iv] - T * dJs[i + 1][iq, iv]));
+      }
+
+      Iv = T * Iv + Jv;
+    }
+  }
+}
+
+void two_level_linear_in_J_step_by_step_full(
+    stokvec_vector &I,
+    std::vector<stokvec_matrix> &dI,
+    const std::vector<muelmat_vector> &Ts,
+    const std::vector<muelmat_vector> &Ls,
+    const std::vector<muelmat_vector> &Pi,
+    const std::vector<muelmat_tensor3> &dTs,
+    const std::vector<muelmat_tensor3> &dLs,
+    const std::vector<stokvec_vector> &Js,
+    const std::vector<stokvec_matrix> &dJs,
+    const stokvec_vector &I0) {
+  const Size nv = I0.size();
+  const Size N  = Ts.size();
+
+  //! FIXME: Add error checks for dLs and Ls sizes
+
+  I = I0;
+  dI.resize(N);
+
+  const Index nq = dJs[0].nrows();
+
+  for (auto &x : dI) {
+    x.resize(nq, nv);
+    x = 0.0;
+  }
+
+  if (N == 0) return;
+
+#pragma omp parallel for if (not arts_omp_in_parallel())
+  for (Size iv = 0; iv < nv; iv++) {
+    stokvec &IM = I[iv];
+    for (Size i = N - 2; i < N; i--) {
+      const stokvec &JN = Js[i + 1][iv];
+      const stokvec &JM = Js[i][iv];
+      const muelmat &T  = Ts[i + 1][iv];
+      const muelmat &L  = Ls[i + 1][iv];
+
+      const stokvec IMmJM = IM - JM;
+      const stokvec JMmJN = JM - JN;
+
+      for (Index iq = 0; iq < nq; iq++) {
+        dI[i][iq, iv] +=
+            Pi[i][iv] * (dJs[i][iq, iv] - L * dJs[i][iq, iv] +
+                         dTs[i][0, iq, iv] * IMmJM + dLs[i][0, iq, iv] * JMmJN);
+
         dI[i + 1][iq, iv] +=
             Pi[i][iv] *
-            (dTs[i + 1][1, iq, iv] * Iv +
-             0.5 * (dJs[i + 1][iq, iv] - Ts[i + 1][iv] * dJs[i + 1][iq, iv]));
+            (dTs[i + 1][1, iq, iv] * IMmJM + dLs[i + 1][1, iq, iv] * JMmJN +
+             L * dJs[i + 1][iq, iv] - T * dJs[i + 1][iq, iv]);
       }
 
-      Iv = Ts[i + 1][iv] * Iv + Jv;
+      IM = T * IMmJM + L * JMmJN + JN;
     }
   }
 }

src/core/rtepack/rtepack_rtestep.cc

@@ -45,6 +45,18 @@ void two_level_linear_emission_step_by_step_full(
     const std::vector<stokvec_matrix> &dJs,
     const stokvec_vector &I0);
 
+void two_level_linear_in_J_step_by_step_full(
+    stokvec_vector &I,
+    std::vector<stokvec_matrix> &dI,
+    const std::vector<muelmat_vector> &Ts,
+    const std::vector<muelmat_vector> &Ls,
+    const std::vector<muelmat_vector> &Pi,
+    const std::vector<muelmat_tensor3> &dTs,
+    const std::vector<muelmat_tensor3> &dLs,
+    const std::vector<stokvec_vector> &Js,
+    const std::vector<stokvec_matrix> &dJs,
+    const stokvec_vector &I0);
+
 void two_level_linear_emission_step_by_step_full(
     std::vector<stokvec_vector> &Is,
     const std::vector<muelmat_vector> &Ts,