transmissionmatrix.h

Go to the documentation of this file.

/* Copyright (C) 2018
 * Richard Larsson <ric.larsson@gmail.com>
 * 
 * This program is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation; either version 2, or (at your option) any
 * later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
 * USA. */
 
#ifndef transmissionmatrix_h
#define transmissionmatrix_h
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wshadow"
#pragma GCC diagnostic ignored "-Wconversion"
#pragma GCC diagnostic ignored "-Wfloat-conversion"
#include <Eigen/Dense>
#include <Eigen/StdVector>
#pragma GCC diagnostic pop
 
#include "jacobian.h"
#include "propagationmatrix.h"
 
struct TransmissionMatrix {
  Index stokes_dim;
  std::vector<Eigen::Matrix4d, Eigen::aligned_allocator<Eigen::Matrix4d>> T4;
  std::vector<Eigen::Matrix3d, Eigen::aligned_allocator<Eigen::Matrix3d>> T3;
  std::vector<Eigen::Matrix2d, Eigen::aligned_allocator<Eigen::Matrix2d>> T2;
  std::vector<Eigen::Matrix<double, 1, 1>,
              Eigen::aligned_allocator<Eigen::Matrix<double, 1, 1>>>
      T1;
 
  TransmissionMatrix(Index nf = 0, Index stokes = 1)
      : stokes_dim(stokes),
        T4(stokes_dim == 4 ? nf : 0, Eigen::Matrix4d::Identity()),
        T3(stokes_dim == 3 ? nf : 0, Eigen::Matrix3d::Identity()),
        T2(stokes_dim == 2 ? nf : 0, Eigen::Matrix2d::Identity()),
        T1(stokes_dim == 1 ? nf : 0, Eigen::Matrix<double, 1, 1>::Identity()) {
    ARTS_ASSERT(stokes_dim < 5 and stokes_dim > 0);
  }
 
  TransmissionMatrix(TransmissionMatrix&& tm) noexcept
      : stokes_dim(tm.stokes_dim),
        T4(std::move(tm.T4)),
        T3(std::move(tm.T3)),
        T2(std::move(tm.T2)),
        T1(std::move(tm.T1)) {}
 
  TransmissionMatrix(const TransmissionMatrix& tm) = default;
 
  TransmissionMatrix& operator=(const TransmissionMatrix& tm) = default;
 
  TransmissionMatrix& operator=(TransmissionMatrix&& tm) noexcept {
    stokes_dim = tm.stokes_dim;
    T4 = std::move(tm.T4);
    T3 = std::move(tm.T3);
    T2 = std::move(tm.T2);
    T1 = std::move(tm.T1);
    return *this;
  }
 
  explicit TransmissionMatrix(const PropagationMatrix& pm, const Numeric& r = 1.0);
 
  operator Tensor3() const {
    Tensor3 T(Frequencies(), stokes_dim, stokes_dim);
    for (size_t i = 0; i < T4.size(); i++)
      for (size_t j = 0; j < 4; j++)
        for (size_t k = 0; k < 4; k++) T(i, j, k) = T4[i](j, k);
    for (size_t i = 0; i < T3.size(); i++)
      for (size_t j = 0; j < 3; j++)
        for (size_t k = 0; k < 3; k++) T(i, j, k) = T3[i](j, k);
    for (size_t i = 0; i < T2.size(); i++)
      for (size_t j = 0; j < 2; j++)
        for (size_t k = 0; k < 2; k++) T(i, j, k) = T2[i](j, k);
    for (size_t i = 0; i < T1.size(); i++) T(i, 0, 0) = T1[i](0, 0);
    return T;
  }
 
  [[nodiscard]] const Eigen::Matrix4d& Mat4(size_t i) const { return T4[i]; }
 
  [[nodiscard]] const Eigen::Matrix3d& Mat3(size_t i) const { return T3[i]; }
 
  [[nodiscard]] const Eigen::Matrix2d& Mat2(size_t i) const { return T2[i]; }
 
  [[nodiscard]] const Eigen::Matrix<double, 1, 1>& Mat1(size_t i) const { return T1[i]; }
 
  #pragma GCC diagnostic push
  #pragma GCC diagnostic ignored "-Wreturn-type"
  [[nodiscard]] Eigen::MatrixXd Mat(size_t i) const {
    switch (stokes_dim) {
      case 1:
        return Mat1(i);
      case 2:
        return Mat2(i);
      case 3:
        return Mat3(i);
      case 4:
        return Mat4(i);
    }
  }
  #pragma GCC diagnostic pop
 
  Eigen::Matrix4d& Mat4(size_t i) { return T4[i]; }
 
