ARTS 2.5.4 (git: 4c0d3b4d)
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#include "constants.h"
31#include "matpack_complex.h"
32
34
35template <int N> constexpr
36Eigen::Matrix<Numeric, N, 1> source_vector(const StokesVector& a,
37 const ConstVectorView& B,
38 const StokesVector& da,
39 const ConstVectorView& dB_dT,
40 const StokesVector& dS,
41 bool dT,
42 Index i) {
43 static_assert(N>0 and N<5, "Bad stokes dimensions");
44
45 if constexpr (N == 1) {
46 if (dT)
47 return Eigen::Matrix<Numeric, 1, 1>(dS.Kjj()[i] + da.Kjj()[i] * B[i] + a.Kjj()[i] * dB_dT[i]);
48 return Eigen::Matrix<Numeric, 1, 1>(dS.Kjj()[i] + da.Kjj()[i] * B[i]);
49 }
50
51 if constexpr (N == 2) {
52 if (dT)
53 return Eigen::Vector2d(dS.Kjj()[i] + da.Kjj()[i] * B[i] + a.Kjj()[i] * dB_dT[i],
54 dS.K12()[i] + da.K12()[i] * B[i] + a.K12()[i] * dB_dT[i]);
55 return Eigen::Vector2d(dS.Kjj()[i] + da.Kjj()[i] * B[i],
56 dS.K12()[i] + da.K12()[i] * B[i]);
57 }
58
59 if constexpr (N == 3) {
60 if (dT)
61 return Eigen::Vector3d(dS.Kjj()[i] + da.Kjj()[i] * B[i] + a.Kjj()[i] * dB_dT[i],
62 dS.K12()[i] + da.K12()[i] * B[i] + a.K12()[i] * dB_dT[i],
63 dS.K13()[i] + da.K13()[i] * B[i] + a.K13()[i] * dB_dT[i]);
64 return Eigen::Vector3d(dS.Kjj()[i] + da.Kjj()[i] * B[i],
65 dS.K12()[i] + da.K12()[i] * B[i],
66 dS.K13()[i] + da.K13()[i] * B[i]);
67 }
68
69 if constexpr (N == 4) {
70 if (dT)
71 return Eigen::Vector4d(dS.Kjj()[i] + da.Kjj()[i] * B[i] + a.Kjj()[i] * dB_dT[i],
72 dS.K12()[i] + da.K12()[i] * B[i] + a.K12()[i] * dB_dT[i],
73 dS.K13()[i] + da.K13()[i] * B[i] + a.K13()[i] * dB_dT[i],
74 dS.K14()[i] + da.K14()[i] * B[i] + a.K14()[i] * dB_dT[i]);
75 return Eigen::Vector4d(dS.Kjj()[i] + da.Kjj()[i] * B[i],
76 dS.K12()[i] + da.K12()[i] * B[i],
77 dS.K13()[i] + da.K13()[i] * B[i],
78 dS.K14()[i] + da.K14()[i] * B[i]);
79 }
80}
81
82template <int N> constexpr
83Eigen::Matrix<Numeric, N, 1> source_vector(const StokesVector& a,
84 const ConstVectorView& B,
85 const StokesVector& da,
86 const ConstVectorView& dB_dT,
87 bool dT,
88 Index i) {
89 static_assert(N>0 and N<5, "Bad stokes dimensions");
90
91 if constexpr (N == 1) {
92 if (dT)
93 return Eigen::Matrix<Numeric, 1, 1>(da.Kjj()[i] * B[i] + a.Kjj()[i] * dB_dT[i]);
94 return Eigen::Matrix<Numeric, 1, 1>(da.Kjj()[i] * B[i]);
95 }
96
97 if constexpr (N == 2) {
98 if (dT)
99 return Eigen::Vector2d(da.Kjj()[i] * B[i] + a.Kjj()[i] * dB_dT[i],
100 da.K12()[i] * B[i] + a.K12()[i] * dB_dT[i]);
101 return Eigen::Vector2d(da.Kjj()[i] * B[i],
102 da.K12()[i] * B[i]);
103 }
104
105 if constexpr (N == 3) {
106 if (dT)
107 return Eigen::Vector3d(da.Kjj()[i] * B[i] + a.Kjj()[i] * dB_dT[i],
108 da.K12()[i] * B[i] + a.K12()[i] * dB_dT[i],
109 da.K13()[i] * B[i] + a.K13()[i] * dB_dT[i]);
110 return Eigen::Vector3d(da.Kjj()[i] * B[i],
111 da.K12()[i] * B[i],
112 da.K13()[i] * B[i]);
113 }
114
115 if constexpr (N == 4) {
116 if (dT)
117 return Eigen::Vector4d(da.Kjj()[i] * B[i] + a.Kjj()[i] * dB_dT[i],
118 da.K12()[i] * B[i] + a.K12()[i] * dB_dT[i],
119 da.K13()[i] * B[i] + a.K13()[i] * dB_dT[i],
120 da.K14()[i] * B[i] + a.K14()[i] * dB_dT[i]);
121 return Eigen::Vector4d(da.Kjj()[i] * B[i],
122 da.K12()[i] * B[i],
123 da.K13()[i] * B[i],
124 da.K14()[i] * B[i]);
125 }
126}
127
129 const PropagationMatrix& K1,
130 const PropagationMatrix& K2,
131 const Numeric& r,
132 const Index iz = 0,
133 const Index ia = 0) noexcept {
134 for (Index i = 0; i < K1.NumberOfFrequencies(); i++)
135 T.Mat1(i)(0, 0) =
136 std::exp(-0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]));
137}
138
140 const PropagationMatrix& K1,
141 const PropagationMatrix& K2,
142 const Numeric& r,
143 const Index iz = 0,
144 const Index ia = 0) noexcept {
145 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
146 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
147 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]);
148 const Numeric exp_a = std::exp(a);
149 const Numeric cb = std::cosh(b), sb = std::sinh(b);
150 T.Mat2(i).noalias() =
151 (Eigen::Matrix2d() << cb, sb, sb, cb).finished() * exp_a;
152 }
153}
154
156 const PropagationMatrix& K1,
157 const PropagationMatrix& K2,
158 const Numeric& r,
159 const Index iz = 0,
160 const Index ia = 0) noexcept {
161 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
162 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
163 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]),
164 c = -0.5 * r * (K1.K13(iz, ia)[i] + K2.K13(iz, ia)[i]),
165 u = -0.5 * r * (K1.K23(iz, ia)[i] + K2.K23(iz, ia)[i]);
166 const Numeric exp_a = std::exp(a);
167
168 if (b == 0. and c == 0. and u == 0.) {
169 T.Mat3(i).noalias() = Eigen::Matrix3d::Identity() * exp_a;
170 } else {
171 const Numeric a2 = a * a, b2 = b * b, c2 = c * c, u2 = u * u;
172 const Numeric Const = b2 + c2 - u2;
173
174 const bool real = Const > 0;
175 const bool imag = Const < 0;
176 const bool either = real or imag;
177
178 const Numeric x =
179 std::sqrt(imag ? -Const : Const); // test to just use real values
180 const Numeric x2 =
181 (real ? 1 : -1) * x * x; // test to change sign if imaginary
182 const Numeric inv_x2 =
183 either ? 1.0 / x2
184 : 1.0; // test so further calculations are replaced as x→0
185
186 const Numeric sx =
187 real ? std::sinh(x) : std::sin(-x); // -i sin(ix) → sinh(x)
188 const Numeric cx =
189 real ? std::cosh(x) : std::cos(+x); // cos(ix) → cosh(x)
190
191 /* Using:
192 * lim x→0 [(cosh(x) - 1) / x^2] → 1/2
193 * lim x→0 [sinh(x) / x] → 1
194 * inv_x2 := 1 for x == 0,
195 * C0, C1, C2 ∝ [1/x^2]
196 */
197 const Numeric C0 =
198 either ? a2 * (cx - 1.0) - a * x * sx + x2 : 1.0 + 0.5 * a2 - a;
199 const Numeric C1 = either ? 2.0 * a * (1.0 - cx) + x * sx : 1.0 - a;
200 const Numeric C2 = either ? cx - 1.0 : 0.5;
201
202 T.Mat3(i).noalias() =
203 exp_a * inv_x2 *
204 (Eigen::Matrix3d() << C0 + C1 * a + C2 * (a2 + b2 + c2),
205 C1 * b + C2 * (2 * a * b - c * u),
206 C1 * c + C2 * (2 * a * c + b * u),
207 C1 * b + C2 * (2 * a * b + c * u),
208 C0 + C1 * a + C2 * (a2 + b2 - u2),
209 C1 * u + C2 * (2 * a * u + b * c),
210 C1 * c + C2 * (2 * a * c - b * u),
211 -C1 * u - C2 * (2 * a * u - b * c),
212 C0 + C1 * a + C2 * (a2 + c2 - u2))
213 .finished();
214 }
215 }
216}
217
219 const PropagationMatrix& K1,
220 const PropagationMatrix& K2,
221 const Numeric& r,
222 const Index iz = 0,
223 const Index ia = 0) noexcept {
224 static constexpr Numeric sqrt_05 = Constant::inv_sqrt_2;
225 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
226 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
227 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]),
228 c = -0.5 * r * (K1.K13(iz, ia)[i] + K2.K13(iz, ia)[i]),
229 d = -0.5 * r * (K1.K14(iz, ia)[i] + K2.K14(iz, ia)[i]),
230 u = -0.5 * r * (K1.K23(iz, ia)[i] + K2.K23(iz, ia)[i]),
231 v = -0.5 * r * (K1.K24(iz, ia)[i] + K2.K24(iz, ia)[i]),
232 w = -0.5 * r * (K1.K34(iz, ia)[i] + K2.K34(iz, ia)[i]);
233 const Numeric exp_a = std::exp(a);
234
235 if (b == 0. and c == 0. and d == 0. and u == 0. and v == 0. and w == 0.)
236 T.Mat4(i).noalias() = Eigen::Matrix4d::Identity() * exp_a;
237 else {
238 const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
239 w2 = w * w;
240
241 const Numeric tmp =
242 w2 * w2 + 2 * (b2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
243 c2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
244 d2 * (d2 * 0.5 + u2 - v2 - w2) +
245 u2 * (u2 * 0.5 + v2 + w2) + v2 * (v2 * 0.5 + w2) +
246 4 * (b * d * u * w - b * c * v * w - c * d * u * v));
247 const Complex Const1 = std::sqrt(Complex(tmp, 0));
248 const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
249
250 const Complex x = std::sqrt(Const2 + Const1) * sqrt_05;
251 const Complex y = std::sqrt(Const2 - Const1) * sqrt_05 * Complex(0, 1);
252 const Complex x2 = x * x;
253 const Complex y2 = y * y;
254 const Complex cy = std::cos(y);
255 const Complex sy = std::sin(y);
256 const Complex cx = std::cosh(x);
257 const Complex sx = std::sinh(x);
258
259 const bool x_zero = std::abs(x) < lower_is_considered_zero_for_sinc_likes;
260 const bool y_zero = std::abs(y) < lower_is_considered_zero_for_sinc_likes;
261 const bool both_zero = y_zero and x_zero;
262 const bool either_zero = y_zero or x_zero;
263
264 /* Using:
265 * lim x→0 [({cosh(x),cos(x)} - 1) / x^2] → 1/2
266 * lim x→0 [{sinh(x),sin(x)} / x] → 1
267 * inv_x2 := 1 for x == 0,
268 * -i sin(ix) → sinh(x)
269 * cos(ix) → cosh(x)
270 * C0, C1, C2 ∝ [1/x^2]
271 */
272 const Complex ix = x_zero ? 0.0 : 1.0 / x;
273 const Complex iy = y_zero ? 0.0 : 1.0 / y;
274 const Complex inv_x2y2 =
275 both_zero
276 ? 1.0
277 : 1.0 /
278 (x2 + y2); // The first "1.0" is the trick for above limits
279
280 const Numeric C0 =
281 either_zero ? 1.0 : ((cy * x2 + cx * y2) * inv_x2y2).real();
282 const Numeric C1 =
283 either_zero ? 1.0 : ((sy * x2 * iy + sx * y2 * ix) * inv_x2y2).real();
284 const Numeric C2 = both_zero ? 0.5 : ((cx - cy) * inv_x2y2).real();
285 const Numeric C3 =
286 both_zero ? 1.0 / 6.0
287 : ((x_zero ? 1.0 - sy * iy
288 : y_zero ? sx * ix - 1.0 : sx * ix - sy * iy) *
289 inv_x2y2)
290 .real();
291 T.Mat4(i).noalias() =
292 exp_a * (Eigen::Matrix4d() << C0 + C2 * (b2 + c2 + d2),
293 C1 * b + C2 * (-c * u - d * v) +
294 C3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) -
295 v * (b * v + c * w)),
296 C1 * c + C2 * (b * u - d * w) +
297 C3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) -
298 w * (b * v + c * w)),
299 C1 * d + C2 * (b * v + c * w) +
300 C3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) +
301 w * (b * u - d * w)),
302
303 C1 * b + C2 * (c * u + d * v) +
304 C3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) +
305 d * (b * d + u * w)),
306 C0 + C2 * (b2 - u2 - v2),
307 C2 * (b * c - v * w) + C1 * u +
308 C3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
309 w * (b * d + u * w)),
310 C2 * (b * d + u * w) + C1 * v +
311 C3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
312 w * (b * c - v * w)),
313
314 C1 * c + C2 * (-b * u + d * w) +
315 C3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) +
316 d * (c * d - u * v)),
317 C2 * (b * c - v * w) - C1 * u +
318 C3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
319 v * (c * d - u * v)),
320 C0 + C2 * (c2 - u2 - w2),
321 C2 * (c * d - u * v) + C1 * w +
322 C3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
323 w * (-c2 + u2 + w2)),
324
325 C1 * d + C2 * (-b * v - c * w) +
326 C3 * (b * (b * d + u * w) + c * (c * d - u * v) -
327 d * (-d2 + v2 + w2)),
328 C2 * (b * d + u * w) - C1 * v +
329 C3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
330 v * (-d2 + v2 + w2)),
331 C2 * (c * d - u * v) - C1 * w +
332 C3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
333 w * (-d2 + v2 + w2)),
334 C0 + C2 * (d2 - v2 - w2))
335 .finished();
336 }
337 }
338}
339
343 const PropagationMatrix& K1,
344 const PropagationMatrix& K2,
345 const ArrayOfPropagationMatrix& dK1,
346 const ArrayOfPropagationMatrix& dK2,
347 const Numeric& r,
348 const Numeric& dr_dT1,
349 const Numeric& dr_dT2,
350 const Index it,
351 const Index iz,
352 const Index ia) noexcept {
353 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
354 T.Mat1(i)(0, 0) =
355 std::exp(-0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]));
356 for (Index j = 0; j < dT1.nelem(); j++) {
357 if (dK1[j].NumberOfFrequencies())
358 dT1[j].Mat1(i)(0, 0) =
359 T.Mat1(i)(0, 0) *
360 (-0.5 *
361 (r * dK1[j].Kjj(iz, ia)[i] +
362 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
363 : 0.0)));
364 if (dK2[j].NumberOfFrequencies())
365 dT2[j].Mat1(i)(0, 0) =
366 T.Mat1(i)(0, 0) *
367 (-0.5 *
368 (r * dK2[j].Kjj(iz, ia)[i] +
369 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
370 : 0.0)));
371 }
372 }
373}
374
378 const PropagationMatrix& K1,
379 const PropagationMatrix& K2,
380 const ArrayOfPropagationMatrix& dK1,
381 const ArrayOfPropagationMatrix& dK2,
382 const Numeric& r,
383 const Numeric& dr_dT1,
384 const Numeric& dr_dT2,
385 const Index it,
386 const Index iz,
387 const Index ia) noexcept {
388 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
389 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
390 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]);
391 const Numeric exp_a = std::exp(a);
392 const Numeric cb = std::cosh(b), sb = std::sinh(b);
393 T.Mat2(i).noalias() =
394 (Eigen::Matrix2d() << cb, sb, sb, cb).finished() * exp_a;
395 for (Index j = 0; j < dT1.nelem(); j++) {
396 if (dK1[j].NumberOfFrequencies()) {
397 const Numeric da = -0.5 * (r * dK1[j].Kjj(iz, ia)[i] +
398 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] +
399 K2.Kjj(iz, ia)[i])
400 : 0.0)),
401 db = -0.5 * (r * dK1[j].K12(iz, ia)[i] +
402 ((j == it) ? dr_dT1 * (K1.K12(iz, ia)[i] +
403 K2.K12(iz, ia)[i])
404 : 0.0));
405 dT1[j].Mat2(i).noalias() =
406 T.Mat2(i) * da +
407 (Eigen::Matrix2d() << sb, cb, cb, sb).finished() * exp_a * db;
408 }
409 if (dK2[j].NumberOfFrequencies()) {
410 const Numeric da = -0.5 * (r * dK2[j].Kjj(iz, ia)[i] +
411 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] +
412 K2.Kjj(iz, ia)[i])
413 : 0.0)),
414 db = -0.5 * (r * dK2[j].K12(iz, ia)[i] +
415 ((j == it) ? dr_dT2 * (K1.K12(iz, ia)[i] +
416 K2.K12(iz, ia)[i])
417 : 0.0));
418 dT2[j].Mat2(i).noalias() =
419 T.Mat2(i) * da +
420 (Eigen::Matrix2d() << sb, cb, cb, sb).finished() * exp_a * db;
421 }
422 }
423 }
424}
425
429 const PropagationMatrix& K1,
430 const PropagationMatrix& K2,
431 const ArrayOfPropagationMatrix& dK1,
432 const ArrayOfPropagationMatrix& dK2,
433 const Numeric& r,
434 const Numeric& dr_dT1,
435 const Numeric& dr_dT2,
436 const Index it,
437 const Index iz,
438 const Index ia) noexcept {
439 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
440 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
441 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]),
442 c = -0.5 * r * (K1.K13(iz, ia)[i] + K2.K13(iz, ia)[i]),
443 u = -0.5 * r * (K1.K23(iz, ia)[i] + K2.K23(iz, ia)[i]);
444 const Numeric exp_a = std::exp(a);
445
446 if (b == 0. and c == 0. and u == 0.) {
447 T.Mat3(i).noalias() = Eigen::Matrix3d::Identity() * exp_a;
448 for (Index j = 0; j < dT1.nelem(); j++) {
449 if (dK1[j].NumberOfFrequencies())
450 dT1[j].Mat3(i).noalias() =
451 T.Mat3(i) *
452 (-0.5 *
453 (r * dK1[j].Kjj(iz, ia)[i] +
454 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
455 : 0.0)));
456 if (dK2[j].NumberOfFrequencies())
457 dT2[j].Mat3(i).noalias() =
458 T.Mat3(i) *
459 (-0.5 *
460 (r * dK2[j].Kjj(iz, ia)[i] +
461 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
462 : 0.0)));
463 }
464 } else {
465 const Numeric a2 = a * a, b2 = b * b, c2 = c * c, u2 = u * u;
466 const Numeric Const = b2 + c2 - u2;
467
468 const bool real = Const > 0;
469 const bool imag = Const < 0;
470 const bool either = real or imag;
471
472 const Numeric x =
473 std::sqrt(imag ? -Const : Const); // test to just use real values
474 const Numeric x2 =
475 (real ? 1 : -1) * x * x; // test to change sign if imaginary
476 const Numeric inv_x =
477 either ? 1.0 / x
478 : 1.0; // test so further calculations are replaced as x→0
479 const Numeric inv_x2 = inv_x * inv_x;
480
481 const Numeric sx =
482 real ? std::sinh(x) : std::sin(-x); // -i sin(ix) → sinh(x)
483 const Numeric cx =
484 real ? std::cosh(x) : std::cos(+x); // cos(ix) → cosh(x)
485
486 /* Using:
487 * lim x→0 [(cosh(x) - 1) / x^2] → 1/2
488 * lim x→0 [sinh(x) / x] → 1
489 * inv_x2 := 1 for x == 0,
490 * C0, C1, C2 ∝ [1/x^2]
491 */
492 const Numeric C0 =
493 either ? a2 * (cx - 1.0) - a * x * sx + x2 : 1.0 + 0.5 * a2 - a;
494 const Numeric C1 = either ? 2.0 * a * (1.0 - cx) + x * sx : 1.0 - a;
495 const Numeric C2 = either ? cx - 1.0 : 0.5;
496
497 T.Mat3(i).noalias() =
498 exp_a * inv_x2 *
499 (Eigen::Matrix3d() << C0 + C1 * a + C2 * (a2 + b2 + c2),
500 C1 * b + C2 * (2 * a * b - c * u),
501 C1 * c + C2 * (2 * a * c + b * u),
502 C1 * b + C2 * (2 * a * b + c * u),
503 C0 + C1 * a + C2 * (a2 + b2 - u2),
504 C1 * u + C2 * (2 * a * u + b * c),
505 C1 * c + C2 * (2 * a * c - b * u),
506 -C1 * u - C2 * (2 * a * u - b * c),
507 C0 + C1 * a + C2 * (a2 + c2 - u2))
508 .finished();
509
510 for (Index j = 0; j < dT1.nelem(); j++) {
511 if (dK1[j].NumberOfFrequencies()) {
512 const Numeric da = -0.5 * (r * dK1[j].Kjj(iz, ia)[i] +
513 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] +
514 K2.Kjj(iz, ia)[i])
515 : 0.0)),
516 db = -0.5 * (r * dK1[j].K12(iz, ia)[i] +
517 ((j == it) ? dr_dT1 * (K1.K12(iz, ia)[i] +
518 K2.K12(iz, ia)[i])
519 : 0.0)),
520 dc = -0.5 * (r * dK1[j].K13(iz, ia)[i] +
521 ((j == it) ? dr_dT1 * (K1.K13(iz, ia)[i] +
522 K2.K13(iz, ia)[i])
523 : 0.0)),
524 du = -0.5 * (r * dK1[j].K23(iz, ia)[i] +
525 ((j == it) ? dr_dT1 * (K1.K23(iz, ia)[i] +
526 K2.K23(iz, ia)[i])
527 : 0.0));
528 const Numeric da2 = 2 * a * da, db2 = 2 * b * db, dc2 = 2 * c * dc,
529 du2 = 2 * u * du;
530 const Numeric dx = either ? ((db2 + dc2 - du2) * inv_x * 0.5) : 0,
531 dx2 = 2 * x * dx;
532 const Numeric dsx = (real ? 1 : -1) * cx * dx;
533 const Numeric dcx = sx * dx;
534
535 const Numeric dC0 = either ? da2 * (cx - 1.0) + da2 * (dcx - 1.0) -
536 da * x * sx - a * dx * sx -
537 a * x * dsx + dx2
538 : 0.5 * da2 - da;
539 const Numeric dC1 =
540 either ? 2.0 * (da * (1.0 - cx) - a * dcx) + dx * sx + x * dsx
541 : -da;
542 const Numeric dC2 = either ? dcx : 0;
543
544 dT1[j].Mat3(i).noalias() =
545 T.Mat3(i) * (da + dx2 * inv_x2) +
546 exp_a * inv_x2 *
547 (Eigen::Matrix3d() << dC0 + dC1 * a + C1 * da +
548 dC2 * (a2 + b2 + c2) +
549 C2 * (da2 + db2 + dc2),
550 dC1 * b + C1 * db + dC2 * (2 * a * b - c * u) +
551 C2 * (2 * da * b - dc * u + 2 * a * db - c * du),
552 dC1 * c + C1 * dc + dC2 * (2 * a * c + b * u) +
553 C2 * (2 * da * c + db * u + 2 * a * dc + b * du),
554 dC1 * b + C1 * db + dC2 * (2 * a * b + c * u) +
555 C2 * (2 * da * b + dc * u + 2 * a * db + c * du),
556 dC0 + dC1 * a + C1 * da + dC2 * (a2 + b2 - u2) +
557 C2 * (da2 + db2 - du2),
558 dC1 * u + C1 * du + dC2 * (2 * a * u + b * c) +
559 C2 * (2 * da * u + db * c + 2 * a * du + b * dc),
560 dC1 * c + C1 * dc + dC2 * (2 * a * c - b * u) +
561 C2 * (2 * da * c - db * u + 2 * a * dc - b * du),
562 -dC1 * u - C1 * du - dC2 * (2 * a * u - b * c) -
563 C2 * (2 * da * u - db * c + 2 * a * du - b * dc),
564 dC0 + dC1 * a + C1 * da + dC2 * (a2 + c2 - u2) +
565 C2 * (da2 + dc2 - du2))
566 .finished();
567 }
568 if (dK2[j].NumberOfFrequencies()) {
569 const Numeric da = -0.5 * (r * dK2[j].Kjj(iz, ia)[i] +
570 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] +
571 K2.Kjj(iz, ia)[i])
572 : 0.0)),
573 db = -0.5 * (r * dK2[j].K12(iz, ia)[i] +
574 ((j == it) ? dr_dT2 * (K1.K12(iz, ia)[i] +
575 K2.K12(iz, ia)[i])
576 : 0.0)),
577 dc = -0.5 * (r * dK2[j].K13(iz, ia)[i] +
578 ((j == it) ? dr_dT2 * (K1.K13(iz, ia)[i] +
579 K2.K13(iz, ia)[i])
580 : 0.0)),
581 du = -0.5 * (r * dK2[j].K23(iz, ia)[i] +
582 ((j == it) ? dr_dT2 * (K1.K23(iz, ia)[i] +
583 K2.K23(iz, ia)[i])
584 : 0.0));
585 const Numeric da2 = 2 * a * da, db2 = 2 * b * db, dc2 = 2 * c * dc,
586 du2 = 2 * u * du;
587 const Numeric dx = either ? ((db2 + dc2 - du2) * inv_x * 0.5) : 0,
588 dx2 = 2 * x * dx;
589 const Numeric dsx = (real ? 1 : -1) * cx * dx;
590 const Numeric dcx = sx * dx;
591
592 const Numeric dC0 = either ? da2 * (cx - 1.0) + da2 * (dcx - 1.0) -
593 da * x * sx - a * dx * sx -
594 a * x * dsx + dx2
595 : 0.5 * da2 - da;
596 const Numeric dC1 =
597 either ? 2.0 * (da * (1.0 - cx) - a * dcx) + dx * sx + x * dsx
598 : -da;
599 const Numeric dC2 = either ? dcx : 0;
600
601 dT2[j].Mat3(i).noalias() =
602 T.Mat3(i) * (da + dx2 * inv_x2) +
603 exp_a * inv_x2 *
604 (Eigen::Matrix3d() << dC0 + dC1 * a + C1 * da +
605 dC2 * (a2 + b2 + c2) +
606 C2 * (da2 + db2 + dc2),
607 dC1 * b + C1 * db + dC2 * (2 * a * b - c * u) +
608 C2 * (2 * da * b - dc * u + 2 * a * db - c * du),
609 dC1 * c + C1 * dc + dC2 * (2 * a * c + b * u) +
610 C2 * (2 * da * c + db * u + 2 * a * dc + b * du),
611 dC1 * b + C1 * db + dC2 * (2 * a * b + c * u) +
612 C2 * (2 * da * b + dc * u + 2 * a * db + c * du),
613 dC0 + dC1 * a + C1 * da + dC2 * (a2 + b2 - u2) +
614 C2 * (da2 + db2 - du2),
615 dC1 * u + C1 * du + dC2 * (2 * a * u + b * c) +
616 C2 * (2 * da * u + db * c + 2 * a * du + b * dc),
617 dC1 * c + C1 * dc + dC2 * (2 * a * c - b * u) +
618 C2 * (2 * da * c - db * u + 2 * a * dc - b * du),
619 -dC1 * u - C1 * du - dC2 * (2 * a * u - b * c) -
620 C2 * (2 * da * u - db * c + 2 * a * du - b * dc),
621 dC0 + dC1 * a + C1 * da + dC2 * (a2 + c2 - u2) +
622 C2 * (da2 + dc2 - du2))
623 .finished();
624 }
625 }
626 }
627 }
628}
629
633 const PropagationMatrix& K1,
634 const PropagationMatrix& K2,
635 const ArrayOfPropagationMatrix& dK1,
636 const ArrayOfPropagationMatrix& dK2,
637 const Numeric& r,
638 const Numeric& dr_dT1,
639 const Numeric& dr_dT2,
640 const Index it,
641 const Index iz,
642 const Index ia) noexcept {
643 static constexpr Numeric sqrt_05 = Constant::inv_sqrt_2;
644 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
645 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
646 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]),
647 c = -0.5 * r * (K1.K13(iz, ia)[i] + K2.K13(iz, ia)[i]),
648 d = -0.5 * r * (K1.K14(iz, ia)[i] + K2.K14(iz, ia)[i]),
649 u = -0.5 * r * (K1.K23(iz, ia)[i] + K2.K23(iz, ia)[i]),
650 v = -0.5 * r * (K1.K24(iz, ia)[i] + K2.K24(iz, ia)[i]),
651 w = -0.5 * r * (K1.K34(iz, ia)[i] + K2.K34(iz, ia)[i]);
652 const Numeric exp_a = std::exp(a);
653
654 if (b == 0. and c == 0. and d == 0. and u == 0. and v == 0. and w == 0.) {
655 T.Mat4(i).noalias() = Eigen::Matrix4d::Identity() * exp_a;
656 for (Index j = 0; j < dK1.nelem(); j++) {
657 if (dK1[j].NumberOfFrequencies())
658 dT1[j].Mat4(i).noalias() =
659 T.Mat4(i) *
660 (-0.5 *
661 (r * dK1[j].Kjj(iz, ia)[i] +
662 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
663 : 0.0)));
664 if (dK2[j].NumberOfFrequencies())
665 dT2[j].Mat4(i).noalias() =
666 T.Mat4(i) *
667 (-0.5 *
668 (r * dK2[j].Kjj(iz, ia)[i] +
669 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
670 : 0.0)));
671 }
672 } else {
673 const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
674 w2 = w * w;
675 const Numeric tmp =
676 w2 * w2 + 2 * (b2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
677 c2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
678 d2 * (d2 * 0.5 + u2 - v2 - w2) +
679 u2 * (u2 * 0.5 + v2 + w2) + v2 * (v2 * 0.5 + w2) +
680 4 * (b * d * u * w - b * c * v * w - c * d * u * v));
681 const Complex Const1 = std::sqrt(Complex(tmp, 0));
682 const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
683 const Complex tmp_x_sqrt = std::sqrt(Const2 + Const1);
684 const Complex tmp_y_sqrt = std::sqrt(Const2 - Const1);
685 const Complex x = tmp_x_sqrt * sqrt_05;
686 const Complex y = tmp_y_sqrt * sqrt_05 * Complex(0, 1);
687 const Complex x2 = x * x;
688 const Complex y2 = y * y;
689 const Complex cy = std::cos(y);
690 const Complex sy = std::sin(y);
691 const Complex cx = std::cosh(x);
692 const Complex sx = std::sinh(x);
693
694 const bool x_zero = std::abs(x) < lower_is_considered_zero_for_sinc_likes;
695 const bool y_zero = std::abs(y) < lower_is_considered_zero_for_sinc_likes;
696 const bool both_zero = y_zero and x_zero;
697 const bool either_zero = y_zero or x_zero;
698
699 /* Using:
700 * lim x→0 [({cosh(x),cos(x)} - 1) / x^2] → 1/2
701 * lim x→0 [{sinh(x),sin(x)} / x] → 1
702 * inv_x2 := 1 for x == 0,
703 * -i sin(ix) → sinh(x)
704 * cos(ix) → cosh(x)
705 * C0, C1, C2 ∝ [1/x^2]
706 */
707 const Complex ix = x_zero ? 0.0 : 1.0 / x;
708 const Complex iy = y_zero ? 0.0 : 1.0 / y;
709 const Complex inv_x2y2 =
710 both_zero
711 ? 1.0
712 : 1.0 /
713 (x2 + y2); // The first "1.0" is the trick for above limits
714 const Complex C0c = either_zero ? 1.0 : (cy * x2 + cx * y2) * inv_x2y2;
715 const Complex C1c =
716 either_zero ? 1.0 : (sy * x2 * iy + sx * y2 * ix) * inv_x2y2;
717 const Complex C2c = both_zero ? 0.5 : (cx - cy) * inv_x2y2;
718 const Complex C3c =
719 both_zero ? 1.0 / 6.0
720 : (x_zero ? 1.0 - sy * iy
721 : y_zero ? sx * ix - 1.0 : sx * ix - sy * iy) *
722 inv_x2y2;
723
724 const Numeric& C0 = real_val(C0c);
725 const Numeric& C1 = real_val(C1c);
726 const Numeric& C2 = real_val(C2c);
727 const Numeric& C3 = real_val(C3c);
728 T.Mat4(i).noalias() =
729 exp_a * (Eigen::Matrix4d() << C0 + C2 * (b2 + c2 + d2),
730 C1 * b + C2 * (-c * u - d * v) +
731 C3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) -
732 v * (b * v + c * w)),
733 C1 * c + C2 * (b * u - d * w) +
734 C3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) -
735 w * (b * v + c * w)),
736 C1 * d + C2 * (b * v + c * w) +
737 C3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) +
738 w * (b * u - d * w)),
739
740 C1 * b + C2 * (c * u + d * v) +
741 C3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) +
742 d * (b * d + u * w)),
743 C0 + C2 * (b2 - u2 - v2),
744 C2 * (b * c - v * w) + C1 * u +
745 C3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
746 w * (b * d + u * w)),
747 C2 * (b * d + u * w) + C1 * v +
748 C3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
749 w * (b * c - v * w)),
750
751 C1 * c + C2 * (-b * u + d * w) +
752 C3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) +
753 d * (c * d - u * v)),
754 C2 * (b * c - v * w) - C1 * u +
755 C3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
756 v * (c * d - u * v)),
757 C0 + C2 * (c2 - u2 - w2),
758 C2 * (c * d - u * v) + C1 * w +
759 C3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
760 w * (-c2 + u2 + w2)),
761
762 C1 * d + C2 * (-b * v - c * w) +
763 C3 * (b * (b * d + u * w) + c * (c * d - u * v) -
764 d * (-d2 + v2 + w2)),
765 C2 * (b * d + u * w) - C1 * v +
766 C3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
767 v * (-d2 + v2 + w2)),
768 C2 * (c * d - u * v) - C1 * w +
769 C3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
770 w * (-d2 + v2 + w2)),
771 C0 + C2 * (d2 - v2 - w2))
772 .finished();
773
774 for (Index j = 0; j < dK1.nelem(); j++) {
775 if (dK1[j].NumberOfFrequencies()) {
776 const Numeric da = -0.5 * (r * dK1[j].Kjj(iz, ia)[i] +
777 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] +
778 K2.Kjj(iz, ia)[i])
779 : 0.0)),
780 db = -0.5 * (r * dK1[j].K12(iz, ia)[i] +
781 ((j == it) ? dr_dT1 * (K1.K12(iz, ia)[i] +
782 K2.K12(iz, ia)[i])
783 : 0.0)),
784 dc = -0.5 * (r * dK1[j].K13(iz, ia)[i] +
785 ((j == it) ? dr_dT1 * (K1.K13(iz, ia)[i] +
786 K2.K13(iz, ia)[i])
787 : 0.0)),
788 dd = -0.5 * (r * dK1[j].K14(iz, ia)[i] +
789 ((j == it) ? dr_dT1 * (K1.K14(iz, ia)[i] +
790 K2.K14(iz, ia)[i])
791 : 0.0)),
792 du = -0.5 * (r * dK1[j].K23(iz, ia)[i] +
793 ((j == it) ? dr_dT1 * (K1.K23(iz, ia)[i] +
794 K2.K23(iz, ia)[i])
795 : 0.0)),
796 dv = -0.5 * (r * dK1[j].K24(iz, ia)[i] +
797 ((j == it) ? dr_dT1 * (K1.K24(iz, ia)[i] +
798 K2.K24(iz, ia)[i])
799 : 0.0)),
800 dw = -0.5 * (r * dK1[j].K34(iz, ia)[i] +
801 ((j == it) ? dr_dT1 * (K1.K34(iz, ia)[i] +
802 K2.K34(iz, ia)[i])
803 : 0.0));
804 const Numeric db2 = 2 * db * b, dc2 = 2 * dc * c, dd2 = 2 * dd * d,
805 du2 = 2 * du * u, dv2 = 2 * dv * v, dw2 = 2 * dw * w;
806 const Numeric dtmp =
807 2 * w2 * dw2 +
808 2 * (db2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
809 b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2) +
810 dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
811 c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2) +
812 dd2 * (d2 * 0.5 + u2 - v2 - w2) +
813 d2 * (dd2 * 0.5 + du2 - dv2 - dw2) +
814 du2 * (u2 * 0.5 + v2 + w2) + u2 * (du2 * 0.5 + dv2 + dw2) +
815 dv2 * (v2 * 0.5 + w2) + v2 * (dv2 * 0.5 + dw2) +
816 4 * (db * d * u * w - db * c * v * w - dc * d * u * v +
817 b * dd * u * w - b * dc * v * w - c * dd * u * v +
818 b * d * du * w - b * c * dv * w - c * d * du * v +
819 b * d * u * dw - b * c * v * dw - c * d * u * dv));
820 const Complex dConst1 = 0.5 * dtmp / Const1;
821 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
822 const Complex dx =
823 x_zero ? 0 : (0.5 * (dConst2 + dConst1) / tmp_x_sqrt) * sqrt_05;
824 const Complex dy = y_zero ? 0
825 : (0.5 * (dConst2 - dConst1) / tmp_y_sqrt) *
826 sqrt_05 * Complex(0, 1);
827 const Complex dx2 = 2 * x * dx;
828 const Complex dy2 = 2 * y * dy;
829 const Complex dcy = -sy * dy;
830 const Complex dsy = cy * dy;
831 const Complex dcx = sx * dx;
832 const Complex dsx = cx * dx;
833 const Complex dix = -dx * ix * ix;
834 const Complex diy = -dy * iy * iy;
835 const Complex dx2dy2 = dx2 + dy2;
836 const Complex dC0c =
837 either_zero
838 ? 0.0
839 : (dcy * x2 + cy * dx2 + dcx * y2 + cx * dy2 - C0c * dx2dy2) *
840 inv_x2y2;
841 const Complex dC1c =
842 either_zero ? 0.0
843 : (dsy * x2 * iy + sy * dx2 * iy + sy * x2 * diy +
844 dsx * y2 * ix + sx * dy2 * ix + sx * y2 * dix -
845 C1c * dx2dy2) *
846 inv_x2y2;
847 const Complex dC2c =
848 both_zero ? 0.0 : (dcx - dcy - C2c * dx2dy2) * inv_x2y2;
849 const Complex dC3c =
850 both_zero ? 0.0
851 : ((x_zero ? -dsy * iy - sy * diy
852 : y_zero ? dsx * ix + sx * dix
853 : dsx * ix + sx * dix - dsy * iy -
854 sy * diy) -
855 C3c * dx2dy2) *
856 inv_x2y2;
857
858 const Numeric& dC0 = real_val(dC0c);
859 const Numeric& dC1 = real_val(dC1c);
860 const Numeric& dC2 = real_val(dC2c);
861 const Numeric& dC3 = real_val(dC3c);
862 dT1[j].Mat4(i).noalias() =
863 T.Mat4(i) * da +
864 exp_a *
865 (Eigen::Matrix4d()
866 << dC0 + dC2 * (b2 + c2 + d2) + C2 * (db2 + dc2 + dd2),
867 db * C1 + b * dC1 + dC2 * (-c * u - d * v) +
868 C2 * (-dc * u - dd * v - c * du - d * dv) +
869 dC3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) -
870 v * (b * v + c * w)) +
871 C3 * (db * (b2 + c2 + d2) - du * (b * u - d * w) -
872 dv * (b * v + c * w) + b * (db2 + dc2 + dd2) -
873 u * (db * u - dd * w) - v * (db * v + dc * w) -
874 u * (b * du - d * dw) - v * (b * dv + c * dw)),
875 dC1 * c + C1 * dc + dC2 * (b * u - d * w) +
876 C2 * (db * u - dd * w + b * du - d * dw) +
877 dC3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) -
878 w * (b * v + c * w)) +
879 C3 * (dc * (b2 + c2 + d2) - du * (c * u + d * v) -
880 dw * (b * v + c * w) + c * (db2 + dc2 + dd2) -
881 u * (dc * u + dd * v) - w * (db * v + dc * w) -
882 u * (c * du + d * dv) - w * (b * dv + c * dw)),
883 dC1 * d + C1 * dd + dC2 * (b * v + c * w) +
884 C2 * (db * v + dc * w + b * dv + c * dw) +
885 dC3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) +
886 w * (b * u - d * w)) +
887 C3 * (dd * (b2 + c2 + d2) - dv * (c * u + d * v) +
888 dw * (b * u - d * w) + d * (db2 + dc2 + dd2) -
889 v * (dc * u + dd * v) + w * (db * u - dd * w) -
890 v * (c * du + d * dv) + w * (b * du - d * dw)),
891
892 db * C1 + b * dC1 + dC2 * (c * u + d * v) +
893 C2 * (dc * u + dd * v + c * du + d * dv) +
894 dC3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) +
895 d * (b * d + u * w)) +
896 C3 * (-db * (-b2 + u2 + v2) + dc * (b * c - v * w) +
897 dd * (b * d + u * w) - b * (-db2 + du2 + dv2) +
898 c * (db * c - dv * w) + d * (db * d + du * w) +
899 c * (b * dc - v * dw) + d * (b * dd + u * dw)),
900 dC0 + dC2 * (b2 - u2 - v2) + C2 * (db2 - du2 - dv2),
901 dC2 * (b * c - v * w) +
902 C2 * (db * c + b * dc - dv * w - v * dw) + dC1 * u +
903 C1 * du +
904 dC3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
905 w * (b * d + u * w)) +
906 C3 * (dc * (c * u + d * v) - du * (-b2 + u2 + v2) -
907 dw * (b * d + u * w) + c * (dc * u + dd * v) -
908 u * (-db2 + du2 + dv2) - w * (db * d + du * w) +
909 c * (c * du + d * dv) - w * (b * dd + u * dw)),
910 dC2 * (b * d + u * w) +
911 C2 * (db * d + b * dd + du * w + u * dw) + dC1 * v +
912 C1 * dv +
913 dC3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
914 w * (b * c - v * w)) +
915 C3 * (dd * (c * u + d * v) - dv * (-b2 + u2 + v2) +
916 dw * (b * c - v * w) + d * (dc * u + dd * v) -
917 v * (-db2 + du2 + dv2) + w * (db * c - dv * w) +
918 d * (c * du + d * dv) + w * (b * dc - v * dw)),
919
920 dC1 * c + C1 * dc + dC2 * (-b * u + d * w) +
921 C2 * (-db * u + dd * w - b * du + d * dw) +
922 dC3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) +
923 d * (c * d - u * v)) +
924 C3 * (db * (b * c - v * w) - dc * (-c2 + u2 + w2) +
925 dd * (c * d - u * v) + b * (db * c - dv * w) -
926 c * (-dc2 + du2 + dw2) + d * (dc * d - du * v) +
927 b * (b * dc - v * dw) + d * (c * dd - u * dv)),
928 dC2 * (b * c - v * w) +
929 C2 * (db * c + b * dc - dv * w - v * dw) - dC1 * u -
930 C1 * du +
931 dC3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
932 v * (c * d - u * v)) +
933 C3 * (-db * (b * u - d * w) + du * (-c2 + u2 + w2) -
934 dv * (c * d - u * v) - b * (db * u - dd * w) +
935 u * (-dc2 + du2 + dw2) - v * (dc * d - du * v) -
936 b * (b * du - d * dw) - v * (c * dd - u * dv)),
937 dC0 + dC2 * (c2 - u2 - w2) + C2 * (dc2 - du2 - dw2),
938 dC2 * (c * d - u * v) +
939 C2 * (dc * d + c * dd - du * v - u * dv) + dC1 * w +
940 C1 * dw +
941 dC3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
942 w * (-c2 + u2 + w2)) +
943 C3 * (-dd * (b * u - d * w) + dv * (b * c - v * w) -
944 dw * (-c2 + u2 + w2) - d * (db * u - dd * w) +
945 v * (db * c - dv * w) - w * (-dc2 + du2 + dw2) -
946 d * (b * du - d * dw) + v * (b * dc - v * dw)),
947
948 dC1 * d + C1 * dd + dC2 * (-b * v - c * w) +
949 C2 * (-db * v - dc * w - b * dv - c * dw) +
950 dC3 * (b * (b * d + u * w) + c * (c * d - u * v) -
951 d * (-d2 + v2 + w2)) +
952 C3 * (db * (b * d + u * w) + dc * (c * d - u * v) -
953 dd * (-d2 + v2 + w2) + b * (db * d + du * w) +
954 c * (dc * d - du * v) - d * (-dd2 + dv2 + dw2) +
955 b * (b * dd + u * dw) + c * (c * dd - u * dv)),
956 dC2 * (b * d + u * w) +
957 C2 * (db * d + b * dd + du * w + u * dw) - dC1 * v -
958 C1 * dv +
959 dC3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
960 v * (-d2 + v2 + w2)) +
961 C3 * (-db * (b * v + c * w) - du * (c * d - u * v) +
962 dv * (-d2 + v2 + w2) - b * (db * v + dc * w) -
963 u * (dc * d - du * v) + v * (-dd2 + dv2 + dw2) -
964 b * (b * dv + c * dw) - u * (c * dd - u * dv)),
965 dC2 * (c * d - u * v) +
966 C2 * (dc * d + c * dd - du * v - u * dv) - dC1 * w -
967 C1 * dw +
968 dC3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
969 w * (-d2 + v2 + w2)) +
970 C3 * (-dc * (b * v + c * w) + du * (b * d + u * w) +
971 dw * (-d2 + v2 + w2) - c * (db * v + dc * w) +
972 u * (db * d + du * w) + w * (-dd2 + dv2 + dw2) -
973 c * (b * dv + c * dw) + u * (b * dd + u * dw)),
974 dC0 + dC2 * (d2 - v2 - w2) + C2 * (dd2 - dv2 - dw2))
975 .finished();
976 }
977 if (dK2[j].NumberOfFrequencies()) {
978 const Numeric da = -0.5 * (r * dK2[j].Kjj(iz, ia)[i] +
979 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] +
980 K2.Kjj(iz, ia)[i])
981 : 0.0)),
982 db = -0.5 * (r * dK2[j].K12(iz, ia)[i] +
983 ((j == it) ? dr_dT2 * (K1.K12(iz, ia)[i] +
984 K2.K12(iz, ia)[i])
985 : 0.0)),
986 dc = -0.5 * (r * dK2[j].K13(iz, ia)[i] +
987 ((j == it) ? dr_dT2 * (K1.K13(iz, ia)[i] +
988 K2.K13(iz, ia)[i])
989 : 0.0)),
990 dd = -0.5 * (r * dK2[j].K14(iz, ia)[i] +
991 ((j == it) ? dr_dT2 * (K1.K14(iz, ia)[i] +
992 K2.K14(iz, ia)[i])
993 : 0.0)),
994 du = -0.5 * (r * dK2[j].K23(iz, ia)[i] +
995 ((j == it) ? dr_dT2 * (K1.K23(iz, ia)[i] +
996 K2.K23(iz, ia)[i])
997 : 0.0)),
998 dv = -0.5 * (r * dK2[j].K24(iz, ia)[i] +
999 ((j == it) ? dr_dT2 * (K1.K24(iz, ia)[i] +
1000 K2.K24(iz, ia)[i])
1001 : 0.0)),
1002 dw = -0.5 * (r * dK2[j].K34(iz, ia)[i] +
1003 ((j == it) ? dr_dT2 * (K1.K34(iz, ia)[i] +
1004 K2.K34(iz, ia)[i])
1005 : 0.0));
1006 const Numeric db2 = 2 * db * b, dc2 = 2 * dc * c, dd2 = 2 * dd * d,
1007 du2 = 2 * du * u, dv2 = 2 * dv * v, dw2 = 2 * dw * w;
1008 const Numeric dtmp =
1009 2 * w2 * dw2 +
1010 2 * (db2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
1011 b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2) +
1012 dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
1013 c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2) +
1014 dd2 * (d2 * 0.5 + u2 - v2 - w2) +
1015 d2 * (dd2 * 0.5 + du2 - dv2 - dw2) +
1016 du2 * (u2 * 0.5 + v2 + w2) + u2 * (du2 * 0.5 + dv2 + dw2) +
1017 dv2 * (v2 * 0.5 + w2) + v2 * (dv2 * 0.5 + dw2) +
1018 4 * (db * d * u * w - db * c * v * w - dc * d * u * v +
1019 b * dd * u * w - b * dc * v * w - c * dd * u * v +
1020 b * d * du * w - b * c * dv * w - c * d * du * v +
1021 b * d * u * dw - b * c * v * dw - c * d * u * dv));
1022 const Complex dConst1 = 0.5 * dtmp / Const1;
1023 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1024 const Complex dx =
1025 x_zero ? 0 : (0.5 * (dConst2 + dConst1) / tmp_x_sqrt) * sqrt_05;
1026 const Complex dy = y_zero ? 0
1027 : (0.5 * (dConst2 - dConst1) / tmp_y_sqrt) *
1028 sqrt_05 * Complex(0, 1);
1029 const Complex dx2 = 2 * x * dx;
1030 const Complex dy2 = 2 * y * dy;
1031 const Complex dcy = -sy * dy;
1032 const Complex dsy = cy * dy;
1033 const Complex dcx = sx * dx;
1034 const Complex dsx = cx * dx;
1035 const Complex dix = -dx * ix * ix;
1036 const Complex diy = -dy * iy * iy;
1037 const Complex dx2dy2 = dx2 + dy2;
1038 const Complex dC0c =
1039 either_zero
1040 ? 0.0
1041 : (dcy * x2 + cy * dx2 + dcx * y2 + cx * dy2 - C0c * dx2dy2) *
1042 inv_x2y2;
1043 const Complex dC1c =
1044 either_zero ? 0.0
1045 : (dsy * x2 * iy + sy * dx2 * iy + sy * x2 * diy +
1046 dsx * y2 * ix + sx * dy2 * ix + sx * y2 * dix -
1047 C1c * dx2dy2) *
1048 inv_x2y2;
1049 const Complex dC2c =
1050 both_zero ? 0.0 : (dcx - dcy - C2c * dx2dy2) * inv_x2y2;
1051 const Complex dC3c =
1052 both_zero ? 0.0
1053 : ((x_zero ? -dsy * iy - sy * diy
1054 : y_zero ? dsx * ix + sx * dix
1055 : dsx * ix + sx * dix - dsy * iy -
1056 sy * diy) -
1057 C3c * dx2dy2) *
1058 inv_x2y2;
1059
1060 const Numeric& dC0 = real_val(dC0c);
1061 const Numeric& dC1 = real_val(dC1c);
1062 const Numeric& dC2 = real_val(dC2c);
1063 const Numeric& dC3 = real_val(dC3c);
1064 dT2[j].Mat4(i).noalias() =
1065 T.Mat4(i) * da +
1066 exp_a *
1067 (Eigen::Matrix4d()
1068 << dC0 + dC2 * (b2 + c2 + d2) + C2 * (db2 + dc2 + dd2),
1069 db * C1 + b * dC1 + dC2 * (-c * u - d * v) +
1070 C2 * (-dc * u - dd * v - c * du - d * dv) +
1071 dC3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) -
1072 v * (b * v + c * w)) +
1073 C3 * (db * (b2 + c2 + d2) - du * (b * u - d * w) -
1074 dv * (b * v + c * w) + b * (db2 + dc2 + dd2) -
1075 u * (db * u - dd * w) - v * (db * v + dc * w) -
1076 u * (b * du - d * dw) - v * (b * dv + c * dw)),
1077 dC1 * c + C1 * dc + dC2 * (b * u - d * w) +
1078 C2 * (db * u - dd * w + b * du - d * dw) +
1079 dC3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) -
1080 w * (b * v + c * w)) +
1081 C3 * (dc * (b2 + c2 + d2) - du * (c * u + d * v) -
1082 dw * (b * v + c * w) + c * (db2 + dc2 + dd2) -
1083 u * (dc * u + dd * v) - w * (db * v + dc * w) -
1084 u * (c * du + d * dv) - w * (b * dv + c * dw)),
1085 dC1 * d + C1 * dd + dC2 * (b * v + c * w) +
1086 C2 * (db * v + dc * w + b * dv + c * dw) +
1087 dC3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) +
1088 w * (b * u - d * w)) +
1089 C3 * (dd * (b2 + c2 + d2) - dv * (c * u + d * v) +
1090 dw * (b * u - d * w) + d * (db2 + dc2 + dd2) -
1091 v * (dc * u + dd * v) + w * (db * u - dd * w) -
1092 v * (c * du + d * dv) + w * (b * du - d * dw)),
1093
1094 db * C1 + b * dC1 + dC2 * (c * u + d * v) +
1095 C2 * (dc * u + dd * v + c * du + d * dv) +
1096 dC3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) +
1097 d * (b * d + u * w)) +
1098 C3 * (-db * (-b2 + u2 + v2) + dc * (b * c - v * w) +
1099 dd * (b * d + u * w) - b * (-db2 + du2 + dv2) +
1100 c * (db * c - dv * w) + d * (db * d + du * w) +
1101 c * (b * dc - v * dw) + d * (b * dd + u * dw)),
1102 dC0 + dC2 * (b2 - u2 - v2) + C2 * (db2 - du2 - dv2),
1103 dC2 * (b * c - v * w) +
1104 C2 * (db * c + b * dc - dv * w - v * dw) + dC1 * u +
1105 C1 * du +
1106 dC3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
1107 w * (b * d + u * w)) +
1108 C3 * (dc * (c * u + d * v) - du * (-b2 + u2 + v2) -
1109 dw * (b * d + u * w) + c * (dc * u + dd * v) -
1110 u * (-db2 + du2 + dv2) - w * (db * d + du * w) +
1111 c * (c * du + d * dv) - w * (b * dd + u * dw)),
1112 dC2 * (b * d + u * w) +
1113 C2 * (db * d + b * dd + du * w + u * dw) + dC1 * v +
1114 C1 * dv +
1115 dC3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
1116 w * (b * c - v * w)) +
1117 C3 * (dd * (c * u + d * v) - dv * (-b2 + u2 + v2) +
1118 dw * (b * c - v * w) + d * (dc * u + dd * v) -
1119 v * (-db2 + du2 + dv2) + w * (db * c - dv * w) +
1120 d * (c * du + d * dv) + w * (b * dc - v * dw)),
1121
1122 dC1 * c + C1 * dc + dC2 * (-b * u + d * w) +
1123 C2 * (-db * u + dd * w - b * du + d * dw) +
1124 dC3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) +
1125 d * (c * d - u * v)) +
1126 C3 * (db * (b * c - v * w) - dc * (-c2 + u2 + w2) +
1127 dd * (c * d - u * v) + b * (db * c - dv * w) -
1128 c * (-dc2 + du2 + dw2) + d * (dc * d - du * v) +
1129 b * (b * dc - v * dw) + d * (c * dd - u * dv)),
1130 dC2 * (b * c - v * w) +
1131 C2 * (db * c + b * dc - dv * w - v * dw) - dC1 * u -
1132 C1 * du +
1133 dC3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
1134 v * (c * d - u * v)) +
1135 C3 * (-db * (b * u - d * w) + du * (-c2 + u2 + w2) -
1136 dv * (c * d - u * v) - b * (db * u - dd * w) +
1137 u * (-dc2 + du2 + dw2) - v * (dc * d - du * v) -
1138 b * (b * du - d * dw) - v * (c * dd - u * dv)),
1139 dC0 + dC2 * (c2 - u2 - w2) + C2 * (dc2 - du2 - dw2),
1140 dC2 * (c * d - u * v) +
1141 C2 * (dc * d + c * dd - du * v - u * dv) + dC1 * w +
1142 C1 * dw +
1143 dC3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
1144 w * (-c2 + u2 + w2)) +
1145 C3 * (-dd * (b * u - d * w) + dv * (b * c - v * w) -
1146 dw * (-c2 + u2 + w2) - d * (db * u - dd * w) +
1147 v * (db * c - dv * w) - w * (-dc2 + du2 + dw2) -
1148 d * (b * du - d * dw) + v * (b * dc - v * dw)),
1149
1150 dC1 * d + C1 * dd + dC2 * (-b * v - c * w) +
1151 C2 * (-db * v - dc * w - b * dv - c * dw) +
1152 dC3 * (b * (b * d + u * w) + c * (c * d - u * v) -
1153 d * (-d2 + v2 + w2)) +
1154 C3 * (db * (b * d + u * w) + dc * (c * d - u * v) -
1155 dd * (-d2 + v2 + w2) + b * (db * d + du * w) +
1156 c * (dc * d - du * v) - d * (-dd2 + dv2 + dw2) +
1157 b * (b * dd + u * dw) + c * (c * dd - u * dv)),
1158 dC2 * (b * d + u * w) +
1159 C2 * (db * d + b * dd + du * w + u * dw) - dC1 * v -
1160 C1 * dv +
1161 dC3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
1162 v * (-d2 + v2 + w2)) +
1163 C3 * (-db * (b * v + c * w) - du * (c * d - u * v) +
1164 dv * (-d2 + v2 + w2) - b * (db * v + dc * w) -
1165 u * (dc * d - du * v) + v * (-dd2 + dv2 + dw2) -
1166 b * (b * dv + c * dw) - u * (c * dd - u * dv)),
1167 dC2 * (c * d - u * v) +
1168 C2 * (dc * d + c * dd - du * v - u * dv) - dC1 * w -
1169 C1 * dw +
1170 dC3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
1171 w * (-d2 + v2 + w2)) +
1172 C3 * (-dc * (b * v + c * w) + du * (b * d + u * w) +
1173 dw * (-d2 + v2 + w2) - c * (db * v + dc * w) +
1174 u * (db * d + du * w) + w * (-dd2 + dv2 + dw2) -
1175 c * (b * dv + c * dw) + u * (b * dd + u * dw)),
1176 dC0 + dC2 * (d2 - v2 - w2) + C2 * (dd2 - dv2 - dw2))
1177 .finished();
1178 }
1179 }
1180 }
1181 }
1182}
1183
1185 const PropagationMatrix& K1,
1186 const PropagationMatrix& K2,
1187 const Numeric& r) noexcept {
1188 switch (K1.StokesDimensions()) {
1189 case 4:
1190 transmat4(T, K1, K2, r);
1191 break;
1192 case 3:
1193 transmat3(T, K1, K2, r);
1194 break;
1195 case 2:
1196 transmat2(T, K1, K2, r);
1197 break;
1198 case 1:
1199 transmat1(T, K1, K2, r);
1200 break;
1201 }
1202}
1203
1207 const PropagationMatrix& K1,
1208 const PropagationMatrix& K2,
1209 const ArrayOfPropagationMatrix& dK1,
1210 const ArrayOfPropagationMatrix& dK2,
1211 const Numeric& r,
1212 const Numeric& dr_dT1 = 0,
1213 const Numeric& dr_dT2 = 0,
1214 const Index it = -1,
1215 const Index iz = 0,
1216 const Index ia = 0) noexcept {
1217 switch (K1.StokesDimensions()) {
1218 case 4:
1219 dtransmat4(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1220 break;
1221 case 3:
1222 dtransmat3(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1223 break;
1224 case 2:
1225 dtransmat2(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1226 break;
1227 case 1:
1228 dtransmat1(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1229 break;
1230 }
1231}
1232
1236 const PropagationMatrix& K1,
1237 const PropagationMatrix& K2,
1238 const ArrayOfPropagationMatrix& dK1,
1239 const ArrayOfPropagationMatrix& dK2,
1240 const Numeric& r,
1241 const Numeric& dr_dtemp1,
1242 const Numeric& dr_dtemp2,
1243 const Index temp_deriv_pos) {
1244 if (not dT1.nelem())
1245 transmat(T, K1, K2, r);
1246 else
1247 dtransmat(
1248 T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dtemp1, dr_dtemp2, temp_deriv_pos);
1249}
1250
1253 const PropagationMatrix& K,
1254 const StokesVector& a,
1255 const StokesVector& S,
1256 const ArrayOfPropagationMatrix& dK,
1257 const ArrayOfStokesVector& da,
1258 const ArrayOfStokesVector& dS,
1259 const ConstVectorView& B,
1260 const ConstVectorView& dB_dT,
1261 const ArrayOfRetrievalQuantity& jacobian_quantities,
1262 const bool& jacobian_do) {
1263 for (Index i = 0; i < K.NumberOfFrequencies(); i++) {
1264 if (K.IsRotational(i)) {
1265 J.SetZero(i);
1266 if (jacobian_do) {
1267 for (Index j = 0; j < jacobian_quantities.nelem(); j++)
1268 if (dJ[j].Frequencies())
1269 dJ[j].SetZero(i);
1270 }
1271 } else {
1272 J.setSource(a, B, S, i);
1273
1274 switch (J.stokes_dim) {
1275 case 4: {
1276 const auto invK = inv_prop_matrix<4>(K.Data()(0, 0, i, joker));
1277 J.Vec4(i) = invK * J.Vec4(i);
1278 if (jacobian_do) {
1279 for (Index j = 0; j < jacobian_quantities.nelem(); j++) {
1280 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1281 dJ[j].Vec4(i).noalias() = 0.5 * invK *
1282 (source_vector<4>(a, B, da[j], dB_dT, dS[j],
1283 jacobian_quantities[j] == Jacobian::Atm::Temperature, i) -
1284 prop_matrix<4>(dK[j].Data()(0, 0, i, joker)) * J.Vec4(i));
1285 }
1286 }
1287 }
1288 } break;
1289 case 3: {
1290 const auto invK = inv_prop_matrix<3>(K.Data()(0, 0, i, joker));
1291 J.Vec3(i) = invK * J.Vec3(i);
1292 if (jacobian_do) {
1293 for (Index j = 0; j < jacobian_quantities.nelem(); j++) {
1294 // Skip others!
1295 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1296 dJ[j].Vec3(i).noalias() = 0.5 * invK *
1297 (source_vector<3>(a, B, da[j], dB_dT, dS[j],
1298 jacobian_quantities[j] == Jacobian::Atm::Temperature, i) -
1299 prop_matrix<3>(dK[j].Data()(0, 0, i, joker)) * J.Vec3(i));
1300 }
1301 }
1302 }
1303 } break;
1304 case 2: {
1305 const auto invK = inv_prop_matrix<2>(K.Data()(0, 0, i, joker));
1306 J.Vec2(i) = invK * J.Vec2(i);
1307 if (jacobian_do) {
1308 for (Index j = 0; j < jacobian_quantities.nelem(); j++) {
1309 // Skip others!
1310 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1311 dJ[j].Vec2(i).noalias() = 0.5 * invK *
1312 (source_vector<2>(a, B, da[j], dB_dT, dS[j],
1313 jacobian_quantities[j] == Jacobian::Atm::Temperature, i) -
1314 prop_matrix<2>(dK[j].Data()(0, 0, i, joker)) * J.Vec2(i));
1315 }
1316 }
1317 }
1318 } break;
1319 default: {
1320 const auto invK = inv_prop_matrix<1>(K.Data()(0, 0, i, joker));
1321 J.Vec1(i) = invK * J.Vec1(i);
1322 if (jacobian_do) {
1323 for (Index j = 0; j < jacobian_quantities.nelem(); j++) {
1324 // Skip others!
1325 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1326 dJ[j].Vec1(i).noalias() = 0.5 * invK *
1327 (source_vector<1>(a, B, da[j], dB_dT, dS[j],
1328 jacobian_quantities[j] == Jacobian::Atm::Temperature, i) -
1329 prop_matrix<1>(dK[j].Data()(0, 0, i, joker)) * J.Vec1(i));
1330 }
1331 }
1332 }
1333 } break;
1334 }
1335 }
1336 }
1337}
1338
1342 const RadiationVector& J1,
1343 const RadiationVector& J2,
1344 const ArrayOfRadiationVector& dJ1,
1345 const ArrayOfRadiationVector& dJ2,
1346 const TransmissionMatrix& T,
1347 const TransmissionMatrix& PiT,
1348 const ArrayOfTransmissionMatrix& dT1,
1349 const ArrayOfTransmissionMatrix& dT2,
1350 [[maybe_unused]] const PropagationMatrix& K1,
1351 [[maybe_unused]] const PropagationMatrix& K2,
1352 [[maybe_unused]] const ArrayOfPropagationMatrix& dK1,
1353 [[maybe_unused]] const ArrayOfPropagationMatrix& dK2,
1354 [[maybe_unused]] const Numeric r,
1355 [[maybe_unused]] const Vector& dr1,
1356 [[maybe_unused]] const Vector& dr2,
1357 [[maybe_unused]] const Index ia,
1358 [[maybe_unused]] const Index iz,
1359 const RadiativeTransferSolver solver) {
1360 switch (solver) {
1362 I.rem_avg(J1, J2);
1363 for (size_t i = 0; i < dI1.size(); i++) {
1364 dI1[i].addDerivEmission(PiT, dT1[i], T, I, dJ1[i]);
1365 dI2[i].addDerivEmission(PiT, dT2[i], T, I, dJ2[i]);
1366 }
1367 I.leftMul(T);
1368 I.add_avg(J1, J2);
1369 } break;
1370
1372 for (size_t i = 0; i < dI1.size(); i++) {
1373 dI1[i].addDerivTransmission(PiT, dT1[i], I);
1374 dI2[i].addDerivTransmission(PiT, dT2[i], I);
1375 }
1376 I.leftMul(T);
1377 } break;
1378
1380 ARTS_USER_ERROR_IF (dI1.size(),
1381 "Cannot support derivatives with current integration method\n");
1382
1383 I.leftMul(T);
1384 I.add_weighted(T, J1, J2, K1.Data()(ia, iz, joker, joker), K2.Data()(ia, iz, joker, joker), r);
1385
1386 } break;
1387 }
1388}
1389
1392 const CumulativeTransmission type) /*[[expects: T.nelem()>0]]*/
1393{
1394 const Index n = T.nelem();
1395 const Index nf = n ? T[0].Frequencies() : 1;
1396 const Index ns = n ? T[0].stokes_dim : 1;
1398 switch (type) {
1400 for (Index i = 1; i < n; i++) PiT[i].mul(PiT[i - 1], T[i]);
1401 } break;
1403 for (Index i = 1; i < n; i++) PiT[i].mul(T[i], PiT[i - 1]);
1404 } break;
1405 }
1406 return PiT; // Note how the output is such that forward transmission is from -1 to 0
1407}
1408
1409// TEST CODE BEGIN
1410
1414 const RadiationVector &I_incoming,
1416 const ArrayOfTransmissionMatrix& PiTf,
1417 const ArrayOfTransmissionMatrix& PiTr,
1422 const BackscatterSolver solver) {
1423 const Index np = I.nelem();
1424 const Index nv = np ? I[0].Frequencies() : 0;
1425 const Index ns = np ? I[0].stokes_dim : 1;
1426 const Index nq = np ? dI[0][0].nelem() : 0;
1427
1428 // For all transmission, the I-vector is the same
1429 for (Index ip = 0; ip < np; ip++)
1430 I[ip].setBackscatterTransmission(I_incoming, PiTr[ip], PiTf[ip], Z[ip]);
1431
1432 for (Index ip = 0; ip < np; ip++) {
1433 for (Index iq = 0; iq < nq; iq++) {
1434 dI[ip][ip][iq].setBackscatterTransmissionDerivative(
1435 I_incoming, PiTr[ip], PiTf[ip], dZ[ip][iq]);
1436 }
1437 }
1438
1439 switch(solver) {
1441 // Forward and backwards transmission derivatives
1442 switch(ns) {
1443 case 1: {
1444 BackscatterSolverCommutativeTransmissionStokesDimOne:
1445 for (Index ip = 0; ip < np; ip++) {
1446 for (Index j = ip; j < np; j++) {
1447 for (Index iq = 0; iq < nq; iq++) {
1448 for (Index iv = 0; iv < nv; iv++) {
1449 dI[ip][j][iq].Vec1(iv).noalias() +=
1450 T[ip].Mat1(iv).inverse() *
1451 (dT1[ip][iq].Mat1(iv) + dT2[ip][iq].Mat1(iv)) *
1452 I[j].Vec1(iv);
1453
1454 if (j < np - 1 and j > ip)
1455 dI[ip][j][iq].Vec1(iv).noalias() +=
1456 T[ip+1].Mat1(iv).inverse() *
1457 (dT1[ip][iq].Mat1(iv) + dT2[ip][iq].Mat1(iv)) *
1458 I[j].Vec1(iv);
1459 }
1460 }
1461 }
1462 }
1463 } break;
1464 case 2: {
1465 for (Index ip = 0; ip < np; ip++) {
1466 for (Index j = ip; j < np; j++) {
1467 for (Index iq = 0; iq < nq; iq++) {
1468 for (Index iv = 0; iv < nv; iv++) {
1469 dI[ip][j][iq].Vec2(iv).noalias() +=
1470 T[ip].Mat2(iv).inverse() *
1471 (dT1[ip][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
1472 I[j].Vec2(iv);
1473
1474 if (j < np - 1 and j > ip)
1475 dI[ip][j][iq].Vec2(iv).noalias() +=
1476 T[ip+1].Mat2(iv).inverse() *
1477 (dT1[ip][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
1478 I[j].Vec2(iv);
1479 }
1480 }
1481 }
1482 }
1483 } break;
1484 case 3: {
1485 for (Index ip = 0; ip < np; ip++) {
1486 for (Index j = ip; j < np; j++) {
1487 for (Index iq = 0; iq < nq; iq++) {
1488 for (Index iv = 0; iv < nv; iv++) {
1489 dI[ip][j][iq].Vec3(iv).noalias() +=
1490 T[ip].Mat3(iv).inverse() *
1491 (dT1[ip][iq].Mat3(iv) + dT2[ip][iq].Mat3(iv)) *
1492 I[j].Vec3(iv);
1493
1494 if (j < np - 1 and j > ip)
1495 dI[ip][j][iq].Vec3(iv).noalias() +=
1496 T[ip+1].Mat3(iv).inverse() *
1497 (dT1[ip][iq].Mat3(iv) + dT2[ip][iq].Mat3(iv)) *
1498 I[j].Vec3(iv);
1499 }
1500 }
1501 }
1502 }
1503 } break;
1504 case 4: {
1505 for (Index ip = 0; ip < np; ip++) {
1506 for (Index j = ip; j < np; j++) {
1507 for (Index iq = 0; iq < nq; iq++) {
1508 for (Index iv = 0; iv < nv; iv++) {
1509 dI[ip][j][iq].Vec4(iv).noalias() +=
1510 T[ip].Mat4(iv).inverse() *
1511 (dT1[ip][iq].Mat4(iv) + dT2[ip][iq].Mat4(iv)) *
1512 I[j].Vec4(iv);
1513
1514 if (j < np - 1 and j > ip)
1515 dI[ip][j][iq].Vec4(iv).noalias() +=
1516 T[ip+1].Mat4(iv).inverse() *
1517 (dT1[ip][iq].Mat4(iv) + dT2[ip][iq].Mat4(iv)) *
1518 I[j].Vec4(iv);
1519 }
1520 }
1521 }
1522 }
1523 } break;
1524 }
1525 } break;
1527 std::runtime_error("Do not activate this code. It is not ready yet");
1528 switch(ns) {
1529 case 1: {
1530 // This is the same as CommutativeTransmission, so use that code
1531 goto BackscatterSolverCommutativeTransmissionStokesDimOne;
1532 } break;
1533 case 2: {
1534 for (Index ip = 0; ip < np; ip++) {
1535 for (Index j = ip; j < np; j++) {
1536 for (Index iq = 0; iq < nq; iq++) {
1537 for (Index iv = 0; iv < nv; iv++) {
1538 if(ip > 1) {
1539 dI[ip][j][iq].Vec2(iv).noalias() +=
1540 PiTr[ip-2].Mat2(iv) *
1541 (dT1[ip-1][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
1542 PiTr[ip-1].Mat2(iv).inverse() * I[j].Vec2(iv);
1543
1544 if(j < np - 1)
1545 dI[ip][j][iq].Vec2(iv).noalias() +=
1546 PiTr[ip].Mat2(iv) * Z[ip].Mat2(iv) * PiTf[ip-2].Mat2(iv) *
1547 (dT1[ip][iq].Mat2(iv) + dT2[ip+1][iq].Mat2(iv)) *
1548 PiTf[ip-1].Mat2(iv).inverse() * PiTf[ip].Mat2(iv) *
1549 I[0].Vec2(iv);
1550 }
1551 else {
1552 dI[ip][j][iq].Vec2(iv).noalias() +=
1553 (dT1[ip-1][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
1554 PiTr[ip-1].Mat2(iv).inverse() * I[j].Vec2(iv);
1555
1556 if(j < np - 1)
1557 dI[ip][j][iq].Vec2(iv).noalias() +=
1558 PiTr[ip].Mat2(iv) * Z[ip].Mat2(iv) *
1559 (dT1[ip][iq].Mat2(iv) + dT2[ip+1][iq].Mat2(iv)) *
1560 PiTf[ip-1].Mat2(iv).inverse() * PiTf[ip].Mat2(iv) *
1561 I[0].Vec2(iv);
1562 }
1563 }
1564 }
1565 }
1566 }
1567 } break;
1568 case 3: {
1569 } break;
1570 case 4: {
1571 } break;
1572 }
1573 } break;
1574 }
1575}
1576
1578 const ConstMatrixView& pnd) {
1579 const Index ns = Pe.ncols();
1580 const Index nv = Pe.npages();
1581 const Index np = Pe.nbooks();
1582 const Index ne = Pe.nshelves();
1583
1585
1586 for (Index ip = 0; ip < np; ip++) {
1587
1588 aotm[ip].setZero();
1589
1590 switch (ns) {
1591 case 4:
1592 for (Index iv = 0; iv < nv; iv++)
1593 for (Index ie = 0; ie < ne; ie++)
1594 aotm[ip].Mat4(iv).noalias() +=
1595 pnd(ie, ip) * prop_matrix<4>(Pe(ie, ip, iv, joker, joker));
1596 break;
1597 case 3:
1598 for (Index iv = 0; iv < nv; iv++)
1599 for (Index ie = 0; ie < ne; ie++)
1600 aotm[ip].Mat3(iv).noalias() +=
1601 pnd(ie, ip) * prop_matrix<3>(Pe(ie, ip, iv, joker, joker));
1602 break;
1603 case 2:
1604 for (Index iv = 0; iv < nv; iv++)
1605 for (Index ie = 0; ie < ne; ie++)
1606 aotm[ip].Mat2(iv).noalias() +=
1607 pnd(ie, ip) * prop_matrix<2>(Pe(ie, ip, iv, joker, joker));
1608 break;
1609 case 1:
1610 for (Index iv = 0; iv < nv; iv++)
1611 for (Index ie = 0; ie < ne; ie++)
1612 aotm[ip].Mat1(iv).noalias() +=
1613 pnd(ie, ip) * prop_matrix<1>(Pe(ie, ip, iv, joker, joker));
1614 break;
1615 }
1616 }
1617 return aotm;
1618}
1619
1621 const ConstTensor5View& Pe, const ArrayOfMatrix& dpnd_dx) {
1622 const Index ns = Pe.ncols();
1623 const Index nv = Pe.npages();
1624 const Index np = Pe.nbooks();
1625 const Index ne = Pe.nshelves();
1626 const Index nq = dpnd_dx.nelem();
1627
1630
1631 for (Index ip = 0; ip < np; ip++) {
1632 for (Index iq = 0; iq < nq; iq++) {
1633
1634 aoaotm[ip][iq].setZero();
1635
1636 if(not dpnd_dx[iq].empty()) {
1637 switch (ns) {
1638 case 4:
1639 for (Index iv = 0; iv < nv; iv++)
1640 for (Index ie = 0; ie < ne; ie++)
1641 aoaotm[ip][iq].Mat4(iv).noalias() +=
1642 dpnd_dx[iq](ie, ip) * prop_matrix<4>(Pe(ie, ip, iv, joker, joker));
1643 break;
1644 case 3:
1645 for (Index iv = 0; iv < nv; iv++)
1646 for (Index ie = 0; ie < ne; ie++)
1647 aoaotm[ip][iq].Mat3(iv).noalias() +=
1648 dpnd_dx[iq](ie, ip) * prop_matrix<3>(Pe(ie, ip, iv, joker, joker));
1649 break;
1650 case 2:
1651 for (Index iv = 0; iv < nv; iv++)
1652 for (Index ie = 0; ie < ne; ie++)
1653 aoaotm[ip][iq].Mat2(iv).noalias() +=
1654 dpnd_dx[iq](ie, ip) * prop_matrix<2>(Pe(ie, ip, iv, joker, joker));
1655 break;
1656 case 1:
1657 for (Index iv = 0; iv < nv; iv++)
1658 for (Index ie = 0; ie < ne; ie++)
1659 aoaotm[ip][iq].Mat1(iv).noalias() +=
1660 dpnd_dx[iq](ie, ip) * prop_matrix<1>(Pe(ie, ip, iv, joker, joker));
1661 break;
1662 }
1663 }
1664 }
1665 }
1666 return aoaotm;
1667}
1668
1669// TEST CODE END
1670
1671std::ostream& operator<<(std::ostream& os, const TransmissionMatrix& tm) {
1672 for (const auto& T : tm.T4) os << T << '\n';
1673 for (const auto& T : tm.T3) os << T << '\n';
1674 for (const auto& T : tm.T2) os << T << '\n';
1675 for (const auto& T : tm.T1) os << T << '\n';
1676 return os;
1677}
1678
1679std::ostream& operator<<(std::ostream& os,
1680 const ArrayOfTransmissionMatrix& atm) {
1681 for (const auto& T : atm) os << T << '\n';
1682 return os;
1683}
1684
1685std::ostream& operator<<(std::ostream& os,
1687 for (const auto& T : aatm) os << T << '\n';
1688 return os;
1689}
1690
1691std::ostream& operator<<(std::ostream& os, const RadiationVector& rv) {
1692 // Write the transpose because it looks better...
1693 for (const auto& R : rv.R4) os << R.transpose() << '\n';
1694 for (const auto& R : rv.R3) os << R.transpose() << '\n';
1695 for (const auto& R : rv.R2) os << R.transpose() << '\n';
1696 for (const auto& R : rv.R1) os << R.transpose() << '\n';
1697 return os;
1698}
1699
1700std::ostream& operator<<(std::ostream& os, const ArrayOfRadiationVector& arv) {
1701 for (const auto& R : arv) os << R << '\n';
1702 return os;
1703}
1704
1705std::ostream& operator<<(std::ostream& os,
1706 const ArrayOfArrayOfRadiationVector& aarv) {
1707 for (const auto& R : aarv) os << R << '\n';
1708 return os;
1709}
1710
1711std::istream& operator>>(std::istream& is, TransmissionMatrix& tm) {
1712 for (auto& T : tm.T4)
1713 is >> T(0, 0) >> T(0, 1) >> T(0, 2) >> T(0, 3) >>
1714 T(1, 0) >> T(1, 1) >> T(1, 2) >> T(1, 3) >>
1715 T(2, 0) >> T(2, 1) >> T(2, 2) >> T(2, 3) >>
1716 T(3, 0) >> T(3, 1) >> T(3, 2) >> T(3, 3);
1717 for (auto& T : tm.T3)
1718 is >> T(0, 0) >> T(0, 1) >> T(0, 2) >>
1719 T(1, 0) >> T(1, 1) >> T(1, 2) >>
1720 T(2, 0) >> T(2, 1) >> T(2, 2);
1721 for (auto& T : tm.T2) is >> T(0, 0) >> T(0, 1) >>
1722 T(1, 0) >> T(1, 1);
1723 for (auto& T : tm.T1) is >> T(0, 0);
1724 return is;
1725}
1726
1727std::istream& operator>>(std::istream& is, RadiationVector& rv) {
1728 for (auto& R : rv.R4) is >> R[0] >> R[1] >> R[2] >> R[3];
1729 for (auto& R : rv.R3) is >> R[0] >> R[1] >> R[2];
1730 for (auto& R : rv.R2) is >> R[0] >> R[1] >> R[2];
1731 for (auto& R : rv.R1) is >> R[0];
1732 return is;
1733}
1734
1736 const Numeric& r) {
1738 transmat(*this, pm, pm, r); // Slower to compute, faster to implement...
1739}
This can be used to make arrays out of anything.
Definition: array.h:108
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:197
A constant view of a Matrix.
Definition: matpackI.h:1050
A constant view of a Tensor5.
Definition: matpackV.h:141
Index ncols() const noexcept
Definition: matpackV.h:153
Index npages() const noexcept
Definition: matpackV.h:151
Index nbooks() const noexcept
Definition: matpackV.h:150
Index nshelves() const noexcept
Definition: matpackV.h:149
A constant view of a Vector.
Definition: matpackI.h:517
Index NumberOfFrequencies() const
The number of frequencies of the propagation matrix.
VectorView K24(const Index iz=0, const Index ia=0)
Vector view to K(1, 3) elements.
Tensor4 & Data()
Get full view to data.
VectorView K13(const Index iz=0, const Index ia=0)
Vector view to K(0, 2) elements.
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
Index StokesDimensions() const
The stokes dimension of the propagation matrix.
VectorView K34(const Index iz=0, const Index ia=0)
Vector view to K(2, 3) elements.
VectorView K12(const Index iz=0, const Index ia=0)
Vector view to K(0, 1) elements.
VectorView K23(const Index iz=0, const Index ia=0)
Vector view to K(1, 2) elements.
VectorView K14(const Index iz=0, const Index ia=0)
Vector view to K(0, 3) elements.
bool IsRotational(const Index iv=0, const Index iz=0, const Index ia=0) const
False if diagonal element is non-zero in internal Matrix representation.
Stokes vector is as Propagation matrix but only has 4 possible values.
The Vector class.
Definition: matpackI.h:908
Constants of physical expressions as constexpr.
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
constexpr std::size_t mul(Inds... inds) noexcept
Definition: grids.h:10
#define abs(x)
#define ns
#define dx
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
Numeric & real_val(Complex &c) noexcept
Return a reference to the real value of c.
constexpr Numeric imag(Complex c) noexcept
imag
const Joker joker
#define a2
std::complex< Numeric > Complex
#define b2
constexpr Numeric real(Complex c) noexcept
real
Temperature
Definition: jacobian.h:58
invlib::MatrixIdentity< Matrix > Identity
Definition: oem.h:38
#define u
#define d
#define v
#define w
#define a
#define c
#define b
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:739
#define N
Definition: rng.cc:164
Radiation Vector for Stokes dimension 1-4.
std::vector< Eigen::Matrix< double, 1, 1 >, Eigen::aligned_allocator< Eigen::Matrix< double, 1, 1 > > > R1
const Eigen::Matrix< double, 1, 1 > & Vec1(size_t i) const
Return Vector at position.
void setSource(const StokesVector &a, const ConstVectorView &B, const StokesVector &S, Index i)
Set this to source vector at position.
const Eigen::Vector4d & Vec4(size_t i) const
Return Vector at position.
void leftMul(const TransmissionMatrix &T)
Multiply radiation vector from the left.
void SetZero(size_t i)
Set Radiation Vector to Zero at position.
void rem_avg(const RadiationVector &O1, const RadiationVector &O2)
Remove the average of two other RadiationVector from *this.
std::vector< Eigen::Vector2d, Eigen::aligned_allocator< Eigen::Vector2d > > R2
void add_avg(const RadiationVector &O1, const RadiationVector &O2)
Add the average of two other RadiationVector to *this.
const Eigen::Vector3d & Vec3(size_t i) const
Return Vector at position.
std::vector< Eigen::Vector4d, Eigen::aligned_allocator< Eigen::Vector4d > > R4
const Eigen::Vector2d & Vec2(size_t i) const
Return Vector at position.
std::vector< Eigen::Vector3d, Eigen::aligned_allocator< Eigen::Vector3d > > R3
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.
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
std::vector< Eigen::Matrix3d, Eigen::aligned_allocator< Eigen::Matrix3d > > T3
const Eigen::Matrix2d & Mat2(size_t i) const
Get Matrix at position.
std::vector< Eigen::Matrix< double, 1, 1 >, Eigen::aligned_allocator< Eigen::Matrix< double, 1, 1 > > > T1
const Eigen::Matrix4d & Mat4(size_t i) const
Get Matrix at position.
const Eigen::Matrix< double, 1, 1 > & Mat1(size_t i) const
Get Matrix at position.
const Eigen::Matrix3d & Mat3(size_t i) const
Get Matrix at position.
std::vector< Eigen::Matrix2d, Eigen::aligned_allocator< Eigen::Matrix2d > > T2
TransmissionMatrix(Index nf=0, Index stokes=1)
Construct a new Transmission Matrix object.
std::vector< Eigen::Matrix4d, Eigen::aligned_allocator< Eigen::Matrix4d > > T4
void transmat3(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r, const Index iz=0, const Index ia=0) noexcept
void stepwise_source(RadiationVector &J, ArrayOfRadiationVector &dJ, 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.
constexpr Numeric lower_is_considered_zero_for_sinc_likes
ArrayOfTransmissionMatrix bulk_backscatter(const ConstTensor5View &Pe, const ConstMatrixView &pnd)
Bulk back-scattering.
void dtransmat(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dT1=0, const Numeric &dr_dT2=0, const Index it=-1, const Index iz=0, const Index ia=0) noexcept
void dtransmat4(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dT1, const Numeric &dr_dT2, const Index it, const Index iz, const Index ia) noexcept
void dtransmat3(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dT1, const Numeric &dr_dT2, const Index it, const Index iz, const Index ia) noexcept
void transmat4(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r, const Index iz=0, const Index ia=0) noexcept
void transmat2(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r, const Index iz=0, const Index ia=0) noexcept
void dtransmat1(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dT1, const Numeric &dr_dT2, const Index it, const Index iz, const Index ia) noexcept
ArrayOfArrayOfTransmissionMatrix bulk_backscatter_derivative(const ConstTensor5View &Pe, const ArrayOfMatrix &dpnd_dx)
Derivatives of bulk back-scattering
constexpr Eigen::Matrix< Numeric, N, 1 > source_vector(const StokesVector &a, const ConstVectorView &B, const StokesVector &da, const ConstVectorView &dB_dT, const StokesVector &dS, bool dT, Index i)
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
void set_backscatter_radiation_vector(ArrayOfRadiationVector &I, ArrayOfArrayOfArrayOfRadiationVector &dI, const RadiationVector &I_incoming, const ArrayOfTransmissionMatrix &T, const ArrayOfTransmissionMatrix &PiTf, const ArrayOfTransmissionMatrix &PiTr, const ArrayOfTransmissionMatrix &Z, const ArrayOfArrayOfTransmissionMatrix &dT1, const ArrayOfArrayOfTransmissionMatrix &dT2, const ArrayOfArrayOfTransmissionMatrix &dZ, const BackscatterSolver solver)
Set the backscatter radiation vector.
std::istream & operator>>(std::istream &is, TransmissionMatrix &tm)
Input operator.
void update_radiation_vector(RadiationVector &I, ArrayOfRadiationVector &dI1, ArrayOfRadiationVector &dI2, const RadiationVector &J1, const RadiationVector &J2, const ArrayOfRadiationVector &dJ1, const ArrayOfRadiationVector &dJ2, const TransmissionMatrix &T, const TransmissionMatrix &PiT, const ArrayOfTransmissionMatrix &dT1, const ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric r, const Vector &dr1, const Vector &dr2, const Index ia, const Index iz, const RadiativeTransferSolver solver)
Update the Radiation Vector.
void transmat(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r) noexcept
std::ostream & operator<<(std::ostream &os, const TransmissionMatrix &tm)
Output operator.
void stepwise_transmission(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dtemp1, const Numeric &dr_dtemp2, const Index temp_deriv_pos)
Set the stepwise transmission matrix.
void dtransmat2(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dT1, const Numeric &dr_dT2, const Index it, const Index iz, const Index ia) noexcept
void transmat1(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r, const Index iz=0, const Index ia=0) noexcept
Stuff related to the transmission matrix.
RadiativeTransferSolver
Intended to hold various forward solvers.
BackscatterSolver
Intended to hold various backscatter solvers.
CumulativeTransmission
Intended to hold various ways to accumulate the transmission matrix.
Array< TransmissionMatrix > ArrayOfTransmissionMatrix