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"
#if defined(__clang__)
#pragma GCC diagnostic ignored "-Wdeprecated-copy-with-dtor"
#pragma GCC diagnostic ignored "-Wdeprecated-copy-with-user-provided-dtor"
#else
#pragma GCC diagnostic ignored "-Wclass-memaccess"
#endif
 
#include <Eigen/Dense>
 
#pragma GCC diagnostic pop
 
#include "jacobian.h"
#include "propagationmatrix.h"
 
template <int N>
Eigen::Matrix<Numeric, N, N> prop_matrix(const ConstVectorView& vec) {
#define a vec[0]
#define b vec[1]
#define c vec[2]
#define d vec[3]
#define u vec[N]
#define v vec[5]
#define w vec[6]
  static_assert (N>0 and N<5, "Bad size N");
  if constexpr (N == 1) {
    return Eigen::Matrix<Numeric, 1, 1>(a);
  } else if constexpr (N==2) {
    return (Eigen::Matrix2d() << a, b, b, a).finished();
  } else if constexpr (N==3) {
    return (Eigen::Matrix3d() << a, b, c, b, a, u, c, -u, a).finished();
  } else if constexpr (N==4) {
    return (Eigen::Matrix4d() << a, b, c, d, b, a, u, v, c, -u, a, w, d, -v, -w, a).finished();
  }
#undef a
#undef b
#undef c
#undef d
#undef u
#undef v
#undef w
}
 
template <int N>
Eigen::Matrix<Numeric, N, N> prop_matrix(const ConstMatrixView& m) {
  static_assert (N>0 and N<5, "Bad size N");
  if constexpr (N == 1) {
    return Eigen::Matrix<Numeric, 1, 1>(m(0, 0));
  } else if constexpr (N == 2) {
    return (Eigen::Matrix2d() << m(0, 0), m(0, 1),
                                 m(1, 0), m(1, 1)).finished();
  } else if constexpr (N == 3) {
    return (Eigen::Matrix3d() << m(0, 0), m(0, 1), m(0, 2),
                                 m(1, 0), m(1, 1), m(1, 2),
                                 m(2, 0), m(2, 1), m(2, 2)) .finished();
  } else if constexpr (N == 4) {
    return (Eigen::Matrix4d() << m(0, 0), m(0, 1), m(0, 2), m(0, 3),
                                 m(1, 0), m(1, 1), m(1, 2), m(1, 3),
                                 m(2, 0), m(2, 1), m(2, 2), m(2, 3),
                                 m(3, 0), m(3, 1), m(3, 2), m(3, 3)) .finished();
  }
}
 