  Eigen::Matrix3d& Mat3(size_t i) { return T3[i]; }
 
  Eigen::Matrix2d& Mat2(size_t i) { return T2[i]; }
 
  Eigen::Matrix<double, 1, 1>& Mat1(size_t i) { return T1[i]; }
 
  void setIdentity() {
    std::fill(T4.begin(), T4.end(), Eigen::Matrix4d::Identity());
    std::fill(T3.begin(), T3.end(), Eigen::Matrix3d::Identity());
    std::fill(T2.begin(), T2.end(), Eigen::Matrix2d::Identity());
    std::fill(T1.begin(), T1.end(), Eigen::Matrix<double, 1, 1>::Identity());
  }
 
  void setZero() {
    std::fill(T4.begin(), T4.end(), Eigen::Matrix4d::Zero());
    std::fill(T3.begin(), T3.end(), Eigen::Matrix3d::Zero());
    std::fill(T2.begin(), T2.end(), Eigen::Matrix2d::Zero());
    std::fill(T1.begin(), T1.end(), Eigen::Matrix<double, 1, 1>::Zero());
  }
 
  void mul(const TransmissionMatrix& A, const TransmissionMatrix& B) {
    for (size_t i = 0; i < T4.size(); i++) T4[i].noalias() = A.T4[i] * B.T4[i];
    for (size_t i = 0; i < T3.size(); i++) T3[i].noalias() = A.T3[i] * B.T3[i];
    for (size_t i = 0; i < T2.size(); i++) T2[i].noalias() = A.T2[i] * B.T2[i];
    for (size_t i = 0; i < T1.size(); i++) T1[i].noalias() = A.T1[i] * B.T1[i];
  }
 
 
  void mul_aliased(const TransmissionMatrix& A, const TransmissionMatrix& B) {
    for (size_t i = 0; i < T4.size(); i++) T4[i] = A.T4[i] * B.T4[i];
    for (size_t i = 0; i < T3.size(); i++) T3[i] = A.T3[i] * B.T3[i];
    for (size_t i = 0; i < T2.size(); i++) T2[i] = A.T2[i] * B.T2[i];
    for (size_t i = 0; i < T1.size(); i++) T1[i] = A.T1[i] * B.T1[i];
  }
 
  Numeric operator()(const Index i, const Index j, const Index k) const {
    switch (stokes_dim) {
      case 4:
        return T4[i](j, k);
      case 3:
        return T3[i](j, k);
      case 2:
        return T2[i](j, k);
      default:
        return T1[i](j, k);
    }
  }
 
  [[nodiscard]] Index Frequencies() const {
    switch (stokes_dim) {
      case 4:
        return Index(T4.size());
      case 3:
        return Index(T3.size());
      case 2:
        return Index(T2.size());
      default:
        return Index(T1.size());
    }
  }
 
  TransmissionMatrix& operator+=(const LazyScale<TransmissionMatrix>& lstm) {
    for (size_t i = 0; i < T4.size(); i++)
      T4[i].noalias() = lstm.scale * lstm.bas.Mat4(i);
    for (size_t i = 0; i < T3.size(); i++)
      T3[i].noalias() = lstm.scale * lstm.bas.Mat3(i);
    for (size_t i = 0; i < T2.size(); i++)
      T2[i].noalias() = lstm.scale * lstm.bas.Mat2(i);
    for (size_t i = 0; i < T1.size(); i++)
      T1[i].noalias() = lstm.scale * lstm.bas.Mat1(i);
    return *this;
  }
 
  TransmissionMatrix& operator*=(const Numeric& scale) {
    std::transform(T4.begin(), T4.end(), T4.begin(), [scale](auto& T){return T*scale;});
    std::transform(T3.begin(), T3.end(), T3.begin(), [scale](auto& T){return T*scale;});
    std::transform(T2.begin(), T2.end(), T2.begin(), [scale](auto& T){return T*scale;});
    std::transform(T1.begin(), T1.end(), T1.begin(), [scale](auto& T){return T*scale;});
    return *this;
  }
 
  TransmissionMatrix& operator=(const LazyScale<TransmissionMatrix>& lstm) {
    operator=(lstm.bas);
    operator*=(lstm.scale);
    return *this;
  }
 
  friend std::ostream& operator<<(std::ostream& os,
                                  const TransmissionMatrix& tm);
 
  friend std::istream& operator>>(std::istream& data, TransmissionMatrix& tm);
  
  template <int N> auto& TraMat(size_t i) noexcept {
    static_assert (N > 0 and N < 5, "Bad size N");
    if constexpr (N == 1) return T1[i];
    else if constexpr (N == 2) return T2[i];
    else if constexpr (N == 3) return T3[i];
    else if constexpr (N == 4) return T4[i];
  }
  
