Radiative-transfer implementation

The notation below follows the linear-source derivation in Radiative Transfer.

Implementation note

It is very important to implement the expression for \(\Lambda_0\) in a numerically stable way. The expression above is not stable for small \(K_0 r\) (i.e., large \(T_0\)) as written. The key instability stems from the subtraction of two nearly equal terms in \(1 - T_0\). The IEEE floating point standard provides a function expm1(x), which computes \(e^x - 1\) in a numerically stable way for small \(x\).

Likewise, the matrix expansion of \(1 - T_0\) might be unstable.

So we use a special solution implementing our own version of the reduced Cayley-Hamilton theorem to compute \(\Lambda_0\) in a numerically stable way for matrices that conform to the propagation matrix notation in ARTS. This makes use of the inversion of \(K\) to remove components from the expansion of the matrix exponential that would otherwise cause numerical instability.