ARTS 2.5.11 (git: 725533f0)
transmissionmatrix.cc
Go to the documentation of this file.
1
11#include "transmissionmatrix.h"
12
13#include "arts_conversions.h"
14#include "double_imanip.h"
15
16TransmissionMatrix::TransmissionMatrix(Index nf, Index stokes)
17 : stokes_dim(stokes),
18 T4(stokes_dim == 4 ? nf : 0, Eigen::Matrix4d::Identity()),
19 T3(stokes_dim == 3 ? nf : 0, Eigen::Matrix3d::Identity()),
20 T2(stokes_dim == 2 ? nf : 0, Eigen::Matrix2d::Identity()),
21 T1(stokes_dim == 1 ? nf : 0, Eigen::Matrix<double, 1, 1>::Identity()) {
22 ARTS_ASSERT(stokes_dim < 5 and stokes_dim > 0);
23}
24
25const Eigen::Matrix4d& TransmissionMatrix::Mat4(size_t i) const { return T4[i]; }
26const Eigen::Matrix3d& TransmissionMatrix::Mat3(size_t i) const { return T3[i]; }
27const Eigen::Matrix2d& TransmissionMatrix::Mat2(size_t i) const { return T2[i]; }
28const Eigen::Matrix<double, 1, 1>& TransmissionMatrix::Mat1(size_t i) const { return T1[i]; }
29
30Eigen::Matrix4d& TransmissionMatrix::Mat4(size_t i) { return T4[i]; }
31Eigen::Matrix3d& TransmissionMatrix::Mat3(size_t i) { return T3[i]; }
32Eigen::Matrix2d& TransmissionMatrix::Mat2(size_t i) { return T2[i]; }
33Eigen::Matrix<double, 1, 1>& TransmissionMatrix::Mat1(size_t i) { return T1[i]; }
34
35TransmissionMatrix& TransmissionMatrix::operator=(
36 const LazyScale<TransmissionMatrix>& lstm) {
37 operator=(lstm.bas);
38 operator*=(lstm.scale);
39 return *this;
40}
41
42TransmissionMatrix::operator Tensor3() const {
43 Tensor3 T(Frequencies(), stokes_dim, stokes_dim);
44 for (size_t i = 0; i < T4.size(); i++)
45 for (size_t j = 0; j < 4; j++)
46 for (size_t k = 0; k < 4; k++) T(i, j, k) = T4[i](j, k);
47 for (size_t i = 0; i < T3.size(); i++)
48 for (size_t j = 0; j < 3; j++)
49 for (size_t k = 0; k < 3; k++) T(i, j, k) = T3[i](j, k);
50 for (size_t i = 0; i < T2.size(); i++)
51 for (size_t j = 0; j < 2; j++)
52 for (size_t k = 0; k < 2; k++) T(i, j, k) = T2[i](j, k);
53 for (size_t i = 0; i < T1.size(); i++) T(i, 0, 0) = T1[i](0, 0);
54 return T;
55}
56
57#pragma GCC diagnostic push
58#pragma GCC diagnostic ignored "-Wreturn-type"
59[[nodiscard]] Eigen::MatrixXd TransmissionMatrix::Mat(size_t i) const {
60 switch (stokes_dim) {
61 case 1:
62 return Mat1(i);
63 case 2:
64 return Mat2(i);
65 case 3:
66 return Mat3(i);
67 case 4:
68 return Mat4(i);
69 }
70}
71#pragma GCC diagnostic pop
72
73void TransmissionMatrix::setIdentity() {
74 std::fill(T4.begin(), T4.end(), Eigen::Matrix4d::Identity());
75 std::fill(T3.begin(), T3.end(), Eigen::Matrix3d::Identity());
76 std::fill(T2.begin(), T2.end(), Eigen::Matrix2d::Identity());
77 std::fill(T1.begin(), T1.end(), Eigen::Matrix<double, 1, 1>::Identity());
78}
79
80void TransmissionMatrix::setZero() {
81 std::fill(T4.begin(), T4.end(), Eigen::Matrix4d::Zero());
82 std::fill(T3.begin(), T3.end(), Eigen::Matrix3d::Zero());
83 std::fill(T2.begin(), T2.end(), Eigen::Matrix2d::Zero());
84 std::fill(T1.begin(), T1.end(), Eigen::Matrix<double, 1, 1>::Zero());
85}
86
87void TransmissionMatrix::mul(const TransmissionMatrix& A,
88 const TransmissionMatrix& B) {
89 for (size_t i = 0; i < T4.size(); i++) T4[i].noalias() = A.T4[i] * B.T4[i];
90 for (size_t i = 0; i < T3.size(); i++) T3[i].noalias() = A.T3[i] * B.T3[i];
91 for (size_t i = 0; i < T2.size(); i++) T2[i].noalias() = A.T2[i] * B.T2[i];
92 for (size_t i = 0; i < T1.size(); i++) T1[i].noalias() = A.T1[i] * B.T1[i];
93}
94
95void TransmissionMatrix::mul_aliased(const TransmissionMatrix& A,
96 const TransmissionMatrix& B) {
97 for (size_t i = 0; i < T4.size(); i++) T4[i] = A.T4[i] * B.T4[i];
98 for (size_t i = 0; i < T3.size(); i++) T3[i] = A.T3[i] * B.T3[i];
99 for (size_t i = 0; i < T2.size(); i++) T2[i] = A.T2[i] * B.T2[i];
100 for (size_t i = 0; i < T1.size(); i++) T1[i] = A.T1[i] * B.T1[i];
101}
102
103Numeric TransmissionMatrix::operator()(const Index i,
104 const Index j,
105 const Index k) const {
106 switch (stokes_dim) {
107 case 4:
108 return T4[i](j, k);
109 case 3:
110 return T3[i](j, k);
111 case 2:
112 return T2[i](j, k);
113 default:
114 return T1[i](j, k);
115 }
116}
117
118Numeric& TransmissionMatrix::operator()(const Index i, const Index j, const Index k) {
119 switch (stokes_dim) {
120 case 4:
121 return T4[i](j, k);
122 case 3:
123 return T3[i](j, k);
124 case 2:
125 return T2[i](j, k);
126 default:
127 return T1[i](j, k);
128 }
129}
130
131[[nodiscard]] Index TransmissionMatrix::Frequencies() const {
132 switch (stokes_dim) {
133 case 4:
134 return Index(T4.size());
135 case 3:
136 return Index(T3.size());
137 case 2:
138 return Index(T2.size());
139 default:
140 return Index(T1.size());
141 }
142}
143
144TransmissionMatrix::TransmissionMatrix(const ConstMatrixView& mat): TransmissionMatrix(1, mat.nrows()){
145 ARTS_ASSERT(mat.nrows() == mat.ncols());
146 for (Index i = 0; i < stokes_dim; i++)
147 for (Index j = 0; j < stokes_dim; j++)
148 operator()(0, i, j) = mat(i,j);
149};
150
151TransmissionMatrix& TransmissionMatrix::operator+=(
152 const LazyScale<TransmissionMatrix>& lstm) {
153 for (size_t i = 0; i < T4.size(); i++)
154 T4[i].noalias() = lstm.scale * lstm.bas.Mat4(i);
155 for (size_t i = 0; i < T3.size(); i++)
156 T3[i].noalias() = lstm.scale * lstm.bas.Mat3(i);
157 for (size_t i = 0; i < T2.size(); i++)
158 T2[i].noalias() = lstm.scale * lstm.bas.Mat2(i);
159 for (size_t i = 0; i < T1.size(); i++)
160 T1[i].noalias() = lstm.scale * lstm.bas.Mat1(i);
161 return *this;
162}
163
164TransmissionMatrix& TransmissionMatrix::operator*=(const Numeric& scale) {
165 std::transform(
166 T4.begin(), T4.end(), T4.begin(), [scale](auto& T) { return T * scale; });
167 std::transform(
168 T3.begin(), T3.end(), T3.begin(), [scale](auto& T) { return T * scale; });
169 std::transform(
170 T2.begin(), T2.end(), T2.begin(), [scale](auto& T) { return T * scale; });
171 std::transform(
172 T1.begin(), T1.end(), T1.begin(), [scale](auto& T) { return T * scale; });
173 return *this;
174}
175
176LazyScale<TransmissionMatrix> operator*(const TransmissionMatrix& tm,
177 const Numeric& x) {
178 return {tm, x};
179}
180
181LazyScale<TransmissionMatrix> operator*(const Numeric& x,
182 const TransmissionMatrix& tm) {
183 return {tm, x};
184}
185
186RadiationVector::RadiationVector(Index nf, Index stokes)
187 : stokes_dim(stokes),
188 R4(stokes_dim == 4 ? nf : 0, Eigen::Vector4d::Zero()),
189 R3(stokes_dim == 3 ? nf : 0, Eigen::Vector3d::Zero()),
190 R2(stokes_dim == 2 ? nf : 0, Eigen::Vector2d::Zero()),
191 R1(stokes_dim == 1 ? nf : 0, Eigen::Matrix<double, 1, 1>::Zero()) {
192 ARTS_ASSERT(stokes_dim < 5 and stokes_dim > 0);
193}
194
195const Eigen::Vector4d& RadiationVector::Vec4(size_t i) const { return R4[i]; }
196const Eigen::Vector3d& RadiationVector::Vec3(size_t i) const { return R3[i]; }
197const Eigen::Vector2d& RadiationVector::Vec2(size_t i) const { return R2[i]; }
198const Eigen::Matrix<double, 1, 1>& RadiationVector::Vec1(size_t i) const {
199 return R1[i];
200}
201
202Eigen::Vector4d& RadiationVector::Vec4(size_t i) { return R4[i]; }
203Eigen::Vector3d& RadiationVector::Vec3(size_t i) { return R3[i]; }
204Eigen::Vector2d& RadiationVector::Vec2(size_t i) { return R2[i]; }
205Eigen::Matrix<double, 1, 1>& RadiationVector::Vec1(size_t i) { return R1[i]; }
206
207Numeric& RadiationVector::operator()(const Index i, const Index j) {
208 switch (stokes_dim) {
209 case 4:
210 return R4[i][j];
211 case 3:
212 return R3[i][j];
213 case 2:
214 return R2[i][j];
215 default:
216 return R1[i][j];
217 }
218}
219
220void RadiationVector::leftMul(const TransmissionMatrix& T) {
221 for (size_t i = 0; i < R4.size(); i++) R4[i] = T.Mat4(i) * R4[i];
222 for (size_t i = 0; i < R3.size(); i++) R3[i] = T.Mat3(i) * R3[i];
223 for (size_t i = 0; i < R2.size(); i++) R2[i] = T.Mat2(i) * R2[i];
224 for (size_t i = 0; i < R1.size(); i++) R1[i] = T.Mat1(i) * R1[i];
225}
226
227void RadiationVector::SetZero(size_t i) {
228 switch (stokes_dim) {
229 case 4:
230 R4[i].noalias() = Eigen::Vector4d::Zero();
231 break;
232 case 3:
233 R3[i].noalias() = Eigen::Vector3d::Zero();
234 break;
235 case 2:
236 R2[i].noalias() = Eigen::Vector2d::Zero();
237 break;
238 case 1:
239 R1[i][0] = 0;
240 break;
241 }
242}
243
244void RadiationVector::SetZero() {
245 for (Index i=0; i<Frequencies(); i++) SetZero(i);
246}
247
248Eigen::VectorXd RadiationVector::Vec(size_t i) const {
249 switch (stokes_dim) {
250 case 4:
251 return Vec4(i);
252 case 3:
253 return Vec3(i);
254 case 2:
255 return Vec2(i);
256 default:
257 return Vec1(i);
258 }
259}
260
261void RadiationVector::rem_avg(const RadiationVector& O1,
262 const RadiationVector& O2) {
263 for (size_t i = 0; i < R4.size(); i++)
264 R4[i].noalias() -= 0.5 * (O1.R4[i] + O2.R4[i]);
265 for (size_t i = 0; i < R3.size(); i++)
266 R3[i].noalias() -= 0.5 * (O1.R3[i] + O2.R3[i]);
267 for (size_t i = 0; i < R2.size(); i++)
268 R2[i].noalias() -= 0.5 * (O1.R2[i] + O2.R2[i]);
269 for (size_t i = 0; i < R1.size(); i++)
270 R1[i].noalias() -= 0.5 * (O1.R1[i] + O2.R1[i]);
271}
272
273void RadiationVector::add_avg(const RadiationVector& O1,
274 const RadiationVector& O2) {
275 for (size_t i = 0; i < R4.size(); i++)
276 R4[i].noalias() += 0.5 * (O1.R4[i] + O2.R4[i]);
277 for (size_t i = 0; i < R3.size(); i++)
278 R3[i].noalias() += 0.5 * (O1.R3[i] + O2.R3[i]);
279 for (size_t i = 0; i < R2.size(); i++)
280 R2[i].noalias() += 0.5 * (O1.R2[i] + O2.R2[i]);
281 for (size_t i = 0; i < R1.size(); i++)
282 R1[i].noalias() += 0.5 * (O1.R1[i] + O2.R1[i]);
283}
284
285void RadiationVector::add_weighted(const TransmissionMatrix& T,
286 const RadiationVector& far,
287 const RadiationVector& close,
288 const ConstMatrixView& Kfar,
289 const ConstMatrixView& Kclose,
290 const Numeric r) {
291 for (size_t i = 0; i < R4.size(); i++) {
292 R4[i].noalias() +=
293 T.second_order_integration_source<4>(T.TraMat<4>(i),
294 far.R4[i],
295 close.R4[i],
296 prop_matrix<4>(Kfar(i, joker)),
297 prop_matrix<4>(Kclose(i, joker)),
298 r);
299 }
300 for (size_t i = 0; i < R3.size(); i++) {
301 R3[i].noalias() +=
302 T.second_order_integration_source<3>(T.TraMat<3>(i),
303 far.R3[i],
304 close.R3[i],
305 prop_matrix<3>(Kfar(i, joker)),
306 prop_matrix<3>(Kclose(i, joker)),
307 r);
308 }
309 for (size_t i = 0; i < R2.size(); i++) {
310 R2[i].noalias() +=
311 T.second_order_integration_source<2>(T.TraMat<2>(i),
312 far.R2[i],
313 close.R2[i],
314 prop_matrix<2>(Kfar(i, joker)),
315 prop_matrix<2>(Kclose(i, joker)),
316 r);
317 }
318 for (size_t i = 0; i < R1.size(); i++) {
319 R1[i].noalias() +=
320 T.second_order_integration_source<1>(T.TraMat<1>(i),
321 far.R1[i],
322 close.R1[i],
323 prop_matrix<1>(Kfar(i, joker)),
324 prop_matrix<1>(Kclose(i, joker)),
325 r);
326 }
327}
328
329void RadiationVector::addDerivEmission(const TransmissionMatrix& PiT,
330 const TransmissionMatrix& dT,
331 const TransmissionMatrix& T,
332 const RadiationVector& ImJ,
333 const RadiationVector& dJ) {
334 for (size_t i = 0; i < R4.size(); i++)
335 R4[i].noalias() += PiT.Mat4(i) * (dT.Mat4(i) * ImJ.R4[i] + dJ.R4[i] -
336 T.Mat4(i) * dJ.R4[i]);
337 for (size_t i = 0; i < R3.size(); i++)
338 R3[i].noalias() += PiT.Mat3(i) * (dT.Mat3(i) * ImJ.R3[i] + dJ.R3[i] -
339 T.Mat3(i) * dJ.R3[i]);
340 for (size_t i = 0; i < R2.size(); i++)
341 R2[i].noalias() += PiT.Mat2(i) * (dT.Mat2(i) * ImJ.R2[i] + dJ.R2[i] -
342 T.Mat2(i) * dJ.R2[i]);
343 for (size_t i = 0; i < R1.size(); i++)
344 R1[i].noalias() += PiT.Mat1(i) * (dT.Mat1(i) * ImJ.R1[i] + dJ.R1[i] -
345 T.Mat1(i) * dJ.R1[i]);
346}
347
348void RadiationVector::addWeightedDerivEmission(const TransmissionMatrix& PiT,
349 const TransmissionMatrix& dT,
350 const TransmissionMatrix& T,
351 const RadiationVector& I,
352 const RadiationVector& far,
353 const RadiationVector& close,
354 const RadiationVector& d,
355 bool isfar) {
356 for (size_t i = 0; i < R4.size(); i++)
357 R4[i].noalias() +=
358 PiT.Mat4(i) * (dT.Mat4(i) * I.R4[i] +
359 T.second_order_integration_dsource<4>(
360 i, dT, far.R4[i], close.R4[i], d.R4[i], isfar));
361 for (size_t i = 0; i < R3.size(); i++)
362 R3[i].noalias() +=
363 PiT.Mat3(i) * (dT.Mat3(i) * I.R3[i] +
364 T.second_order_integration_dsource<3>(
365 i, dT, far.R3[i], close.R3[i], d.R3[i], isfar));
366 for (size_t i = 0; i < R2.size(); i++)
367 R2[i].noalias() +=
368 PiT.Mat2(i) * (dT.Mat2(i) * I.R2[i] +
369 T.second_order_integration_dsource<2>(
370 i, dT, far.R2[i], close.R2[i], d.R2[i], isfar));
371 for (size_t i = 0; i < R1.size(); i++)
372 R1[i].noalias() +=
373 PiT.Mat1(i) * (dT.Mat1(i) * I.R1[i] +
374 T.second_order_integration_dsource<1>(
375 i, dT, far.R1[i], close.R1[i], d.R1[i], isfar));
376}
377
378void RadiationVector::addDerivTransmission(const TransmissionMatrix& PiT,
379 const TransmissionMatrix& dT,
380 const RadiationVector& I) {
381 for (size_t i = 0; i < R4.size(); i++)
382 R4[i].noalias() += PiT.Mat4(i) * dT.Mat4(i) * I.R4[i];
383 for (size_t i = 0; i < R3.size(); i++)
384 R3[i].noalias() += PiT.Mat3(i) * dT.Mat3(i) * I.R3[i];
385 for (size_t i = 0; i < R2.size(); i++)
386 R2[i].noalias() += PiT.Mat2(i) * dT.Mat2(i) * I.R2[i];
387 for (size_t i = 0; i < R1.size(); i++)
388 R1[i].noalias() += PiT.Mat1(i) * dT.Mat1(i) * I.R1[i];
389}
390
391void RadiationVector::addMultiplied(const TransmissionMatrix& A,
392 const RadiationVector& x) {
393 for (size_t i = 0; i < R4.size(); i++) R4[i].noalias() += A.Mat4(i) * x.R4[i];
394 for (size_t i = 0; i < R3.size(); i++) R3[i].noalias() += A.Mat3(i) * x.R3[i];
395 for (size_t i = 0; i < R2.size(); i++) R2[i].noalias() += A.Mat2(i) * x.R2[i];
396 for (size_t i = 0; i < R1.size(); i++) R1[i].noalias() += A.Mat1(i) * x.R1[i];
397}
398
399void RadiationVector::setDerivReflection(const RadiationVector& I,
400 const TransmissionMatrix& PiT,
401 const TransmissionMatrix& Z,
402 const TransmissionMatrix& dZ) {
403 for (size_t i = 0; i < R4.size(); i++)
404 R4[i] = PiT.Mat4(i) * (Z.Mat4(i) * R4[i] + dZ.Mat4(i) * I.R4[i]);
405 for (size_t i = 0; i < R3.size(); i++)
406 R3[i] = PiT.Mat3(i) * (Z.Mat3(i) * R3[i] + dZ.Mat3(i) * I.R3[i]);
407 for (size_t i = 0; i < R2.size(); i++)
408 R2[i] = PiT.Mat2(i) * (Z.Mat2(i) * R2[i] + dZ.Mat2(i) * I.R2[i]);
409 for (size_t i = 0; i < R1.size(); i++)
410 R4[i][0] =
411 PiT.Mat1(i)[0] * (Z.Mat1(i)[0] * R1[i][0] + dZ.Mat1(i)[0] * I.R1[i][0]);
412}
413
414void RadiationVector::setBackscatterTransmission(const RadiationVector& I0,
415 const TransmissionMatrix& Tr,
416 const TransmissionMatrix& Tf,
417 const TransmissionMatrix& Z) {
418 for (size_t i = 0; i < R4.size(); i++)
419 R4[i].noalias() = Tr.Mat4(i) * Z.Mat4(i) * Tf.Mat4(i) * I0.R4[i];
420 for (size_t i = 0; i < R3.size(); i++)
421 R3[i].noalias() = Tr.Mat3(i) * Z.Mat3(i) * Tf.Mat3(i) * I0.R3[i];
422 for (size_t i = 0; i < R2.size(); i++)
423 R2[i].noalias() = Tr.Mat2(i) * Z.Mat2(i) * Tf.Mat2(i) * I0.R2[i];
424 for (size_t i = 0; i < R1.size(); i++)
425 R1[i].noalias() = Tr.Mat1(i) * Z.Mat1(i) * Tf.Mat1(i) * I0.R1[i];
426}
427
428void RadiationVector::setBackscatterTransmissionDerivative(
429 const RadiationVector& I0,
430 const TransmissionMatrix& Tr,
431 const TransmissionMatrix& Tf,
432 const TransmissionMatrix& dZ) {
433 for (size_t i = 0; i < R4.size(); i++)
434 R4[i].noalias() += Tr.Mat4(i) * dZ.Mat4(i) * Tf.Mat4(i) * I0.R4[i];
435 for (size_t i = 0; i < R3.size(); i++)
436 R3[i].noalias() += Tr.Mat3(i) * dZ.Mat3(i) * Tf.Mat3(i) * I0.R3[i];
437 for (size_t i = 0; i < R2.size(); i++)
438 R2[i].noalias() += Tr.Mat2(i) * dZ.Mat2(i) * Tf.Mat2(i) * I0.R2[i];
439 for (size_t i = 0; i < R1.size(); i++)
440 R1[i].noalias() += Tr.Mat1(i) * dZ.Mat1(i) * Tf.Mat1(i) * I0.R1[i];
441}
442
443RadiationVector& RadiationVector::operator=(const ConstMatrixView& M) {
444 ARTS_ASSERT(M.ncols() == stokes_dim and M.nrows() == Frequencies());
445 for (size_t i = 0; i < R4.size(); i++) {
446 R4[i][0] = M(i, 0);
447 R4[i][1] = M(i, 1);
448 R4[i][2] = M(i, 2);
449 R4[i][3] = M(i, 3);
450 }
451 for (size_t i = 0; i < R3.size(); i++) {
452 R3[i][0] = M(i, 0);
453 R3[i][1] = M(i, 1);
454 R3[i][2] = M(i, 2);
455 }
456 for (size_t i = 0; i < R2.size(); i++) {
457 R2[i][0] = M(i, 0);
458 R2[i][1] = M(i, 1);
459 }
460 for (size_t i = 0; i < R1.size(); i++) {
461 R1[i][0] = M(i, 0);
462 }
463 return *this;
464}
465
466RadiationVector& RadiationVector::operator+=(const RadiationVector& rv) {
467 for (size_t i = 0; i < R4.size(); i++)
468 R4[i].noalias() += rv.R4[i];
469 for (size_t i = 0; i < R3.size(); i++)
470 R3[i].noalias() += rv.R3[i];
471 for (size_t i = 0; i < R2.size(); i++)
472 R2[i].noalias() += rv.R2[i];
473 for (size_t i = 0; i < R1.size(); i++)
474 R1[i].noalias() += rv.R1[i];
475
476 return *this;
477}
478
479
480const Numeric& RadiationVector::operator()(const Index i, const Index j) const {
481 switch (stokes_dim) {
482 case 4:
483 return R4[i][j];
484 case 3:
485 return R3[i][j];
486 case 2:
487 return R2[i][j];
488 default:
489 return R1[i][j];
490 }
491}
492
493RadiationVector::operator Matrix() const {
494 Matrix M(Frequencies(), stokes_dim);
495 for (size_t i = 0; i < R4.size(); i++)
496 for (size_t j = 0; j < 4; j++) M(i, j) = R4[i](j);
497 for (size_t i = 0; i < R3.size(); i++)
498 for (size_t j = 0; j < 3; j++) M(i, j) = R3[i](j);
499 for (size_t i = 0; i < R2.size(); i++)
500 for (size_t j = 0; j < 2; j++) M(i, j) = R2[i](j);
501 for (size_t i = 0; i < R1.size(); i++) M(i, 0) = R1[i](0);
502 return M;
503}
504
505void RadiationVector::setSource(const StokesVector& a,
506 const ConstVectorView& B,
507 const StokesVector& S,
508 Index i) {
509 ARTS_ASSERT(a.NumberOfAzimuthAngles() == 1);
510 ARTS_ASSERT(a.NumberOfZenithAngles() == 1);
511 ARTS_ASSERT(S.NumberOfAzimuthAngles() == 1);
512 ARTS_ASSERT(S.NumberOfZenithAngles() == 1);
513 switch (stokes_dim) {
514 case 4:
515 if (not S.IsEmpty())
516 R4[i].noalias() =
517 Eigen::Vector4d(a.Kjj()[i], a.K12()[i], a.K13()[i], a.K14()[i]) *
518 B[i] +
519 Eigen::Vector4d(S.Kjj()[i], S.K12()[i], S.K13()[i], S.K14()[i]);
520 else
521 R4[i].noalias() =
522 Eigen::Vector4d(a.Kjj()[i], a.K12()[i], a.K13()[i], a.K14()[i]) *
523 B[i];
524 break;
525 case 3:
526 if (not S.IsEmpty())
527 R3[i].noalias() =
528 Eigen::Vector3d(a.Kjj()[i], a.K12()[i], a.K13()[i]) * B[i] +
529 Eigen::Vector3d(S.Kjj()[i], S.K12()[i], S.K13()[i]);
530 else
531 R3[i].noalias() =
532 Eigen::Vector3d(a.Kjj()[i], a.K12()[i], a.K13()[i]) * B[i];
533 break;
534 case 2:
535 if (not S.IsEmpty())
536 R2[i].noalias() = Eigen::Vector2d(a.Kjj()[i], a.K12()[i]) * B[i] +
537 Eigen::Vector2d(S.Kjj()[i], S.K12()[i]);
538 else
539 R2[i].noalias() = Eigen::Vector2d(a.Kjj()[i], a.K12()[i]) * B[i];
540 break;
541 default:
542 if (not S.IsEmpty())
543 R1[i][0] = a.Kjj()[i] * B[i] + S.Kjj()[i];
544 else
545 R1[i][0] = a.Kjj()[i] * B[i];
546 }
547}
548
549Index RadiationVector::Frequencies() const {
550 switch (stokes_dim) {
551 case 4:
552 return Index(R4.size());
553 case 3:
554 return Index(R3.size());
555 case 2:
556 return Index(R2.size());
557 default:
558 return Index(R1.size());
559 }
560}
561
562constexpr Numeric lower_is_considered_zero_for_sinc_likes = 1e-4;
563
564template <int N>
565constexpr Eigen::Matrix<Numeric, N, 1> source_vector(
566 const StokesVector& a,
567 const ConstVectorView& B,
568 const StokesVector& da,
569 const ConstVectorView& dB_dT,
570 const StokesVector& dS,
571 bool dT,
572 Index i) {
573 static_assert(N > 0 and N < 5, "Bad stokes dimensions");
574
575 if constexpr (N == 1) {
576 if (dT)
577 return Eigen::Matrix<Numeric, 1, 1>(dS.Kjj()[i] + da.Kjj()[i] * B[i] +
578 a.Kjj()[i] * dB_dT[i]);
579 return Eigen::Matrix<Numeric, 1, 1>(dS.Kjj()[i] + da.Kjj()[i] * B[i]);
580 }
581
582 if constexpr (N == 2) {
583 if (dT)
584 return Eigen::Vector2d(
585 dS.Kjj()[i] + da.Kjj()[i] * B[i] + a.Kjj()[i] * dB_dT[i],
586 dS.K12()[i] + da.K12()[i] * B[i] + a.K12()[i] * dB_dT[i]);
587 return Eigen::Vector2d(dS.Kjj()[i] + da.Kjj()[i] * B[i],
588 dS.K12()[i] + da.K12()[i] * B[i]);
589 }
590
591 if constexpr (N == 3) {
592 if (dT)
593 return Eigen::Vector3d(
594 dS.Kjj()[i] + da.Kjj()[i] * B[i] + a.Kjj()[i] * dB_dT[i],
595 dS.K12()[i] + da.K12()[i] * B[i] + a.K12()[i] * dB_dT[i],
596 dS.K13()[i] + da.K13()[i] * B[i] + a.K13()[i] * dB_dT[i]);
597 return Eigen::Vector3d(dS.Kjj()[i] + da.Kjj()[i] * B[i],
598 dS.K12()[i] + da.K12()[i] * B[i],
599 dS.K13()[i] + da.K13()[i] * B[i]);
600 }
601
602 if constexpr (N == 4) {
603 if (dT)
604 return Eigen::Vector4d(
605 dS.Kjj()[i] + da.Kjj()[i] * B[i] + a.Kjj()[i] * dB_dT[i],
606 dS.K12()[i] + da.K12()[i] * B[i] + a.K12()[i] * dB_dT[i],
607 dS.K13()[i] + da.K13()[i] * B[i] + a.K13()[i] * dB_dT[i],
608 dS.K14()[i] + da.K14()[i] * B[i] + a.K14()[i] * dB_dT[i]);
609 return Eigen::Vector4d(dS.Kjj()[i] + da.Kjj()[i] * B[i],
610 dS.K12()[i] + da.K12()[i] * B[i],
611 dS.K13()[i] + da.K13()[i] * B[i],
612 dS.K14()[i] + da.K14()[i] * B[i]);
613 }
614}
615
616template <int N>
617constexpr Eigen::Matrix<Numeric, N, 1> source_vector(
618 const StokesVector& a,
619 const ConstVectorView& B,
620 const StokesVector& da,
621 const ConstVectorView& dB_dT,
622 bool dT,
623 Index i) {
624 static_assert(N > 0 and N < 5, "Bad stokes dimensions");
625
626 if constexpr (N == 1) {
627 if (dT)
628 return Eigen::Matrix<Numeric, 1, 1>(da.Kjj()[i] * B[i] +
629 a.Kjj()[i] * dB_dT[i]);
630 return Eigen::Matrix<Numeric, 1, 1>(da.Kjj()[i] * B[i]);
631 }
632
633 if constexpr (N == 2) {
634 if (dT)
635 return Eigen::Vector2d(da.Kjj()[i] * B[i] + a.Kjj()[i] * dB_dT[i],
636 da.K12()[i] * B[i] + a.K12()[i] * dB_dT[i]);
637 return Eigen::Vector2d(da.Kjj()[i] * B[i], da.K12()[i] * B[i]);
638 }
639
640 if constexpr (N == 3) {
641 if (dT)
642 return Eigen::Vector3d(da.Kjj()[i] * B[i] + a.Kjj()[i] * dB_dT[i],
643 da.K12()[i] * B[i] + a.K12()[i] * dB_dT[i],
644 da.K13()[i] * B[i] + a.K13()[i] * dB_dT[i]);
645 return Eigen::Vector3d(
646 da.Kjj()[i] * B[i], da.K12()[i] * B[i], da.K13()[i] * B[i]);
647 }
648
649 if constexpr (N == 4) {
650 if (dT)
651 return Eigen::Vector4d(da.Kjj()[i] * B[i] + a.Kjj()[i] * dB_dT[i],
652 da.K12()[i] * B[i] + a.K12()[i] * dB_dT[i],
653 da.K13()[i] * B[i] + a.K13()[i] * dB_dT[i],
654 da.K14()[i] * B[i] + a.K14()[i] * dB_dT[i]);
655 return Eigen::Vector4d(da.Kjj()[i] * B[i],
656 da.K12()[i] * B[i],
657 da.K13()[i] * B[i],
658 da.K14()[i] * B[i]);
659 }
660}
661
662inline void transmat1(TransmissionMatrix& T,
663 const PropagationMatrix& K1,
664 const PropagationMatrix& K2,
665 const Numeric& r,
666 const Index iz = 0,
667 const Index ia = 0) noexcept {
668 for (Index i = 0; i < K1.NumberOfFrequencies(); i++)
669 T.Mat1(i)(0, 0) =
670 std::exp(-0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]));
671}
672
673inline void transmat2(TransmissionMatrix& T,
674 const PropagationMatrix& K1,
675 const PropagationMatrix& K2,
676 const Numeric& r,
677 const Index iz = 0,
678 const Index ia = 0) noexcept {
679 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
680 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
681 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]);
682 const Numeric exp_a = std::exp(a);
683 const Numeric cb = std::cosh(b), sb = std::sinh(b);
684 T.Mat2(i).noalias() =
685 (Eigen::Matrix2d() << cb, sb, sb, cb).finished() * exp_a;
686 }
687}
688
689inline void transmat3(TransmissionMatrix& T,
690 const PropagationMatrix& K1,
691 const PropagationMatrix& K2,
692 const Numeric& r,
693 const Index iz = 0,
694 const Index ia = 0) noexcept {
695 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
696 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
697 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]),
698 c = -0.5 * r * (K1.K13(iz, ia)[i] + K2.K13(iz, ia)[i]),
699 u = -0.5 * r * (K1.K23(iz, ia)[i] + K2.K23(iz, ia)[i]);
700 const Numeric exp_a = std::exp(a);
701
702 if (b == 0. and c == 0. and u == 0.) {
703 T.Mat3(i).noalias() = Eigen::Matrix3d::Identity() * exp_a;
704 } else {
705 const Numeric a2 = a * a, b2 = b * b, c2 = c * c, u2 = u * u;
706 const Numeric Const = b2 + c2 - u2;
707
708 const bool real = Const > 0;
709 const bool imag = Const < 0;
710 const bool either = real or imag;
711
712 const Numeric x =
713 std::sqrt(imag ? -Const : Const); // test to just use real values
714 const Numeric x2 =
715 (real ? 1 : -1) * x * x; // test to change sign if imaginary
716 const Numeric inv_x2 =
717 either ? 1.0 / x2
718 : 1.0; // test so further calculations are replaced as x→0
719
720 const Numeric sx =
721 real ? std::sinh(x) : std::sin(-x); // -i sin(ix) → sinh(x)
722 const Numeric cx =
723 real ? std::cosh(x) : std::cos(+x); // cos(ix) → cosh(x)
724
725 /* Using:
726 * lim x→0 [(cosh(x) - 1) / x^2] → 1/2
727 * lim x→0 [sinh(x) / x] → 1
728 * inv_x2 := 1 for x == 0,
729 * C0, C1, C2 ∝ [1/x^2]
730 */
731 const Numeric C0 =
732 either ? a2 * (cx - 1.0) - a * x * sx + x2 : 1.0 + 0.5 * a2 - a;
733 const Numeric C1 = either ? 2.0 * a * (1.0 - cx) + x * sx : 1.0 - a;
734 const Numeric C2 = either ? cx - 1.0 : 0.5;
735
736 T.Mat3(i).noalias() =
737 exp_a * inv_x2 *
738 (Eigen::Matrix3d() << C0 + C1 * a + C2 * (a2 + b2 + c2),
739 C1 * b + C2 * (2 * a * b - c * u),
740 C1 * c + C2 * (2 * a * c + b * u),
741 C1 * b + C2 * (2 * a * b + c * u),
742 C0 + C1 * a + C2 * (a2 + b2 - u2),
743 C1 * u + C2 * (2 * a * u + b * c),
744 C1 * c + C2 * (2 * a * c - b * u),
745 -C1 * u - C2 * (2 * a * u - b * c),
746 C0 + C1 * a + C2 * (a2 + c2 - u2))
747 .finished();
748 }
749 }
750}
751
752inline void transmat4(TransmissionMatrix& T,
753 const PropagationMatrix& K1,
754 const PropagationMatrix& K2,
755 const Numeric& r,
756 const Index iz = 0,
757 const Index ia = 0) noexcept {
758 static constexpr Numeric sqrt_05 = Constant::inv_sqrt_2;
759 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
760 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
761 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]),
762 c = -0.5 * r * (K1.K13(iz, ia)[i] + K2.K13(iz, ia)[i]),
763 d = -0.5 * r * (K1.K14(iz, ia)[i] + K2.K14(iz, ia)[i]),
764 u = -0.5 * r * (K1.K23(iz, ia)[i] + K2.K23(iz, ia)[i]),
765 v = -0.5 * r * (K1.K24(iz, ia)[i] + K2.K24(iz, ia)[i]),
766 w = -0.5 * r * (K1.K34(iz, ia)[i] + K2.K34(iz, ia)[i]);
767 const Numeric exp_a = std::exp(a);
768
769 if (b == 0. and c == 0. and d == 0. and u == 0. and v == 0. and w == 0.)
770 T.Mat4(i).noalias() = Eigen::Matrix4d::Identity() * exp_a;
771 else {
772 const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
773 w2 = w * w;
774
775 const Numeric tmp =
776 w2 * w2 + 2 * (b2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
777 c2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
778 d2 * (d2 * 0.5 + u2 - v2 - w2) +
779 u2 * (u2 * 0.5 + v2 + w2) + v2 * (v2 * 0.5 + w2) +
780 4 * (b * d * u * w - b * c * v * w - c * d * u * v));
781 const Complex Const1 = std::sqrt(Complex(tmp, 0));
782 const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
783
784 const Complex x = std::sqrt(Const2 + Const1) * sqrt_05;
785 const Complex y = std::sqrt(Const2 - Const1) * sqrt_05 * Complex(0, 1);
786 const Complex x2 = x * x;
787 const Complex y2 = y * y;
788 const Complex cy = std::cos(y);
789 const Complex sy = std::sin(y);
790 const Complex cx = std::cosh(x);
791 const Complex sx = std::sinh(x);
792
793 const bool x_zero = std::abs(x) < lower_is_considered_zero_for_sinc_likes;
794 const bool y_zero = std::abs(y) < lower_is_considered_zero_for_sinc_likes;
795 const bool both_zero = y_zero and x_zero;
796 const bool either_zero = y_zero or x_zero;
797
798 /* Using:
799 * lim x→0 [({cosh(x),cos(x)} - 1) / x^2] → 1/2
800 * lim x→0 [{sinh(x),sin(x)} / x] → 1
801 * inv_x2 := 1 for x == 0,
802 * -i sin(ix) → sinh(x)
803 * cos(ix) → cosh(x)
804 * C0, C1, C2 ∝ [1/x^2]
805 */
806 const Complex ix = x_zero ? 0.0 : 1.0 / x;
807 const Complex iy = y_zero ? 0.0 : 1.0 / y;
808 const Complex inv_x2y2 =
809 both_zero
810 ? 1.0
811 : 1.0 /
812 (x2 + y2); // The first "1.0" is the trick for above limits
813
814 const Numeric C0 =
815 either_zero ? 1.0 : ((cy * x2 + cx * y2) * inv_x2y2).real();
816 const Numeric C1 =
817 either_zero ? 1.0 : ((sy * x2 * iy + sx * y2 * ix) * inv_x2y2).real();
818 const Numeric C2 = both_zero ? 0.5 : ((cx - cy) * inv_x2y2).real();
819 const Numeric C3 = both_zero ? 1.0 / 6.0
820 : ((x_zero ? 1.0 - sy * iy
821 : y_zero ? sx * ix - 1.0
822 : sx * ix - sy * iy) *
823 inv_x2y2)
824 .real();
825 T.Mat4(i).noalias() =
826 exp_a * (Eigen::Matrix4d() << C0 + C2 * (b2 + c2 + d2),
827 C1 * b + C2 * (-c * u - d * v) +
828 C3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) -
829 v * (b * v + c * w)),
830 C1 * c + C2 * (b * u - d * w) +
831 C3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) -
832 w * (b * v + c * w)),
833 C1 * d + C2 * (b * v + c * w) +
834 C3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) +
835 w * (b * u - d * w)),
836
837 C1 * b + C2 * (c * u + d * v) +
838 C3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) +
839 d * (b * d + u * w)),
840 C0 + C2 * (b2 - u2 - v2),
841 C2 * (b * c - v * w) + C1 * u +
842 C3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
843 w * (b * d + u * w)),
844 C2 * (b * d + u * w) + C1 * v +
845 C3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
846 w * (b * c - v * w)),
847
848 C1 * c + C2 * (-b * u + d * w) +
849 C3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) +
850 d * (c * d - u * v)),
851 C2 * (b * c - v * w) - C1 * u +
852 C3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
853 v * (c * d - u * v)),
854 C0 + C2 * (c2 - u2 - w2),
855 C2 * (c * d - u * v) + C1 * w +
856 C3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
857 w * (-c2 + u2 + w2)),
858
859 C1 * d + C2 * (-b * v - c * w) +
860 C3 * (b * (b * d + u * w) + c * (c * d - u * v) -
861 d * (-d2 + v2 + w2)),
862 C2 * (b * d + u * w) - C1 * v +
863 C3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
864 v * (-d2 + v2 + w2)),
865 C2 * (c * d - u * v) - C1 * w +
866 C3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
867 w * (-d2 + v2 + w2)),
868 C0 + C2 * (d2 - v2 - w2))
869 .finished();
870 }
871 }
872}
873
874inline void dtransmat1(TransmissionMatrix& T,
875 ArrayOfTransmissionMatrix& dT1,
876 ArrayOfTransmissionMatrix& dT2,
877 const PropagationMatrix& K1,
878 const PropagationMatrix& K2,
879 const ArrayOfPropagationMatrix& dK1,
880 const ArrayOfPropagationMatrix& dK2,
881 const Numeric& r,
882 const Numeric& dr_dT1,
883 const Numeric& dr_dT2,
884 const Index it,
885 const Index iz,
886 const Index ia) noexcept {
887 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
888 T.Mat1(i)(0, 0) =
889 std::exp(-0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]));
890 for (Index j = 0; j < dT1.nelem(); j++) {
891 if (dK1[j].NumberOfFrequencies())
892 dT1[j].Mat1(i)(0, 0) =
893 T.Mat1(i)(0, 0) *
894 (-0.5 *
895 (r * dK1[j].Kjj(iz, ia)[i] +
896 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
897 : 0.0)));
898 if (dK2[j].NumberOfFrequencies())
899 dT2[j].Mat1(i)(0, 0) =
900 T.Mat1(i)(0, 0) *
901 (-0.5 *
902 (r * dK2[j].Kjj(iz, ia)[i] +
903 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
904 : 0.0)));
905 }
906 }
907}
908
909inline void dtransmat2(TransmissionMatrix& T,
910 ArrayOfTransmissionMatrix& dT1,
911 ArrayOfTransmissionMatrix& dT2,
912 const PropagationMatrix& K1,
913 const PropagationMatrix& K2,
914 const ArrayOfPropagationMatrix& dK1,
915 const ArrayOfPropagationMatrix& dK2,
916 const Numeric& r,
917 const Numeric& dr_dT1,
918 const Numeric& dr_dT2,
919 const Index it,
920 const Index iz,
921 const Index ia) noexcept {
922 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
923 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
924 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]);
925 const Numeric exp_a = std::exp(a);
926 const Numeric cb = std::cosh(b), sb = std::sinh(b);
927 T.Mat2(i).noalias() =
928 (Eigen::Matrix2d() << cb, sb, sb, cb).finished() * exp_a;
929 for (Index j = 0; j < dT1.nelem(); j++) {
930 if (dK1[j].NumberOfFrequencies()) {
931 const Numeric da = -0.5 * (r * dK1[j].Kjj(iz, ia)[i] +
932 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] +
933 K2.Kjj(iz, ia)[i])
934 : 0.0)),
935 db = -0.5 * (r * dK1[j].K12(iz, ia)[i] +
936 ((j == it) ? dr_dT1 * (K1.K12(iz, ia)[i] +
937 K2.K12(iz, ia)[i])
938 : 0.0));
939 dT1[j].Mat2(i).noalias() =
940 T.Mat2(i) * da +
941 (Eigen::Matrix2d() << sb, cb, cb, sb).finished() * exp_a * db;
942 }
943 if (dK2[j].NumberOfFrequencies()) {
944 const Numeric da = -0.5 * (r * dK2[j].Kjj(iz, ia)[i] +
945 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] +
946 K2.Kjj(iz, ia)[i])
947 : 0.0)),
948 db = -0.5 * (r * dK2[j].K12(iz, ia)[i] +
949 ((j == it) ? dr_dT2 * (K1.K12(iz, ia)[i] +
950 K2.K12(iz, ia)[i])
951 : 0.0));
952 dT2[j].Mat2(i).noalias() =
953 T.Mat2(i) * da +
954 (Eigen::Matrix2d() << sb, cb, cb, sb).finished() * exp_a * db;
955 }
956 }
957 }
958}
959
960inline void dtransmat3(TransmissionMatrix& T,
961 ArrayOfTransmissionMatrix& dT1,
962 ArrayOfTransmissionMatrix& dT2,
963 const PropagationMatrix& K1,
964 const PropagationMatrix& K2,
965 const ArrayOfPropagationMatrix& dK1,
966 const ArrayOfPropagationMatrix& dK2,
967 const Numeric& r,
968 const Numeric& dr_dT1,
969 const Numeric& dr_dT2,
970 const Index it,
971 const Index iz,
972 const Index ia) noexcept {
973 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
974 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
975 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]),
976 c = -0.5 * r * (K1.K13(iz, ia)[i] + K2.K13(iz, ia)[i]),
977 u = -0.5 * r * (K1.K23(iz, ia)[i] + K2.K23(iz, ia)[i]);
978 const Numeric exp_a = std::exp(a);
979
980 if (b == 0. and c == 0. and u == 0.) {
981 T.Mat3(i).noalias() = Eigen::Matrix3d::Identity() * exp_a;
982 for (Index j = 0; j < dT1.nelem(); j++) {
983 if (dK1[j].NumberOfFrequencies())
984 dT1[j].Mat3(i).noalias() =
985 T.Mat3(i) *
986 (-0.5 *
987 (r * dK1[j].Kjj(iz, ia)[i] +
988 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
989 : 0.0)));
990 if (dK2[j].NumberOfFrequencies())
991 dT2[j].Mat3(i).noalias() =
992 T.Mat3(i) *
993 (-0.5 *
994 (r * dK2[j].Kjj(iz, ia)[i] +
995 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
996 : 0.0)));
997 }
998 } else {
999 const Numeric a2 = a * a, b2 = b * b, c2 = c * c, u2 = u * u;
1000 const Numeric Const = b2 + c2 - u2;
1001
1002 const bool real = Const > 0;
1003 const bool imag = Const < 0;
1004 const bool either = real or imag;
1005
1006 const Numeric x =
1007 std::sqrt(imag ? -Const : Const); // test to just use real values
1008 const Numeric x2 =
1009 (real ? 1 : -1) * x * x; // test to change sign if imaginary
1010 const Numeric inv_x =
1011 either ? 1.0 / x
1012 : 1.0; // test so further calculations are replaced as x→0
1013 const Numeric inv_x2 = inv_x * inv_x;
1014
1015 const Numeric sx =
1016 real ? std::sinh(x) : std::sin(-x); // -i sin(ix) → sinh(x)
1017 const Numeric cx =
1018 real ? std::cosh(x) : std::cos(+x); // cos(ix) → cosh(x)
1019
1020 /* Using:
1021 * lim x→0 [(cosh(x) - 1) / x^2] → 1/2
1022 * lim x→0 [sinh(x) / x] → 1
1023 * inv_x2 := 1 for x == 0,
1024 * C0, C1, C2 ∝ [1/x^2]
1025 */
1026 const Numeric C0 =
1027 either ? a2 * (cx - 1.0) - a * x * sx + x2 : 1.0 + 0.5 * a2 - a;
1028 const Numeric C1 = either ? 2.0 * a * (1.0 - cx) + x * sx : 1.0 - a;
1029 const Numeric C2 = either ? cx - 1.0 : 0.5;
1030
1031 T.Mat3(i).noalias() =
1032 exp_a * inv_x2 *
1033 (Eigen::Matrix3d() << C0 + C1 * a + C2 * (a2 + b2 + c2),
1034 C1 * b + C2 * (2 * a * b - c * u),
1035 C1 * c + C2 * (2 * a * c + b * u),
1036 C1 * b + C2 * (2 * a * b + c * u),
1037 C0 + C1 * a + C2 * (a2 + b2 - u2),
1038 C1 * u + C2 * (2 * a * u + b * c),
1039 C1 * c + C2 * (2 * a * c - b * u),
1040 -C1 * u - C2 * (2 * a * u - b * c),
1041 C0 + C1 * a + C2 * (a2 + c2 - u2))
1042 .finished();
1043
1044 for (Index j = 0; j < dT1.nelem(); j++) {
1045 if (dK1[j].NumberOfFrequencies()) {
1046 const Numeric da = -0.5 * (r * dK1[j].Kjj(iz, ia)[i] +
1047 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] +
1048 K2.Kjj(iz, ia)[i])
1049 : 0.0)),
1050 db = -0.5 * (r * dK1[j].K12(iz, ia)[i] +
1051 ((j == it) ? dr_dT1 * (K1.K12(iz, ia)[i] +
1052 K2.K12(iz, ia)[i])
1053 : 0.0)),
1054 dc = -0.5 * (r * dK1[j].K13(iz, ia)[i] +
1055 ((j == it) ? dr_dT1 * (K1.K13(iz, ia)[i] +
1056 K2.K13(iz, ia)[i])
1057 : 0.0)),
1058 du = -0.5 * (r * dK1[j].K23(iz, ia)[i] +
1059 ((j == it) ? dr_dT1 * (K1.K23(iz, ia)[i] +
1060 K2.K23(iz, ia)[i])
1061 : 0.0));
1062 const Numeric da2 = 2 * a * da, db2 = 2 * b * db, dc2 = 2 * c * dc,
1063 du2 = 2 * u * du;
1064 const Numeric dx = either ? ((db2 + dc2 - du2) * inv_x * 0.5) : 0,
1065 dx2 = 2 * x * dx;
1066 const Numeric dsx = (real ? 1 : -1) * cx * dx;
1067 const Numeric dcx = sx * dx;
1068
1069 const Numeric dC0 = either ? da2 * (cx - 1.0) + da2 * (dcx - 1.0) -
1070 da * x * sx - a * dx * sx -
1071 a * x * dsx + dx2
1072 : 0.5 * da2 - da;
1073 const Numeric dC1 =
1074 either ? 2.0 * (da * (1.0 - cx) - a * dcx) + dx * sx + x * dsx
1075 : -da;
1076 const Numeric dC2 = either ? dcx : 0;
1077
1078 dT1[j].Mat3(i).noalias() =
1079 T.Mat3(i) * (da + dx2 * inv_x2) +
1080 exp_a * inv_x2 *
1081 (Eigen::Matrix3d() << dC0 + dC1 * a + C1 * da +
1082 dC2 * (a2 + b2 + c2) +
1083 C2 * (da2 + db2 + dc2),
1084 dC1 * b + C1 * db + dC2 * (2 * a * b - c * u) +
1085 C2 * (2 * da * b - dc * u + 2 * a * db - c * du),
1086 dC1 * c + C1 * dc + dC2 * (2 * a * c + b * u) +
1087 C2 * (2 * da * c + db * u + 2 * a * dc + b * du),
1088 dC1 * b + C1 * db + dC2 * (2 * a * b + c * u) +
1089 C2 * (2 * da * b + dc * u + 2 * a * db + c * du),
1090 dC0 + dC1 * a + C1 * da + dC2 * (a2 + b2 - u2) +
1091 C2 * (da2 + db2 - du2),
1092 dC1 * u + C1 * du + dC2 * (2 * a * u + b * c) +
1093 C2 * (2 * da * u + db * c + 2 * a * du + b * dc),
1094 dC1 * c + C1 * dc + dC2 * (2 * a * c - b * u) +
1095 C2 * (2 * da * c - db * u + 2 * a * dc - b * du),
1096 -dC1 * u - C1 * du - dC2 * (2 * a * u - b * c) -
1097 C2 * (2 * da * u - db * c + 2 * a * du - b * dc),
1098 dC0 + dC1 * a + C1 * da + dC2 * (a2 + c2 - u2) +
1099 C2 * (da2 + dc2 - du2))
1100 .finished();
1101 }
1102 if (dK2[j].NumberOfFrequencies()) {
1103 const Numeric da = -0.5 * (r * dK2[j].Kjj(iz, ia)[i] +
1104 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] +
1105 K2.Kjj(iz, ia)[i])
1106 : 0.0)),
1107 db = -0.5 * (r * dK2[j].K12(iz, ia)[i] +
1108 ((j == it) ? dr_dT2 * (K1.K12(iz, ia)[i] +
1109 K2.K12(iz, ia)[i])
1110 : 0.0)),
1111 dc = -0.5 * (r * dK2[j].K13(iz, ia)[i] +
1112 ((j == it) ? dr_dT2 * (K1.K13(iz, ia)[i] +
1113 K2.K13(iz, ia)[i])
1114 : 0.0)),
1115 du = -0.5 * (r * dK2[j].K23(iz, ia)[i] +
1116 ((j == it) ? dr_dT2 * (K1.K23(iz, ia)[i] +
1117 K2.K23(iz, ia)[i])
1118 : 0.0));
1119 const Numeric da2 = 2 * a * da, db2 = 2 * b * db, dc2 = 2 * c * dc,
1120 du2 = 2 * u * du;
1121 const Numeric dx = either ? ((db2 + dc2 - du2) * inv_x * 0.5) : 0,
1122 dx2 = 2 * x * dx;
1123 const Numeric dsx = (real ? 1 : -1) * cx * dx;
1124 const Numeric dcx = sx * dx;
1125
1126 const Numeric dC0 = either ? da2 * (cx - 1.0) + da2 * (dcx - 1.0) -
1127 da * x * sx - a * dx * sx -
1128 a * x * dsx + dx2
1129 : 0.5 * da2 - da;
1130 const Numeric dC1 =
1131 either ? 2.0 * (da * (1.0 - cx) - a * dcx) + dx * sx + x * dsx
1132 : -da;
1133 const Numeric dC2 = either ? dcx : 0;
1134
1135 dT2[j].Mat3(i).noalias() =
1136 T.Mat3(i) * (da + dx2 * inv_x2) +
1137 exp_a * inv_x2 *
1138 (Eigen::Matrix3d() << dC0 + dC1 * a + C1 * da +
1139 dC2 * (a2 + b2 + c2) +
1140 C2 * (da2 + db2 + dc2),
1141 dC1 * b + C1 * db + dC2 * (2 * a * b - c * u) +
1142 C2 * (2 * da * b - dc * u + 2 * a * db - c * du),
1143 dC1 * c + C1 * dc + dC2 * (2 * a * c + b * u) +
1144 C2 * (2 * da * c + db * u + 2 * a * dc + b * du),
1145 dC1 * b + C1 * db + dC2 * (2 * a * b + c * u) +
1146 C2 * (2 * da * b + dc * u + 2 * a * db + c * du),
1147 dC0 + dC1 * a + C1 * da + dC2 * (a2 + b2 - u2) +
1148 C2 * (da2 + db2 - du2),
1149 dC1 * u + C1 * du + dC2 * (2 * a * u + b * c) +
1150 C2 * (2 * da * u + db * c + 2 * a * du + b * dc),
1151 dC1 * c + C1 * dc + dC2 * (2 * a * c - b * u) +
1152 C2 * (2 * da * c - db * u + 2 * a * dc - b * du),
1153 -dC1 * u - C1 * du - dC2 * (2 * a * u - b * c) -
1154 C2 * (2 * da * u - db * c + 2 * a * du - b * dc),
1155 dC0 + dC1 * a + C1 * da + dC2 * (a2 + c2 - u2) +
1156 C2 * (da2 + dc2 - du2))
1157 .finished();
1158 }
1159 }
1160 }
1161 }
1162}
1163
1164inline void dtransmat4(TransmissionMatrix& T,
1165 ArrayOfTransmissionMatrix& dT1,
1166 ArrayOfTransmissionMatrix& dT2,
1167 const PropagationMatrix& K1,
1168 const PropagationMatrix& K2,
1169 const ArrayOfPropagationMatrix& dK1,
1170 const ArrayOfPropagationMatrix& dK2,
1171 const Numeric& r,
1172 const Numeric& dr_dT1,
1173 const Numeric& dr_dT2,
1174 const Index it,
1175 const Index iz,
1176 const Index ia) noexcept {
1177 static constexpr Numeric sqrt_05 = Constant::inv_sqrt_2;
1178 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
1179 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
1180 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]),
1181 c = -0.5 * r * (K1.K13(iz, ia)[i] + K2.K13(iz, ia)[i]),
1182 d = -0.5 * r * (K1.K14(iz, ia)[i] + K2.K14(iz, ia)[i]),
1183 u = -0.5 * r * (K1.K23(iz, ia)[i] + K2.K23(iz, ia)[i]),
1184 v = -0.5 * r * (K1.K24(iz, ia)[i] + K2.K24(iz, ia)[i]),
1185 w = -0.5 * r * (K1.K34(iz, ia)[i] + K2.K34(iz, ia)[i]);
1186 const Numeric exp_a = std::exp(a);
1187
1188 if (b == 0. and c == 0. and d == 0. and u == 0. and v == 0. and w == 0.) {
1189 T.Mat4(i).noalias() = Eigen::Matrix4d::Identity() * exp_a;
1190 for (Index j = 0; j < dK1.nelem(); j++) {
1191 if (dK1[j].NumberOfFrequencies())
1192 dT1[j].Mat4(i).noalias() =
1193 T.Mat4(i) *
1194 (-0.5 *
1195 (r * dK1[j].Kjj(iz, ia)[i] +
1196 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
1197 : 0.0)));
1198 if (dK2[j].NumberOfFrequencies())
1199 dT2[j].Mat4(i).noalias() =
1200 T.Mat4(i) *
1201 (-0.5 *
1202 (r * dK2[j].Kjj(iz, ia)[i] +
1203 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
1204 : 0.0)));
1205 }
1206 } else {
1207 const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
1208 w2 = w * w;
1209 const Numeric tmp =
1210 w2 * w2 + 2 * (b2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
1211 c2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
1212 d2 * (d2 * 0.5 + u2 - v2 - w2) +
1213 u2 * (u2 * 0.5 + v2 + w2) + v2 * (v2 * 0.5 + w2) +
1214 4 * (b * d * u * w - b * c * v * w - c * d * u * v));
1215 const Complex Const1 = std::sqrt(Complex(tmp, 0));
1216 const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
1217 const Complex tmp_x_sqrt = std::sqrt(Const2 + Const1);
1218 const Complex tmp_y_sqrt = std::sqrt(Const2 - Const1);
1219 const Complex x = tmp_x_sqrt * sqrt_05;
1220 const Complex y = tmp_y_sqrt * sqrt_05 * Complex(0, 1);
1221 const Complex x2 = x * x;
1222 const Complex y2 = y * y;
1223 const Complex cy = std::cos(y);
1224 const Complex sy = std::sin(y);
1225 const Complex cx = std::cosh(x);
1226 const Complex sx = std::sinh(x);
1227
1228 const bool x_zero = std::abs(x) < lower_is_considered_zero_for_sinc_likes;
1229 const bool y_zero = std::abs(y) < lower_is_considered_zero_for_sinc_likes;
1230 const bool both_zero = y_zero and x_zero;
1231 const bool either_zero = y_zero or x_zero;
1232
1233 /* Using:
1234 * lim x→0 [({cosh(x),cos(x)} - 1) / x^2] → 1/2
1235 * lim x→0 [{sinh(x),sin(x)} / x] → 1
1236 * inv_x2 := 1 for x == 0,
1237 * -i sin(ix) → sinh(x)
1238 * cos(ix) → cosh(x)
1239 * C0, C1, C2 ∝ [1/x^2]
1240 */
1241 const Complex ix = x_zero ? 0.0 : 1.0 / x;
1242 const Complex iy = y_zero ? 0.0 : 1.0 / y;
1243 const Complex inv_x2y2 =
1244 both_zero
1245 ? 1.0
1246 : 1.0 /
1247 (x2 + y2); // The first "1.0" is the trick for above limits
1248 const Complex C0c = either_zero ? 1.0 : (cy * x2 + cx * y2) * inv_x2y2;
1249 const Complex C1c =
1250 either_zero ? 1.0 : (sy * x2 * iy + sx * y2 * ix) * inv_x2y2;
1251 const Complex C2c = both_zero ? 0.5 : (cx - cy) * inv_x2y2;
1252 const Complex C3c = both_zero ? 1.0 / 6.0
1253 : (x_zero ? 1.0 - sy * iy
1254 : y_zero ? sx * ix - 1.0
1255 : sx * ix - sy * iy) *
1256 inv_x2y2;
1257
1258 const Numeric& C0 = real_val(C0c);
1259 const Numeric& C1 = real_val(C1c);
1260 const Numeric& C2 = real_val(C2c);
1261 const Numeric& C3 = real_val(C3c);
1262 T.Mat4(i).noalias() =
1263 exp_a * (Eigen::Matrix4d() << C0 + C2 * (b2 + c2 + d2),
1264 C1 * b + C2 * (-c * u - d * v) +
1265 C3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) -
1266 v * (b * v + c * w)),
1267 C1 * c + C2 * (b * u - d * w) +
1268 C3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) -
1269 w * (b * v + c * w)),
1270 C1 * d + C2 * (b * v + c * w) +
1271 C3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) +
1272 w * (b * u - d * w)),
1273
1274 C1 * b + C2 * (c * u + d * v) +
1275 C3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) +
1276 d * (b * d + u * w)),
1277 C0 + C2 * (b2 - u2 - v2),
1278 C2 * (b * c - v * w) + C1 * u +
1279 C3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
1280 w * (b * d + u * w)),
1281 C2 * (b * d + u * w) + C1 * v +
1282 C3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
1283 w * (b * c - v * w)),
1284
1285 C1 * c + C2 * (-b * u + d * w) +
1286 C3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) +
1287 d * (c * d - u * v)),
1288 C2 * (b * c - v * w) - C1 * u +
1289 C3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
1290 v * (c * d - u * v)),
1291 C0 + C2 * (c2 - u2 - w2),
1292 C2 * (c * d - u * v) + C1 * w +
1293 C3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
1294 w * (-c2 + u2 + w2)),
1295
1296 C1 * d + C2 * (-b * v - c * w) +
1297 C3 * (b * (b * d + u * w) + c * (c * d - u * v) -
1298 d * (-d2 + v2 + w2)),
1299 C2 * (b * d + u * w) - C1 * v +
1300 C3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
1301 v * (-d2 + v2 + w2)),
1302 C2 * (c * d - u * v) - C1 * w +
1303 C3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
1304 w * (-d2 + v2 + w2)),
1305 C0 + C2 * (d2 - v2 - w2))
1306 .finished();
1307
1308 for (Index j = 0; j < dK1.nelem(); j++) {
1309 if (dK1[j].NumberOfFrequencies()) {
1310 const Numeric da = -0.5 * (r * dK1[j].Kjj(iz, ia)[i] +
1311 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] +
1312 K2.Kjj(iz, ia)[i])
1313 : 0.0)),
1314 db = -0.5 * (r * dK1[j].K12(iz, ia)[i] +
1315 ((j == it) ? dr_dT1 * (K1.K12(iz, ia)[i] +
1316 K2.K12(iz, ia)[i])
1317 : 0.0)),
1318 dc = -0.5 * (r * dK1[j].K13(iz, ia)[i] +
1319 ((j == it) ? dr_dT1 * (K1.K13(iz, ia)[i] +
1320 K2.K13(iz, ia)[i])
1321 : 0.0)),
1322 dd = -0.5 * (r * dK1[j].K14(iz, ia)[i] +
1323 ((j == it) ? dr_dT1 * (K1.K14(iz, ia)[i] +
1324 K2.K14(iz, ia)[i])
1325 : 0.0)),
1326 du = -0.5 * (r * dK1[j].K23(iz, ia)[i] +
1327 ((j == it) ? dr_dT1 * (K1.K23(iz, ia)[i] +
1328 K2.K23(iz, ia)[i])
1329 : 0.0)),
1330 dv = -0.5 * (r * dK1[j].K24(iz, ia)[i] +
1331 ((j == it) ? dr_dT1 * (K1.K24(iz, ia)[i] +
1332 K2.K24(iz, ia)[i])
1333 : 0.0)),
1334 dw = -0.5 * (r * dK1[j].K34(iz, ia)[i] +
1335 ((j == it) ? dr_dT1 * (K1.K34(iz, ia)[i] +
1336 K2.K34(iz, ia)[i])
1337 : 0.0));
1338 const Numeric db2 = 2 * db * b, dc2 = 2 * dc * c, dd2 = 2 * dd * d,
1339 du2 = 2 * du * u, dv2 = 2 * dv * v, dw2 = 2 * dw * w;
1340 const Numeric dtmp =
1341 2 * w2 * dw2 +
1342 2 * (db2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
1343 b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2) +
1344 dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
1345 c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2) +
1346 dd2 * (d2 * 0.5 + u2 - v2 - w2) +
1347 d2 * (dd2 * 0.5 + du2 - dv2 - dw2) +
1348 du2 * (u2 * 0.5 + v2 + w2) + u2 * (du2 * 0.5 + dv2 + dw2) +
1349 dv2 * (v2 * 0.5 + w2) + v2 * (dv2 * 0.5 + dw2) +
1350 4 * (db * d * u * w - db * c * v * w - dc * d * u * v +
1351 b * dd * u * w - b * dc * v * w - c * dd * u * v +
1352 b * d * du * w - b * c * dv * w - c * d * du * v +
1353 b * d * u * dw - b * c * v * dw - c * d * u * dv));
1354 const Complex dConst1 = 0.5 * dtmp / Const1;
1355 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1356 const Complex dx =
1357 x_zero ? 0 : (0.5 * (dConst2 + dConst1) / tmp_x_sqrt) * sqrt_05;
1358 const Complex dy = y_zero ? 0
1359 : (0.5 * (dConst2 - dConst1) / tmp_y_sqrt) *
1360 sqrt_05 * Complex(0, 1);
1361 const Complex dx2 = 2 * x * dx;
1362 const Complex dy2 = 2 * y * dy;
1363 const Complex dcy = -sy * dy;
1364 const Complex dsy = cy * dy;
1365 const Complex dcx = sx * dx;
1366 const Complex dsx = cx * dx;
1367 const Complex dix = -dx * ix * ix;
1368 const Complex diy = -dy * iy * iy;
1369 const Complex dx2dy2 = dx2 + dy2;
1370 const Complex dC0c =
1371 either_zero
1372 ? 0.0
1373 : (dcy * x2 + cy * dx2 + dcx * y2 + cx * dy2 - C0c * dx2dy2) *
1374 inv_x2y2;
1375 const Complex dC1c =
1376 either_zero ? 0.0
1377 : (dsy * x2 * iy + sy * dx2 * iy + sy * x2 * diy +
1378 dsx * y2 * ix + sx * dy2 * ix + sx * y2 * dix -
1379 C1c * dx2dy2) *
1380 inv_x2y2;
1381 const Complex dC2c =
1382 both_zero ? 0.0 : (dcx - dcy - C2c * dx2dy2) * inv_x2y2;
1383 const Complex dC3c =
1384 both_zero
1385 ? 0.0
1386 : ((x_zero ? -dsy * iy - sy * diy
1387 : y_zero ? dsx * ix + sx * dix
1388 : dsx * ix + sx * dix - dsy * iy - sy * diy) -
1389 C3c * dx2dy2) *
1390 inv_x2y2;
1391
1392 const Numeric& dC0 = real_val(dC0c);
1393 const Numeric& dC1 = real_val(dC1c);
1394 const Numeric& dC2 = real_val(dC2c);
1395 const Numeric& dC3 = real_val(dC3c);
1396 dT1[j].Mat4(i).noalias() =
1397 T.Mat4(i) * da +
1398 exp_a *
1399 (Eigen::Matrix4d()
1400 << dC0 + dC2 * (b2 + c2 + d2) + C2 * (db2 + dc2 + dd2),
1401 db * C1 + b * dC1 + dC2 * (-c * u - d * v) +
1402 C2 * (-dc * u - dd * v - c * du - d * dv) +
1403 dC3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) -
1404 v * (b * v + c * w)) +
1405 C3 * (db * (b2 + c2 + d2) - du * (b * u - d * w) -
1406 dv * (b * v + c * w) + b * (db2 + dc2 + dd2) -
1407 u * (db * u - dd * w) - v * (db * v + dc * w) -
1408 u * (b * du - d * dw) - v * (b * dv + c * dw)),
1409 dC1 * c + C1 * dc + dC2 * (b * u - d * w) +
1410 C2 * (db * u - dd * w + b * du - d * dw) +
1411 dC3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) -
1412 w * (b * v + c * w)) +
1413 C3 * (dc * (b2 + c2 + d2) - du * (c * u + d * v) -
1414 dw * (b * v + c * w) + c * (db2 + dc2 + dd2) -
1415 u * (dc * u + dd * v) - w * (db * v + dc * w) -
1416 u * (c * du + d * dv) - w * (b * dv + c * dw)),
1417 dC1 * d + C1 * dd + dC2 * (b * v + c * w) +
1418 C2 * (db * v + dc * w + b * dv + c * dw) +
1419 dC3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) +
1420 w * (b * u - d * w)) +
1421 C3 * (dd * (b2 + c2 + d2) - dv * (c * u + d * v) +
1422 dw * (b * u - d * w) + d * (db2 + dc2 + dd2) -
1423 v * (dc * u + dd * v) + w * (db * u - dd * w) -
1424 v * (c * du + d * dv) + w * (b * du - d * dw)),
1425
1426 db * C1 + b * dC1 + dC2 * (c * u + d * v) +
1427 C2 * (dc * u + dd * v + c * du + d * dv) +
1428 dC3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) +
1429 d * (b * d + u * w)) +
1430 C3 * (-db * (-b2 + u2 + v2) + dc * (b * c - v * w) +
1431 dd * (b * d + u * w) - b * (-db2 + du2 + dv2) +
1432 c * (db * c - dv * w) + d * (db * d + du * w) +
1433 c * (b * dc - v * dw) + d * (b * dd + u * dw)),
1434 dC0 + dC2 * (b2 - u2 - v2) + C2 * (db2 - du2 - dv2),
1435 dC2 * (b * c - v * w) +
1436 C2 * (db * c + b * dc - dv * w - v * dw) + dC1 * u +
1437 C1 * du +
1438 dC3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
1439 w * (b * d + u * w)) +
1440 C3 * (dc * (c * u + d * v) - du * (-b2 + u2 + v2) -
1441 dw * (b * d + u * w) + c * (dc * u + dd * v) -
1442 u * (-db2 + du2 + dv2) - w * (db * d + du * w) +
1443 c * (c * du + d * dv) - w * (b * dd + u * dw)),
1444 dC2 * (b * d + u * w) +
1445 C2 * (db * d + b * dd + du * w + u * dw) + dC1 * v +
1446 C1 * dv +
1447 dC3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
1448 w * (b * c - v * w)) +
1449 C3 * (dd * (c * u + d * v) - dv * (-b2 + u2 + v2) +
1450 dw * (b * c - v * w) + d * (dc * u + dd * v) -
1451 v * (-db2 + du2 + dv2) + w * (db * c - dv * w) +
1452 d * (c * du + d * dv) + w * (b * dc - v * dw)),
1453
1454 dC1 * c + C1 * dc + dC2 * (-b * u + d * w) +
1455 C2 * (-db * u + dd * w - b * du + d * dw) +
1456 dC3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) +
1457 d * (c * d - u * v)) +
1458 C3 * (db * (b * c - v * w) - dc * (-c2 + u2 + w2) +
1459 dd * (c * d - u * v) + b * (db * c - dv * w) -
1460 c * (-dc2 + du2 + dw2) + d * (dc * d - du * v) +
1461 b * (b * dc - v * dw) + d * (c * dd - u * dv)),
1462 dC2 * (b * c - v * w) +
1463 C2 * (db * c + b * dc - dv * w - v * dw) - dC1 * u -
1464 C1 * du +
1465 dC3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
1466 v * (c * d - u * v)) +
1467 C3 * (-db * (b * u - d * w) + du * (-c2 + u2 + w2) -
1468 dv * (c * d - u * v) - b * (db * u - dd * w) +
1469 u * (-dc2 + du2 + dw2) - v * (dc * d - du * v) -
1470 b * (b * du - d * dw) - v * (c * dd - u * dv)),
1471 dC0 + dC2 * (c2 - u2 - w2) + C2 * (dc2 - du2 - dw2),
1472 dC2 * (c * d - u * v) +
1473 C2 * (dc * d + c * dd - du * v - u * dv) + dC1 * w +
1474 C1 * dw +
1475 dC3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
1476 w * (-c2 + u2 + w2)) +
1477 C3 * (-dd * (b * u - d * w) + dv * (b * c - v * w) -
1478 dw * (-c2 + u2 + w2) - d * (db * u - dd * w) +
1479 v * (db * c - dv * w) - w * (-dc2 + du2 + dw2) -
1480 d * (b * du - d * dw) + v * (b * dc - v * dw)),
1481
1482 dC1 * d + C1 * dd + dC2 * (-b * v - c * w) +
1483 C2 * (-db * v - dc * w - b * dv - c * dw) +
1484 dC3 * (b * (b * d + u * w) + c * (c * d - u * v) -
1485 d * (-d2 + v2 + w2)) +
1486 C3 * (db * (b * d + u * w) + dc * (c * d - u * v) -
1487 dd * (-d2 + v2 + w2) + b * (db * d + du * w) +
1488 c * (dc * d - du * v) - d * (-dd2 + dv2 + dw2) +
1489 b * (b * dd + u * dw) + c * (c * dd - u * dv)),
1490 dC2 * (b * d + u * w) +
1491 C2 * (db * d + b * dd + du * w + u * dw) - dC1 * v -
1492 C1 * dv +
1493 dC3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
1494 v * (-d2 + v2 + w2)) +
1495 C3 * (-db * (b * v + c * w) - du * (c * d - u * v) +
1496 dv * (-d2 + v2 + w2) - b * (db * v + dc * w) -
1497 u * (dc * d - du * v) + v * (-dd2 + dv2 + dw2) -
1498 b * (b * dv + c * dw) - u * (c * dd - u * dv)),
1499 dC2 * (c * d - u * v) +
1500 C2 * (dc * d + c * dd - du * v - u * dv) - dC1 * w -
1501 C1 * dw +
1502 dC3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
1503 w * (-d2 + v2 + w2)) +
1504 C3 * (-dc * (b * v + c * w) + du * (b * d + u * w) +
1505 dw * (-d2 + v2 + w2) - c * (db * v + dc * w) +
1506 u * (db * d + du * w) + w * (-dd2 + dv2 + dw2) -
1507 c * (b * dv + c * dw) + u * (b * dd + u * dw)),
1508 dC0 + dC2 * (d2 - v2 - w2) + C2 * (dd2 - dv2 - dw2))
1509 .finished();
1510 }
1511 if (dK2[j].NumberOfFrequencies()) {
1512 const Numeric da = -0.5 * (r * dK2[j].Kjj(iz, ia)[i] +
1513 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] +
1514 K2.Kjj(iz, ia)[i])
1515 : 0.0)),
1516 db = -0.5 * (r * dK2[j].K12(iz, ia)[i] +
1517 ((j == it) ? dr_dT2 * (K1.K12(iz, ia)[i] +
1518 K2.K12(iz, ia)[i])
1519 : 0.0)),
1520 dc = -0.5 * (r * dK2[j].K13(iz, ia)[i] +
1521 ((j == it) ? dr_dT2 * (K1.K13(iz, ia)[i] +
1522 K2.K13(iz, ia)[i])
1523 : 0.0)),
1524 dd = -0.5 * (r * dK2[j].K14(iz, ia)[i] +
1525 ((j == it) ? dr_dT2 * (K1.K14(iz, ia)[i] +
1526 K2.K14(iz, ia)[i])
1527 : 0.0)),
1528 du = -0.5 * (r * dK2[j].K23(iz, ia)[i] +
1529 ((j == it) ? dr_dT2 * (K1.K23(iz, ia)[i] +
1530 K2.K23(iz, ia)[i])
1531 : 0.0)),
1532 dv = -0.5 * (r * dK2[j].K24(iz, ia)[i] +
1533 ((j == it) ? dr_dT2 * (K1.K24(iz, ia)[i] +
1534 K2.K24(iz, ia)[i])
1535 : 0.0)),
1536 dw = -0.5 * (r * dK2[j].K34(iz, ia)[i] +
1537 ((j == it) ? dr_dT2 * (K1.K34(iz, ia)[i] +
1538 K2.K34(iz, ia)[i])
1539 : 0.0));
1540 const Numeric db2 = 2 * db * b, dc2 = 2 * dc * c, dd2 = 2 * dd * d,
1541 du2 = 2 * du * u, dv2 = 2 * dv * v, dw2 = 2 * dw * w;
1542 const Numeric dtmp =
1543 2 * w2 * dw2 +
1544 2 * (db2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
1545 b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2) +
1546 dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
1547 c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2) +
1548 dd2 * (d2 * 0.5 + u2 - v2 - w2) +
1549 d2 * (dd2 * 0.5 + du2 - dv2 - dw2) +
1550 du2 * (u2 * 0.5 + v2 + w2) + u2 * (du2 * 0.5 + dv2 + dw2) +
1551 dv2 * (v2 * 0.5 + w2) + v2 * (dv2 * 0.5 + dw2) +
1552 4 * (db * d * u * w - db * c * v * w - dc * d * u * v +
1553 b * dd * u * w - b * dc * v * w - c * dd * u * v +
1554 b * d * du * w - b * c * dv * w - c * d * du * v +
1555 b * d * u * dw - b * c * v * dw - c * d * u * dv));
1556 const Complex dConst1 = 0.5 * dtmp / Const1;
1557 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1558 const Complex dx =
1559 x_zero ? 0 : (0.5 * (dConst2 + dConst1) / tmp_x_sqrt) * sqrt_05;
1560 const Complex dy = y_zero ? 0
1561 : (0.5 * (dConst2 - dConst1) / tmp_y_sqrt) *
1562 sqrt_05 * Complex(0, 1);
1563 const Complex dx2 = 2 * x * dx;
1564 const Complex dy2 = 2 * y * dy;
1565 const Complex dcy = -sy * dy;
1566 const Complex dsy = cy * dy;
1567 const Complex dcx = sx * dx;
1568 const Complex dsx = cx * dx;
1569 const Complex dix = -dx * ix * ix;
1570 const Complex diy = -dy * iy * iy;
1571 const Complex dx2dy2 = dx2 + dy2;
1572 const Complex dC0c =
1573 either_zero
1574 ? 0.0
1575 : (dcy * x2 + cy * dx2 + dcx * y2 + cx * dy2 - C0c * dx2dy2) *
1576 inv_x2y2;
1577 const Complex dC1c =
1578 either_zero ? 0.0
1579 : (dsy * x2 * iy + sy * dx2 * iy + sy * x2 * diy +
1580 dsx * y2 * ix + sx * dy2 * ix + sx * y2 * dix -
1581 C1c * dx2dy2) *
1582 inv_x2y2;
1583 const Complex dC2c =
1584 both_zero ? 0.0 : (dcx - dcy - C2c * dx2dy2) * inv_x2y2;
1585 const Complex dC3c =
1586 both_zero
1587 ? 0.0
1588 : ((x_zero ? -dsy * iy - sy * diy
1589 : y_zero ? dsx * ix + sx * dix
1590 : dsx * ix + sx * dix - dsy * iy - sy * diy) -
1591 C3c * dx2dy2) *
1592 inv_x2y2;
1593
1594 const Numeric& dC0 = real_val(dC0c);
1595 const Numeric& dC1 = real_val(dC1c);
1596 const Numeric& dC2 = real_val(dC2c);
1597 const Numeric& dC3 = real_val(dC3c);
1598 dT2[j].Mat4(i).noalias() =
1599 T.Mat4(i) * da +
1600 exp_a *
1601 (Eigen::Matrix4d()
1602 << dC0 + dC2 * (b2 + c2 + d2) + C2 * (db2 + dc2 + dd2),
1603 db * C1 + b * dC1 + dC2 * (-c * u - d * v) +
1604 C2 * (-dc * u - dd * v - c * du - d * dv) +
1605 dC3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) -
1606 v * (b * v + c * w)) +
1607 C3 * (db * (b2 + c2 + d2) - du * (b * u - d * w) -
1608 dv * (b * v + c * w) + b * (db2 + dc2 + dd2) -
1609 u * (db * u - dd * w) - v * (db * v + dc * w) -
1610 u * (b * du - d * dw) - v * (b * dv + c * dw)),
1611 dC1 * c + C1 * dc + dC2 * (b * u - d * w) +
1612 C2 * (db * u - dd * w + b * du - d * dw) +
1613 dC3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) -
1614 w * (b * v + c * w)) +
1615 C3 * (dc * (b2 + c2 + d2) - du * (c * u + d * v) -
1616 dw * (b * v + c * w) + c * (db2 + dc2 + dd2) -
1617 u * (dc * u + dd * v) - w * (db * v + dc * w) -
1618 u * (c * du + d * dv) - w * (b * dv + c * dw)),
1619 dC1 * d + C1 * dd + dC2 * (b * v + c * w) +
1620 C2 * (db * v + dc * w + b * dv + c * dw) +
1621 dC3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) +
1622 w * (b * u - d * w)) +
1623 C3 * (dd * (b2 + c2 + d2) - dv * (c * u + d * v) +
1624 dw * (b * u - d * w) + d * (db2 + dc2 + dd2) -
1625 v * (dc * u + dd * v) + w * (db * u - dd * w) -
1626 v * (c * du + d * dv) + w * (b * du - d * dw)),
1627
1628 db * C1 + b * dC1 + dC2 * (c * u + d * v) +
1629 C2 * (dc * u + dd * v + c * du + d * dv) +
1630 dC3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) +
1631 d * (b * d + u * w)) +
1632 C3 * (-db * (-b2 + u2 + v2) + dc * (b * c - v * w) +
1633 dd * (b * d + u * w) - b * (-db2 + du2 + dv2) +
1634 c * (db * c - dv * w) + d * (db * d + du * w) +
1635 c * (b * dc - v * dw) + d * (b * dd + u * dw)),
1636 dC0 + dC2 * (b2 - u2 - v2) + C2 * (db2 - du2 - dv2),
1637 dC2 * (b * c - v * w) +
1638 C2 * (db * c + b * dc - dv * w - v * dw) + dC1 * u +
1639 C1 * du +
1640 dC3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
1641 w * (b * d + u * w)) +
1642 C3 * (dc * (c * u + d * v) - du * (-b2 + u2 + v2) -
1643 dw * (b * d + u * w) + c * (dc * u + dd * v) -
1644 u * (-db2 + du2 + dv2) - w * (db * d + du * w) +
1645 c * (c * du + d * dv) - w * (b * dd + u * dw)),
1646 dC2 * (b * d + u * w) +
1647 C2 * (db * d + b * dd + du * w + u * dw) + dC1 * v +
1648 C1 * dv +
1649 dC3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
1650 w * (b * c - v * w)) +
1651 C3 * (dd * (c * u + d * v) - dv * (-b2 + u2 + v2) +
1652 dw * (b * c - v * w) + d * (dc * u + dd * v) -
1653 v * (-db2 + du2 + dv2) + w * (db * c - dv * w) +
1654 d * (c * du + d * dv) + w * (b * dc - v * dw)),
1655
1656 dC1 * c + C1 * dc + dC2 * (-b * u + d * w) +
1657 C2 * (-db * u + dd * w - b * du + d * dw) +
1658 dC3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) +
1659 d * (c * d - u * v)) +
1660 C3 * (db * (b * c - v * w) - dc * (-c2 + u2 + w2) +
1661 dd * (c * d - u * v) + b * (db * c - dv * w) -
1662 c * (-dc2 + du2 + dw2) + d * (dc * d - du * v) +
1663 b * (b * dc - v * dw) + d * (c * dd - u * dv)),
1664 dC2 * (b * c - v * w) +
1665 C2 * (db * c + b * dc - dv * w - v * dw) - dC1 * u -
1666 C1 * du +
1667 dC3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
1668 v * (c * d - u * v)) +
1669 C3 * (-db * (b * u - d * w) + du * (-c2 + u2 + w2) -
1670 dv * (c * d - u * v) - b * (db * u - dd * w) +
1671 u * (-dc2 + du2 + dw2) - v * (dc * d - du * v) -
1672 b * (b * du - d * dw) - v * (c * dd - u * dv)),
1673 dC0 + dC2 * (c2 - u2 - w2) + C2 * (dc2 - du2 - dw2),
1674 dC2 * (c * d - u * v) +
1675 C2 * (dc * d + c * dd - du * v - u * dv) + dC1 * w +
1676 C1 * dw +
1677 dC3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
1678 w * (-c2 + u2 + w2)) +
1679 C3 * (-dd * (b * u - d * w) + dv * (b * c - v * w) -
1680 dw * (-c2 + u2 + w2) - d * (db * u - dd * w) +
1681 v * (db * c - dv * w) - w * (-dc2 + du2 + dw2) -
1682 d * (b * du - d * dw) + v * (b * dc - v * dw)),
1683
1684 dC1 * d + C1 * dd + dC2 * (-b * v - c * w) +
1685 C2 * (-db * v - dc * w - b * dv - c * dw) +
1686 dC3 * (b * (b * d + u * w) + c * (c * d - u * v) -
1687 d * (-d2 + v2 + w2)) +
1688 C3 * (db * (b * d + u * w) + dc * (c * d - u * v) -
1689 dd * (-d2 + v2 + w2) + b * (db * d + du * w) +
1690 c * (dc * d - du * v) - d * (-dd2 + dv2 + dw2) +
1691 b * (b * dd + u * dw) + c * (c * dd - u * dv)),
1692 dC2 * (b * d + u * w) +
1693 C2 * (db * d + b * dd + du * w + u * dw) - dC1 * v -
1694 C1 * dv +
1695 dC3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
1696 v * (-d2 + v2 + w2)) +
1697 C3 * (-db * (b * v + c * w) - du * (c * d - u * v) +
1698 dv * (-d2 + v2 + w2) - b * (db * v + dc * w) -
1699 u * (dc * d - du * v) + v * (-dd2 + dv2 + dw2) -
1700 b * (b * dv + c * dw) - u * (c * dd - u * dv)),
1701 dC2 * (c * d - u * v) +
1702 C2 * (dc * d + c * dd - du * v - u * dv) - dC1 * w -
1703 C1 * dw +
1704 dC3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
1705 w * (-d2 + v2 + w2)) +
1706 C3 * (-dc * (b * v + c * w) + du * (b * d + u * w) +
1707 dw * (-d2 + v2 + w2) - c * (db * v + dc * w) +
1708 u * (db * d + du * w) + w * (-dd2 + dv2 + dw2) -
1709 c * (b * dv + c * dw) + u * (b * dd + u * dw)),
1710 dC0 + dC2 * (d2 - v2 - w2) + C2 * (dd2 - dv2 - dw2))
1711 .finished();
1712 }
1713 }
1714 }
1715 }
1716}
1717
1718inline void transmat(TransmissionMatrix& T,
1719 const PropagationMatrix& K1,
1720 const PropagationMatrix& K2,
1721 const Numeric& r) noexcept {
1722 switch (K1.StokesDimensions()) {
1723 case 4:
1724 transmat4(T, K1, K2, r);
1725 break;
1726 case 3:
1727 transmat3(T, K1, K2, r);
1728 break;
1729 case 2:
1730 transmat2(T, K1, K2, r);
1731 break;
1732 case 1:
1733 transmat1(T, K1, K2, r);
1734 break;
1735 }
1736}
1737
1738inline void dtransmat(TransmissionMatrix& T,
1739 ArrayOfTransmissionMatrix& dT1,
1740 ArrayOfTransmissionMatrix& dT2,
1741 const PropagationMatrix& K1,
1742 const PropagationMatrix& K2,
1743 const ArrayOfPropagationMatrix& dK1,
1744 const ArrayOfPropagationMatrix& dK2,
1745 const Numeric& r,
1746 const Numeric& dr_dT1 = 0,
1747 const Numeric& dr_dT2 = 0,
1748 const Index it = -1,
1749 const Index iz = 0,
1750 const Index ia = 0) noexcept {
1751 switch (K1.StokesDimensions()) {
1752 case 4:
1753 dtransmat4(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1754 break;
1755 case 3:
1756 dtransmat3(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1757 break;
1758 case 2:
1759 dtransmat2(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1760 break;
1761 case 1:
1762 dtransmat1(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1763 break;
1764 }
1765}
1766
1767void stepwise_transmission(TransmissionMatrix& T,
1768 ArrayOfTransmissionMatrix& dT1,
1769 ArrayOfTransmissionMatrix& dT2,
1770 const PropagationMatrix& K1,
1771 const PropagationMatrix& K2,
1772 const ArrayOfPropagationMatrix& dK1,
1773 const ArrayOfPropagationMatrix& dK2,
1774 const Numeric& r,
1775 const Numeric& dr_dtemp1,
1776 const Numeric& dr_dtemp2,
1777 const Index temp_deriv_pos) {
1778 if (not dT1.nelem())
1779 transmat(T, K1, K2, r);
1780 else
1781 dtransmat(
1782 T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dtemp1, dr_dtemp2, temp_deriv_pos);
1783}
1784
1785void stepwise_source(RadiationVector& J,
1786 ArrayOfRadiationVector& dJ,
1787 RadiationVector& J_add,
1788 const PropagationMatrix& K,
1789 const StokesVector& a,
1790 const StokesVector& S,
1791 const ArrayOfPropagationMatrix& dK,
1792 const ArrayOfStokesVector& da,
1793 const ArrayOfStokesVector& dS,
1794 const ConstVectorView& B,
1795 const ConstVectorView& dB_dT,
1796 const ArrayOfRetrievalQuantity& jacobian_quantities,
1797 const bool& jacobian_do) {
1798 for (Index i = 0; i < K.NumberOfFrequencies(); i++) {
1799 if (K.IsRotational(i)) {
1800 J.SetZero(i);
1801 if (jacobian_do) {
1802 for (Index j = 0; j < jacobian_quantities.nelem(); j++)
1803 if (dJ[j].Frequencies()) dJ[j].SetZero(i);
1804 }
1805 } else {
1806 J.setSource(a, B, S, i);
1807
1808 switch (J.stokes_dim) {
1809 case 4: {
1810 const auto invK = inv_prop_matrix<4>(K.Data()(0, 0, i, joker));
1811 J.Vec4(i) = invK * J.Vec4(i);
1812 if (jacobian_do) {
1813 for (Index j = 0; j < jacobian_quantities.nelem(); j++) {
1814 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and
1815 dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1816 dJ[j].Vec4(i).noalias() =
1817 0.5 * invK *
1818 (source_vector<4>(
1819 a,
1820 B,
1821 da[j],
1822 dB_dT,
1823 dS[j],
1824 jacobian_quantities[j] == Jacobian::Atm::Temperature,
1825 i) -
1826 prop_matrix<4>(dK[j].Data()(0, 0, i, joker)) * J.Vec4(i));
1827 }
1828 }
1829 }
1830 if (J_add.Frequencies()) {
1831 J_add.Vec4(i) = invK * J_add.Vec4(i);
1832 //TODO: Add jacobians dJ_add of additional source
1833 }
1834 } break;
1835 case 3: {
1836 const auto invK = inv_prop_matrix<3>(K.Data()(0, 0, i, joker));
1837 J.Vec3(i) = invK * J.Vec3(i);
1838 if (jacobian_do) {
1839 for (Index j = 0; j < jacobian_quantities.nelem(); j++) {
1840 // Skip others!
1841 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and
1842 dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1843 dJ[j].Vec3(i).noalias() =
1844 0.5 * invK *
1845 (source_vector<3>(
1846 a,
1847 B,
1848 da[j],
1849 dB_dT,
1850 dS[j],
1851 jacobian_quantities[j] == Jacobian::Atm::Temperature,
1852 i) -
1853 prop_matrix<3>(dK[j].Data()(0, 0, i, joker)) * J.Vec3(i));
1854 }
1855 }
1856 }
1857 if (J_add.Frequencies()) {
1858 J_add.Vec3(i) = invK * J_add.Vec3(i);
1859 //TODO: Add jacobians dJ_add of additional source
1860 }
1861 } break;
1862 case 2: {
1863 const auto invK = inv_prop_matrix<2>(K.Data()(0, 0, i, joker));
1864 J.Vec2(i) = invK * J.Vec2(i);
1865 if (jacobian_do) {
1866 for (Index j = 0; j < jacobian_quantities.nelem(); j++) {
1867 // Skip others!
1868 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and
1869 dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1870 dJ[j].Vec2(i).noalias() =
1871 0.5 * invK *
1872 (source_vector<2>(
1873 a,
1874 B,
1875 da[j],
1876 dB_dT,
1877 dS[j],
1878 jacobian_quantities[j] == Jacobian::Atm::Temperature,
1879 i) -
1880 prop_matrix<2>(dK[j].Data()(0, 0, i, joker)) * J.Vec2(i));
1881 }
1882 }
1883 }
1884 if (J_add.Frequencies()) {
1885 J_add.Vec2(i) = invK * J_add.Vec2(i);
1886 //TODO: Add jacobians dJ_add of additional source
1887 }
1888 } break;
1889 default: {
1890 const auto invK = inv_prop_matrix<1>(K.Data()(0, 0, i, joker));
1891 J.Vec1(i) = invK * J.Vec1(i);
1892 if (jacobian_do) {
1893 for (Index j = 0; j < jacobian_quantities.nelem(); j++) {
1894 // Skip others!
1895 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and
1896 dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1897 dJ[j].Vec1(i).noalias() =
1898 0.5 * invK *
1899 (source_vector<1>(
1900 a,
1901 B,
1902 da[j],
1903 dB_dT,
1904 dS[j],
1905 jacobian_quantities[j] == Jacobian::Atm::Temperature,
1906 i) -
1907 prop_matrix<1>(dK[j].Data()(0, 0, i, joker)) * J.Vec1(i));
1908 }
1909 }
1910 }
1911 if (J_add.Frequencies()) {
1912 J_add.Vec1(i) = invK * J_add.Vec1(i);
1913 //TODO: Add jacobians dJ_add of additional source
1914 }
1915 } break;
1916 }
1917 }
1918 }
1919
1920 if (J_add.Frequencies()) {
1921 J += J_add;
1922 }
1923}
1924
1925void update_radiation_vector(
1926 RadiationVector& I,
1927 ArrayOfRadiationVector& dI1,
1928 ArrayOfRadiationVector& dI2,
1929 const RadiationVector& J1,
1930 const RadiationVector& J2,
1931 const ArrayOfRadiationVector& dJ1,
1932 const ArrayOfRadiationVector& dJ2,
1933 const TransmissionMatrix& T,
1934 const TransmissionMatrix& PiT,
1935 const ArrayOfTransmissionMatrix& dT1,
1936 const ArrayOfTransmissionMatrix& dT2,
1937 [[maybe_unused]] const PropagationMatrix& K1,
1938 [[maybe_unused]] const PropagationMatrix& K2,
1939 [[maybe_unused]] const ArrayOfPropagationMatrix& dK1,
1940 [[maybe_unused]] const ArrayOfPropagationMatrix& dK2,
1941 [[maybe_unused]] const Numeric r,
1942 [[maybe_unused]] const Vector& dr1,
1943 [[maybe_unused]] const Vector& dr2,
1944 [[maybe_unused]] const Index ia,
1945 [[maybe_unused]] const Index iz,
1946 const RadiativeTransferSolver solver) {
1947 switch (solver) {
1948 case RadiativeTransferSolver::Emission: {
1949 I.rem_avg(J1, J2);
1950 for (size_t i = 0; i < dI1.size(); i++) {
1951 dI1[i].addDerivEmission(PiT, dT1[i], T, I, dJ1[i]);
1952 dI2[i].addDerivEmission(PiT, dT2[i], T, I, dJ2[i]);
1953 }
1954 I.leftMul(T);
1955 I.add_avg(J1, J2);
1956 } break;
1957
1958 case RadiativeTransferSolver::Transmission: {
1959 for (size_t i = 0; i < dI1.size(); i++) {
1960 dI1[i].addDerivTransmission(PiT, dT1[i], I);
1961 dI2[i].addDerivTransmission(PiT, dT2[i], I);
1962 }
1963 I.leftMul(T);
1964 } break;
1965
1966 case RadiativeTransferSolver::LinearWeightedEmission: {
1967 ARTS_USER_ERROR_IF(
1968 dI1.size(),
1969 "Cannot support derivatives with current integration method\n");
1970
1971 I.leftMul(T);
1972 I.add_weighted(T,
1973 J1,
1974 J2,
1975 K1.Data()(ia, iz, joker, joker),
1976 K2.Data()(ia, iz, joker, joker),
1977 r);
1978
1979 } break;
1980 }
1981}
1982
1983ArrayOfTransmissionMatrix cumulative_transmission(
1984 const ArrayOfTransmissionMatrix& T,
1985 const CumulativeTransmission type) /*[[expects: T.nelem()>0]]*/
1986{
1987 const Index n = T.nelem();
1988 const Index nf = n ? T[0].Frequencies() : 1;
1989 const Index ns = n ? T[0].stokes_dim : 1;
1990 ArrayOfTransmissionMatrix PiT(n, TransmissionMatrix(nf, ns));
1991 switch (type) {
1992 case CumulativeTransmission::Forward: {
1993 for (Index i = 1; i < n; i++) PiT[i].mul(PiT[i - 1], T[i]);
1994 } break;
1995 case CumulativeTransmission::Reverse: {
1996 for (Index i = 1; i < n; i++) PiT[i].mul(T[i], PiT[i - 1]);
1997 } break;
1998 }
1999 return PiT; // Note how the output is such that forward transmission is from -1 to 0
2000}
2001
2002// TEST CODE BEGIN
2003
2004void set_backscatter_radiation_vector(
2005 ArrayOfRadiationVector& I,
2006 ArrayOfArrayOfArrayOfRadiationVector& dI,
2007 const RadiationVector& I_incoming,
2008 const ArrayOfTransmissionMatrix& T,
2009 const ArrayOfTransmissionMatrix& PiTf,
2010 const ArrayOfTransmissionMatrix& PiTr,
2011 const ArrayOfTransmissionMatrix& Z,
2012 const ArrayOfArrayOfTransmissionMatrix& dT1,
2013 const ArrayOfArrayOfTransmissionMatrix& dT2,
2014 const ArrayOfArrayOfTransmissionMatrix& dZ,
2015 const BackscatterSolver solver) {
2016 const Index np = I.nelem();
2017 const Index nv = np ? I[0].Frequencies() : 0;
2018 const Index ns = np ? I[0].stokes_dim : 1;
2019 const Index nq = np ? dI[0][0].nelem() : 0;
2020
2021 // For all transmission, the I-vector is the same
2022 for (Index ip = 0; ip < np; ip++)
2023 I[ip].setBackscatterTransmission(I_incoming, PiTr[ip], PiTf[ip], Z[ip]);
2024
2025 for (Index ip = 0; ip < np; ip++) {
2026 for (Index iq = 0; iq < nq; iq++) {
2027 dI[ip][ip][iq].setBackscatterTransmissionDerivative(
2028 I_incoming, PiTr[ip], PiTf[ip], dZ[ip][iq]);
2029 }
2030 }
2031
2032 switch (solver) {
2033 case BackscatterSolver::CommutativeTransmission: {
2034 // Forward and backwards transmission derivatives
2035 switch (ns) {
2036 case 1: {
2037 BackscatterSolverCommutativeTransmissionStokesDimOne:
2038 for (Index ip = 0; ip < np; ip++) {
2039 for (Index j = ip; j < np; j++) {
2040 for (Index iq = 0; iq < nq; iq++) {
2041 for (Index iv = 0; iv < nv; iv++) {
2042 dI[ip][j][iq].Vec1(iv).noalias() +=
2043 T[ip].Mat1(iv).inverse() *
2044 (dT1[ip][iq].Mat1(iv) + dT2[ip][iq].Mat1(iv)) *
2045 I[j].Vec1(iv);
2046
2047 if (j < np - 1 and j > ip)
2048 dI[ip][j][iq].Vec1(iv).noalias() +=
2049 T[ip + 1].Mat1(iv).inverse() *
2050 (dT1[ip][iq].Mat1(iv) + dT2[ip][iq].Mat1(iv)) *
2051 I[j].Vec1(iv);
2052 }
2053 }
2054 }
2055 }
2056 } break;
2057 case 2: {
2058 for (Index ip = 0; ip < np; ip++) {
2059 for (Index j = ip; j < np; j++) {
2060 for (Index iq = 0; iq < nq; iq++) {
2061 for (Index iv = 0; iv < nv; iv++) {
2062 dI[ip][j][iq].Vec2(iv).noalias() +=
2063 T[ip].Mat2(iv).inverse() *
2064 (dT1[ip][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
2065 I[j].Vec2(iv);
2066
2067 if (j < np - 1 and j > ip)
2068 dI[ip][j][iq].Vec2(iv).noalias() +=
2069 T[ip + 1].Mat2(iv).inverse() *
2070 (dT1[ip][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
2071 I[j].Vec2(iv);
2072 }
2073 }
2074 }
2075 }
2076 } break;
2077 case 3: {
2078 for (Index ip = 0; ip < np; ip++) {
2079 for (Index j = ip; j < np; j++) {
2080 for (Index iq = 0; iq < nq; iq++) {
2081 for (Index iv = 0; iv < nv; iv++) {
2082 dI[ip][j][iq].Vec3(iv).noalias() +=
2083 T[ip].Mat3(iv).inverse() *
2084 (dT1[ip][iq].Mat3(iv) + dT2[ip][iq].Mat3(iv)) *
2085 I[j].Vec3(iv);
2086
2087 if (j < np - 1 and j > ip)
2088 dI[ip][j][iq].Vec3(iv).noalias() +=
2089 T[ip + 1].Mat3(iv).inverse() *
2090 (dT1[ip][iq].Mat3(iv) + dT2[ip][iq].Mat3(iv)) *
2091 I[j].Vec3(iv);
2092 }
2093 }
2094 }
2095 }
2096 } break;
2097 case 4: {
2098 for (Index ip = 0; ip < np; ip++) {
2099 for (Index j = ip; j < np; j++) {
2100 for (Index iq = 0; iq < nq; iq++) {
2101 for (Index iv = 0; iv < nv; iv++) {
2102 dI[ip][j][iq].Vec4(iv).noalias() +=
2103 T[ip].Mat4(iv).inverse() *
2104 (dT1[ip][iq].Mat4(iv) + dT2[ip][iq].Mat4(iv)) *
2105 I[j].Vec4(iv);
2106
2107 if (j < np - 1 and j > ip)
2108 dI[ip][j][iq].Vec4(iv).noalias() +=
2109 T[ip + 1].Mat4(iv).inverse() *
2110 (dT1[ip][iq].Mat4(iv) + dT2[ip][iq].Mat4(iv)) *
2111 I[j].Vec4(iv);
2112 }
2113 }
2114 }
2115 }
2116 } break;
2117 }
2118 } break;
2119 case BackscatterSolver::FullTransmission: {
2120 std::runtime_error("Do not activate this code. It is not ready yet");
2121 switch (ns) {
2122 case 1: {
2123 // This is the same as CommutativeTransmission, so use that code
2124 goto BackscatterSolverCommutativeTransmissionStokesDimOne;
2125 } break;
2126 case 2: {
2127 for (Index ip = 0; ip < np; ip++) {
2128 for (Index j = ip; j < np; j++) {
2129 for (Index iq = 0; iq < nq; iq++) {
2130 for (Index iv = 0; iv < nv; iv++) {
2131 if (ip > 1) {
2132 dI[ip][j][iq].Vec2(iv).noalias() +=
2133 PiTr[ip - 2].Mat2(iv) *
2134 (dT1[ip - 1][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
2135 PiTr[ip - 1].Mat2(iv).inverse() * I[j].Vec2(iv);
2136
2137 if (j < np - 1)
2138 dI[ip][j][iq].Vec2(iv).noalias() +=
2139 PiTr[ip].Mat2(iv) * Z[ip].Mat2(iv) *
2140 PiTf[ip - 2].Mat2(iv) *
2141 (dT1[ip][iq].Mat2(iv) + dT2[ip + 1][iq].Mat2(iv)) *
2142 PiTf[ip - 1].Mat2(iv).inverse() * PiTf[ip].Mat2(iv) *
2143 I[0].Vec2(iv);
2144 } else {
2145 dI[ip][j][iq].Vec2(iv).noalias() +=
2146 (dT1[ip - 1][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
2147 PiTr[ip - 1].Mat2(iv).inverse() * I[j].Vec2(iv);
2148
2149 if (j < np - 1)
2150 dI[ip][j][iq].Vec2(iv).noalias() +=
2151 PiTr[ip].Mat2(iv) * Z[ip].Mat2(iv) *
2152 (dT1[ip][iq].Mat2(iv) + dT2[ip + 1][iq].Mat2(iv)) *
2153 PiTf[ip - 1].Mat2(iv).inverse() * PiTf[ip].Mat2(iv) *
2154 I[0].Vec2(iv);
2155 }
2156 }
2157 }
2158 }
2159 }
2160 } break;
2161 case 3: {
2162 } break;
2163 case 4: {
2164 } break;
2165 }
2166 } break;
2167 }
2168}
2169
2170ArrayOfTransmissionMatrix bulk_backscatter(const ConstTensor5View& Pe,
2171 const ConstMatrixView& pnd) {
2172 const Index ns = Pe.ncols();
2173 const Index nv = Pe.npages();
2174 const Index np = Pe.nbooks();
2175 const Index ne = Pe.nshelves();
2176
2177 ArrayOfTransmissionMatrix aotm(np, TransmissionMatrix(nv, ns));
2178
2179 for (Index ip = 0; ip < np; ip++) {
2180 aotm[ip].setZero();
2181
2182 switch (ns) {
2183 case 4:
2184 for (Index iv = 0; iv < nv; iv++)
2185 for (Index ie = 0; ie < ne; ie++)
2186 aotm[ip].Mat4(iv).noalias() +=
2187 pnd(ie, ip) * prop_matrix<4>(Pe(ie, ip, iv, joker, joker));
2188 break;
2189 case 3:
2190 for (Index iv = 0; iv < nv; iv++)
2191 for (Index ie = 0; ie < ne; ie++)
2192 aotm[ip].Mat3(iv).noalias() +=
2193 pnd(ie, ip) * prop_matrix<3>(Pe(ie, ip, iv, joker, joker));
2194 break;
2195 case 2:
2196 for (Index iv = 0; iv < nv; iv++)
2197 for (Index ie = 0; ie < ne; ie++)
2198 aotm[ip].Mat2(iv).noalias() +=
2199 pnd(ie, ip) * prop_matrix<2>(Pe(ie, ip, iv, joker, joker));
2200 break;
2201 case 1:
2202 for (Index iv = 0; iv < nv; iv++)
2203 for (Index ie = 0; ie < ne; ie++)
2204 aotm[ip].Mat1(iv).noalias() +=
2205 pnd(ie, ip) * prop_matrix<1>(Pe(ie, ip, iv, joker, joker));
2206 break;
2207 }
2208 }
2209 return aotm;
2210}
2211
2212ArrayOfArrayOfTransmissionMatrix bulk_backscatter_derivative(
2213 const ConstTensor5View& Pe, const ArrayOfMatrix& dpnd_dx) {
2214 const Index ns = Pe.ncols();
2215 const Index nv = Pe.npages();
2216 const Index np = Pe.nbooks();
2217 const Index ne = Pe.nshelves();
2218 const Index nq = dpnd_dx.nelem();
2219
2220 ArrayOfArrayOfTransmissionMatrix aoaotm(
2221 np, ArrayOfTransmissionMatrix(nq, TransmissionMatrix(nv, ns)));
2222
2223 for (Index ip = 0; ip < np; ip++) {
2224 for (Index iq = 0; iq < nq; iq++) {
2225 aoaotm[ip][iq].setZero();
2226
2227 if (not dpnd_dx[iq].empty()) {
2228 switch (ns) {
2229 case 4:
2230 for (Index iv = 0; iv < nv; iv++)
2231 for (Index ie = 0; ie < ne; ie++)
2232 aoaotm[ip][iq].Mat4(iv).noalias() +=
2233 dpnd_dx[iq](ie, ip) *
2234 prop_matrix<4>(Pe(ie, ip, iv, joker, joker));
2235 break;
2236 case 3:
2237 for (Index iv = 0; iv < nv; iv++)
2238 for (Index ie = 0; ie < ne; ie++)
2239 aoaotm[ip][iq].Mat3(iv).noalias() +=
2240 dpnd_dx[iq](ie, ip) *
2241 prop_matrix<3>(Pe(ie, ip, iv, joker, joker));
2242 break;
2243 case 2:
2244 for (Index iv = 0; iv < nv; iv++)
2245 for (Index ie = 0; ie < ne; ie++)
2246 aoaotm[ip][iq].Mat2(iv).noalias() +=
2247 dpnd_dx[iq](ie, ip) *
2248 prop_matrix<2>(Pe(ie, ip, iv, joker, joker));
2249 break;
2250 case 1:
2251 for (Index iv = 0; iv < nv; iv++)
2252 for (Index ie = 0; ie < ne; ie++)
2253 aoaotm[ip][iq].Mat1(iv).noalias() +=
2254 dpnd_dx[iq](ie, ip) *
2255 prop_matrix<1>(Pe(ie, ip, iv, joker, joker));
2256 break;
2257 }
2258 }
2259 }
2260 }
2261 return aoaotm;
2262}
2263
2264// TEST CODE END
2265
2266std::ostream& operator<<(std::ostream& os, const TransmissionMatrix& tm) {
2267 for (const auto& T : tm.T4) os << T << '\n';
2268 for (const auto& T : tm.T3) os << T << '\n';
2269 for (const auto& T : tm.T2) os << T << '\n';
2270 for (const auto& T : tm.T1) os << T << '\n';
2271 return os;
2272}
2273
2274std::ostream& operator<<(std::ostream& os, const RadiationVector& rv) {
2275 // Write the transpose because it looks better...
2276 for (const auto& R : rv.R4) os << R.transpose() << '\n';
2277 for (const auto& R : rv.R3) os << R.transpose() << '\n';
2278 for (const auto& R : rv.R2) os << R.transpose() << '\n';
2279 for (const auto& R : rv.R1) os << R.transpose() << '\n';
2280 return os;
2281}
2282
2283std::istream& operator>>(std::istream& is, TransmissionMatrix& tm) {
2284 for (auto& T : tm.T4)
2285 is >> double_imanip() >> T(0, 0) >> T(0, 1) >> T(0, 2) >> T(0, 3) >>
2286 T(1, 0) >> T(1, 1) >> T(1, 2) >> T(1, 3) >> T(2, 0) >> T(2, 1) >>
2287 T(2, 2) >> T(2, 3) >> T(3, 0) >> T(3, 1) >> T(3, 2) >> T(3, 3);
2288 for (auto& T : tm.T3)
2289 is >> double_imanip() >> T(0, 0) >> T(0, 1) >> T(0, 2) >> T(1, 0) >>
2290 T(1, 1) >> T(1, 2) >> T(2, 0) >> T(2, 1) >> T(2, 2);
2291 for (auto& T : tm.T2)
2292 is >> double_imanip() >> T(0, 0) >> T(0, 1) >> T(1, 0) >> T(1, 1);
2293 for (auto& T : tm.T1) is >> double_imanip() >> T(0, 0);
2294 return is;
2295}
2296
2297std::istream& operator>>(std::istream& is, RadiationVector& rv) {
2298 for (auto& R : rv.R4) is >> double_imanip() >> R[0] >> R[1] >> R[2] >> R[3];
2299 for (auto& R : rv.R3) is >> double_imanip() >> R[0] >> R[1] >> R[2];
2300 for (auto& R : rv.R2) is >> double_imanip() >> R[0] >> R[1] >> R[2];
2301 for (auto& R : rv.R1) is >> double_imanip() >> R[0];
2302 return is;
2303}
2304
2305TransmissionMatrix::TransmissionMatrix(const PropagationMatrix& pm,
2306 const Numeric& r) {
2307 *this = TransmissionMatrix(pm.NumberOfFrequencies(), pm.StokesDimensions());
2308 transmat(*this, pm, pm, r); // Slower to compute, faster to implement...
2309}
arts_conversions.h
Common ARTS conversions.
Array
This can be used to make arrays out of anything.
Definition array.h:31
Array::nelem
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
double_imanip
Input manipulator class for doubles to enable nan and inf parsing.
Definition double_imanip.h:25
ARTS_ASSERT
#define ARTS_ASSERT(condition,...)
Definition debug.h:86
ARTS_USER_ERROR_IF
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
double_imanip.h
Fast double input stream with support for parsing nan and inf.
mul
constexpr std::size_t mul(Inds... inds) noexcept
Definition grids.h:10
Constant::inv_sqrt_2
constexpr Numeric inv_sqrt_2
Inverse of the square root of 2.
Definition arts_constants.h:90
RadiationVector
Radiation Vector for Stokes dimension 1-4.
Definition transmissionmatrix.h:480
RadiationVector::R2
std::vector< Eigen::Vector2d > R2
Definition transmissionmatrix.h:484
RadiationVector::RadiationVector
RadiationVector(Index nf=0, Index stokes=1)
Construct a new Radiation Vector object.
Definition transmissionmatrix.cc:186
RadiationVector::Frequencies
Index Frequencies() const
Get frequency count.
Definition transmissionmatrix.cc:549
RadiationVector::Vec1
const Eigen::Matrix< double, 1, 1 > & Vec1(size_t i) const
Return Vector at position.
Definition transmissionmatrix.cc:198
RadiationVector::setSource
void setSource(const StokesVector &a, const ConstVectorView &B, const StokesVector &S, Index i)
Set this to source vector at position.
Definition transmissionmatrix.cc:505
RadiationVector::addMultiplied
void addMultiplied(const TransmissionMatrix &A, const RadiationVector &x)
Add multiply.
Definition transmissionmatrix.cc:391
RadiationVector::Vec4
const Eigen::Vector4d & Vec4(size_t i) const
Return Vector at position.
Definition transmissionmatrix.cc:195
RadiationVector::leftMul
void leftMul(const TransmissionMatrix &T)
Multiply radiation vector from the left.
Definition transmissionmatrix.cc:220
RadiationVector::operator()
const Numeric & operator()(const Index i, const Index j) const
Access operator.
Definition transmissionmatrix.cc:480
RadiationVector::SetZero
void SetZero(size_t i)
Set Radiation Vector to Zero at position.
Definition transmissionmatrix.cc:227
RadiationVector::rem_avg
void rem_avg(const RadiationVector &O1, const RadiationVector &O2)
Remove the average of two other RadiationVector from *this.
Definition transmissionmatrix.cc:261
RadiationVector::R1
std::vector< Eigen::Matrix< double, 1, 1 > > R1
Definition transmissionmatrix.h:485
RadiationVector::add_avg
void add_avg(const RadiationVector &O1, const RadiationVector &O2)
Add the average of two other RadiationVector to *this.
Definition transmissionmatrix.cc:273
RadiationVector::R4
std::vector< Eigen::Vector4d > R4
Definition transmissionmatrix.h:482
RadiationVector::operator+=
RadiationVector & operator+=(const RadiationVector &rv)
Addition operator.
Definition transmissionmatrix.cc:466
RadiationVector::Vec
Eigen::VectorXd Vec(size_t i) const
Return Vector at position by copy.
Definition transmissionmatrix.cc:248
RadiationVector::addDerivTransmission
void addDerivTransmission(const TransmissionMatrix &PiT, const TransmissionMatrix &dT, const RadiationVector &I)
Add the transmission derivative to this.
Definition transmissionmatrix.cc:378
RadiationVector::setDerivReflection
void setDerivReflection(const RadiationVector &I, const TransmissionMatrix &PiT, const TransmissionMatrix &Z, const TransmissionMatrix &dZ)
Sets *this to the reflection derivative.
Definition transmissionmatrix.cc:399
RadiationVector::Vec3
const Eigen::Vector3d & Vec3(size_t i) const
Return Vector at position.
Definition transmissionmatrix.cc:196
RadiationVector::Vec2
const Eigen::Vector2d & Vec2(size_t i) const
Return Vector at position.
Definition transmissionmatrix.cc:197
RadiationVector::R3
std::vector< Eigen::Vector3d > R3
Definition transmissionmatrix.h:483
RadiationVector::addWeightedDerivEmission
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.
Definition transmissionmatrix.cc:348
RadiationVector::setBackscatterTransmission
void setBackscatterTransmission(const RadiationVector &I0, const TransmissionMatrix &Tr, const TransmissionMatrix &Tf, const TransmissionMatrix &Z)
Set this to backscatter transmission.
Definition transmissionmatrix.cc:414
RadiationVector::addDerivEmission
void addDerivEmission(const TransmissionMatrix &PiT, const TransmissionMatrix &dT, const TransmissionMatrix &T, const RadiationVector &ImJ, const RadiationVector &dJ)
Add the emission derivative to this.
Definition transmissionmatrix.cc:329
RadiationVector::operator=
RadiationVector & operator=(const RadiationVector &rv)=default
Assign old radiation vector to this.
RadiationVector::add_weighted
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.
Definition transmissionmatrix.cc:285
RadiationVector::setBackscatterTransmissionDerivative
void setBackscatterTransmissionDerivative(const RadiationVector &I0, const TransmissionMatrix &Tr, const TransmissionMatrix &Tf, const TransmissionMatrix &dZ)
Set this to backscatter transmission scatter derivative.
Definition transmissionmatrix.cc:428
TransmissionMatrix
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
Definition transmissionmatrix.h:148
TransmissionMatrix::mul_aliased
void mul_aliased(const TransmissionMatrix &A, const TransmissionMatrix &B)
Set this to a multiple of A by B.
Definition transmissionmatrix.cc:95
TransmissionMatrix::T1
std::vector< Eigen::Matrix< double, 1, 1 > > T1
Definition transmissionmatrix.h:153
TransmissionMatrix::T4
std::vector< Eigen::Matrix4d > T4
Definition transmissionmatrix.h:150
TransmissionMatrix::operator*=
TransmissionMatrix & operator*=(const Numeric &scale)
Scale self.
Definition transmissionmatrix.cc:164
TransmissionMatrix::Mat2
const Eigen::Matrix2d & Mat2(size_t i) const
Get Matrix at position.
Definition transmissionmatrix.cc:27
TransmissionMatrix::Mat
Eigen::MatrixXd Mat(size_t i) const
Get Matrix at position by copy.
Definition transmissionmatrix.cc:59
TransmissionMatrix::TraMat
auto & TraMat(size_t i) noexcept
Simple template access for the transmission.
Definition transmissionmatrix.h:340
TransmissionMatrix::Frequencies
Index Frequencies() const
Number of frequencies.
Definition transmissionmatrix.cc:131
TransmissionMatrix::second_order_integration_dsource
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
Definition transmissionmatrix.h:425
TransmissionMatrix::mul
void mul(const TransmissionMatrix &A, const TransmissionMatrix &B)
Set this to a multiple of A by B.
Definition transmissionmatrix.cc:87
TransmissionMatrix::T3
std::vector< Eigen::Matrix3d > T3
Definition transmissionmatrix.h:151
TransmissionMatrix::Mat4
const Eigen::Matrix4d & Mat4(size_t i) const
Get Matrix at position.
Definition transmissionmatrix.cc:25
TransmissionMatrix::second_order_integration_source
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
Definition transmissionmatrix.h:392
TransmissionMatrix::setIdentity
void setIdentity()
Set to identity matrix.
Definition transmissionmatrix.cc:73
TransmissionMatrix::Mat1
const Eigen::Matrix< double, 1, 1 > & Mat1(size_t i) const
Get Matrix at position.
Definition transmissionmatrix.cc:28
TransmissionMatrix::T2
std::vector< Eigen::Matrix2d > T2
Definition transmissionmatrix.h:152
TransmissionMatrix::Mat3
const Eigen::Matrix3d & Mat3(size_t i) const
Get Matrix at position.
Definition transmissionmatrix.cc:26
TransmissionMatrix::operator()
Numeric operator()(const Index i, const Index j, const Index k) const
Access value in matrix.
Definition transmissionmatrix.cc:103
TransmissionMatrix::operator=
TransmissionMatrix & operator=(const TransmissionMatrix &tm)=default
Assignment operator.
TransmissionMatrix::operator+=
TransmissionMatrix & operator+=(const LazyScale< TransmissionMatrix > &lstm)
Assign to *this lazily.
Definition transmissionmatrix.cc:151
TransmissionMatrix::setZero
void setZero()
Set to zero matrix.
Definition transmissionmatrix.cc:80
TransmissionMatrix::TransmissionMatrix
TransmissionMatrix(Index nf=0, Index stokes=1)
Construct a new Transmission Matrix object.
Definition transmissionmatrix.cc:16
transmat3
void transmat3(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r, const Index iz=0, const Index ia=0) noexcept
Definition transmissionmatrix.cc:689
lower_is_considered_zero_for_sinc_likes
constexpr Numeric lower_is_considered_zero_for_sinc_likes
Definition transmissionmatrix.cc:562
bulk_backscatter
ArrayOfTransmissionMatrix bulk_backscatter(const ConstTensor5View &Pe, const ConstMatrixView &pnd)
Bulk back-scattering.
Definition transmissionmatrix.cc:2170
dtransmat
void dtransmat(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dT1=0, const Numeric &dr_dT2=0, const Index it=-1, const Index iz=0, const Index ia=0) noexcept
Definition transmissionmatrix.cc:1738
dtransmat4
void dtransmat4(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dT1, const Numeric &dr_dT2, const Index it, const Index iz, const Index ia) noexcept
Definition transmissionmatrix.cc:1164
dtransmat3
void dtransmat3(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dT1, const Numeric &dr_dT2, const Index it, const Index iz, const Index ia) noexcept
Definition transmissionmatrix.cc:960
transmat4
void transmat4(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r, const Index iz=0, const Index ia=0) noexcept
Definition transmissionmatrix.cc:752
transmat2
void transmat2(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r, const Index iz=0, const Index ia=0) noexcept
Definition transmissionmatrix.cc:673
dtransmat1
void dtransmat1(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dT1, const Numeric &dr_dT2, const Index it, const Index iz, const Index ia) noexcept
Definition transmissionmatrix.cc:874
bulk_backscatter_derivative
ArrayOfArrayOfTransmissionMatrix bulk_backscatter_derivative(const ConstTensor5View &Pe, const ArrayOfMatrix &dpnd_dx)
Derivatives of bulk back-scattering
Definition transmissionmatrix.cc:2212
source_vector
constexpr Eigen::Matrix< Numeric, N, 1 > source_vector(const StokesVector &a, const ConstVectorView &B, const StokesVector &da, const ConstVectorView &dB_dT, const StokesVector &dS, bool dT, Index i)
Definition transmissionmatrix.cc:565
operator*
LazyScale< TransmissionMatrix > operator*(const TransmissionMatrix &tm, const Numeric &x)
Lazy scale of Transmission Matrix.
Definition transmissionmatrix.cc:176
cumulative_transmission
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
Definition transmissionmatrix.cc:1983
set_backscatter_radiation_vector
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.
Definition transmissionmatrix.cc:2004
stepwise_source
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.
Definition transmissionmatrix.cc:1785
operator>>
std::istream & operator>>(std::istream &is, TransmissionMatrix &tm)
Definition transmissionmatrix.cc:2283
update_radiation_vector
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.
Definition transmissionmatrix.cc:1925
transmat
void transmat(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r) noexcept
Definition transmissionmatrix.cc:1718
operator<<
std::ostream & operator<<(std::ostream &os, const TransmissionMatrix &tm)
Definition transmissionmatrix.cc:2266
stepwise_transmission
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.
Definition transmissionmatrix.cc:1767
dtransmat2
void dtransmat2(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dT1, const Numeric &dr_dT2, const Index it, const Index iz, const Index ia) noexcept
Definition transmissionmatrix.cc:909
transmat1
void transmat1(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r, const Index iz=0, const Index ia=0) noexcept
Definition transmissionmatrix.cc:662
transmissionmatrix.h
Stuff related to the transmission matrix.
RadiativeTransferSolver
RadiativeTransferSolver
Intended to hold various forward solvers.
Definition transmissionmatrix.h:791
u
#define u
d
#define d
BackscatterSolver
BackscatterSolver
Intended to hold various backscatter solvers.
Definition transmissionmatrix.h:779
CumulativeTransmission
CumulativeTransmission
Intended to hold various ways to accumulate the transmission matrix.
Definition transmissionmatrix.h:785
v
#define v
w
#define w
a
#define a
ArrayOfTransmissionMatrix
Array< TransmissionMatrix > ArrayOfTransmissionMatrix
Definition transmissionmatrix.h:769
c
#define c
b
#define b