  template <int N> [[nodiscard]] auto& TraMat(size_t i) const noexcept {
    static_assert (N > 0 and N < 5, "Bad size N");
    if constexpr (N == 1) return T1[i];
    else if constexpr (N == 2) return T2[i];
    else if constexpr (N == 3) return T3[i];
    else if constexpr (N == 4) return T4[i];
  }
  
  template <int N> [[nodiscard]] Eigen::Matrix<Numeric, N, N> OptDepth(size_t i) const noexcept {
    return (-TraMat<N>(i).diagonal().array().log().matrix()).asDiagonal();
  }
  
  template <int N>
  [[nodiscard]] Eigen::Matrix<Numeric, N, 1> second_order_integration_source(const Eigen::Matrix<Numeric, N, N> T,
                                                                             const Eigen::Matrix<Numeric, N, 1> far,
                                                                             const Eigen::Matrix<Numeric, N, 1> close,
                                                                             const Eigen::Matrix<Numeric, N, N> Kfar,
                                                                             const Eigen::Matrix<Numeric, N, N> Kclose,
                                                                             const Numeric r)
  const noexcept {
    static_assert (N>0 and N<5, "Bad size N");
    
    const auto I = Eigen::Matrix<Numeric, N, N>::Identity();
    if (T(0, 0) < 0.99) {
      const auto od = 0.5 * r * (Kfar + Kclose);
      Eigen::Matrix<Numeric, N, 1> second = od.inverse() * ((I - (I + od) * T) * far + (od - I + T) * close);
      if ((far[0] < close[0] and Kfar(0, 0) > Kclose(0, 0)) or (far[0] > close[0] and Kfar(0, 0) < Kclose(0, 0))) {
        Eigen::Matrix<Numeric, N, 1> second_limit = 0.5 * (I - T) * (far + close);
        
        // FIXME: This is the equation given in the source material...
        // const Eigen::Matrix<Numeric, N, 1> second_limit = 0.25 * r * (Kfar * far + Kclose * close);
        
        if (second_limit[0] > second[0])
          return second;
        return second_limit;
      }
 
      return second;
    }
 
    return 0.5 * (I - T) * (far + close);
  }
  
  
  template <int N>
  [[nodiscard]] Eigen::Matrix<Numeric, N, 1> second_order_integration_dsource(size_t i,
                                                                              const TransmissionMatrix& dx,
                                                                              const Eigen::Matrix<Numeric, N, 1> far,
                                                                              const Eigen::Matrix<Numeric, N, 1> close,
                                                                              const Eigen::Matrix<Numeric, N, 1> d,
                                                                              bool isfar) const noexcept {
    static_assert (N > 0 and N < 5, "Bad size N");
    
    const auto I = Eigen::Matrix<Numeric, N, N>::Identity();
    const auto T = TraMat<N>(i);
    const auto& dTdx = dx.TraMat<N>(i);
    const Eigen::Matrix<Numeric, N, 1> first = 0.5 * (I - T) * (far + close);
    if (T(0, 0) < 0.99) {
      const auto od = OptDepth<N>(i);
      const auto doddx = dx.OptDepth<N>(i);
      const Eigen::Matrix<Numeric, N, 1> second = od.inverse() * ((I - (I + od) * T) * far + (od - I + T) * close);
      if (first[0] > second[0]) {
        if (isfar) {
          return od.inverse() * ((I - (I + od) * T) * d
                              - (I + od) * dTdx * far
                              + (doddx + T) * close
                              - doddx * second);
        }
 
        return od.inverse() * (
        - (I + od) * dTdx * far
        + (od - I + T) * d
        + (doddx + T) * close
        - doddx * second);
      }
 
      return 0.5 * ((I - T) * d - dTdx * (far + close));
    }
 
    return 0.5 * ((I - T) * d - dTdx * (far + close));
  }
};
 
inline LazyScale<TransmissionMatrix> operator*(const TransmissionMatrix& tm,
                                               const Numeric& x) {
  return LazyScale<TransmissionMatrix>(tm, x);
}
 
inline LazyScale<TransmissionMatrix> operator*(const Numeric& x,
                                               const TransmissionMatrix& tm) {
  return LazyScale<TransmissionMatrix>(tm, x);
}
 
