ARTS 2.5.10 (git: 2f1c442c)
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
34TransmissionMatrix::TransmissionMatrix(Index nf, Index stokes)
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
53TransmissionMatrix& TransmissionMatrix::operator=(
54 const LazyScale<TransmissionMatrix>& lstm) {
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
91void TransmissionMatrix::setIdentity() {
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
98void TransmissionMatrix::setZero() {
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
105void TransmissionMatrix::mul(const TransmissionMatrix& A,
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
113void TransmissionMatrix::mul_aliased(const TransmissionMatrix& A,
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
121Numeric TransmissionMatrix::operator()(const Index i,
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
136Numeric& TransmissionMatrix::operator()(const Index i, const Index j, const Index k) {
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
149[[nodiscard]] Index TransmissionMatrix::Frequencies() const {
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
162TransmissionMatrix::TransmissionMatrix(const ConstMatrixView& mat): TransmissionMatrix(1, mat.nrows()){
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
169TransmissionMatrix& TransmissionMatrix::operator+=(
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
182TransmissionMatrix& TransmissionMatrix::operator*=(const Numeric& scale) {
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
194LazyScale<TransmissionMatrix> operator*(const TransmissionMatrix& tm,
195 const Numeric& x) {
196 return {tm, x};
197}
198
199LazyScale<TransmissionMatrix> operator*(const Numeric& x,
200 const TransmissionMatrix& tm) {
201 return {tm, x};
202}
203
204RadiationVector::RadiationVector(Index nf, Index stokes)
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
225Numeric& RadiationVector::operator()(const Index i, const Index j) {
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
238void RadiationVector::leftMul(const TransmissionMatrix& T) {
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
245void RadiationVector::SetZero(size_t i) {
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
262void RadiationVector::SetZero() {
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
279void RadiationVector::rem_avg(const RadiationVector& O1,
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
291void RadiationVector::add_avg(const RadiationVector& O1,
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
303void RadiationVector::add_weighted(const TransmissionMatrix& T,
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() +=
311 T.second_order_integration_source<4>(T.TraMat<4>(i),
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() +=
320 T.second_order_integration_source<3>(T.TraMat<3>(i),
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() +=
329 T.second_order_integration_source<2>(T.TraMat<2>(i),
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() +=
338 T.second_order_integration_source<1>(T.TraMat<1>(i),
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
347void RadiationVector::addDerivEmission(const TransmissionMatrix& PiT,
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
366void RadiationVector::addWeightedDerivEmission(const TransmissionMatrix& PiT,
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] +
377 T.second_order_integration_dsource<4>(
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] +
382 T.second_order_integration_dsource<3>(
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] +
387 T.second_order_integration_dsource<2>(
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] +
392 T.second_order_integration_dsource<1>(
393 i, dT, far.R1[i], close.R1[i], d.R1[i], isfar));
394}
395
396void RadiationVector::addDerivTransmission(const TransmissionMatrix& PiT,
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
409void RadiationVector::addMultiplied(const TransmissionMatrix& A,
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
417void RadiationVector::setDerivReflection(const RadiationVector& I,
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
432void RadiationVector::setBackscatterTransmission(const RadiationVector& I0,
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
446void RadiationVector::setBackscatterTransmissionDerivative(
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
461RadiationVector& RadiationVector::operator=(const ConstMatrixView& M) {
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
484RadiationVector& RadiationVector::operator+=(const RadiationVector& rv) {
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
523void RadiationVector::setSource(const StokesVector& a,
524 const ConstVectorView& B,
525 const StokesVector& S,
526 Index i) {
527 ARTS_ASSERT(a.NumberOfAzimuthAngles() == 1);
528 ARTS_ASSERT(a.NumberOfZenithAngles() == 1);
529 ARTS_ASSERT(S.NumberOfAzimuthAngles() == 1);
530 ARTS_ASSERT(S.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
567Index RadiationVector::Frequencies() const {
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
580constexpr Numeric lower_is_considered_zero_for_sinc_likes = 1e-4;
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
680inline void transmat1(TransmissionMatrix& T,
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
691inline void transmat2(TransmissionMatrix& T,
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
707inline void transmat3(TransmissionMatrix& T,
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
770inline void transmat4(TransmissionMatrix& T,
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
892inline void dtransmat1(TransmissionMatrix& T,
893 ArrayOfTransmissionMatrix& dT1,
894 ArrayOfTransmissionMatrix& dT2,
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
927inline void dtransmat2(TransmissionMatrix& T,
928 ArrayOfTransmissionMatrix& dT1,
929 ArrayOfTransmissionMatrix& dT2,
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
978inline void dtransmat3(TransmissionMatrix& T,
979 ArrayOfTransmissionMatrix& dT1,
980 ArrayOfTransmissionMatrix& dT2,
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
1182inline void dtransmat4(TransmissionMatrix& T,
1183 ArrayOfTransmissionMatrix& dT1,
1184 ArrayOfTransmissionMatrix& dT2,
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
1736inline void transmat(TransmissionMatrix& T,
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
1756inline void dtransmat(TransmissionMatrix& T,
1757 ArrayOfTransmissionMatrix& dT1,
1758 ArrayOfTransmissionMatrix& dT2,
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
1785void stepwise_transmission(TransmissionMatrix& T,
1786 ArrayOfTransmissionMatrix& dT1,
1787 ArrayOfTransmissionMatrix& dT2,
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
1803void stepwise_source(RadiationVector& J,
1804 ArrayOfRadiationVector& dJ,
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
1943void update_radiation_vector(
1944 RadiationVector& I,
1945 ArrayOfRadiationVector& dI1,
1946 ArrayOfRadiationVector& dI2,
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) {
1966 case RadiativeTransferSolver::Emission: {
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
1976 case RadiativeTransferSolver::Transmission: {
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
1984 case RadiativeTransferSolver::LinearWeightedEmission: {
1985 ARTS_USER_ERROR_IF(
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
2001ArrayOfTransmissionMatrix cumulative_transmission(
2002 const ArrayOfTransmissionMatrix& T,
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;
2008 ArrayOfTransmissionMatrix PiT(n, TransmissionMatrix(nf, ns));
2009 switch (type) {
2010 case CumulativeTransmission::Forward: {
2011 for (Index i = 1; i < n; i++) PiT[i].mul(PiT[i - 1], T[i]);
2012 } break;
2013 case CumulativeTransmission::Reverse: {
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
2022void set_backscatter_radiation_vector(
2023 ArrayOfRadiationVector& I,
2024 ArrayOfArrayOfArrayOfRadiationVector& dI,
2025 const RadiationVector& I_incoming,
2026 const ArrayOfTransmissionMatrix& T,
2027 const ArrayOfTransmissionMatrix& PiTf,
2028 const ArrayOfTransmissionMatrix& PiTr,
2029 const ArrayOfTransmissionMatrix& Z,
2030 const ArrayOfArrayOfTransmissionMatrix& dT1,
2031 const ArrayOfArrayOfTransmissionMatrix& dT2,
2032 const ArrayOfArrayOfTransmissionMatrix& dZ,
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) {
2051 case BackscatterSolver::CommutativeTransmission: {
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;
2137 case BackscatterSolver::FullTransmission: {
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
2188ArrayOfTransmissionMatrix bulk_backscatter(const ConstTensor5View& Pe,
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
2195 ArrayOfTransmissionMatrix aotm(np, TransmissionMatrix(nv, ns));
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
2230ArrayOfArrayOfTransmissionMatrix bulk_backscatter_derivative(
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
2238 ArrayOfArrayOfTransmissionMatrix aoaotm(
2239 np, ArrayOfTransmissionMatrix(nq, TransmissionMatrix(nv, ns)));
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
2323TransmissionMatrix::TransmissionMatrix(const PropagationMatrix& pm,
2324 const Numeric& r) {
2325 *this = TransmissionMatrix(pm.NumberOfFrequencies(), pm.StokesDimensions());
2326 transmat(*this, pm, pm, r); // Slower to compute, faster to implement...
2327}
arts_conversions.h
Common ARTS conversions.
Array
This can be used to make arrays out of anything.
Definition: array.h:48
Array::nelem
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:1065
ConstTensor5View
A constant view of a Tensor5.
Definition: matpackV.h:144
ConstTensor5View::ncols
Index ncols() const noexcept
Definition: matpackV.h:158
ConstTensor5View::npages
Index npages() const noexcept
Definition: matpackV.h:156
ConstTensor5View::nbooks
Index nbooks() const noexcept
Definition: matpackV.h:155
ConstTensor5View::nshelves
Index nshelves() const noexcept
Definition: matpackV.h:154
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:521
LazyScale
Class to help with hidden temporary variables for operations of type Numeric times Class.
Definition: propagationmatrix.h:46
LazyScale::bas
const base & bas
A const reference to a value.
Definition: propagationmatrix.h:57
LazyScale::scale
const Numeric & scale
A const reference to a Numeric.
Definition: propagationmatrix.h:60
Matrix
The Matrix class.
Definition: matpackI.h:1285
PropagationMatrix::NumberOfFrequencies
Index NumberOfFrequencies() const
The number of frequencies of the propagation matrix.
Definition: propagationmatrix.h:224
PropagationMatrix::K24
VectorView K24(const Index iz=0, const Index ia=0)
Vector view to K(1, 3) elements.
Definition: propagationmatrix.h:840
PropagationMatrix::Data
Tensor4 & Data()
Get full view to data.
Definition: propagationmatrix.h:928
PropagationMatrix::K13
VectorView K13(const Index iz=0, const Index ia=0)
Vector view to K(0, 2) elements.
Definition: propagationmatrix.h:810
PropagationMatrix::Kjj
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
Definition: propagationmatrix.h:790
PropagationMatrix::StokesDimensions
Index StokesDimensions() const
The stokes dimension of the propagation matrix.
Definition: propagationmatrix.h:221
PropagationMatrix::K34
VectorView K34(const Index iz=0, const Index ia=0)
Vector view to K(2, 3) elements.
Definition: propagationmatrix.h:850
PropagationMatrix::NumberOfAzimuthAngles
Index NumberOfAzimuthAngles() const
The number of azimuth angles of the propagation matrix.
Definition: propagationmatrix.h:230
PropagationMatrix::K12
VectorView K12(const Index iz=0, const Index ia=0)
Vector view to K(0, 1) elements.
Definition: propagationmatrix.h:800
PropagationMatrix::K23
VectorView K23(const Index iz=0, const Index ia=0)
Vector view to K(1, 2) elements.
Definition: propagationmatrix.h:830
PropagationMatrix::K14
VectorView K14(const Index iz=0, const Index ia=0)
Vector view to K(0, 3) elements.
Definition: propagationmatrix.h:820
PropagationMatrix::IsEmpty
bool IsEmpty() const
Asks if the class is empty.
Definition: propagationmatrix.h:235
PropagationMatrix::NumberOfZenithAngles
Index NumberOfZenithAngles() const
The number of zenith angles of the propagation matrix.
Definition: propagationmatrix.h:227
PropagationMatrix::IsRotational
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.
Definition: propagationmatrix.h:261
StokesVector
Stokes vector is as Propagation matrix but only has 4 possible values.
Definition: propagationmatrix.h:1061
Tensor3
The Tensor3 class.
Definition: matpackIII.h:352
Vector
The Vector class.
Definition: matpackI.h:910
double_imanip
Input manipulator class for doubles to enable nan and inf parsing.
Definition: double_imanip.h:42
ARTS_ASSERT
#define ARTS_ASSERT(condition,...)
Definition: debug.h:102
ARTS_USER_ERROR_IF
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:153
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
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
real_val
Numeric & real_val(Complex &c) noexcept
Return a reference to the real value of c.
Definition: matpack_complex.h:40
imag
constexpr Numeric imag(Complex c) noexcept
imag
Definition: matpack_complex.h:87
joker
const Joker joker
a2
#define a2
Definition: matpack_complex.h:74
Complex
std::complex< Numeric > Complex
Definition: matpack_complex.h:36
b2
#define b2
Definition: matpack_complex.h:75
real
constexpr Numeric real(Complex c) noexcept
real
Definition: matpack_complex.h:84
Constant::inv_sqrt_2
constexpr Numeric inv_sqrt_2
Inverse of the square root of 2.
Definition: arts_constants.h:108
N
#define N
Definition: rng.cc:164
M
#define M
Definition: rng.cc:165
RadiationVector
Radiation Vector for Stokes dimension 1-4.
Definition: transmissionmatrix.h:498
RadiationVector::R2
std::vector< Eigen::Vector2d > R2
Definition: transmissionmatrix.h:502
RadiationVector::RadiationVector
RadiationVector(Index nf=0, Index stokes=1)
Construct a new Radiation Vector object.
Definition: transmissionmatrix.cc:204
RadiationVector::Frequencies
Index Frequencies() const
Get frequency count.
Definition: transmissionmatrix.cc:567
RadiationVector::Vec1
const Eigen::Matrix< double, 1, 1 > & Vec1(size_t i) const
Return Vector at position.
Definition: transmissionmatrix.cc:216
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:523
RadiationVector::addMultiplied
void addMultiplied(const TransmissionMatrix &A, const RadiationVector &x)
Add multiply.
Definition: transmissionmatrix.cc:409
RadiationVector::Vec4
const Eigen::Vector4d & Vec4(size_t i) const
Return Vector at position.
Definition: transmissionmatrix.cc:213
RadiationVector::leftMul
void leftMul(const TransmissionMatrix &T)
Multiply radiation vector from the left.
Definition: transmissionmatrix.cc:238
RadiationVector::operator()
const Numeric & operator()(const Index i, const Index j) const
Access operator.
Definition: transmissionmatrix.cc:498
RadiationVector::SetZero
void SetZero(size_t i)
Set Radiation Vector to Zero at position.
Definition: transmissionmatrix.cc:245
RadiationVector::rem_avg
void rem_avg(const RadiationVector &O1, const RadiationVector &O2)
Remove the average of two other RadiationVector from *this.
Definition: transmissionmatrix.cc:279
RadiationVector::R1
std::vector< Eigen::Matrix< double, 1, 1 > > R1
Definition: transmissionmatrix.h:503
RadiationVector::add_avg
void add_avg(const RadiationVector &O1, const RadiationVector &O2)
Add the average of two other RadiationVector to *this.
Definition: transmissionmatrix.cc:291
RadiationVector::R4
std::vector< Eigen::Vector4d > R4
Definition: transmissionmatrix.h:500
RadiationVector::operator+=
RadiationVector & operator+=(const RadiationVector &rv)
Addition operator.
Definition: transmissionmatrix.cc:484
RadiationVector::Vec
Eigen::VectorXd Vec(size_t i) const
Return Vector at position by copy.
Definition: transmissionmatrix.cc:266
RadiationVector::addDerivTransmission
void addDerivTransmission(const TransmissionMatrix &PiT, const TransmissionMatrix &dT, const RadiationVector &I)
Add the transmission derivative to this.
Definition: transmissionmatrix.cc:396
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:417
RadiationVector::Vec3
const Eigen::Vector3d & Vec3(size_t i) const
Return Vector at position.
Definition: transmissionmatrix.cc:214
RadiationVector::Vec2
const Eigen::Vector2d & Vec2(size_t i) const
Return Vector at position.
Definition: transmissionmatrix.cc:215
RadiationVector::R3
std::vector< Eigen::Vector3d > R3
Definition: transmissionmatrix.h:501
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:366
RadiationVector::setBackscatterTransmission
void setBackscatterTransmission(const RadiationVector &I0, const TransmissionMatrix &Tr, const TransmissionMatrix &Tf, const TransmissionMatrix &Z)
Set this to backscatter transmission.
Definition: transmissionmatrix.cc:432
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:347
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:303
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:446
TransmissionMatrix
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
Definition: transmissionmatrix.h:166
TransmissionMatrix::mul_aliased
void mul_aliased(const TransmissionMatrix &A, const TransmissionMatrix &B)
Set this to a multiple of A by B.
Definition: transmissionmatrix.cc:113
TransmissionMatrix::T1
std::vector< Eigen::Matrix< double, 1, 1 > > T1
Definition: transmissionmatrix.h:171
TransmissionMatrix::T4
std::vector< Eigen::Matrix4d > T4
Definition: transmissionmatrix.h:168
TransmissionMatrix::operator*=
TransmissionMatrix & operator*=(const Numeric &scale)
Scale self.
Definition: transmissionmatrix.cc:182
TransmissionMatrix::Mat2
const Eigen::Matrix2d & Mat2(size_t i) const
Get Matrix at position.
Definition: transmissionmatrix.cc:45
TransmissionMatrix::Mat
Eigen::MatrixXd Mat(size_t i) const
Get Matrix at position by copy.
Definition: transmissionmatrix.cc:77
TransmissionMatrix::TraMat
auto & TraMat(size_t i) noexcept
Simple template access for the transmission.
Definition: transmissionmatrix.h:358
TransmissionMatrix::Frequencies
Index Frequencies() const
Number of frequencies.
Definition: transmissionmatrix.cc:149
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:443
TransmissionMatrix::mul
void mul(const TransmissionMatrix &A, const TransmissionMatrix &B)
Set this to a multiple of A by B.
Definition: transmissionmatrix.cc:105
TransmissionMatrix::T3
std::vector< Eigen::Matrix3d > T3
Definition: transmissionmatrix.h:169
TransmissionMatrix::Mat4
const Eigen::Matrix4d & Mat4(size_t i) const
Get Matrix at position.
Definition: transmissionmatrix.cc:43
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:410
TransmissionMatrix::setIdentity
void setIdentity()
Set to identity matrix.
Definition: transmissionmatrix.cc:91
TransmissionMatrix::Mat1
const Eigen::Matrix< double, 1, 1 > & Mat1(size_t i) const
Get Matrix at position.
Definition: transmissionmatrix.cc:46
TransmissionMatrix::T2
std::vector< Eigen::Matrix2d > T2
Definition: transmissionmatrix.h:170
TransmissionMatrix::Mat3
const Eigen::Matrix3d & Mat3(size_t i) const
Get Matrix at position.
Definition: transmissionmatrix.cc:44
TransmissionMatrix::operator()
Numeric operator()(const Index i, const Index j, const Index k) const
Access value in matrix.
Definition: transmissionmatrix.cc:121
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:169
TransmissionMatrix::setZero
void setZero()
Set to zero matrix.
Definition: transmissionmatrix.cc:98
TransmissionMatrix::TransmissionMatrix
TransmissionMatrix(Index nf=0, Index stokes=1)
Construct a new Transmission Matrix object.
Definition: transmissionmatrix.cc:34
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:707
lower_is_considered_zero_for_sinc_likes
constexpr Numeric lower_is_considered_zero_for_sinc_likes
Definition: transmissionmatrix.cc:580
bulk_backscatter
ArrayOfTransmissionMatrix bulk_backscatter(const ConstTensor5View &Pe, const ConstMatrixView &pnd)
Bulk back-scattering.
Definition: transmissionmatrix.cc:2188
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:1756
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:1182
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:978
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:770
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:691
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:892
bulk_backscatter_derivative
ArrayOfArrayOfTransmissionMatrix bulk_backscatter_derivative(const ConstTensor5View &Pe, const ArrayOfMatrix &dpnd_dx)
Derivatives of bulk back-scattering
Definition: transmissionmatrix.cc:2230
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:583
operator*
LazyScale< TransmissionMatrix > operator*(const TransmissionMatrix &tm, const Numeric &x)
Lazy scale of Transmission Matrix.
Definition: transmissionmatrix.cc:194
cumulative_transmission
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
Definition: transmissionmatrix.cc:2001
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:2022
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:1803
operator>>
std::istream & operator>>(std::istream &is, TransmissionMatrix &tm)
Definition: transmissionmatrix.cc:2301
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:1943
transmat
void transmat(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r) noexcept
Definition: transmissionmatrix.cc:1736
operator<<
std::ostream & operator<<(std::ostream &os, const TransmissionMatrix &tm)
Definition: transmissionmatrix.cc:2284
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:1785
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:927
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:680
transmissionmatrix.h
Stuff related to the transmission matrix.
RadiativeTransferSolver
RadiativeTransferSolver
Intended to hold various forward solvers.
Definition: transmissionmatrix.h:809
u
#define u
d
#define d
BackscatterSolver
BackscatterSolver
Intended to hold various backscatter solvers.
Definition: transmissionmatrix.h:797
CumulativeTransmission
CumulativeTransmission
Intended to hold various ways to accumulate the transmission matrix.
Definition: transmissionmatrix.h:803
v
#define v
w
#define w
a
#define a
ArrayOfTransmissionMatrix
Array< TransmissionMatrix > ArrayOfTransmissionMatrix
Definition: transmissionmatrix.h:787
c
#define c
b
#define b