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