struct RadiationVector {
  Index stokes_dim;
  std::vector<Eigen::Vector4d, Eigen::aligned_allocator<Eigen::Vector4d>> R4;
  std::vector<Eigen::Vector3d, Eigen::aligned_allocator<Eigen::Vector3d>> R3;
  std::vector<Eigen::Vector2d, Eigen::aligned_allocator<Eigen::Vector2d>> R2;
  std::vector<Eigen::Matrix<double, 1, 1>,
              Eigen::aligned_allocator<Eigen::Matrix<double, 1, 1>>>
      R1;
 
  RadiationVector(Index nf = 0, Index stokes = 1)
      : stokes_dim(stokes),
        R4(stokes_dim == 4 ? nf : 0, Eigen::Vector4d::Zero()),
        R3(stokes_dim == 3 ? nf : 0, Eigen::Vector3d::Zero()),
        R2(stokes_dim == 2 ? nf : 0, Eigen::Vector2d::Zero()),
        R1(stokes_dim == 1 ? nf : 0, Eigen::Matrix<double, 1, 1>::Zero()) {
    ARTS_ASSERT(stokes_dim < 5 and stokes_dim > 0);
  }
 
  RadiationVector(RadiationVector&& rv) noexcept
      : stokes_dim(rv.stokes_dim),
        R4(std::move(rv.R4)),
        R3(std::move(rv.R3)),
        R2(std::move(rv.R2)),
        R1(std::move(rv.R1)) {}
 
  RadiationVector(const RadiationVector& rv) = default;
 
  RadiationVector& operator=(const RadiationVector& rv) = default;
 
  RadiationVector& operator=(RadiationVector&& rv) noexcept {
    stokes_dim = rv.stokes_dim;
    R4 = std::move(rv.R4);
    R3 = std::move(rv.R3);
    R2 = std::move(rv.R2);
    R1 = std::move(rv.R1);
    return *this;
  }
 
  void leftMul(const TransmissionMatrix& T) {
    for (size_t i = 0; i < R4.size(); i++) R4[i] = T.Mat4(i) * R4[i];
    for (size_t i = 0; i < R3.size(); i++) R3[i] = T.Mat3(i) * R3[i];
    for (size_t i = 0; i < R2.size(); i++) R2[i] = T.Mat2(i) * R2[i];
    for (size_t i = 0; i < R1.size(); i++) R1[i] = T.Mat1(i) * R1[i];
  }
 
  void SetZero(size_t i) {
    switch (stokes_dim) {
      case 4:
        R4[i].noalias() = Eigen::Vector4d::Zero();
        break;
      case 3:
        R3[i].noalias() = Eigen::Vector3d::Zero();
        break;
      case 2:
        R2[i].noalias() = Eigen::Vector2d::Zero();
        break;
      case 1:
        R1[i][0] = 0;
        break;
    }
  }
 
  [[nodiscard]] const Eigen::Vector4d& Vec4(size_t i) const { return R4[i]; }
 
  [[nodiscard]] const Eigen::Vector3d& Vec3(size_t i) const { return R3[i]; }
 
  [[nodiscard]] const Eigen::Vector2d& Vec2(size_t i) const { return R2[i]; }
 
  [[nodiscard]] const Eigen::Matrix<double, 1, 1>& Vec1(size_t i) const { return R1[i]; }
  
  [[nodiscard]] Eigen::VectorXd Vec(size_t i) const { 
    switch (stokes_dim) {
      case 4:
        return Vec4(i);
      case 3:
        return Vec3(i);
      case 2:
        return Vec2(i);
      default:
        return Vec1(i);
    }
  }
 
 
  Eigen::Vector4d& Vec4(size_t i) { return R4[i]; }
 
  Eigen::Vector3d& Vec3(size_t i) { return R3[i]; }
 
  Eigen::Vector2d& Vec2(size_t i) { return R2[i]; }
 
  Eigen::Matrix<double, 1, 1>& Vec1(size_t i) { return R1[i]; }
 
  void rem_avg(const RadiationVector& O1, const RadiationVector& O2) {
    for (size_t i = 0; i < R4.size(); i++)
      R4[i].noalias() -= 0.5 * (O1.R4[i] + O2.R4[i]);
    for (size_t i = 0; i < R3.size(); i++)
      R3[i].noalias() -= 0.5 * (O1.R3[i] + O2.R3[i]);
    for (size_t i = 0; i < R2.size(); i++)
      R2[i].noalias() -= 0.5 * (O1.R2[i] + O2.R2[i]);
    for (size_t i = 0; i < R1.size(); i++)
      R1[i].noalias() -= 0.5 * (O1.R1[i] + O2.R1[i]);
  }
 
 
  void add_avg(const RadiationVector& O1, const RadiationVector& O2) {
    for (size_t i = 0; i < R4.size(); i++)
      R4[i].noalias() += 0.5 * (O1.R4[i] + O2.R4[i]);
    for (size_t i = 0; i < R3.size(); i++)
      R3[i].noalias() += 0.5 * (O1.R3[i] + O2.R3[i]);
    for (size_t i = 0; i < R2.size(); i++)
      R2[i].noalias() += 0.5 * (O1.R2[i] + O2.R2[i]);
    for (size_t i = 0; i < R1.size(); i++)
      R1[i].noalias() += 0.5 * (O1.R1[i] + O2.R1[i]);
  }
  