template <int N>
Eigen::Matrix<Numeric, N, N> inv_prop_matrix(const ConstVectorView& vec) {
#define a vec[0]
#define b vec[1]
#define c vec[2]
#define d vec[3]
#define u vec[N]
#define v vec[5]
#define w vec[6]
  static_assert (N>0 and N<5, "Bad size N");
  if constexpr (N == 1) {
    return Eigen::Matrix<double, 1, 1>(1 / a);
  } else if constexpr (N==2) {
    return (Eigen::Matrix2d() << a, -b, -b, a).finished() / (a * a - b * b);
  } else if constexpr (N==3) {
    return (Eigen::Matrix3d() << a * a + u * u,
          -a * b - c * u,
          -a * c + b * u,
          -a * b + c * u,
          a * a - c * c,
          -a * u + b * c,
          -a * c - b * u,
          a * u + b * c,
          a * a - b * b)
             .finished() /
         (a * a * a - a * b * b - a * c * c + a * u * u);
  } else if constexpr (N==4) {
    return (Eigen::Matrix4d() << a * a * a + a * u * u + a * v * v + a * w * w,
          -a * a * b - a * c * u - a * d * v - b * w * w + c * v * w -
              d * u * w,
          -a * a * c + a * b * u - a * d * w + b * v * w - c * v * v +
              d * u * v,
          -a * a * d + a * b * v + a * c * w - b * u * w + c * u * v -
              d * u * u,
          -a * a * b + a * c * u + a * d * v - b * w * w + c * v * w -
              d * u * w,
          a * a * a - a * c * c - a * d * d + a * w * w,
          -a * a * u + a * b * c - a * v * w + b * d * w - c * d * v +
              d * d * u,
          -a * a * v + a * b * d + a * u * w - b * c * w + c * c * v -
              c * d * u,
          -a * a * c - a * b * u + a * d * w + b * v * w - c * v * v +
              d * u * v,
          a * a * u + a * b * c - a * v * w - b * d * w + c * d * v - d * d * u,
          a * a * a - a * b * b - a * d * d + a * v * v,
          -a * a * w + a * c * d - a * u * v + b * b * w - b * c * v +
              b * d * u,
          -a * a * d - a * b * v - a * c * w - b * u * w + c * u * v -
              d * u * u,
          a * a * v + a * b * d + a * u * w + b * c * w - c * c * v + c * d * u,
          a * a * w + a * c * d - a * u * v - b * b * w + b * c * v - b * d * u,
          a * a * a - a * b * b - a * c * c + a * u * u)
             .finished() /
         (a * a * a * a - a * a * b * b - a * a * c * c - a * a * d * d +
          a * a * u * u + a * a * v * v + a * a * w * w - b * b * w * w +
          2 * b * c * v * w - 2 * b * d * u * w - c * c * v * v +
          2 * c * d * u * v - d * d * u * u);
  }
#undef a
#undef b
#undef c
#undef d
#undef u
#undef v
#undef w
}
 
struct TransmissionMatrix {
  Index stokes_dim;
  std::vector<Eigen::Matrix4d> T4;
  std::vector<Eigen::Matrix3d> T3;
  std::vector<Eigen::Matrix2d> T2;
  std::vector<Eigen::Matrix<double, 1, 1>> T1;
 
  TransmissionMatrix(Index nf = 0, Index stokes = 1);
 
  TransmissionMatrix(TransmissionMatrix&& tm) noexcept = default;
 
  TransmissionMatrix(const TransmissionMatrix& tm) = default;
 
  TransmissionMatrix& operator=(const TransmissionMatrix& tm) = default;
 
  TransmissionMatrix& operator=(TransmissionMatrix&& tm) noexcept = default;
 
  explicit TransmissionMatrix(const PropagationMatrix& pm,
                              const Numeric& r = 1.0);
 
  operator Tensor3() const;
 
  explicit TransmissionMatrix(const ConstMatrixView& mat);
 
  [[nodiscard]] const Eigen::Matrix4d& Mat4(size_t i) const;
 
  [[nodiscard]] const Eigen::Matrix3d& Mat3(size_t i) const;
 
  [[nodiscard]] const Eigen::Matrix2d& Mat2(size_t i) const;
 
  [[nodiscard]] const Eigen::Matrix<double, 1, 1>& Mat1(size_t i) const;
 
  [[nodiscard]] Eigen::MatrixXd Mat(size_t i) const;
 
  Eigen::Matrix4d& Mat4(size_t i);
 
  Eigen::Matrix3d& Mat3(size_t i);
 
  Eigen::Matrix2d& Mat2(size_t i);
 
  Eigen::Matrix<double, 1, 1>& Mat1(size_t i);
 
  void setIdentity();
 
  void setZero();
 
  void mul(const TransmissionMatrix& A, const TransmissionMatrix& B);
 
  void mul_aliased(const TransmissionMatrix& A, const TransmissionMatrix& B);
 
  [[nodiscard]] Numeric operator()(const Index i,
                                   const Index j,
                                   const Index k) const;
 
  Numeric& operator()(const Index i, const Index j, const Index k);
 
  [[nodiscard]] Index Frequencies() const;
 
  TransmissionMatrix& operator+=(const LazyScale<TransmissionMatrix>& lstm);
 
  TransmissionMatrix& operator*=(const Numeric& scale);
 
