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

#include <Eigen/Dense>
#include <Eigen/StdVector>
#include "jacobian.h"
#include "propagationmatrix.h"

class TransmissionMatrix {
 private:
  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;

 public:

  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()) {
    assert(stokes_dim < 5 and stokes_dim > 0);
  }

  TransmissionMatrix(TransmissionMatrix&& tm) noexcept
      : stokes_dim(std::move(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 = std::move(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;
  }

  const Eigen::Matrix4d& Mat4(size_t i) const { return T4[i]; }

  const Eigen::Matrix3d& Mat3(size_t i) const { return T3[i]; }

  const Eigen::Matrix2d& Mat2(size_t i) const { return T2[i]; }

  const Eigen::Matrix<double, 1, 1>& Mat1(size_t i) const { return T1[i]; }

  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);
    }
    std::terminate();
  }

  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);
    }
  }

  Index StokesDim() const { return stokes_dim; }

  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 <size_t N>
  struct Weights {Eigen::Matrix<Numeric, int(N), int(N)> far, close;};
  
  template <size_t N>
  Weights<N> linear_in_tau_weights(size_t) const noexcept {
    // FIXME: Make the complete implementation here by using "if constexpr" in C++17
    static_assert (N not_eq 0 and N < 5, "Bad size N");
  }
};

template <> inline TransmissionMatrix::Weights<1> TransmissionMatrix::linear_in_tau_weights(size_t i) const noexcept {
  const Numeric od = std::log(this -> operator()(i, 0, 0));
  return TransmissionMatrix::Weights<1>{(Eigen::Matrix<Numeric, 1, 1>::Identity() - (1 + od) * T1[i])/od,
    (od * Eigen::Matrix<Numeric, 1, 1>::Identity() - Eigen::Matrix<Numeric, 1, 1>::Identity() + T1[i])/od};
}

template <> inline TransmissionMatrix::Weights<2> TransmissionMatrix::linear_in_tau_weights(size_t i) const noexcept {
  const Numeric od = std::log(this -> operator()(i, 0, 0));
  return TransmissionMatrix::Weights<2>{(Eigen::Matrix<Numeric, 2, 2>::Identity() - (1 + od) * T2[i])/od,
    (od * Eigen::Matrix<Numeric, 2, 2>::Identity() - Eigen::Matrix<Numeric, 2, 2>::Identity() + T2[i])/od};
}

template <> inline TransmissionMatrix::Weights<3> TransmissionMatrix::linear_in_tau_weights(size_t i) const noexcept {
  const Numeric od = std::log(this -> operator()(i, 0, 0));
  return TransmissionMatrix::Weights<3>{(Eigen::Matrix<Numeric, 3, 3>::Identity() - (1 + od) * T3[i])/od,
    (od * Eigen::Matrix<Numeric, 3, 3>::Identity() - Eigen::Matrix<Numeric, 3, 3>::Identity() + T3[i])/od};
}

template <> inline TransmissionMatrix::Weights<4> TransmissionMatrix::linear_in_tau_weights(size_t i) const noexcept {
  const Numeric od = std::log(this -> operator()(i, 0, 0));
  return TransmissionMatrix::Weights<4>{(Eigen::Matrix<Numeric, 4, 4>::Identity() - (1 + od) * T4[i])/od,
    (od * Eigen::Matrix<Numeric, 4, 4>::Identity() - Eigen::Matrix<Numeric, 4, 4>::Identity() + T4[i])/od};
}

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);
}

class RadiationVector {
 private:
  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;

 public:

  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()) {
    assert(stokes_dim < 5 and stokes_dim > 0);
  }

  RadiationVector(RadiationVector&& rv) noexcept
      : stokes_dim(std::move(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 = std::move(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;
    }
  }

  const Eigen::Vector4d& Vec4(size_t i) const { return R4[i]; }

  const Eigen::Vector3d& Vec3(size_t i) const { return R3[i]; }

  const Eigen::Vector2d& Vec2(size_t i) const { return R2[i]; }

  const Eigen::Matrix<double, 1, 1>& Vec1(size_t i) const { return R1[i]; }
  
  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& close, const RadiationVector& far) {
    for (size_t i = 0; i < R4.size(); i++) {
      const auto w = T.linear_in_tau_weights<4>(i);
      R4[i].noalias() += w.far * far.R4[i] + w.close * close.R4[i];
    }
    for (size_t i = 0; i < R3.size(); i++) {
      const auto w = T.linear_in_tau_weights<3>(i);
      R3[i].noalias() += w.far * far.R3[i] + w.close * close.R3[i];
    }
    for (size_t i = 0; i < R2.size(); i++) {
      const auto w = T.linear_in_tau_weights<2>(i);
      R2[i].noalias() += w.far * far.R2[i] + w.close * close.R2[i];
    }
    for (size_t i = 0; i < R1.size(); i++) {
      const auto w = T.linear_in_tau_weights<1>(i);
      R1[i].noalias() += w.far * far.R1[i] + w.close * close.R1[i];
    }
  }

  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 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) {
    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) {
    assert(a.NumberOfAzimuthAngles() == 1);
    assert(a.NumberOfZenithAngles() == 1);
    assert(S.NumberOfAzimuthAngles() == 1);
    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];
    }
  }

  Index StokesDim() const { return stokes_dim; }

  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);
};

typedef Array<TransmissionMatrix> ArrayOfTransmissionMatrix;
typedef Array<ArrayOfTransmissionMatrix> ArrayOfArrayOfTransmissionMatrix;
typedef Array<ArrayOfArrayOfTransmissionMatrix> ArrayOfArrayOfArrayOfTransmissionMatrix;
typedef Array<RadiationVector> ArrayOfRadiationVector;
typedef Array<ArrayOfRadiationVector> ArrayOfArrayOfRadiationVector;
typedef Array<ArrayOfArrayOfRadiationVector> ArrayOfArrayOfArrayOfRadiationVector;

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 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 cumulative_backscatter(ConstTensor5View t,
                                                 ConstMatrixView m);

ArrayOfArrayOfTransmissionMatrix cumulative_backscatter_derivative(
    ConstTensor5View t, const ArrayOfMatrix& aom);

#endif  // transmissionmatrix_h