11#ifndef transmissionmatrix_h
12#define transmissionmatrix_h
14#pragma GCC diagnostic push
15#pragma GCC diagnostic ignored "-Wshadow"
16#pragma GCC diagnostic ignored "-Wconversion"
17#pragma GCC diagnostic ignored "-Wfloat-conversion"
19#pragma GCC diagnostic ignored "-Wdeprecated-copy-with-dtor"
20#pragma GCC diagnostic ignored "-Wdeprecated-copy-with-user-provided-dtor"
22#pragma GCC diagnostic ignored "-Wclass-memaccess"
27#pragma GCC diagnostic pop
30#include "propagationmatrix.h"
33Eigen::Matrix<Numeric, N, N>
prop_matrix(
const ConstVectorView& vec) {
41 static_assert (N>0 and N<5,
"Bad size N");
42 if constexpr (N == 1) {
43 return Eigen::Matrix<Numeric, 1, 1>(
a);
44 }
else if constexpr (N==2) {
45 return (Eigen::Matrix2d() <<
a,
b,
b,
a).finished();
46 }
else if constexpr (N==3) {
47 return (Eigen::Matrix3d() <<
a,
b,
c,
b,
a,
u,
c, -
u,
a).finished();
48 }
else if constexpr (N==4) {
49 return (Eigen::Matrix4d() <<
a,
b,
c,
d,
b,
a,
u,
v,
c, -
u,
a,
w,
d, -
v, -
w,
a).finished();
61Eigen::Matrix<Numeric, N, N>
prop_matrix(
const ConstMatrixView& m) {
62 static_assert (N>0 and N<5,
"Bad size N");
63 if constexpr (N == 1) {
64 return Eigen::Matrix<Numeric, 1, 1>(m(0, 0));
65 }
else if constexpr (N == 2) {
66 return (Eigen::Matrix2d() << m(0, 0), m(0, 1),
67 m(1, 0), m(1, 1)).finished();
68 }
else if constexpr (N == 3) {
69 return (Eigen::Matrix3d() << m(0, 0), m(0, 1), m(0, 2),
70 m(1, 0), m(1, 1), m(1, 2),
71 m(2, 0), m(2, 1), m(2, 2)) .finished();
72 }
else if constexpr (N == 4) {
73 return (Eigen::Matrix4d() << m(0, 0), m(0, 1), m(0, 2), m(0, 3),
74 m(1, 0), m(1, 1), m(1, 2), m(1, 3),
75 m(2, 0), m(2, 1), m(2, 2), m(2, 3),
76 m(3, 0), m(3, 1), m(3, 2), m(3, 3)) .finished();
89 static_assert (N>0 and N<5,
"Bad size N");
90 if constexpr (N == 1) {
91 return Eigen::Matrix<double, 1, 1>(1 /
a);
92 }
else if constexpr (N==2) {
93 return (Eigen::Matrix2d() <<
a, -
b, -
b,
a).finished() / (
a *
a -
b *
b);
94 }
else if constexpr (N==3) {
95 return (Eigen::Matrix3d() <<
a *
a +
u *
u,
106 }
else if constexpr (N==4) {
107 return (Eigen::Matrix4d() <<
a *
a *
a +
a *
u *
u +
a *
v *
v +
a *
w *
w,
150 std::vector<Eigen::Matrix4d>
T4;
151 std::vector<Eigen::Matrix3d>
T3;
152 std::vector<Eigen::Matrix2d>
T2;
153 std::vector<Eigen::Matrix<double, 1, 1>>
T1;
194 const Numeric& r = 1.0);
196 operator Tensor3()
const;
205 [[nodiscard]]
const Eigen::Matrix4d&
Mat4(
size_t i)
const;
212 [[nodiscard]]
const Eigen::Matrix3d&
Mat3(
size_t i)
const;
219 [[nodiscard]]
const Eigen::Matrix2d&
Mat2(
size_t i)
const;
226 [[nodiscard]]
const Eigen::Matrix<double, 1, 1>&
Mat1(
size_t i)
const;
233 [[nodiscard]] Eigen::MatrixXd
Mat(
size_t i)
const;
240 Eigen::Matrix4d&
Mat4(
size_t i);
247 Eigen::Matrix3d&
Mat3(
size_t i);
254 Eigen::Matrix2d&
Mat2(
size_t i);
261 Eigen::Matrix<double, 1, 1>&
Mat1(
size_t i);
294 [[nodiscard]] Numeric
operator()(
const Index i,
296 const Index k)
const;
305 Numeric&
operator()(
const Index i,
const Index j,
const Index k);
332 friend std::ostream&
operator<<(std::ostream& os,
341 static_assert(N > 0 and N < 5,
"Bad size N");
342 if constexpr (N == 1)
344 else if constexpr (N == 2)
346 else if constexpr (N == 3)
348 else if constexpr (N == 4)
354 [[nodiscard]]
auto&
TraMat(
size_t i)
const noexcept {
355 static_assert(N > 0 and N < 5,
"Bad size N");
356 if constexpr (N == 1)
358 else if constexpr (N == 2)
360 else if constexpr (N == 3)
362 else if constexpr (N == 4)
368 [[nodiscard]] Eigen::Matrix<Numeric, N, N>
OptDepth(
size_t i)
const noexcept {
369 return (-TraMat<N>(i).diagonal().array().log().matrix()).asDiagonal();
393 const Eigen::Matrix<Numeric, N, N> T,
394 const Eigen::Matrix<Numeric, N, 1> far,
395 const Eigen::Matrix<Numeric, N, 1> close,
396 const Eigen::Matrix<Numeric, N, N> Kfar,
397 const Eigen::Matrix<Numeric, N, N> Kclose,
398 const Numeric r)
const noexcept {
399 static_assert(N > 0 and N < 5,
"Bad size N");
401 const auto I = Eigen::Matrix<Numeric, N, N>::Identity();
402 if (T(0, 0) < 0.99) {
403 const auto od = 0.5 * r * (Kfar + Kclose);
404 Eigen::Matrix<Numeric, N, 1> second =
405 od.inverse() * ((I - (I + od) * T) * far + (od - I + T) * close);
406 if ((far[0] < close[0] and Kfar(0, 0) > Kclose(0, 0)) or
407 (far[0] > close[0] and Kfar(0, 0) < Kclose(0, 0))) {
408 Eigen::Matrix<Numeric, N, 1> second_limit =
409 0.5 * (I - T) * (far + close);
414 if (second_limit[0] > second[0])
return second;
421 return 0.5 * (I - T) * (far + close);
428 const Eigen::Matrix<Numeric, N, 1> far,
429 const Eigen::Matrix<Numeric, N, 1> close,
430 const Eigen::Matrix<Numeric, N, 1>
d,
431 bool isfar)
const noexcept {
432 static_assert(N > 0 and N < 5,
"Bad size N");
434 const auto I = Eigen::Matrix<Numeric, N, N>::Identity();
435 const auto T = TraMat<N>(i);
436 const auto& dTdx = dx.TraMat<N>(i);
437 const Eigen::Matrix<Numeric, N, 1> first = 0.5 * (I - T) * (far + close);
438 if (T(0, 0) < 0.99) {
439 const auto od = OptDepth<N>(i);
440 const auto doddx = dx.OptDepth<N>(i);
441 const Eigen::Matrix<Numeric, N, 1> second =
442 od.inverse() * ((I - (I + od) * T) * far + (od - I + T) * close);
443 if (first[0] > second[0]) {
445 return od.inverse() *
446 ((I - (I + od) * T) *
d - (I + od) * dTdx * far +
447 (doddx + T) * close - doddx * second);
450 return od.inverse() * (-(I + od) * dTdx * far + (od - I + T) *
d +
451 (doddx + T) * close - doddx * second);
454 return 0.5 * ((I - T) *
d - dTdx * (far + close));
457 return 0.5 * ((I - T) *
d - dTdx * (far + close));
467[[nodiscard]] LazyScale<TransmissionMatrix>
operator*(
476[[nodiscard]] LazyScale<TransmissionMatrix>
operator*(
482 std::vector<Eigen::Vector4d>
R4;
483 std::vector<Eigen::Vector3d>
R3;
484 std::vector<Eigen::Vector2d>
R2;
485 std::vector<Eigen::Matrix<double, 1, 1>>
R1;
546 [[nodiscard]]
const Eigen::Vector4d&
Vec4(
size_t i)
const;
553 [[nodiscard]]
const Eigen::Vector3d&
Vec3(
size_t i)
const;
560 [[nodiscard]]
const Eigen::Vector2d&
Vec2(
size_t i)
const;
567 [[nodiscard]]
const Eigen::Matrix<double, 1, 1>&
Vec1(
size_t i)
const;
574 [[nodiscard]] Eigen::VectorXd
Vec(
size_t i)
const;
581 Eigen::Vector4d&
Vec4(
size_t i);
588 Eigen::Vector3d&
Vec3(
size_t i);
595 Eigen::Vector2d&
Vec2(
size_t i);
602 Eigen::Matrix<double, 1, 1>&
Vec1(
size_t i);
627 const ConstMatrixView& Kfar,
628 const ConstMatrixView& Kclose,
732 const Numeric&
operator()(
const Index i,
const Index j)
const;
740 [[nodiscard]] Numeric&
operator()(
const Index i,
const Index j);
745 operator Matrix()
const;
755 const ConstVectorView& B,
756 const StokesVector& S,
823 const PropagationMatrix& K1,
824 const PropagationMatrix& K2,
825 const ArrayOfPropagationMatrix& dK1,
826 const ArrayOfPropagationMatrix& dK2,
853 const PropagationMatrix& K,
854 const StokesVector&
a,
855 const StokesVector& S,
856 const ArrayOfPropagationMatrix& dK,
857 const ArrayOfStokesVector& da,
858 const ArrayOfStokesVector& dS,
859 const ConstVectorView& B,
860 const ConstVectorView& dB_dT,
862 const bool& jacobian_do);
881 const PropagationMatrix& K1,
882 const PropagationMatrix& K2,
883 const ArrayOfPropagationMatrix& dK1,
884 const ArrayOfPropagationMatrix& dK2,
886 const Numeric& dr_dtemp1,
887 const Numeric& dr_dtemp2,
888 const Index temp_deriv_pos);
939 const ConstMatrixView& pnd);
951 const ConstTensor5View& Pe,
const ArrayOfMatrix& dpnd_dx);
This can be used to make arrays out of anything.
Routines for setting up the jacobian.
Radiation Vector for Stokes dimension 1-4.
std::vector< Eigen::Vector2d > R2
Index Frequencies() const
Get frequency count.
RadiationVector(RadiationVector &&rv) noexcept=default
Construct a new Radiation Vector object.
const Eigen::Matrix< double, 1, 1 > & Vec1(size_t i) const
Return Vector at position.
void setSource(const StokesVector &a, const ConstVectorView &B, const StokesVector &S, Index i)
Set this to source vector at position.
void addMultiplied(const TransmissionMatrix &A, const RadiationVector &x)
Add multiply.
const Eigen::Vector4d & Vec4(size_t i) const
Return Vector at position.
void leftMul(const TransmissionMatrix &T)
Multiply radiation vector from the left.
const Numeric & operator()(const Index i, const Index j) const
Access operator.
void rem_avg(const RadiationVector &O1, const RadiationVector &O2)
Remove the average of two other RadiationVector from *this.
std::vector< Eigen::Matrix< double, 1, 1 > > R1
void add_avg(const RadiationVector &O1, const RadiationVector &O2)
Add the average of two other RadiationVector to *this.
std::vector< Eigen::Vector4d > R4
RadiationVector & operator+=(const RadiationVector &rv)
Addition operator.
Eigen::VectorXd Vec(size_t i) const
Return Vector at position by copy.
void addDerivTransmission(const TransmissionMatrix &PiT, const TransmissionMatrix &dT, const RadiationVector &I)
Add the transmission derivative to this.
friend std::ostream & operator<<(std::ostream &os, const RadiationVector &rv)
Output operator.
void setDerivReflection(const RadiationVector &I, const TransmissionMatrix &PiT, const TransmissionMatrix &Z, const TransmissionMatrix &dZ)
Sets *this to the reflection derivative.
const Eigen::Vector3d & Vec3(size_t i) const
Return Vector at position.
const Eigen::Vector2d & Vec2(size_t i) const
Return Vector at position.
std::vector< Eigen::Vector3d > R3
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)
Add the emission derivative to this.
void setBackscatterTransmission(const RadiationVector &I0, const TransmissionMatrix &Tr, const TransmissionMatrix &Tf, const TransmissionMatrix &Z)
Set this to backscatter transmission.
void addDerivEmission(const TransmissionMatrix &PiT, const TransmissionMatrix &dT, const TransmissionMatrix &T, const RadiationVector &ImJ, const RadiationVector &dJ)
Add the emission derivative to this.
RadiationVector & operator=(const RadiationVector &rv)=default
Assign old radiation vector to this.
void add_weighted(const TransmissionMatrix &T, const RadiationVector &far, const RadiationVector &close, const ConstMatrixView &Kfar, const ConstMatrixView &Kclose, const Numeric r)
Add the weighted source of two RadiationVector to *this.
friend std::istream & operator>>(std::istream &data, RadiationVector &rv)
Input operator.
RadiationVector & operator=(RadiationVector &&rv) noexcept=default
Assign old radiation vector to this.
void setBackscatterTransmissionDerivative(const RadiationVector &I0, const TransmissionMatrix &Tr, const TransmissionMatrix &Tf, const TransmissionMatrix &dZ)
Set this to backscatter transmission scatter derivative.
RadiationVector(const RadiationVector &rv)=default
Construct a new Radiation Vector object.
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
TransmissionMatrix & operator=(TransmissionMatrix &&tm) noexcept=default
Move operator.
void mul_aliased(const TransmissionMatrix &A, const TransmissionMatrix &B)
Set this to a multiple of A by B.
std::vector< Eigen::Matrix< double, 1, 1 > > T1
std::vector< Eigen::Matrix4d > T4
Eigen::Matrix< Numeric, N, N > OptDepth(size_t i) const noexcept
Simple template access for the optical depth.
TransmissionMatrix & operator*=(const Numeric &scale)
Scale self.
const Eigen::Matrix2d & Mat2(size_t i) const
Get Matrix at position.
TransmissionMatrix(const TransmissionMatrix &tm)=default
Construct a new Transmission Matrix object.
Eigen::MatrixXd Mat(size_t i) const
Get Matrix at position by copy.
auto & TraMat(size_t i) noexcept
Simple template access for the transmission.
Index Frequencies() const
Number of frequencies.
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
auto & TraMat(size_t i) const noexcept
Simple template access for the transmission.
void mul(const TransmissionMatrix &A, const TransmissionMatrix &B)
Set this to a multiple of A by B.
std::vector< Eigen::Matrix3d > T3
const Eigen::Matrix4d & Mat4(size_t i) const
Get Matrix at position.
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
void setIdentity()
Set to identity matrix.
const Eigen::Matrix< double, 1, 1 > & Mat1(size_t i) const
Get Matrix at position.
std::vector< Eigen::Matrix2d > T2
friend std::istream & operator>>(std::istream &data, TransmissionMatrix &tm)
Input operator.
const Eigen::Matrix3d & Mat3(size_t i) const
Get Matrix at position.
TransmissionMatrix(TransmissionMatrix &&tm) noexcept=default
Construct a new Transmission Matrix object.
friend std::ostream & operator<<(std::ostream &os, const TransmissionMatrix &tm)
Output operator.
Numeric operator()(const Index i, const Index j, const Index k) const
Access value in matrix.
TransmissionMatrix & operator=(const TransmissionMatrix &tm)=default
Assignment operator.
TransmissionMatrix & operator+=(const LazyScale< TransmissionMatrix > &lstm)
Assign to *this lazily.
void setZero()
Set to zero matrix.
RadiativeTransferSolver
Intended to hold various forward solvers.
ArrayOfTransmissionMatrix bulk_backscatter(const ConstTensor5View &Pe, const ConstMatrixView &pnd)
Bulk back-scattering.
BackscatterSolver
Intended to hold various backscatter solvers.
@ CommutativeTransmission
CumulativeTransmission
Intended to hold various ways to accumulate the transmission matrix.
ArrayOfArrayOfTransmissionMatrix bulk_backscatter_derivative(const ConstTensor5View &Pe, const ArrayOfMatrix &dpnd_dx)
Derivatives of bulk back-scattering
LazyScale< TransmissionMatrix > operator*(const TransmissionMatrix &tm, const Numeric &x)
Lazy scale of Transmission Matrix.
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
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)
Set the backscatter radiation vector.
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)
Set the stepwise source.
Eigen::Matrix< Numeric, N, N > inv_prop_matrix(const ConstVectorView &vec)
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)
Update the Radiation Vector.
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)
Set the stepwise transmission matrix.
Eigen::Matrix< Numeric, N, N > prop_matrix(const ConstVectorView &vec)