Linear in propagation and source (#1077)

This introduce a RTE forward option that assumes both the source
function and the propagation matrix changes linearly along the path.
There is a gross simplification involved with the use of a Dawson
function call that integrates the exponential. This is implemented
elementwise, so it is likely it will cause problems for strong
polarization.

Author: riclarsson

doc/arts/concept.radiative_transfer.rst

@@ -460,3 +460,54 @@ The derivative contribution is then given by
       &+& \frac{\partial T_{N-1}} {\partial \vec{x}_{N}} \left(I_{N-1} - J_{N-1}\right)
       &+& \frac{\partial \Lambda_{N-1}} {\partial \vec{x}_{N}} \left(J_{N-1} - J_{N}\right)
     &\Bigr]
+
+Linear propagation matrix and source function
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This relaxes the assumption of constant propagation matrix within
+the layer but is otherwise similar to linear in source.
+The indexing is the same for the
+transmittance and incoming background radiation as above - in layers.  However, the source
+function and propagation matrix indexing below is now shifted to represent levels - the same indexing
+as we use for the spectral radiance term.
+
+Linear in source means that we assume :math:`J = J_0 + (J_1 - J_0) \frac{r}{r_x}`,
+where :math:`r` is the distance through the layer and :math:`r_x` is the total
+thickness of the layer.  With a linearly varying propagation matrix :math:`K = K_0 + (K_1 - K_0) \frac{r}{r_x}`, we define
+:math:`T_0 = \exp\left(- \overline{K_{0, 1}} r_x\right)`, with the overline indicate
+the arithmetic mean, and get
+
+.. math::
+  I_1 = J_1 + T_0 \left(I_0 - J_0\right) + \left(\Lambda_0^{(0)} + \Lambda_0^{(1)}\right) \left(J_0 - J_1\right).
+
+This is a classical solution to the radiative transfer equation to improve
+numerical stability when the propagation matrix :math:`K` is
+large.  Going through the motion of extending this to multiple steps as for the
+constant source case is then just a matter of adding 1 extra term
+per layer in the dot products above.  If we define
+
+.. math::
+  \Lambda_i = \Lambda_i^{(0)} + \Lambda_i^{(1)},
+
+where the first term is identical to the linear source function case,
+and the second term is a second order correction for the linear
+variation of the propagation matrix, we get
+
+
+.. math::
+
+  \Lambda_i^{(0)} &=& \frac{1}{r_i} \overline{K_{i, i+1}}^{-1} &\left(1 - T_i\right), \\
+  \Lambda_i^{(1)} &=& \frac{1}{r_i^2} \overline{K_{i, i+1}}^{-1} &D\left(\left[K_{i+1} - K_{i}\right] r_i \middle/ 2\right),
+
+where :math:`D` is the Dawson function of a matrix.
+
+The rest of the steps that follows are the same as for the linear source function
+case above, just replacing :math:`\Lambda_i` with the new definition.
+
+.. caution::
+
+  The Dawson function of a matrix is not implemented in ARTS yet.  Instead,
+  the element-wise Dawson function is used as an approximation.  This means
+  that it might be better to use the linear source function method.
+
+  Indeed, we do not know how well this approximation performs in practice.

examples/2-clearsky-radiative-transfer/2-zeeman/10-zeeman-linear-prop.py