  TransmissionMatrix& operator=(const LazyScale<TransmissionMatrix>& lstm);
 
  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));
  }
};
 
[[nodiscard]] LazyScale<TransmissionMatrix> operator*(
    const TransmissionMatrix& tm, const Numeric& x);
 
[[nodiscard]] LazyScale<TransmissionMatrix> operator*(
    const Numeric& x, const TransmissionMatrix& tm);
 
struct RadiationVector {
  Index stokes_dim;
  std::vector<Eigen::Vector4d> R4;
  std::vector<Eigen::Vector3d> R3;
  std::vector<Eigen::Vector2d> R2;
  std::vector<Eigen::Matrix<double, 1, 1>> R1;
 
  RadiationVector(Index nf = 0, Index stokes = 1);
 
  RadiationVector(RadiationVector&& rv) noexcept = default;
 
  RadiationVector(const RadiationVector& rv) = default;
 
  RadiationVector& operator=(const RadiationVector& rv) = default;
 
  RadiationVector& operator=(RadiationVector&& rv) noexcept = default;
 
  RadiationVector& operator+=(const RadiationVector& rv);
 
  void leftMul(const TransmissionMatrix& T);
 
  void SetZero(size_t i);
 
  void SetZero();
 
  [[nodiscard]] const Eigen::Vector4d& Vec4(size_t i) const;
 
  [[nodiscard]] const Eigen::Vector3d& Vec3(size_t i) const;
 
  [[nodiscard]] const Eigen::Vector2d& Vec2(size_t i) const;
 
  [[nodiscard]] const Eigen::Matrix<double, 1, 1>& Vec1(size_t i) const;
 
  [[nodiscard]] Eigen::VectorXd Vec(size_t i) const;
 
  Eigen::Vector4d& Vec4(size_t i);
 
  Eigen::Vector3d& Vec3(size_t i);
 
  Eigen::Vector2d& Vec2(size_t i);
 
  Eigen::Matrix<double, 1, 1>& Vec1(size_t i);
 
  void rem_avg(const RadiationVector& O1, const RadiationVector& O2);
 
  void add_avg(const RadiationVector& O1, const RadiationVector& O2);
 
  void add_weighted(const TransmissionMatrix& T,
                    const RadiationVector& far,
                    const RadiationVector& close,
                    const ConstMatrixView& Kfar,
                    const ConstMatrixView& Kclose,
                    const Numeric r);
 
  void addDerivEmission(const TransmissionMatrix& PiT,
                        const TransmissionMatrix& dT,
                        const TransmissionMatrix& T,
                        const RadiationVector& ImJ,
                        const RadiationVector& dJ);
 
  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);
 
  void addDerivTransmission(const TransmissionMatrix& PiT,
                            const TransmissionMatrix& dT,
                            const RadiationVector& I);
 
  void addMultiplied(const TransmissionMatrix& A, const RadiationVector& x);
 
  void setDerivReflection(const RadiationVector& I,
                          const TransmissionMatrix& PiT,
                          const TransmissionMatrix& Z,
                          const TransmissionMatrix& dZ);
 
  void setBackscatterTransmission(const RadiationVector& I0,
                                  const TransmissionMatrix& Tr,
                                  const TransmissionMatrix& Tf,
                                  const TransmissionMatrix& Z);
 
  void setBackscatterTransmissionDerivative(const RadiationVector& I0,
                                            const TransmissionMatrix& Tr,
                                            const TransmissionMatrix& Tf,
                                            const TransmissionMatrix& dZ);
 
  RadiationVector& operator=(const ConstMatrixView& M);
 
  const Numeric& operator()(const Index i, const Index j) const;
 
  [[nodiscard]] Numeric& operator()(const Index i, const Index j);
  operator Matrix() const;
 
  void setSource(const StokesVector& a,
                 const ConstVectorView& B,
                 const StokesVector& S,
                 Index i);
 
  [[nodiscard]] Index Frequencies() const;
 
  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>;
 
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,
                     RadiationVector& J_add,
                     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