  void add_weighted(const TransmissionMatrix& T, const RadiationVector& far, const RadiationVector& close, 
                    const ConstMatrixView& Kfar, const ConstMatrixView& Kclose, const Numeric r) {
    for (size_t i = 0; i < R4.size(); i++) {
      R4[i].noalias() += T.second_order_integration_source<4>(T.TraMat<4>(i), far.R4[i], close.R4[i], prop_matrix<4>(Kfar(i, joker)), prop_matrix<4>(Kclose(i, joker)), r);
    }
    for (size_t i = 0; i < R3.size(); i++) {
      R3[i].noalias() += T.second_order_integration_source<3>(T.TraMat<3>(i), far.R3[i], close.R3[i], prop_matrix<3>(Kfar(i, joker)), prop_matrix<3>(Kclose(i, joker)), r);
    }
    for (size_t i = 0; i < R2.size(); i++) {
      R2[i].noalias() += T.second_order_integration_source<2>(T.TraMat<2>(i), far.R2[i], close.R2[i], prop_matrix<2>(Kfar(i, joker)), prop_matrix<2>(Kclose(i, joker)), r);
    }
    for (size_t i = 0; i < R1.size(); i++) {
      R1[i].noalias() += T.second_order_integration_source<1>(T.TraMat<1>(i), far.R1[i], close.R1[i], prop_matrix<1>(Kfar(i, joker)), prop_matrix<1>(Kclose(i, joker)), r);
    }
  }
 
  void addDerivEmission(const TransmissionMatrix& PiT,
                        const TransmissionMatrix& dT,
                        const TransmissionMatrix& T,
                        const RadiationVector& ImJ,
                        const RadiationVector& dJ) {
    for (size_t i = 0; i < R4.size(); i++)
      R4[i].noalias() += PiT.Mat4(i) * (dT.Mat4(i) * ImJ.R4[i] + dJ.R4[i] -
                                        T.Mat4(i) * dJ.R4[i]);
    for (size_t i = 0; i < R3.size(); i++)
      R3[i].noalias() += PiT.Mat3(i) * (dT.Mat3(i) * ImJ.R3[i] + dJ.R3[i] -
                                        T.Mat3(i) * dJ.R3[i]);
    for (size_t i = 0; i < R2.size(); i++)
      R2[i].noalias() += PiT.Mat2(i) * (dT.Mat2(i) * ImJ.R2[i] + dJ.R2[i] -
                                        T.Mat2(i) * dJ.R2[i]);
    for (size_t i = 0; i < R1.size(); i++)
      R1[i].noalias() += PiT.Mat1(i) * (dT.Mat1(i) * ImJ.R1[i] + dJ.R1[i] -
                                        T.Mat1(i) * dJ.R1[i]);
  }
 