@@ -0,0 +1,71 @@
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pyarts3 as pyarts
+
+# Download catalogs
+pyarts.data.download()
+
+ws = pyarts.workspace.Workspace()
+
+# %% Sampled frequency range
+line_f0 = 118750348044.712
+ws.freq_grid = np.linspace(-50e6, 50e6, 1001) + line_f0
+
+# %% Species and line absorption
+ws.abs_speciesSet(species=["O2-66"])
+ws.ReadCatalogData()
+ws.abs_bandsSelectFrequencyByLine(fmin=40e9, fmax=120e9)
+ws.abs_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9)
+ws.WignerInit()
+
+# %% Use the automatic agenda setter for propagation matrix calculations
+ws.spectral_propmat_agendaAuto()
+
+# %% Grids and planet
+ws.surf_fieldPlanet(option="Earth")
+ws.surf_field[pyarts.arts.SurfaceKey("t")] = 295.0
+ws.atm_fieldRead(
+    toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
+)
+ws.atm_fieldSchmidthFieldFromIGRF(time="2000-03-11 14:39:37")
+
+# %% Checks and settings
+ws.spectral_rad_transform_operatorSet(option="Tb")
+
+# %% Core calculations
+pos = [100e3, 0, 0]
+los = [180.0, 0.0]
+ws.ray_pathGeometric(pos=pos, los=los, max_stepsize=1000.0)
+ws.spectral_radClearskyLinearInTauAndPropEmission()
+ws.spectral_radApplyUnitFromSpectralRadiance()
+
+# %% Show results
+fig, ax = pyarts.plot(ws.spectral_rad, freqs=(
+    ws.freq_grid - line_f0) / 1e6)
+[a.set_xlabel("Frequency offset [MHz]") for a in ax.flatten()]
+[a.set_ylabel("Spectral radiance [K]") for a in ax.flatten()]
+fig.suptitle(f"Zeeman effect of {round(line_f0 / 1e6)} MHz O$_2$ line")
+
+if "ARTS_HEADLESS" not in os.environ:
+    plt.show()
+
+# %% Test
+
+assert np.allclose(
+    ws.spectral_rad[::100],
+    np.array(
+        [[ 2.27830141e+02,  4.25489128e-04,  1.02599187e-04,  5.68881918e-02],
+         [ 2.30907190e+02,  6.59070807e-04,  1.59097448e-04,  7.03735321e-02],
+         [ 2.34847968e+02,  1.16702319e-03,  2.82212718e-04,  9.33486760e-02],
+         [ 2.40400744e+02,  2.61388422e-03,  6.34254900e-04,  1.40005198e-01],
+         [ 2.49821088e+02,  9.73795166e-03,  2.38733210e-03,  2.69700244e-01],
+         [ 2.08949873e+02,  2.37747404e+01,  1.64966784e+00,  5.47066239e-06],
+         [ 2.49820521e+02,  9.74149272e-03,  2.38822964e-03, -2.69777785e-01],
+         [ 2.40399650e+02,  2.61568960e-03,  6.34700330e-04, -1.40081856e-01],
+         [ 2.34846390e+02,  1.16822346e-03,  2.82506203e-04, -9.34247505e-02],
+         [ 2.30905136e+02,  6.59975316e-04,  1.59317615e-04, -7.04496271e-02],
+         [ 2.27827622e+02,  4.26220137e-04,  1.02776641e-04, -5.69642694e-02]]
+    ),
+), "Values have drifted from expected results in spectral radiance"

examples/2-clearsky-radiative-transfer/2-zeeman/9-zeeman-linear-src.py

@@ -56,16 +56,16 @@
 assert np.allclose(
     ws.spectral_rad[::100],
     np.array(
-        [[2.27742531e+02,  4.26846049e-04,  1.02930910e-04, 5.69747534e-02],
-         [2.30824032e+02,  6.60877182e-04,  1.59541966e-04, 7.04658183e-02],
-         [2.34769234e+02,  1.16969415e-03,  2.82878204e-04, 9.34339220e-02],
-         [2.40325915e+02,  2.61895652e-03,  6.35550519e-04, 1.40077284e-01],
-         [2.49750254e+02,  9.75738899e-03,  2.39258319e-03, 2.69811017e-01],
-         [2.09105759e+02,  2.38744101e+01,  1.64300556e+00, 5.42840543e-06],
-         [2.49749687e+02,  9.76093079e-03,  2.39348131e-03, -2.69888391e-01],
-         [2.40324823e+02,  2.62076265e-03,  6.35996235e-04, -1.40153790e-01],
-         [2.34767658e+02,  1.17089544e-03,  2.83171984e-04, -9.35098578e-02],
-         [2.30821982e+02,  6.61782801e-04,  1.59762430e-04, -7.05417543e-02],
-         [2.27740017e+02,  4.27578039e-04,  1.03108621e-04, -5.70506086e-02]]
+        [[ 2.27827140e+02,  4.25383147e-04,  1.02572526e-04,  5.68785875e-02],
+         [ 2.30903679e+02,  6.58907241e-04,  1.59056068e-04,  7.03617869e-02],
+         [ 2.34843814e+02,  1.16673604e-03,  2.82139468e-04,  9.33332773e-02],
+         [ 2.40395699e+02,  2.61325506e-03,  6.34092170e-04,  1.39982571e-01],
+         [ 2.49814592e+02,  9.73578709e-03,  2.38675820e-03,  2.69659250e-01],
+         [ 2.08948687e+02,  2.37725593e+01,  1.64969081e+00,  5.47121870e-06],
+         [ 2.49814025e+02,  9.73932755e-03,  2.38765555e-03, -2.69736783e-01],
+         [ 2.40394605e+02,  2.61506011e-03,  6.34537500e-04, -1.40059220e-01],
+         [ 2.34842236e+02,  1.16793609e-03,  2.82432888e-04, -9.34093435e-02],
+         [ 2.30901627e+02,  6.59811580e-04,  1.59276187e-04, -7.04378734e-02],
+         [ 2.27824621e+02,  4.26114017e-04,  1.02749941e-04, -5.69546563e-02]]
     ),
 ), "Values have drifted from expected results in spectral radiance"

src/core/rtepack/CMakeLists.txt

