ARTS 2.5.9 (git: 825fa5f2)
transmissionmatrix.h
Go to the documentation of this file.
1/* Copyright (C) 2018
2 * Richard Larsson <ric.larsson@gmail.com>
3 *
4 * This program is free software; you can redistribute it and/or modify it
5 * under the terms of the GNU General Public License as published by the
6 * Free Software Foundation; either version 2, or (at your option) any
7 * later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program; if not, write to the Free Software
16 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17 * USA. */
18
29#ifndef transmissionmatrix_h
30#define transmissionmatrix_h
31
32#pragma GCC diagnostic push
33#pragma GCC diagnostic ignored "-Wshadow"
34#pragma GCC diagnostic ignored "-Wconversion"
35#pragma GCC diagnostic ignored "-Wfloat-conversion"
36#if defined(__clang__)
37#pragma GCC diagnostic ignored "-Wdeprecated-copy-with-dtor"
38#pragma GCC diagnostic ignored "-Wdeprecated-copy-with-user-provided-dtor"
39#else
40#pragma GCC diagnostic ignored "-Wclass-memaccess"
41#endif
42
43#include <Eigen/Dense>
44
45#pragma GCC diagnostic pop
46
47#include "jacobian.h"
48#include "propagationmatrix.h"
49
50template <int N>
51Eigen::Matrix<Numeric, N, N> prop_matrix(const ConstVectorView& vec) {
52#define a vec[0]
53#define b vec[1]
54#define c vec[2]
55#define d vec[3]
56#define u vec[N]
57#define v vec[5]
58#define w vec[6]
59 static_assert (N>0 and N<5, "Bad size N");
60 if constexpr (N == 1) {
61 return Eigen::Matrix<Numeric, 1, 1>(a);
62 } else if constexpr (N==2) {
63 return (Eigen::Matrix2d() << a, b, b, a).finished();
64 } else if constexpr (N==3) {
65 return (Eigen::Matrix3d() << a, b, c, b, a, u, c, -u, a).finished();
66 } else if constexpr (N==4) {
67 return (Eigen::Matrix4d() << a, b, c, d, b, a, u, v, c, -u, a, w, d, -v, -w, a).finished();
68 }
69#undef a
70#undef b
71#undef c
72#undef d
73#undef u
74#undef v
75#undef w
76}
77
78template <int N>
79Eigen::Matrix<Numeric, N, N> prop_matrix(const ConstMatrixView& m) {
80 static_assert (N>0 and N<5, "Bad size N");
81 if constexpr (N == 1) {
82 return Eigen::Matrix<Numeric, 1, 1>(m(0, 0));
83 } else if constexpr (N == 2) {
84 return (Eigen::Matrix2d() << m(0, 0), m(0, 1),
85 m(1, 0), m(1, 1)).finished();
86 } else if constexpr (N == 3) {
87 return (Eigen::Matrix3d() << m(0, 0), m(0, 1), m(0, 2),
88 m(1, 0), m(1, 1), m(1, 2),
89 m(2, 0), m(2, 1), m(2, 2)) .finished();
90 } else if constexpr (N == 4) {
91 return (Eigen::Matrix4d() << m(0, 0), m(0, 1), m(0, 2), m(0, 3),
92 m(1, 0), m(1, 1), m(1, 2), m(1, 3),
93 m(2, 0), m(2, 1), m(2, 2), m(2, 3),
94 m(3, 0), m(3, 1), m(3, 2), m(3, 3)) .finished();
95 }
96}
97
98template <int N>
99Eigen::Matrix<Numeric, N, N> inv_prop_matrix(const ConstVectorView& vec) {
100#define a vec[0]
101#define b vec[1]
102#define c vec[2]
103#define d vec[3]
104#define u vec[N]
105#define v vec[5]
106#define w vec[6]
107 static_assert (N>0 and N<5, "Bad size N");
108 if constexpr (N == 1) {
109 return Eigen::Matrix<double, 1, 1>(1 / a);
110 } else if constexpr (N==2) {
111 return (Eigen::Matrix2d() << a, -b, -b, a).finished() / (a * a - b * b);
112 } else if constexpr (N==3) {
113 return (Eigen::Matrix3d() << a * a + u * u,
114 -a * b - c * u,
115 -a * c + b * u,
116 -a * b + c * u,
117 a * a - c * c,
118 -a * u + b * c,
119 -a * c - b * u,
120 a * u + b * c,
121 a * a - b * b)
122 .finished() /
123 (a * a * a - a * b * b - a * c * c + a * u * u);
124 } else if constexpr (N==4) {
125 return (Eigen::Matrix4d() << a * a * a + a * u * u + a * v * v + a * w * w,
126 -a * a * b - a * c * u - a * d * v - b * w * w + c * v * w -
127 d * u * w,
128 -a * a * c + a * b * u - a * d * w + b * v * w - c * v * v +
129 d * u * v,
130 -a * a * d + a * b * v + a * c * w - b * u * w + c * u * v -
131 d * u * u,
132 -a * a * b + a * c * u + a * d * v - b * w * w + c * v * w -
133 d * u * w,
134 a * a * a - a * c * c - a * d * d + a * w * w,
135 -a * a * u + a * b * c - a * v * w + b * d * w - c * d * v +
136 d * d * u,
137 -a * a * v + a * b * d + a * u * w - b * c * w + c * c * v -
138 c * d * u,
139 -a * a * c - a * b * u + a * d * w + b * v * w - c * v * v +
140 d * u * v,
141 a * a * u + a * b * c - a * v * w - b * d * w + c * d * v - d * d * u,
142 a * a * a - a * b * b - a * d * d + a * v * v,
143 -a * a * w + a * c * d - a * u * v + b * b * w - b * c * v +
144 b * d * u,
145 -a * a * d - a * b * v - a * c * w - b * u * w + c * u * v -
146 d * u * u,
147 a * a * v + a * b * d + a * u * w + b * c * w - c * c * v + c * d * u,
148 a * a * w + a * c * d - a * u * v - b * b * w + b * c * v - b * d * u,
149 a * a * a - a * b * b - a * c * c + a * u * u)
150 .finished() /
151 (a * a * a * a - a * a * b * b - a * a * c * c - a * a * d * d +
152 a * a * u * u + a * a * v * v + a * a * w * w - b * b * w * w +
153 2 * b * c * v * w - 2 * b * d * u * w - c * c * v * v +
154 2 * c * d * u * v - d * d * u * u);
155 }
156#undef a
157#undef b
158#undef c
159#undef d
160#undef u
161#undef v
162#undef w
163}
164
168 std::vector<Eigen::Matrix4d> T4;
169 std::vector<Eigen::Matrix3d> T3;
170 std::vector<Eigen::Matrix2d> T2;
171 std::vector<Eigen::Matrix<double, 1, 1>> T1;
172
178 TransmissionMatrix(Index nf = 0, Index stokes = 1);
179
184 TransmissionMatrix(TransmissionMatrix&& tm) noexcept = default;
185
191
198
205
211 explicit TransmissionMatrix(const PropagationMatrix& pm,
212 const Numeric& r = 1.0);
213
214 operator Tensor3() const;
215
216 explicit TransmissionMatrix(const ConstMatrixView& mat);
217
223 [[nodiscard]] const Eigen::Matrix4d& Mat4(size_t i) const;
224
230 [[nodiscard]] const Eigen::Matrix3d& Mat3(size_t i) const;
231
237 [[nodiscard]] const Eigen::Matrix2d& Mat2(size_t i) const;
238
244 [[nodiscard]] const Eigen::Matrix<double, 1, 1>& Mat1(size_t i) const;
245
251 [[nodiscard]] Eigen::MatrixXd Mat(size_t i) const;
252
258 Eigen::Matrix4d& Mat4(size_t i);
259
265 Eigen::Matrix3d& Mat3(size_t i);
266
272 Eigen::Matrix2d& Mat2(size_t i);
273
279 Eigen::Matrix<double, 1, 1>& Mat1(size_t i);
280
282 void setIdentity();
283
285 void setZero();
286
294 void mul(const TransmissionMatrix& A, const TransmissionMatrix& B);
295
303 void mul_aliased(const TransmissionMatrix& A, const TransmissionMatrix& B);
304
312 [[nodiscard]] Numeric operator()(const Index i,
313 const Index j,
314 const Index k) const;
315
323 Numeric& operator()(const Index i, const Index j, const Index k);
324
326 [[nodiscard]] Index Frequencies() const;
327
334
341
348
350 friend std::ostream& operator<<(std::ostream& os,
351 const TransmissionMatrix& tm);
352
354 friend std::istream& operator>>(std::istream& data, TransmissionMatrix& tm);
355
357 template <int N>
358 auto& TraMat(size_t i) noexcept {
359 static_assert(N > 0 and N < 5, "Bad size N");
360 if constexpr (N == 1)
361 return T1[i];
362 else if constexpr (N == 2)
363 return T2[i];
364 else if constexpr (N == 3)
365 return T3[i];
366 else if constexpr (N == 4)
367 return T4[i];
368 }
369
371 template <int N>
372 [[nodiscard]] auto& TraMat(size_t i) const noexcept {
373 static_assert(N > 0 and N < 5, "Bad size N");
374 if constexpr (N == 1)
375 return T1[i];
376 else if constexpr (N == 2)
377 return T2[i];
378 else if constexpr (N == 3)
379 return T3[i];
380 else if constexpr (N == 4)
381 return T4[i];
382 }
383
385 template <int N>
386 [[nodiscard]] Eigen::Matrix<Numeric, N, N> OptDepth(size_t i) const noexcept {
387 return (-TraMat<N>(i).diagonal().array().log().matrix()).asDiagonal();
388 }
389
409 template <int N>
410 [[nodiscard]] Eigen::Matrix<Numeric, N, 1> second_order_integration_source(
411 const Eigen::Matrix<Numeric, N, N> T,
412 const Eigen::Matrix<Numeric, N, 1> far,
413 const Eigen::Matrix<Numeric, N, 1> close,
414 const Eigen::Matrix<Numeric, N, N> Kfar,
415 const Eigen::Matrix<Numeric, N, N> Kclose,
416 const Numeric r) const noexcept {
417 static_assert(N > 0 and N < 5, "Bad size N");
418
419 const auto I = Eigen::Matrix<Numeric, N, N>::Identity();
420 if (T(0, 0) < 0.99) {
421 const auto od = 0.5 * r * (Kfar + Kclose);
422 Eigen::Matrix<Numeric, N, 1> second =
423 od.inverse() * ((I - (I + od) * T) * far + (od - I + T) * close);
424 if ((far[0] < close[0] and Kfar(0, 0) > Kclose(0, 0)) or
425 (far[0] > close[0] and Kfar(0, 0) < Kclose(0, 0))) {
426 Eigen::Matrix<Numeric, N, 1> second_limit =
427 0.5 * (I - T) * (far + close);
428
429 // FIXME: This is the equation given in the source material...
430 // const Eigen::Matrix<Numeric, N, 1> second_limit = 0.25 * r * (Kfar * far + Kclose * close);
431
432 if (second_limit[0] > second[0]) return second;
433 return second_limit;
434 }
435
436 return second;
437 }
438
439 return 0.5 * (I - T) * (far + close);
440 }
441
442 template <int N>
443 [[nodiscard]] Eigen::Matrix<Numeric, N, 1> second_order_integration_dsource(
444 size_t i,
445 const TransmissionMatrix& dx,
446 const Eigen::Matrix<Numeric, N, 1> far,
447 const Eigen::Matrix<Numeric, N, 1> close,
448 const Eigen::Matrix<Numeric, N, 1> d,
449 bool isfar) const noexcept {
450 static_assert(N > 0 and N < 5, "Bad size N");
451
452 const auto I = Eigen::Matrix<Numeric, N, N>::Identity();
453 const auto T = TraMat<N>(i);
454 const auto& dTdx = dx.TraMat<N>(i);
455 const Eigen::Matrix<Numeric, N, 1> first = 0.5 * (I - T) * (far + close);
456 if (T(0, 0) < 0.99) {
457 const auto od = OptDepth<N>(i);
458 const auto doddx = dx.OptDepth<N>(i);
459 const Eigen::Matrix<Numeric, N, 1> second =
460 od.inverse() * ((I - (I + od) * T) * far + (od - I + T) * close);
461 if (first[0] > second[0]) {
462 if (isfar) {
463 return od.inverse() *
464 ((I - (I + od) * T) * d - (I + od) * dTdx * far +
465 (doddx + T) * close - doddx * second);
466 }
467
468 return od.inverse() * (-(I + od) * dTdx * far + (od - I + T) * d +
469 (doddx + T) * close - doddx * second);
470 }
471
472 return 0.5 * ((I - T) * d - dTdx * (far + close));
473 }
474
475 return 0.5 * ((I - T) * d - dTdx * (far + close));
476 }
477};
478
486 const TransmissionMatrix& tm, const Numeric& x);
487
495 const Numeric& x, const TransmissionMatrix& tm);
496
500 std::vector<Eigen::Vector4d> R4;
501 std::vector<Eigen::Vector3d> R3;
502 std::vector<Eigen::Vector2d> R2;
503 std::vector<Eigen::Matrix<double, 1, 1>> R1;
504
510 RadiationVector(Index nf = 0, Index stokes = 1);
511
516 RadiationVector(RadiationVector&& rv) noexcept = default;
517
522 RadiationVector(const RadiationVector& rv) = default;
523
530
536 RadiationVector& operator=(RadiationVector&& rv) noexcept = default;
537
544
549 void leftMul(const TransmissionMatrix& T);
550
555 void SetZero(size_t i);
556
557 void SetZero();
558
564 [[nodiscard]] const Eigen::Vector4d& Vec4(size_t i) const;
565
571 [[nodiscard]] const Eigen::Vector3d& Vec3(size_t i) const;
572
578 [[nodiscard]] const Eigen::Vector2d& Vec2(size_t i) const;
579
585 [[nodiscard]] const Eigen::Matrix<double, 1, 1>& Vec1(size_t i) const;
586
592 [[nodiscard]] Eigen::VectorXd Vec(size_t i) const;
593
599 Eigen::Vector4d& Vec4(size_t i);
600
606 Eigen::Vector3d& Vec3(size_t i);
607
613 Eigen::Vector2d& Vec2(size_t i);
614
620 Eigen::Matrix<double, 1, 1>& Vec1(size_t i);
621
627 void rem_avg(const RadiationVector& O1, const RadiationVector& O2);
628
634 void add_avg(const RadiationVector& O1, const RadiationVector& O2);
635
642 void add_weighted(const TransmissionMatrix& T,
643 const RadiationVector& far,
644 const RadiationVector& close,
645 const ConstMatrixView& Kfar,
646 const ConstMatrixView& Kclose,
647 const Numeric r);
648
657 void addDerivEmission(const TransmissionMatrix& PiT,
658 const TransmissionMatrix& dT,
659 const TransmissionMatrix& T,
660 const RadiationVector& ImJ,
661 const RadiationVector& dJ);
662
672 const TransmissionMatrix& dT,
673 const TransmissionMatrix& T,
674 const RadiationVector& I,
675 const RadiationVector& far,
676 const RadiationVector& close,
677 const RadiationVector& d,
678 bool isfar);
679
687 const TransmissionMatrix& dT,
688 const RadiationVector& I);
689
699 void addMultiplied(const TransmissionMatrix& A, const RadiationVector& x);
700
709 const TransmissionMatrix& PiT,
710 const TransmissionMatrix& Z,
711 const TransmissionMatrix& dZ);
712
721 const TransmissionMatrix& Tr,
722 const TransmissionMatrix& Tf,
723 const TransmissionMatrix& Z);
724
733 const TransmissionMatrix& Tr,
734 const TransmissionMatrix& Tf,
735 const TransmissionMatrix& dZ);
736
743
750 const Numeric& operator()(const Index i, const Index j) const;
751
758 [[nodiscard]] Numeric& operator()(const Index i, const Index j);
763 operator Matrix() const;
764
772 void setSource(const StokesVector& a,
773 const ConstVectorView& B,
774 const StokesVector& S,
775 Index i);
776
778 [[nodiscard]] Index Frequencies() const;
779
781 friend std::ostream& operator<<(std::ostream& os, const RadiationVector& rv);
782
784 friend std::istream& operator>>(std::istream& data, RadiationVector& rv);
785};
786
795
800};
801
804 Forward,
805 Reverse,
806};
807
810 Emission,
813};
814
833 const RadiationVector& J1,
834 const RadiationVector& J2,
835 const ArrayOfRadiationVector& dJ1,
836 const ArrayOfRadiationVector& dJ2,
837 const TransmissionMatrix& T,
838 const TransmissionMatrix& PiT,
839 const ArrayOfTransmissionMatrix& dT1,
840 const ArrayOfTransmissionMatrix& dT2,
841 const PropagationMatrix& K1,
842 const PropagationMatrix& K2,
843 const ArrayOfPropagationMatrix& dK1,
844 const ArrayOfPropagationMatrix& dK2,
845 const Numeric r,
846 const Vector& dr1,
847 const Vector& dr2,
848 const Index ia,
849 const Index iz,
850 const RadiativeTransferSolver solver);
851
870 RadiationVector& J_add,
871 const PropagationMatrix& K,
872 const StokesVector& a,
873 const StokesVector& S,
874 const ArrayOfPropagationMatrix& dK,
875 const ArrayOfStokesVector& da,
876 const ArrayOfStokesVector& dS,
877 const ConstVectorView& B,
878 const ConstVectorView& dB_dT,
879 const ArrayOfRetrievalQuantity& jacobian_quantities,
880 const bool& jacobian_do);
881
899 const PropagationMatrix& K1,
900 const PropagationMatrix& K2,
901 const ArrayOfPropagationMatrix& dK1,
902 const ArrayOfPropagationMatrix& dK2,
903 const Numeric& r,
904 const Numeric& dr_dtemp1,
905 const Numeric& dr_dtemp2,
906 const Index temp_deriv_pos);
907
916 const CumulativeTransmission type) /*[[expects: T.nelem()>0]]*/;
917
935 const RadiationVector& I_incoming,
937 const ArrayOfTransmissionMatrix& PiTf,
938 const ArrayOfTransmissionMatrix& PiTr,
943 const BackscatterSolver solver);
944
957 const ConstMatrixView& pnd);
958
969 const ConstTensor5View& Pe, const ArrayOfMatrix& dpnd_dx);
970
971#endif // transmissionmatrix_h
This can be used to make arrays out of anything.
Definition: array.h:48
A constant view of a Matrix.
Definition: matpackI.h:1065
A constant view of a Tensor5.
Definition: matpackV.h:144
A constant view of a Vector.
Definition: matpackI.h:521
Class to help with hidden temporary variables for operations of type Numeric times Class.
The Matrix class.
Definition: matpackI.h:1285
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3 class.
Definition: matpackIII.h:352
The Vector class.
Definition: matpackI.h:910
Routines for setting up the jacobian.
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
Stuff related to the propagation matrix.
#define N
Definition: rng.cc:164
#define M
Definition: rng.cc:165
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.
#define u
ArrayOfTransmissionMatrix bulk_backscatter(const ConstTensor5View &Pe, const ConstMatrixView &pnd)
Bulk back-scattering.
#define d
BackscatterSolver
Intended to hold various backscatter solvers.
CumulativeTransmission
Intended to hold various ways to accumulate the transmission matrix.
#define v
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.
#define w
#define a
Eigen::Matrix< Numeric, N, N > inv_prop_matrix(const ConstVectorView &vec)
#define c
#define b
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)