  void addWeightedDerivEmission(const TransmissionMatrix& PiT,
                        const TransmissionMatrix& dT,
                        const TransmissionMatrix& T,
                        const RadiationVector& I,
                        const RadiationVector& far,
                        const RadiationVector& close,
                        const RadiationVector& d,
                        bool isfar) {
    for (size_t i = 0; i < R4.size(); i++)
      R4[i].noalias() += PiT.Mat4(i) * (dT.Mat4(i) * I.R4[i] + T.second_order_integration_dsource<4>(i, dT, far.R4[i], close.R4[i], d.R4[i], isfar));
    for (size_t i = 0; i < R3.size(); i++)
      R3[i].noalias() += PiT.Mat3(i) * (dT.Mat3(i) * I.R3[i] + T.second_order_integration_dsource<3>(i, dT, far.R3[i], close.R3[i], d.R3[i], isfar));
    for (size_t i = 0; i < R2.size(); i++)
      R2[i].noalias() += PiT.Mat2(i) * (dT.Mat2(i) * I.R2[i] + T.second_order_integration_dsource<2>(i, dT, far.R2[i], close.R2[i], d.R2[i], isfar));
    for (size_t i = 0; i < R1.size(); i++)
      R1[i].noalias() += PiT.Mat1(i) * (dT.Mat1(i) * I.R1[i] + T.second_order_integration_dsource<1>(i, dT, far.R1[i], close.R1[i], d.R1[i], isfar));
  }
 
 
  void addDerivTransmission(const TransmissionMatrix& PiT,
                            const TransmissionMatrix& dT,
                            const RadiationVector& I) {
    for (size_t i = 0; i < R4.size(); i++)
      R4[i].noalias() += PiT.Mat4(i) * dT.Mat4(i) * I.R4[i];
    for (size_t i = 0; i < R3.size(); i++)
      R3[i].noalias() += PiT.Mat3(i) * dT.Mat3(i) * I.R3[i];
    for (size_t i = 0; i < R2.size(); i++)
      R2[i].noalias() += PiT.Mat2(i) * dT.Mat2(i) * I.R2[i];
    for (size_t i = 0; i < R1.size(); i++)
      R1[i].noalias() += PiT.Mat1(i) * dT.Mat1(i) * I.R1[i];
  }
 
 
  void addMultiplied(const TransmissionMatrix& A,
                     const RadiationVector& x) {
    for (size_t i = 0; i < R4.size(); i++)
      R4[i].noalias() += A.Mat4(i) * x.R4[i];
    for (size_t i = 0; i < R3.size(); i++)
      R3[i].noalias() += A.Mat3(i) * x.R3[i];
    for (size_t i = 0; i < R2.size(); i++)
      R2[i].noalias() += A.Mat2(i) * x.R2[i];
    for (size_t i = 0; i < R1.size(); i++)
      R1[i].noalias() += A.Mat1(i) * x.R1[i];
  }
 
  void setDerivReflection(const RadiationVector& I,
                          const TransmissionMatrix& PiT,
                          const TransmissionMatrix& Z,
                          const TransmissionMatrix& dZ) {
    for (size_t i = 0; i < R4.size(); i++)
      R4[i] = PiT.Mat4(i) * (Z.Mat4(i) * R4[i] + dZ.Mat4(i) * I.R4[i]);
    for (size_t i = 0; i < R3.size(); i++)
      R3[i] = PiT.Mat3(i) * (Z.Mat3(i) * R3[i] + dZ.Mat3(i) * I.R3[i]);
    for (size_t i = 0; i < R2.size(); i++)
      R2[i] = PiT.Mat2(i) * (Z.Mat2(i) * R2[i] + dZ.Mat2(i) * I.R2[i]);
    for (size_t i = 0; i < R1.size(); i++)
      R4[i][0] = PiT.Mat1(i)[0] *
                 (Z.Mat1(i)[0] * R1[i][0] + dZ.Mat1(i)[0] * I.R1[i][0]);
  }
 
  void setBackscatterTransmission(const RadiationVector& I0,
                                  const TransmissionMatrix& Tr,
                                  const TransmissionMatrix& Tf,
                                  const TransmissionMatrix& Z) {
    for (size_t i = 0; i < R4.size(); i++)
      R4[i].noalias() = Tr.Mat4(i) * Z.Mat4(i) * Tf.Mat4(i) * I0.R4[i];
    for (size_t i = 0; i < R3.size(); i++)
      R3[i].noalias() = Tr.Mat3(i) * Z.Mat3(i) * Tf.Mat3(i) * I0.R3[i];
    for (size_t i = 0; i < R2.size(); i++)
      R2[i].noalias() = Tr.Mat2(i) * Z.Mat2(i) * Tf.Mat2(i) * I0.R2[i];
    for (size_t i = 0; i < R1.size(); i++)
      R1[i].noalias() = Tr.Mat1(i) * Z.Mat1(i) * Tf.Mat1(i) * I0.R1[i];
  }
 
  void setBackscatterTransmissionDerivative(const RadiationVector& I0,
                                            const TransmissionMatrix& Tr,
                                            const TransmissionMatrix& Tf,
                                            const TransmissionMatrix& dZ) {
    for (size_t i = 0; i < R4.size(); i++)
      R4[i].noalias() += Tr.Mat4(i) * dZ.Mat4(i) * Tf.Mat4(i) * I0.R4[i];
    for (size_t i = 0; i < R3.size(); i++)
      R3[i].noalias() += Tr.Mat3(i) * dZ.Mat3(i) * Tf.Mat3(i) * I0.R3[i];
    for (size_t i = 0; i < R2.size(); i++)
      R2[i].noalias() += Tr.Mat2(i) * dZ.Mat2(i) * Tf.Mat2(i) * I0.R2[i];
    for (size_t i = 0; i < R1.size(); i++)
      R1[i].noalias() += Tr.Mat1(i) * dZ.Mat1(i) * Tf.Mat1(i) * I0.R1[i];
  }
 