@@ -11,4 +11,4 @@ add_library(rtepack STATIC
   rtepack_spectral_matrix.cc
 )
 target_include_directories(rtepack PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
-target_link_libraries(rtepack PUBLIC matpack physics)
+target_link_libraries(rtepack PUBLIC matpack physics Faddeeva)

src/core/rtepack/rtepack_mueller_matrix.h

@@ -12,23 +12,38 @@ struct muelmat final : mat44 {
   constexpr muelmat(Numeric tau = 1.0) noexcept
       : mat44{tau, 0, 0, 0, 0, tau, 0, 0, 0, 0, tau, 0, 0, 0, 0, tau} {}
 
-  constexpr muelmat(Numeric a,
-                    Numeric b,
-                    Numeric c,
-                    Numeric d,
-                    Numeric e,
-                    Numeric f,
-                    Numeric g,
-                    Numeric h,
-                    Numeric i,
-                    Numeric j,
-                    Numeric k,
-                    Numeric l,
-                    Numeric m,
-                    Numeric n,
-                    Numeric o,
-                    Numeric p) noexcept
-      : mat44{a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p} {}
+  constexpr muelmat(Numeric m00,
+                    Numeric m01,
+                    Numeric m02,
+                    Numeric m03,
+                    Numeric m10,
+                    Numeric m11,
+                    Numeric m12,
+                    Numeric m13,
+                    Numeric m20,
+                    Numeric m21,
+                    Numeric m22,
+                    Numeric m23,
+                    Numeric m30,
+                    Numeric m31,
+                    Numeric m32,
+                    Numeric m33) noexcept
+      : mat44{m00,
+              m01,
+              m02,
+              m03,
+              m10,
+              m11,
+              m12,
+              m13,
+              m20,
+              m21,
+              m22,
+              m23,
+              m30,
+              m31,
+              m32,
+              m33} {}
 
   constexpr muelmat(std::array<Numeric, 16> data) noexcept : mat44{data} {}
 

src/core/rtepack/rtepack_multitype.h

@@ -128,6 +128,26 @@ constexpr stokvec operator*(const propmat &k, const stokvec s) {
                  a * s4 + d * s1 - s2 * v - s3 * w};
 }
 
+//! Element-wise product between a propmat and a muelmat
+constexpr muelmat elem_prod(const propmat &a, const muelmat &b) {
+  return muelmat{a.A() * b[0, 0],
+                 a.B() * b[0, 1],
+                 a.C() * b[0, 2],
+                 a.D() * b[0, 3],
+                 a.B() * b[1, 0],
+                 a.A() * b[1, 1],
+                 a.U() * b[1, 2],
+                 a.V() * b[1, 3],
+                 a.C() * b[2, 0],
+                 -a.U() * b[2, 1],
+                 a.A() * b[2, 2],
+                 a.W() * b[2, 3],
+                 a.D() * b[3, 0],
+                 -a.V() * b[3, 1],
+                 -a.W() * b[3, 2],
+                 a.A() * b[3, 3]};
+}
+
 //! Transform a matrix of shape (N, 4) to a list of Stokes (absorption) vectors
 stokvec_vector to_stokvec_vector(const ConstMatrixView &v);
 

src/core/rtepack/rtepack_propagation_matrix.cc

@@ -1,5 +1,7 @@
 #include "rtepack_propagation_matrix.h"
 
+#include <Faddeeva.hh>
+
 namespace rtepack {
 propmat_vector operator*(Numeric x, const propmat_vector_const_view &y) {
   propmat_vector z(y.size());
@@ -8,6 +10,16 @@ propmat_vector operator*(Numeric x, const propmat_vector_const_view &y) {
   }
   return z;
 }
+
+propmat dawson(const propmat &pm) {
+  return {Faddeeva::Dawson(pm.A()),
+          Faddeeva::Dawson(pm.B()),
+          Faddeeva::Dawson(pm.C()),
+          Faddeeva::Dawson(pm.D()),
+          Faddeeva::Dawson(pm.U()),
+          Faddeeva::Dawson(pm.V()),
+          Faddeeva::Dawson(pm.W())};
+}
 }  // namespace rtepack
 
 void xml_io_stream<rtepack::propmat>::write(std::ostream &os,

src/core/rtepack/rtepack_propagation_matrix.h

@@ -103,6 +103,9 @@ constexpr propmat avg(const propmat &a, const propmat &b) {
           std::midpoint(a.W(), b.W())};
 }
 
+//! Element-wise dawson function of a propmat matrix (FIXME: implement as matrix?)
+propmat dawson(const propmat &pm);
+
 using propmat_vector            = matpack::data_t<propmat, 1>;
 using propmat_vector_view       = matpack::view_t<propmat, 1>;
 using propmat_vector_const_view = matpack::view_t<const propmat, 1>;