  RadiationVector& operator=(const ConstMatrixView& M) {
    ARTS_ASSERT(M.ncols() == stokes_dim and M.nrows() == Frequencies());
    for (size_t i = 0; i < R4.size(); i++) {
      R4[i][0] = M(i, 0);
      R4[i][1] = M(i, 1);
      R4[i][2] = M(i, 2);
      R4[i][3] = M(i, 3);
    }
    for (size_t i = 0; i < R3.size(); i++) {
      R3[i][0] = M(i, 0);
      R3[i][1] = M(i, 1);
      R3[i][2] = M(i, 2);
    }
    for (size_t i = 0; i < R2.size(); i++) {
      R2[i][0] = M(i, 0);
      R2[i][1] = M(i, 1);
    }
    for (size_t i = 0; i < R1.size(); i++) {
      R1[i][0] = M(i, 0);
    }
    return *this;
  }
 
  const Numeric& operator()(const Index i, const Index j) const {
    switch (stokes_dim) {
      case 4:
        return R4[i][j];
      case 3:
        return R3[i][j];
      case 2:
        return R2[i][j];
      default:
        return R1[i][j];
    }
  }
 
  operator Matrix() const {
    Matrix M(Frequencies(), stokes_dim);
    for (size_t i = 0; i < R4.size(); i++)
      for (size_t j = 0; j < 4; j++) M(i, j) = R4[i](j);
    for (size_t i = 0; i < R3.size(); i++)
      for (size_t j = 0; j < 3; j++) M(i, j) = R3[i](j);
    for (size_t i = 0; i < R2.size(); i++)
      for (size_t j = 0; j < 2; j++) M(i, j) = R2[i](j);
    for (size_t i = 0; i < R1.size(); i++) M(i, 0) = R1[i](0);
    return M;
  }
 
  void setSource(const StokesVector& a,
                 const ConstVectorView& B,
                 const StokesVector& S,
                 Index i) {
    ARTS_ASSERT(a.NumberOfAzimuthAngles() == 1);
    ARTS_ASSERT(a.NumberOfZenithAngles() == 1);
    ARTS_ASSERT(S.NumberOfAzimuthAngles() == 1);
    ARTS_ASSERT(S.NumberOfZenithAngles() == 1);
    switch (stokes_dim) {
      case 4:
        if (not S.IsEmpty())
          R4[i].noalias() =
              Eigen::Vector4d(a.Kjj()[i], a.K12()[i], a.K13()[i], a.K14()[i]) *
                  B[i] +
              Eigen::Vector4d(S.Kjj()[i], S.K12()[i], S.K13()[i], S.K14()[i]);
        else
          R4[i].noalias() =
              Eigen::Vector4d(a.Kjj()[i], a.K12()[i], a.K13()[i], a.K14()[i]) *
              B[i];
        break;
      case 3:
        if (not S.IsEmpty())
          R3[i].noalias() =
              Eigen::Vector3d(a.Kjj()[i], a.K12()[i], a.K13()[i]) * B[i] +
              Eigen::Vector3d(S.Kjj()[i], S.K12()[i], S.K13()[i]);
        else
          R3[i].noalias() =
              Eigen::Vector3d(a.Kjj()[i], a.K12()[i], a.K13()[i]) * B[i];
        break;
      case 2:
        if (not S.IsEmpty())
          R2[i].noalias() = Eigen::Vector2d(a.Kjj()[i], a.K12()[i]) * B[i] +
                            Eigen::Vector2d(S.Kjj()[i], S.K12()[i]);
        else
          R2[i].noalias() = Eigen::Vector2d(a.Kjj()[i], a.K12()[i]) * B[i];
        break;
      default:
        if (not S.IsEmpty())
          R1[i][0] = a.Kjj()[i] * B[i] + S.Kjj()[i];
        else
          R1[i][0] = a.Kjj()[i] * B[i];
    }
  }
 
  [[nodiscard]] Index Frequencies() const {
    switch (stokes_dim) {
      case 4:
        return Index(R4.size());
      case 3:
        return Index(R3.size());
      case 2:
        return Index(R2.size());
      default:
        return Index(R1.size());
    }
  }
 
  friend std::ostream& operator<<(std::ostream& os, const RadiationVector& rv);
 
  friend std::istream& operator>>(std::istream& data, RadiationVector& rv);
};
 
using ArrayOfTransmissionMatrix = Array<TransmissionMatrix>;
using ArrayOfArrayOfTransmissionMatrix = Array<ArrayOfTransmissionMatrix>;
using ArrayOfArrayOfArrayOfTransmissionMatrix = Array<ArrayOfArrayOfTransmissionMatrix>;
using ArrayOfRadiationVector = Array<RadiationVector>;
using ArrayOfArrayOfRadiationVector = Array<ArrayOfRadiationVector>;
using ArrayOfArrayOfArrayOfRadiationVector = Array<ArrayOfArrayOfRadiationVector>;
 
std::ostream& operator<<(std::ostream& os, const TransmissionMatrix& tm);
 
std::ostream& operator<<(std::ostream& os,
                         const ArrayOfTransmissionMatrix& atm);
 
std::ostream& operator<<(std::ostream& os,
                         const ArrayOfArrayOfTransmissionMatrix& aatm);
 
std::ostream& operator<<(std::ostream& os, const RadiationVector& rv);
 
std::ostream& operator<<(std::ostream& os, const ArrayOfRadiationVector& arv);
 
std::ostream& operator<<(std::ostream& os,
                         const ArrayOfArrayOfRadiationVector& aarv);
 
std::istream& operator>>(std::istream& is, TransmissionMatrix& tm);
 
std::istream& operator>>(std::istream& is, RadiationVector& rv);
 
enum class BackscatterSolver {
  CommutativeTransmission,
  FullTransmission,
};
 
enum class CumulativeTransmission {
  Forward,
  Reverse,
};
 
enum class RadiativeTransferSolver {
  Emission,
  Transmission,
  LinearWeightedEmission,
};
 
void update_radiation_vector(RadiationVector& I,
                             ArrayOfRadiationVector& dI1,
                             ArrayOfRadiationVector& dI2,
                             const RadiationVector& J1,
                             const RadiationVector& J2,
                             const ArrayOfRadiationVector& dJ1,
                             const ArrayOfRadiationVector& dJ2,
                             const TransmissionMatrix& T,
                             const TransmissionMatrix& PiT,
                             const ArrayOfTransmissionMatrix& dT1,
                             const ArrayOfTransmissionMatrix& dT2,
                             const PropagationMatrix& K1,
                             const PropagationMatrix& K2,
                             const ArrayOfPropagationMatrix& dK1,
                             const ArrayOfPropagationMatrix& dK2,
                             const Numeric r,
                             const Vector& dr1,
                             const Vector& dr2,
                             const Index ia,
                             const Index iz,
                             const RadiativeTransferSolver solver);
 
void stepwise_source(RadiationVector& J,
                     ArrayOfRadiationVector& dJ,
                     const PropagationMatrix& K,
                     const StokesVector& a,
                     const StokesVector& S,
                     const ArrayOfPropagationMatrix& dK,
                     const ArrayOfStokesVector& da,
                     const ArrayOfStokesVector& dS,
                     const ConstVectorView& B,
                     const ConstVectorView& dB_dT,
                     const ArrayOfRetrievalQuantity& jacobian_quantities,
                     const bool& jacobian_do);
 
void stepwise_transmission(TransmissionMatrix& T,
                           ArrayOfTransmissionMatrix& dT1,
                           ArrayOfTransmissionMatrix& dT2,
                           const PropagationMatrix& K1,
                           const PropagationMatrix& K2,
                           const ArrayOfPropagationMatrix& dK1,
                           const ArrayOfPropagationMatrix& dK2,
                           const Numeric& r,
                           const Numeric& dr_dtemp1,
                           const Numeric& dr_dtemp2,
                           const Index temp_deriv_pos);
 
ArrayOfTransmissionMatrix cumulative_transmission(
    const ArrayOfTransmissionMatrix& T,
    const CumulativeTransmission type) /*[[expects: T.nelem()>0]]*/;
 
void set_backscatter_radiation_vector(
    ArrayOfRadiationVector& I,
    ArrayOfArrayOfArrayOfRadiationVector& dI,
    const RadiationVector & I_incoming,
    const ArrayOfTransmissionMatrix& T,
    const ArrayOfTransmissionMatrix& PiTf,
    const ArrayOfTransmissionMatrix& PiTr,
    const ArrayOfTransmissionMatrix& Z,
    const ArrayOfArrayOfTransmissionMatrix& dT1,
    const ArrayOfArrayOfTransmissionMatrix& dT2,
    const ArrayOfArrayOfTransmissionMatrix& dZ,
    const BackscatterSolver solver);
 
ArrayOfTransmissionMatrix bulk_backscatter(const ConstTensor5View& Pe,
                                           const ConstMatrixView& pnd);
 
ArrayOfArrayOfTransmissionMatrix bulk_backscatter_derivative(
    const ConstTensor5View& Pe, const ArrayOfMatrix& dpnd_dx);
 
#endif  // transmissionmatrix_h