ARTS 2.5.4 (git: 4c0d3b4d)
propagationmatrix.cc
Go to the documentation of this file.
1/* Copyright (C) 2017
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 "propagationmatrix.h"
30#include "arts_omp.h"
31#include "lin_alg.h"
32
33void compute_transmission_matrix(Tensor3View T,
34 const Numeric& r,
35 const PropagationMatrix& upper_level,
36 const PropagationMatrix& lower_level,
37 const Index iz,
38 const Index ia) {
39 const Index mstokes_dim = upper_level.StokesDimensions();
40 const Index mfreqs = upper_level.NumberOfFrequencies();
41
42 if (mstokes_dim == 1) {
43 for (Index i = 0; i < mfreqs; i++)
44 T(i, 0, 0) = exp(
45 -0.5 * r * (upper_level.Kjj(iz, ia)[i] + lower_level.Kjj(iz, ia)[i]));
46 } else if (mstokes_dim == 2) {
47 for (Index i = 0; i < mfreqs; i++) {
48 MatrixView F = T(i, joker, joker);
49
50 const Numeric a = -0.5 * r *
51 (upper_level.Kjj(iz, ia)[i] +
52 lower_level.Kjj(iz, ia)[i]),
53 b = -0.5 * r *
54 (upper_level.K12(iz, ia)[i] +
55 lower_level.K12(iz, ia)[i]);
56
57 const Numeric exp_a = exp(a);
58
59 if (b == 0.) {
60 F(0, 1) = F(1, 0) = 0.;
61 F(0, 0) = F(1, 1) = exp_a;
62 continue;
63 }
64
65 const Numeric C0 = (b * cosh(b) - a * sinh(b)) / b;
66 const Numeric C1 = sinh(b) / b;
67
68 F(0, 0) = F(1, 1) = C0 + C1 * a;
69 F(0, 1) = F(1, 0) = C1 * b;
70
71 F *= exp_a;
72 }
73 } else if (mstokes_dim == 3) {
74 for (Index i = 0; i < mfreqs; i++) {
75 MatrixView F = T(i, joker, joker);
76
77 const Numeric
78 a = -0.5 * r *
79 (upper_level.Kjj(iz, ia)[i] + lower_level.Kjj(iz, ia)[i]),
80 b = -0.5 * r *
81 (upper_level.K12(iz, ia)[i] + lower_level.K12(iz, ia)[i]),
82 c = -0.5 * r *
83 (upper_level.K13(iz, ia)[i] + lower_level.K13(iz, ia)[i]),
84 u = -0.5 * r *
85 (upper_level.K23(iz, ia)[i] + lower_level.K23(iz, ia)[i]);
86
87 const Numeric exp_a = exp(a);
88
89 if (b == 0. and c == 0. and u == 0.) {
90 F = 0.;
91 F(0, 0) = F(1, 1) = F(2, 2) = exp_a;
92 continue;
93 }
94
95 const Numeric a2 = a * a, b2 = b * b, c2 = c * c, u2 = u * u;
96
97 const Numeric x = sqrt(b2 + c2 - u2), x2 = x * x, inv_x2 = 1.0 / x2;
98 const Numeric sinh_x = sinh(x), cosh_x = cosh(x);
99
100 const Numeric C0 = (a2 * (cosh_x - 1) - a * x * sinh_x + x2) *
101 inv_x2; // approaches (1-a)*exp_a for low x
102 const Numeric C1 = (2 * a * (1 - cosh_x) + x * sinh_x) *
103 inv_x2; // approaches (exp_a) for low_x
104 const Numeric C2 =
105 (cosh_x - 1) * inv_x2; // Approaches infinity for low x
106
107 F(0, 0) = F(1, 1) = F(2, 2) = C0 + C1 * a;
108 F(0, 0) += C2 * (a2 + b2 + c2);
109 F(1, 1) += C2 * (a2 + b2 - u2);
110 F(2, 2) += C2 * (a2 + c2 - u2);
111
112 F(0, 1) = F(1, 0) = C1 * b;
113 F(0, 1) += C2 * (2 * a * b - c * u);
114 F(1, 0) += C2 * (2 * a * b + c * u);
115
116 F(0, 2) = F(2, 0) = C1 * c;
117 F(0, 2) += C2 * (2 * a * c + b * u);
118 F(2, 0) += C2 * (2 * a * c - b * u);
119
120 F(1, 2) = C1 * u + C2 * (2 * a * u + b * c);
121 F(2, 1) = -C1 * u - C2 * (2 * a * u - b * c);
122
123 F *= exp_a;
124 }
125 } else if (mstokes_dim == 4) {
126 static const Numeric sqrt_05 = sqrt(0.5);
127 for (Index i = 0; i < mfreqs; i++) {
128 MatrixView F = T(i, joker, joker);
129
130 const Numeric
131 a = -0.5 * r *
132 (upper_level.Kjj(iz, ia)[i] + lower_level.Kjj(iz, ia)[i]),
133 b = -0.5 * r *
134 (upper_level.K12(iz, ia)[i] + lower_level.K12(iz, ia)[i]),
135 c = -0.5 * r *
136 (upper_level.K13(iz, ia)[i] + lower_level.K13(iz, ia)[i]),
137 d = -0.5 * r *
138 (upper_level.K14(iz, ia)[i] + lower_level.K14(iz, ia)[i]),
139 u = -0.5 * r *
140 (upper_level.K23(iz, ia)[i] + lower_level.K23(iz, ia)[i]),
141 v = -0.5 * r *
142 (upper_level.K24(iz, ia)[i] + lower_level.K24(iz, ia)[i]),
143 w = -0.5 * r *
144 (upper_level.K34(iz, ia)[i] + lower_level.K34(iz, ia)[i]);
145
146 const Numeric exp_a = exp(a);
147
148 if (b == 0. and c == 0. and d == 0. and u == 0. and v == 0. and w == 0.) {
149 F = 0.;
150 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
151 continue;
152 }
153
154 const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
155 w2 = w * w;
156
157 const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
158
159 Numeric Const1;
160 Const1 = b2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
161 c2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
162 d2 * (d2 * 0.5 + u2 - v2 - w2) + u2 * (u2 * 0.5 + v2 + w2) +
163 v2 * (v2 * 0.5 + w2) +
164 4 * (b * d * u * w - b * c * v * w - c * d * u * v);
165 Const1 *= 2;
166 Const1 += w2 * w2;
167
168 if (Const1 > 0.0)
169 Const1 = sqrt(Const1);
170 else
171 Const1 = 0.0;
172
173 const Complex sqrt_BpA = sqrt(Complex(Const2 + Const1, 0.0));
174 const Complex sqrt_BmA = sqrt(Complex(Const2 - Const1, 0.0));
175 const Numeric x = sqrt_BpA.real() * sqrt_05;
176 const Numeric y = sqrt_BmA.imag() * sqrt_05;
177 const Numeric x2 = x * x;
178 const Numeric y2 = y * y;
179 const Numeric cos_y = cos(y);
180 const Numeric sin_y = sin(y);
181 const Numeric cosh_x = cosh(x);
182 const Numeric sinh_x = sinh(x);
183 const Numeric x2y2 = x2 + y2;
184 const Numeric inv_x2y2 = 1.0 / x2y2;
185
186 Numeric C0, C1, C2, C3;
187 Numeric inv_y = 0.0, inv_x = 0.0; // Init'd to remove warnings
188
189 // X and Y cannot both be zero
190 if (x == 0.0) {
191 inv_y = 1.0 / y;
192 C0 = 1.0;
193 C1 = 1.0;
194 C2 = (1.0 - cos_y) * inv_x2y2;
195 C3 = (1.0 - sin_y * inv_y) * inv_x2y2;
196 } else if (y == 0.0) {
197 inv_x = 1.0 / x;
198 C0 = 1.0;
199 C1 = 1.0;
200 C2 = (cosh_x - 1.0) * inv_x2y2;
201 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
202 } else {
203 inv_x = 1.0 / x;
204 inv_y = 1.0 / y;
205
206 C0 = (cos_y * x2 + cosh_x * y2) * inv_x2y2;
207 C1 = (sin_y * x2 * inv_y + sinh_x * y2 * inv_x) * inv_x2y2;
208 C2 = (cosh_x - cos_y) * inv_x2y2;
209 C3 = (sinh_x * inv_x - sin_y * inv_y) * inv_x2y2;
210 }
211
212 // Diagonal Elements
213 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = C0;
214 F(0, 0) += C2 * (b2 + c2 + d2);
215 F(1, 1) += C2 * (b2 - u2 - v2);
216 F(2, 2) += C2 * (c2 - u2 - w2);
217 F(3, 3) += C2 * (d2 - v2 - w2);
218
219 // Linear main-axis polarization
220 F(0, 1) = F(1, 0) = C1 * b;
221 F(0, 1) +=
222 C2 * (-c * u - d * v) +
223 C3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) - v * (b * v + c * w));
224 F(1, 0) += C2 * (c * u + d * v) +
225 C3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) +
226 d * (b * d + u * w));
227
228 // Linear off-axis polarization
229 F(0, 2) = F(2, 0) = C1 * c;
230 F(0, 2) +=
231 C2 * (b * u - d * w) +
232 C3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) - w * (b * v + c * w));
233 F(2, 0) += C2 * (-b * u + d * w) +
234 C3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) +
235 d * (c * d - u * v));
236
237 // Circular polarization
238 F(0, 3) = F(3, 0) = C1 * d;
239 F(0, 3) +=
240 C2 * (b * v + c * w) +
241 C3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) + w * (b * u - d * w));
242 F(3, 0) += C2 * (-b * v - c * w) +
243 C3 * (b * (b * d + u * w) + c * (c * d - u * v) -
244 d * (-d2 + v2 + w2));
245
246 // Circular polarization rotation
247 F(1, 2) = F(2, 1) = C2 * (b * c - v * w);
248 F(1, 2) += C1 * u + C3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
249 w * (b * d + u * w));
250 F(2, 1) += -C1 * u + C3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
251 v * (c * d - u * v));
252
253 // Linear off-axis polarization rotation
254 F(1, 3) = F(3, 1) = C2 * (b * d + u * w);
255 F(1, 3) += C1 * v + C3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
256 w * (b * c - v * w));
257 F(3, 1) += -C1 * v + C3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
258 v * (-d2 + v2 + w2));
259
260 // Linear main-axis polarization rotation
261 F(2, 3) = F(3, 2) = C2 * (c * d - u * v);
262 F(2, 3) += C1 * w + C3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
263 w * (-c2 + u2 + w2));
264 F(3, 2) += -C1 * w + C3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
265 w * (-d2 + v2 + w2));
266
267 F *= exp_a;
268 }
269 }
270}
271
272void compute_transmission_matrix_from_averaged_matrix_at_frequency(
273 MatrixView T,
274 const Numeric& r,
275 const PropagationMatrix& averaged_propagation_matrix,
276 const Index iv,
277 const Index iz,
278 const Index ia) {
279 static const Numeric sqrt_05 = sqrt(0.5);
280 const Index mstokes_dim = averaged_propagation_matrix.StokesDimensions();
281
282 if (mstokes_dim == 1) {
283 T(0, 0) = exp(-r * averaged_propagation_matrix.Kjj(iz, ia)[iv]);
284 } else if (mstokes_dim == 2) {
285 const Numeric a = -r * averaged_propagation_matrix.Kjj(iz, ia)[iv],
286 b = -r * averaged_propagation_matrix.K12(iz, ia)[iv];
287
288 const Numeric exp_a = exp(a);
289
290 if (b == 0.) {
291 T(0, 1) = T(1, 0) = 0.;
292 T(0, 0) = T(1, 1) = exp_a;
293 } else {
294 const Numeric C0 = (b * cosh(b) - a * sinh(b)) / b;
295 const Numeric C1 = sinh(b) / b;
296
297 T(0, 0) = T(1, 1) = C0 + C1 * a;
298 T(0, 1) = T(1, 0) = C1 * b;
299
300 T *= exp_a;
301 }
302 } else if (mstokes_dim == 3) {
303 const Numeric a = -r * averaged_propagation_matrix.Kjj(iz, ia)[iv],
304 b = -r * averaged_propagation_matrix.K12(iz, ia)[iv],
305 c = -r * averaged_propagation_matrix.K13(iz, ia)[iv],
306 u = -r * averaged_propagation_matrix.K23(iz, ia)[iv];
307
308 const Numeric exp_a = exp(a);
309
310 if (b == 0. and c == 0. and u == 0.) {
311 T = 0.;
312 T(0, 0) = T(1, 1) = T(2, 2) = exp_a;
313 } else {
314 const Numeric a2 = a * a, b2 = b * b, c2 = c * c, u2 = u * u;
315
316 const Numeric x = sqrt(b2 + c2 - u2), x2 = x * x, inv_x2 = 1.0 / x2;
317 const Numeric sinh_x = sinh(x), cosh_x = cosh(x);
318
319 const Numeric C0 = (a2 * (cosh_x - 1) - a * x * sinh_x + x2) *
320 inv_x2; // approaches (1-a)*exp_a for low x
321 const Numeric C1 = (2 * a * (1 - cosh_x) + x * sinh_x) *
322 inv_x2; // approaches (exp_a) for low_x
323 const Numeric C2 =
324 (cosh_x - 1) * inv_x2; // Approaches infinity for low x
325
326 T(0, 0) = T(1, 1) = T(2, 2) = C0 + C1 * a;
327 T(0, 0) += C2 * (a2 + b2 + c2);
328 T(1, 1) += C2 * (a2 + b2 - u2);
329 T(2, 2) += C2 * (a2 + c2 - u2);
330
331 T(0, 1) = T(1, 0) = C1 * b;
332 T(0, 1) += C2 * (2 * a * b - c * u);
333 T(1, 0) += C2 * (2 * a * b + c * u);
334
335 T(0, 2) = T(2, 0) = C1 * c;
336 T(0, 2) += C2 * (2 * a * c + b * u);
337 T(2, 0) += C2 * (2 * a * c - b * u);
338
339 T(1, 2) = C1 * u + C2 * (2 * a * u + b * c);
340 T(2, 1) = -C1 * u - C2 * (2 * a * u - b * c);
341
342 T *= exp_a;
343 }
344 } else if (mstokes_dim == 4) {
345 const Numeric a = -r * averaged_propagation_matrix.Kjj(iz, ia)[iv],
346 b = -r * averaged_propagation_matrix.K12(iz, ia)[iv],
347 c = -r * averaged_propagation_matrix.K13(iz, ia)[iv],
348 d = -r * averaged_propagation_matrix.K14(iz, ia)[iv],
349 u = -r * averaged_propagation_matrix.K23(iz, ia)[iv],
350 v = -r * averaged_propagation_matrix.K24(iz, ia)[iv],
351 w = -r * averaged_propagation_matrix.K34(iz, ia)[iv];
352
353 const Numeric exp_a = exp(a);
354
355 if (b == 0. and c == 0. and d == 0. and u == 0. and v == 0. and w == 0.) {
356 T = 0.;
357 T(0, 0) = T(1, 1) = T(2, 2) = T(3, 3) = exp_a;
358 } else {
359 const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
360 w2 = w * w;
361
362 const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
363
364 Numeric Const1;
365 Const1 = b2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
366 Const1 += c2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
367 Const1 += d2 * (d2 * 0.5 + u2 - v2 - w2);
368 Const1 += u2 * (u2 * 0.5 + v2 + w2);
369 Const1 += v2 * (v2 * 0.5 + w2);
370 Const1 *= 2;
371 Const1 += 8 * (b * d * u * w - b * c * v * w - c * d * u * v);
372 Const1 += w2 * w2;
373
374 if (Const1 > 0.0)
375 Const1 = sqrt(Const1);
376 else
377 Const1 = 0.0;
378
379 const Complex sqrt_BpA = sqrt(Complex(Const2 + Const1, 0.0));
380 const Complex sqrt_BmA = sqrt(Complex(Const2 - Const1, 0.0));
381 const Numeric x = sqrt_BpA.real() * sqrt_05;
382 const Numeric y = sqrt_BmA.imag() * sqrt_05;
383 const Numeric x2 = x * x;
384 const Numeric y2 = y * y;
385 const Numeric cos_y = cos(y);
386 const Numeric sin_y = sin(y);
387 const Numeric cosh_x = cosh(x);
388 const Numeric sinh_x = sinh(x);
389 const Numeric x2y2 = x2 + y2;
390 const Numeric inv_x2y2 = 1.0 / x2y2;
391
392 Numeric C0, C1, C2, C3;
393 Numeric inv_y = 0.0, inv_x = 0.0; // Init'd to remove warnings
394
395 // X and Y cannot both be zero
396 if (x == 0.0) {
397 inv_y = 1.0 / y;
398 C0 = 1.0;
399 C1 = 1.0;
400 C2 = (1.0 - cos_y) * inv_x2y2;
401 C3 = (1.0 - sin_y * inv_y) * inv_x2y2;
402 } else if (y == 0.0) {
403 inv_x = 1.0 / x;
404 C0 = 1.0;
405 C1 = 1.0;
406 C2 = (cosh_x - 1.0) * inv_x2y2;
407 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
408 } else {
409 inv_x = 1.0 / x;
410 inv_y = 1.0 / y;
411
412 C0 = (cos_y * x2 + cosh_x * y2) * inv_x2y2;
413 C1 = (sin_y * x2 * inv_y + sinh_x * y2 * inv_x) * inv_x2y2;
414 C2 = (cosh_x - cos_y) * inv_x2y2;
415 C3 = (sinh_x * inv_x - sin_y * inv_y) * inv_x2y2;
416 }
417
418 // Diagonal Elements
419 T(0, 0) = T(1, 1) = T(2, 2) = T(3, 3) = C0;
420 T(0, 0) += C2 * (b2 + c2 + d2);
421 T(1, 1) += C2 * (b2 - u2 - v2);
422 T(2, 2) += C2 * (c2 - u2 - w2);
423 T(3, 3) += C2 * (d2 - v2 - w2);
424
425 // Linear main-axis polarization
426 T(0, 1) = T(1, 0) = C1 * b;
427 T(0, 1) +=
428 C2 * (-c * u - d * v) +
429 C3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) - v * (b * v + c * w));
430 T(1, 0) += C2 * (c * u + d * v) +
431 C3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) +
432 d * (b * d + u * w));
433
434 // Linear off-axis polarization
435 T(0, 2) = T(2, 0) = C1 * c;
436 T(0, 2) +=
437 C2 * (b * u - d * w) +
438 C3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) - w * (b * v + c * w));
439 T(2, 0) += C2 * (-b * u + d * w) +
440 C3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) +
441 d * (c * d - u * v));
442
443 // Circular polarization
444 T(0, 3) = T(3, 0) = C1 * d;
445 T(0, 3) +=
446 C2 * (b * v + c * w) +
447 C3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) + w * (b * u - d * w));
448 T(3, 0) += C2 * (-b * v - c * w) +
449 C3 * (b * (b * d + u * w) + c * (c * d - u * v) -
450 d * (-d2 + v2 + w2));
451
452 // Circular polarization rotation
453 T(1, 2) = T(2, 1) = C2 * (b * c - v * w);
454 T(1, 2) += C1 * u + C3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
455 w * (b * d + u * w));
456 T(2, 1) += -C1 * u + C3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
457 v * (c * d - u * v));
458
459 // Linear off-axis polarization rotation
460 T(1, 3) = T(3, 1) = C2 * (b * d + u * w);
461 T(1, 3) += C1 * v + C3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
462 w * (b * c - v * w));
463 T(3, 1) += -C1 * v + C3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
464 v * (-d2 + v2 + w2));
465
466 // Linear main-axis polarization rotation
467 T(2, 3) = T(3, 2) = C2 * (c * d - u * v);
468 T(2, 3) += C1 * w + C3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
469 w * (-c2 + u2 + w2));
470 T(3, 2) += -C1 * w + C3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
471 w * (-d2 + v2 + w2));
472
473 T *= exp_a;
474 }
475 }
476}
477
478void compute_transmission_matrix_and_derivative(
479 Tensor3View T,
480 Tensor4View dT_dx_upper_level,
481 Tensor4View dT_dx_lower_level,
482 const Numeric& r,
483 const PropagationMatrix& upper_level,
484 const PropagationMatrix& lower_level,
485 const ArrayOfPropagationMatrix& dupper_level_dx,
486 const ArrayOfPropagationMatrix& dlower_level_dx,
487 const Numeric& dr_dTu,
488 const Numeric& dr_dTl,
489 const Index it,
490 const Index iz,
491 const Index ia) {
492 const Index mstokes_dim = upper_level.StokesDimensions();
493 const Index mfreqs = upper_level.NumberOfFrequencies();
494 const Index nppd = dupper_level_dx.nelem();
495
496 if (mstokes_dim == 1) {
497 for (Index i = 0; i < mfreqs; i++) {
498 T(i, 0, 0) = exp(
499 -0.5 * r * (upper_level.Kjj(iz, ia)[i] + lower_level.Kjj(iz, ia)[i]));
500 for (Index j = 0; j < nppd; j++) {
501 if (dupper_level_dx[j].NumberOfFrequencies()) {
502 const Numeric da =
503 -0.5 * (r * dupper_level_dx[j].Kjj(iz, ia)[i] +
504 ((j == it) ? (dr_dTu * (upper_level.Kjj(iz, ia)[i] +
505 lower_level.Kjj(iz, ia)[i]))
506 : 0.0));
507 dT_dx_upper_level(j, i, 0, 0) = T(i, 0, 0) * da;
508 }
509
510 if (dlower_level_dx[j].NumberOfFrequencies()) {
511 const Numeric da =
512 -0.5 * (r * dlower_level_dx[j].Kjj(iz, ia)[i] +
513 ((j == it) ? (dr_dTl * (upper_level.Kjj(iz, ia)[i] +
514 lower_level.Kjj(iz, ia)[i]))
515 : 0.0));
516 dT_dx_lower_level(j, i, 0, 0) = T(i, 0, 0) * da;
517 }
518 }
519 }
520 } else if (mstokes_dim == 2) {
521 for (Index i = 0; i < mfreqs; i++) {
522 MatrixView F = T(i, joker, joker);
523
524 const Numeric a = -0.5 * r *
525 (upper_level.Kjj(iz, ia)[i] +
526 lower_level.Kjj(iz, ia)[i]),
527 b = -0.5 * r *
528 (upper_level.K12(iz, ia)[i] +
529 lower_level.K12(iz, ia)[i]);
530
531 const Numeric exp_a = exp(a);
532
533 if (b == 0.) {
534 F(0, 1) = F(1, 0) = 0.;
535 F(0, 0) = F(1, 1) = exp_a;
536 for (Index j = 0; j < nppd; j++) {
537 if (dupper_level_dx[j].NumberOfFrequencies()) {
538 const Numeric da =
539 -0.5 * (r * dupper_level_dx[j].Kjj(iz, ia)[i] +
540 ((j == it) ? dr_dTu * (upper_level.Kjj(iz, ia)[i] +
541 lower_level.Kjj(iz, ia)[i])
542 : 0.0));
543 dT_dx_upper_level(j, i, joker, joker) = F;
544 dT_dx_upper_level(j, i, 0, 0) = dT_dx_upper_level(j, i, 1, 1) *= da;
545 }
546
547 if (dlower_level_dx[j].NumberOfFrequencies()) {
548 const Numeric da =
549 -0.5 * (r * dlower_level_dx[j].Kjj(iz, ia)[i] +
550 ((j == it) ? dr_dTl * (upper_level.Kjj(iz, ia)[i] +
551 lower_level.Kjj(iz, ia)[i])
552 : 0.0));
553 dT_dx_lower_level(j, i, joker, joker) = F;
554 dT_dx_lower_level(j, i, 0, 0) = dT_dx_lower_level(j, i, 1, 1) *= da;
555 }
556 }
557 continue;
558 }
559
560 const Numeric C0 = cosh(b) - a * sinh(b) / b;
561 const Numeric C1 = sinh(b) / b;
562
563 F(0, 0) = F(1, 1) = C0 + C1 * a;
564 F(0, 1) = F(1, 0) = C1 * b;
565
566 F *= exp_a;
567
568 for (Index j = 0; j < nppd; j++) {
569 if (not dlower_level_dx[j].NumberOfFrequencies()) continue;
570 MatrixView dF = dT_dx_lower_level(j, i, joker, joker);
571
572 const Numeric da = -0.5 *
573 (r * dlower_level_dx[j].Kjj(iz, ia)[i] +
574 ((j == it) ? dr_dTl * (upper_level.Kjj(iz, ia)[i] +
575 lower_level.Kjj(iz, ia)[i])
576 : 0.0)),
577 db = -0.5 *
578 (r * dlower_level_dx[j].K12(iz, ia)[i] +
579 ((j == it) ? dr_dTl * (upper_level.K12(iz, ia)[i] +
580 lower_level.K12(iz, ia)[i])
581 : 0.0));
582
583 const Numeric dC0 = -a * cosh(b) * db / b + a * sinh(b) * db / b / b +
584 sinh(b) * db - sinh(b) * da / b;
585 const Numeric dC1 = (cosh(b) - C1) * db / b;
586
587 dF(0, 0) = dF(1, 1) = (dC0 + C1 * da + dC1 * a) * exp_a + F(0, 0) * da;
588 dF(0, 1) = dF(1, 0) = (C1 * db + dC1 * b) * exp_a + F(0, 1) * da;
589 }
590
591 for (Index j = 0; j < nppd; j++) {
592 if (not dupper_level_dx[j].NumberOfFrequencies()) continue;
593
594 MatrixView dF = dT_dx_upper_level(j, i, joker, joker);
595
596 const Numeric da = -0.5 *
597 (r * dupper_level_dx[j].Kjj(iz, ia)[i] +
598 ((j == it) ? dr_dTu * (upper_level.Kjj(iz, ia)[i] +
599 lower_level.Kjj(iz, ia)[i])
600 : 0.0)),
601 db = -0.5 *
602 (r * dupper_level_dx[j].K12(iz, ia)[i] +
603 ((j == it) ? dr_dTu * (upper_level.K12(iz, ia)[i] +
604 lower_level.K12(iz, ia)[i])
605 : 0.0));
606
607 const Numeric dC0 = -a * cosh(b) * db / b + a * sinh(b) * db / b / b +
608 sinh(b) * db - sinh(b) * da / b;
609 const Numeric dC1 = (cosh(b) - C1) * db / b;
610
611 dF(0, 0) = dF(1, 1) = (dC0 + C1 * da + dC1 * a) * exp_a + F(0, 0) * da;
612 dF(0, 1) = dF(1, 0) = (C1 * db + dC1 * b) * exp_a + F(0, 1) * da;
613 }
614 }
615 } else if (mstokes_dim == 3) {
616 for (Index i = 0; i < mfreqs; i++) {
617 MatrixView F = T(i, joker, joker);
618
619 const Numeric
620 a = -0.5 * r *
621 (upper_level.Kjj(iz, ia)[i] + lower_level.Kjj(iz, ia)[i]),
622 b = -0.5 * r *
623 (upper_level.K12(iz, ia)[i] + lower_level.K12(iz, ia)[i]),
624 c = -0.5 * r *
625 (upper_level.K13(iz, ia)[i] + lower_level.K13(iz, ia)[i]),
626 u = -0.5 * r *
627 (upper_level.K23(iz, ia)[i] + lower_level.K23(iz, ia)[i]);
628
629 const Numeric exp_a = exp(a);
630
631 if (b == 0. and c == 0. and u == 0.) {
632 F(0, 1) = F(1, 0) = F(2, 0) = F(0, 2) = F(2, 1) = F(1, 2) = 0.;
633 F(0, 0) = F(1, 1) = F(2, 2) = exp_a;
634 for (Index j = 0; j < nppd; j++) {
635 if (dupper_level_dx[j].NumberOfFrequencies()) {
636 const Numeric da =
637 -0.5 * (r * dupper_level_dx[j].Kjj(iz, ia)[i] +
638 ((j == it) ? dr_dTu * (upper_level.Kjj(iz, ia)[i] +
639 lower_level.Kjj(iz, ia)[i])
640 : 0.0));
641 dT_dx_upper_level(j, i, joker, joker) = F;
642 dT_dx_upper_level(j, i, 0, 0) = dT_dx_upper_level(j, i, 1, 1) =
643 dT_dx_upper_level(j, i, 2, 2) *= da;
644 }
645
646 if (dlower_level_dx[j].NumberOfFrequencies()) {
647 const Numeric da =
648 -0.5 * (r * dlower_level_dx[j].Kjj(iz, ia)[i] +
649 ((j == it) ? dr_dTl * (upper_level.Kjj(iz, ia)[i] +
650 lower_level.Kjj(iz, ia)[i])
651 : 0.0));
652 dT_dx_lower_level(j, i, joker, joker) = F;
653 dT_dx_lower_level(j, i, 0, 0) = dT_dx_lower_level(j, i, 1, 1) =
654 dT_dx_lower_level(j, i, 2, 2) *= da;
655 }
656 }
657 continue;
658 }
659
660 const Numeric a2 = a * a, b2 = b * b, c2 = c * c, u2 = u * u;
661
662 const Numeric x = sqrt(b2 + c2 - u2), x2 = x * x, inv_x2 = 1.0 / x2;
663 const Numeric sinh_x = sinh(x), cosh_x = cosh(x);
664
665 const Numeric C0 = (a2 * (cosh_x - 1) - a * x * sinh_x) * inv_x2 +
666 1; // approaches (1-a)*exp_a for low x
667 const Numeric C1 = (2 * a * (1 - cosh_x) + x * sinh_x) *
668 inv_x2; // approaches (exp_a) for low_x
669 const Numeric C2 =
670 (cosh_x - 1) * inv_x2; // Approaches infinity for low x
671
672 F(0, 0) = F(1, 1) = F(2, 2) = C0 + C1 * a;
673 F(0, 0) += C2 * (a2 + b2 + c2);
674 F(1, 1) += C2 * (a2 + b2 - u2);
675 F(2, 2) += C2 * (a2 + c2 - u2);
676
677 F(0, 1) = F(1, 0) = C1 * b;
678 F(0, 1) += C2 * (2 * a * b - c * u);
679 F(1, 0) += C2 * (2 * a * b + c * u);
680
681 F(0, 2) = F(2, 0) = C1 * c;
682 F(0, 2) += C2 * (2 * a * c + b * u);
683 F(2, 0) += C2 * (2 * a * c - b * u);
684
685 F(1, 2) = C1 * u + C2 * (2 * a * u + b * c);
686 F(2, 1) = -C1 * u - C2 * (2 * a * u - b * c);
687
688 F *= exp_a;
689
690 for (Index j = 0; j < nppd; j++) {
691 if (not dlower_level_dx[j].NumberOfFrequencies()) continue;
692
693 MatrixView dF = dT_dx_lower_level(j, i, joker, joker);
694
695 const Numeric da = -0.5 *
696 (r * dlower_level_dx[j].Kjj(iz, ia)[i] +
697 ((j == it) ? dr_dTl * (upper_level.Kjj(iz, ia)[i] +
698 lower_level.Kjj(iz, ia)[i])
699 : 0.0)),
700 db = -0.5 *
701 (r * dlower_level_dx[j].K12(iz, ia)[i] +
702 ((j == it) ? dr_dTl * (upper_level.K12(iz, ia)[i] +
703 lower_level.K12(iz, ia)[i])
704 : 0.0)),
705 dc = -0.5 *
706 (r * dlower_level_dx[j].K13(iz, ia)[i] +
707 ((j == it) ? dr_dTl * (upper_level.K13(iz, ia)[i] +
708 lower_level.K13(iz, ia)[i])
709 : 0.0)),
710 du = -0.5 *
711 (r * dlower_level_dx[j].K23(iz, ia)[i] +
712 ((j == it) ? dr_dTl * (upper_level.K23(iz, ia)[i] +
713 lower_level.K23(iz, ia)[i])
714 : 0.0));
715
716 const Numeric da2 = 2 * a * da, db2 = 2 * b * db, dc2 = 2 * c * dc,
717 du2 = 2 * u * du;
718
719 const Numeric dx = (db2 + dc2 - du2) / x / 2;
720
721 const Numeric dC0 =
722 -2 * (C0 - 1) * dx / x +
723 (da2 * (cosh_x - 1) + a2 * sinh_x * dx - a * b * cosh_x * dx -
724 a * sinh_x * dx - b * sinh_x * da) *
725 inv_x2;
726 const Numeric dC1 =
727 -2 * C1 * dx / x + (2 * da * (1 - cosh_x) - 2 * a * sinh_x * dx +
728 x * cosh_x * dx + sinh_x * dx) *
729 inv_x2;
730 const Numeric dC2 = (sinh_x / x - 2 * C2) * dx / x;
731
732 dF(0, 0) = dF(1, 1) = dF(2, 2) = dC0 + dC1 * a + C1 * da;
733 dF(0, 0) += dC2 * (a2 + b2 + c2) + C2 * (da2 + db2 + dc2);
734 dF(1, 1) += dC2 * (a2 + b2 - u2) + C2 * (da2 + db2 - du2);
735 dF(2, 2) += dC2 * (a2 + c2 - u2) + C2 * (da2 + dc2 - du2);
736
737 dF(0, 1) = dF(1, 0) = dC1 * b + C1 * db;
738 dF(0, 1) += dC2 * (2 * a * b - c * u) +
739 C2 * (2 * da * b + 2 * a * db - dc * u - c * du);
740 dF(1, 0) += dC2 * (2 * a * b + c * u) +
741 C2 * (2 * da * b + 2 * a * db + dc * u + c * du);
742
743 dF(0, 2) = dF(2, 0) = dC1 * c + C1 * dc;
744 dF(0, 2) += dC2 * (2 * a * c + b * u) +
745 C2 * (2 * da * c + 2 * a * dc + db * u + b * du);
746 dF(2, 0) += dC2 * (2 * a * c - b * u) +
747 C2 * (2 * da * c + 2 * a * dc - db * u - b * du);
748
749 dF(1, 2) = dC1 * u + C1 * du + dC2 * (2 * a * u + b * c) +
750 C2 * (2 * da * u + 2 * a * du + db * c + b * dc);
751 dF(2, 1) = -dC1 * u - C1 * du - dC2 * (2 * a * u - b * c) -
752 C2 * (2 * da * u + 2 * a * du - db * c - b * dc);
753
754 dF *= exp_a;
755 for (int s1 = 0; s1 < 3; s1++)
756 for (int s2 = 0; s2 < 3; s2++) dF(s1, s2) += F(s1, s2) * da;
757 }
758
759 for (Index j = 0; j < nppd; j++) {
760 if (not dupper_level_dx[j].NumberOfFrequencies()) continue;
761
762 MatrixView dF = dT_dx_upper_level(j, i, joker, joker);
763
764 const Numeric da = -0.5 *
765 (r * dupper_level_dx[j].Kjj(iz, ia)[i] +
766 ((j == it) ? dr_dTu * (upper_level.Kjj(iz, ia)[i] +
767 lower_level.Kjj(iz, ia)[i])
768 : 0.0)),
769 db = -0.5 *
770 (r * dupper_level_dx[j].K12(iz, ia)[i] +
771 ((j == it) ? dr_dTu * (upper_level.K12(iz, ia)[i] +
772 lower_level.K12(iz, ia)[i])
773 : 0.0)),
774 dc = -0.5 *
775 (r * dupper_level_dx[j].K13(iz, ia)[i] +
776 ((j == it) ? dr_dTu * (upper_level.K13(iz, ia)[i] +
777 lower_level.K13(iz, ia)[i])
778 : 0.0)),
779 du = -0.5 *
780 (r * dupper_level_dx[j].K23(iz, ia)[i] +
781 ((j == it) ? dr_dTu * (upper_level.K23(iz, ia)[i] +
782 lower_level.K23(iz, ia)[i])
783 : 0.0));
784
785 const Numeric da2 = 2 * a * da, db2 = 2 * b * db, dc2 = 2 * c * dc,
786 du2 = 2 * u * du;
787
788 const Numeric dx = (db2 + dc2 - du2) / x / 2;
789
790 const Numeric dC0 =
791 -2 * (C0 - 1) * dx / x +
792 (da2 * (cosh_x - 1) + a2 * sinh_x * dx - a * b * cosh_x * dx -
793 a * sinh_x * dx - b * sinh_x * da) *
794 inv_x2;
795 const Numeric dC1 =
796 -2 * C1 * dx / x + (2 * da * (1 - cosh_x) - 2 * a * sinh_x * dx +
797 x * cosh_x * dx + sinh_x * dx) *
798 inv_x2;
799 const Numeric dC2 = (sinh_x / x - 2 * C2) * dx / x;
800
801 dF(0, 0) = dF(1, 1) = dF(2, 2) = dC0 + dC1 * a + C1 * da;
802 dF(0, 0) += dC2 * (a2 + b2 + c2) + C2 * (da2 + db2 + dc2);
803 dF(1, 1) += dC2 * (a2 + b2 - u2) + C2 * (da2 + db2 - du2);
804 dF(2, 2) += dC2 * (a2 + c2 - u2) + C2 * (da2 + dc2 - du2);
805
806 dF(0, 1) = dF(1, 0) = dC1 * b + C1 * db;
807 dF(0, 1) += dC2 * (2 * a * b - c * u) +
808 C2 * (2 * da * b + 2 * a * db - dc * u - c * du);
809 dF(1, 0) += dC2 * (2 * a * b + c * u) +
810 C2 * (2 * da * b + 2 * a * db + dc * u + c * du);
811
812 dF(0, 2) = dF(2, 0) = dC1 * c + C1 * dc;
813 dF(0, 2) += dC2 * (2 * a * c + b * u) +
814 C2 * (2 * da * c + 2 * a * dc + db * u + b * du);
815 dF(2, 0) += dC2 * (2 * a * c - b * u) +
816 C2 * (2 * da * c + 2 * a * dc - db * u - b * du);
817
818 dF(1, 2) = dC1 * u + C1 * du + dC2 * (2 * a * u + b * c) +
819 C2 * (2 * da * u + 2 * a * du + db * c + b * dc);
820 dF(2, 1) = -dC1 * u - C1 * du - dC2 * (2 * a * u - b * c) -
821 C2 * (2 * da * u + 2 * a * du - db * c - b * dc);
822
823 dF *= exp_a;
824 for (int s1 = 0; s1 < 3; s1++)
825 for (int s2 = 0; s2 < 3; s2++) dF(s1, s2) += F(s1, s2) * da;
826 }
827 }
828 } else if (mstokes_dim == 4) {
829 static const Numeric sqrt_05 = sqrt(0.5);
830 for (Index i = 0; i < mfreqs; i++) {
831 MatrixView F = T(i, joker, joker);
832
833 const Numeric
834 a = -0.5 * r *
835 (upper_level.Kjj(iz, ia)[i] + lower_level.Kjj(iz, ia)[i]),
836 b = -0.5 * r *
837 (upper_level.K12(iz, ia)[i] + lower_level.K12(iz, ia)[i]),
838 c = -0.5 * r *
839 (upper_level.K13(iz, ia)[i] + lower_level.K13(iz, ia)[i]),
840 d = -0.5 * r *
841 (upper_level.K14(iz, ia)[i] + lower_level.K14(iz, ia)[i]),
842 u = -0.5 * r *
843 (upper_level.K23(iz, ia)[i] + lower_level.K23(iz, ia)[i]),
844 v = -0.5 * r *
845 (upper_level.K24(iz, ia)[i] + lower_level.K24(iz, ia)[i]),
846 w = -0.5 * r *
847 (upper_level.K34(iz, ia)[i] + lower_level.K34(iz, ia)[i]);
848
849 const Numeric exp_a = exp(a);
850
851 if (b == 0. and c == 0. and d == 0. and u == 0. and v == 0. and w == 0.) {
852 F(0, 1) = F(0, 2) = F(0, 3) = F(1, 0) = F(1, 2) = F(1, 3) = F(2, 0) =
853 F(2, 1) = F(2, 3) = F(3, 0) = F(3, 1) = F(3, 2) = 0.;
854 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
855 for (Index j = 0; j < nppd; j++) {
856 if (dupper_level_dx[j].NumberOfFrequencies()) {
857 const Numeric da =
858 -0.5 * (r * dupper_level_dx[j].Kjj(iz, ia)[i] +
859 ((j == it) ? dr_dTu * (upper_level.Kjj(iz, ia)[i] +
860 lower_level.Kjj(iz, ia)[i])
861 : 0.0));
862 dT_dx_upper_level(j, i, joker, joker) = F;
863 dT_dx_upper_level(j, i, 0, 0) = dT_dx_upper_level(j, i, 1, 1) =
864 dT_dx_upper_level(j, i, 2, 2) = dT_dx_upper_level(j, i, 3, 3) *=
865 da;
866 }
867
868 if (dlower_level_dx[j].NumberOfFrequencies()) {
869 const Numeric da =
870 -0.5 * (r * dlower_level_dx[j].Kjj(iz, ia)[i] +
871 ((j == it) ? dr_dTl * (upper_level.Kjj(iz, ia)[i] +
872 lower_level.Kjj(iz, ia)[i])
873 : 0.0));
874 dT_dx_lower_level(j, i, joker, joker) = F;
875 dT_dx_lower_level(j, i, 0, 0) = dT_dx_lower_level(j, i, 1, 1) =
876 dT_dx_lower_level(j, i, 2, 2) = dT_dx_lower_level(j, i, 3, 3) *=
877 da;
878 }
879 }
880 continue;
881 }
882
883 const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
884 w2 = w * w;
885
886 const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
887
888 Numeric Const1;
889 Const1 = b2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
890 Const1 += c2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
891 Const1 += d2 * (d2 * 0.5 + u2 - v2 - w2);
892 Const1 += u2 * (u2 * 0.5 + v2 + w2);
893 Const1 += v2 * (v2 * 0.5 + w2);
894 Const1 *= 2;
895 Const1 += 8 * (b * d * u * w - b * c * v * w - c * d * u * v);
896 Const1 += w2 * w2;
897
898 if (Const1 > 0.0)
899 Const1 = sqrt(Const1);
900 else
901 Const1 = 0.0;
902
903 const Complex sqrt_BpA = sqrt(Complex(Const2 + Const1, 0.0));
904 const Complex sqrt_BmA = sqrt(Complex(Const2 - Const1, 0.0));
905 const Numeric x = sqrt_BpA.real() * sqrt_05;
906 const Numeric y = sqrt_BmA.imag() * sqrt_05;
907 const Numeric x2 = x * x;
908 const Numeric y2 = y * y;
909 const Numeric cos_y = cos(y);
910 const Numeric sin_y = sin(y);
911 const Numeric cosh_x = cosh(x);
912 const Numeric sinh_x = sinh(x);
913 const Numeric x2y2 = x2 + y2;
914 const Numeric inv_x2y2 = 1.0 / x2y2;
915
916 Numeric C0, C1, C2, C3;
917 Numeric inv_y = 0.0, inv_x = 0.0; // Init'd to remove warnings
918
919 // X and Y cannot both be zero
920 if (x == 0.0) {
921 inv_y = 1.0 / y;
922 C0 = 1.0;
923 C1 = 1.0;
924 C2 = (1.0 - cos_y) * inv_x2y2;
925 C3 = (1.0 - sin_y * inv_y) * inv_x2y2;
926 } else if (y == 0.0) {
927 inv_x = 1.0 / x;
928 C0 = 1.0;
929 C1 = 1.0;
930 C2 = (cosh_x - 1.0) * inv_x2y2;
931 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
932 } else {
933 inv_x = 1.0 / x;
934 inv_y = 1.0 / y;
935
936 C0 = (cos_y * x2 + cosh_x * y2) * inv_x2y2;
937 C1 = (sin_y * x2 * inv_y + sinh_x * y2 * inv_x) * inv_x2y2;
938 C2 = (cosh_x - cos_y) * inv_x2y2;
939 C3 = (sinh_x * inv_x - sin_y * inv_y) * inv_x2y2;
940 }
941
942 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = C0;
943 F(0, 0) += C2 * (b2 + c2 + d2);
944 F(1, 1) += C2 * (b2 - u2 - v2);
945 F(2, 2) += C2 * (c2 - u2 - w2);
946 F(3, 3) += C2 * (d2 - v2 - w2);
947
948 F(0, 1) = F(1, 0) = C1 * b;
949 F(0, 1) +=
950 C2 * (-c * u - d * v) +
951 C3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) - v * (b * v + c * w));
952 F(1, 0) += C2 * (c * u + d * v) +
953 C3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) +
954 d * (b * d + u * w));
955
956 F(0, 2) = F(2, 0) = C1 * c;
957 F(0, 2) +=
958 C2 * (b * u - d * w) +
959 C3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) - w * (b * v + c * w));
960 F(2, 0) += C2 * (-b * u + d * w) +
961 C3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) +
962 d * (c * d - u * v));
963
964 F(0, 3) = F(3, 0) = C1 * d;
965 F(0, 3) +=
966 C2 * (b * v + c * w) +
967 C3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) + w * (b * u - d * w));
968 F(3, 0) += C2 * (-b * v - c * w) +
969 C3 * (b * (b * d + u * w) + c * (c * d - u * v) -
970 d * (-d2 + v2 + w2));
971
972 F(1, 2) = F(2, 1) = C2 * (b * c - v * w);
973 F(1, 2) += C1 * u + C3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
974 w * (b * d + u * w));
975 F(2, 1) += -C1 * u + C3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
976 v * (c * d - u * v));
977
978 F(1, 3) = F(3, 1) = C2 * (b * d + u * w);
979 F(1, 3) += C1 * v + C3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
980 w * (b * c - v * w));
981 F(3, 1) += -C1 * v + C3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
982 v * (-d2 + v2 + w2));
983
984 F(2, 3) = F(3, 2) = C2 * (c * d - u * v);
985 F(2, 3) += C1 * w + C3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
986 w * (-c2 + u2 + w2));
987 F(3, 2) += -C1 * w + C3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
988 w * (-d2 + v2 + w2));
989
990 F *= exp_a;
991
992 if (nppd) {
993 const Numeric inv_x2 = inv_x * inv_x;
994 const Numeric inv_y2 = inv_y * inv_y;
995
996 for (Index j = 0; j < nppd; j++) {
997 if (not dupper_level_dx[j].NumberOfFrequencies()) continue;
998
999 MatrixView dF = dT_dx_upper_level(j, i, joker, joker);
1000
1001 const Numeric
1002 da = -0.5 * (r * dupper_level_dx[j].Kjj(iz, ia)[i] +
1003 ((j == it) ? dr_dTu * (upper_level.Kjj(iz, ia)[i] +
1004 lower_level.Kjj(iz, ia)[i])
1005 : 0.0)),
1006 db = -0.5 * (r * dupper_level_dx[j].K12(iz, ia)[i] +
1007 ((j == it) ? dr_dTu * (upper_level.K12(iz, ia)[i] +
1008 lower_level.K12(iz, ia)[i])
1009 : 0.0)),
1010 dc = -0.5 * (r * dupper_level_dx[j].K13(iz, ia)[i] +
1011 ((j == it) ? dr_dTu * (upper_level.K13(iz, ia)[i] +
1012 lower_level.K13(iz, ia)[i])
1013 : 0.0)),
1014 dd = -0.5 * (r * dupper_level_dx[j].K14(iz, ia)[i] +
1015 ((j == it) ? dr_dTu * (upper_level.K14(iz, ia)[i] +
1016 lower_level.K14(iz, ia)[i])
1017 : 0.0)),
1018 du = -0.5 * (r * dupper_level_dx[j].K23(iz, ia)[i] +
1019 ((j == it) ? dr_dTu * (upper_level.K23(iz, ia)[i] +
1020 lower_level.K23(iz, ia)[i])
1021 : 0.0)),
1022 dv = -0.5 * (r * dupper_level_dx[j].K24(iz, ia)[i] +
1023 ((j == it) ? dr_dTu * (upper_level.K24(iz, ia)[i] +
1024 lower_level.K24(iz, ia)[i])
1025 : 0.0)),
1026 dw = -0.5 * (r * dupper_level_dx[j].K34(iz, ia)[i] +
1027 ((j == it) ? dr_dTu * (upper_level.K34(iz, ia)[i] +
1028 lower_level.K34(iz, ia)[i])
1029 : 0.0));
1030
1031 const Numeric db2 = 2 * db * b, dc2 = 2 * dc * c, dd2 = 2 * dd * d,
1032 du2 = 2 * du * u, dv2 = 2 * dv * v, dw2 = 2 * dw * w;
1033
1034 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1035
1036 Numeric dConst1;
1037 if (Const1 > 0.) {
1038 dConst1 = db2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
1039 dConst1 += b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2);
1040
1041 dConst1 += dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
1042 dConst1 += c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2);
1043
1044 dConst1 += dd2 * (d2 * 0.5 + u2 - v2 - w2);
1045 dConst1 += d2 * (dd2 * 0.5 + du2 - dv2 - dw2);
1046
1047 dConst1 += du2 * (u2 * 0.5 + v2 + w2);
1048 dConst1 += u2 * (du2 * 0.5 + dv2 + dw2);
1049
1050 dConst1 += dv2 * (v2 * 0.5 + w2);
1051 dConst1 += v2 * (dv2 * 0.5 + dw2);
1052
1053 dConst1 += 4 * ((db * d * u * w - db * c * v * w - dc * d * u * v +
1054 b * dd * u * w - b * dc * v * w - c * dd * u * v +
1055 b * d * du * w - b * c * dv * w - c * d * du * v +
1056 b * d * u * dw - b * c * v * dw - c * d * u * dv));
1057 dConst1 += dw2 * w2;
1058 dConst1 /= Const1;
1059 } else
1060 dConst1 = 0.0;
1061
1062 Numeric dC0, dC1, dC2, dC3;
1063 if (x == 0.0) {
1064 const Numeric dy =
1065 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1066
1067 dC0 = 0.0;
1068 dC1 = 0.0;
1069 dC2 = -2 * y * dy * C2 * inv_x2y2 + dy * sin_y * inv_x2y2;
1070 dC3 = -2 * y * dy * C3 * inv_x2y2 +
1071 (dy * sin_y * inv_y2 - cos_y * dy * inv_y) * inv_x2y2;
1072 ;
1073 } else if (y == 0.0) {
1074 const Numeric dx =
1075 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1076
1077 dC0 = 0.0;
1078 dC1 = 0.0;
1079 dC2 = -2 * x * dx * C2 * inv_x2y2 + dx * sinh_x * inv_x2y2;
1080 dC3 = -2 * x * dx * C3 * inv_x2y2 +
1081 (cosh_x * dx * inv_x - dx * sinh_x * inv_x2) * inv_x2y2;
1082 } else {
1083 const Numeric dx =
1084 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1085 const Numeric dy =
1086 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1087 const Numeric dy2 = 2 * y * dy;
1088 const Numeric dx2 = 2 * x * dx;
1089 const Numeric dx2dy2 = dx2 + dy2;
1090
1091 dC0 = -dx2dy2 * C0 * inv_x2y2 +
1092 (2 * cos_y * dx * x + 2 * cosh_x * dy * y + dx * sinh_x * y2 -
1093 dy * sin_y * x2) *
1094 inv_x2y2;
1095
1096 dC1 = -dx2dy2 * C1 * inv_x2y2 +
1097 (cos_y * dy * x2 * inv_y + dx2 * sin_y * inv_y -
1098 dy * sin_y * x2 * inv_y2 - dx * sinh_x * y2 * inv_x2 +
1099 cosh_x * dx * y2 * inv_x + dy2 * sinh_x * inv_x) *
1100 inv_x2y2;
1101
1102 dC2 =
1103 -dx2dy2 * C2 * inv_x2y2 + (dx * sinh_x + dy * sin_y) * inv_x2y2;
1104
1105 dC3 = -dx2dy2 * C3 * inv_x2y2 +
1106 (dy * sin_y * inv_y2 - cos_y * dy * inv_y +
1107 cosh_x * dx * inv_x - dx * sinh_x * inv_x2) *
1108 inv_x2y2;
1109 }
1110
1111 dF(0, 0) = dF(1, 1) = dF(2, 2) = dF(3, 3) = dC0;
1112 dF(0, 0) += dC2 * (b2 + c2 + d2) + C2 * (db2 + dc2 + dd2);
1113 dF(1, 1) += dC2 * (b2 - u2 - v2) + C2 * (db2 - du2 - dv2);
1114 dF(2, 2) += dC2 * (c2 - u2 - w2) + C2 * (dc2 - du2 - dw2);
1115 dF(3, 3) += dC2 * (d2 - v2 - w2) + C2 * (dd2 - dv2 - dw2);
1116
1117 dF(0, 1) = dF(1, 0) = db * C1 + b * dC1;
1118
1119 dF(0, 1) += dC2 * (-c * u - d * v) +
1120 C2 * (-dc * u - dd * v - c * du - d * dv) +
1121 dC3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) -
1122 v * (b * v + c * w)) +
1123 C3 * (db * (b2 + c2 + d2) - du * (b * u - d * w) -
1124 dv * (b * v + c * w) + b * (db2 + dc2 + dd2) -
1125 u * (db * u - dd * w) - v * (db * v + dc * w) -
1126 u * (b * du - d * dw) - v * (b * dv + c * dw));
1127 dF(1, 0) += dC2 * (c * u + d * v) +
1128 C2 * (dc * u + dd * v + c * du + d * dv) +
1129 dC3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) +
1130 d * (b * d + u * w)) +
1131 C3 * (-db * (-b2 + u2 + v2) + dc * (b * c - v * w) +
1132 dd * (b * d + u * w) - b * (-db2 + du2 + dv2) +
1133 c * (db * c - dv * w) + d * (db * d + du * w) +
1134 c * (b * dc - v * dw) + d * (b * dd + u * dw));
1135
1136 dF(0, 2) = dF(2, 0) = dC1 * c + C1 * dc;
1137
1138 dF(0, 2) += dC2 * (b * u - d * w) +
1139 C2 * (db * u - dd * w + b * du - d * dw) +
1140 dC3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) -
1141 w * (b * v + c * w)) +
1142 C3 * (dc * (b2 + c2 + d2) - du * (c * u + d * v) -
1143 dw * (b * v + c * w) + c * (db2 + dc2 + dd2) -
1144 u * (dc * u + dd * v) - w * (db * v + dc * w) -
1145 u * (c * du + d * dv) - w * (b * dv + c * dw));
1146 dF(2, 0) += dC2 * (-b * u + d * w) +
1147 C2 * (-db * u + dd * w - b * du + d * dw) +
1148 dC3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) +
1149 d * (c * d - u * v)) +
1150 C3 * (db * (b * c - v * w) - dc * (-c2 + u2 + w2) +
1151 dd * (c * d - u * v) + b * (db * c - dv * w) -
1152 c * (-dc2 + du2 + dw2) + d * (dc * d - du * v) +
1153 b * (b * dc - v * dw) + d * (c * dd - u * dv));
1154
1155 dF(0, 3) = dF(3, 0) = dC1 * d + C1 * dd;
1156
1157 dF(0, 3) += dC2 * (b * v + c * w) +
1158 C2 * (db * v + dc * w + b * dv + c * dw) +
1159 dC3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) +
1160 w * (b * u - d * w)) +
1161 C3 * (dd * (b2 + c2 + d2) - dv * (c * u + d * v) +
1162 dw * (b * u - d * w) + d * (db2 + dc2 + dd2) -
1163 v * (dc * u + dd * v) + w * (db * u - dd * w) -
1164 v * (c * du + d * dv) + w * (b * du - d * dw));
1165 dF(3, 0) += dC2 * (-b * v - c * w) +
1166 C2 * (-db * v - dc * w - b * dv - c * dw) +
1167 dC3 * (b * (b * d + u * w) + c * (c * d - u * v) -
1168 d * (-d2 + v2 + w2)) +
1169 C3 * (db * (b * d + u * w) + dc * (c * d - u * v) -
1170 dd * (-d2 + v2 + w2) + b * (db * d + du * w) +
1171 c * (dc * d - du * v) - d * (-dd2 + dv2 + dw2) +
1172 b * (b * dd + u * dw) + c * (c * dd - u * dv));
1173
1174 dF(1, 2) = dF(2, 1) =
1175 dC2 * (b * c - v * w) + C2 * (db * c + b * dc - dv * w - v * dw);
1176
1177 dF(1, 2) += dC1 * u + C1 * du +
1178 dC3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
1179 w * (b * d + u * w)) +
1180 C3 * (dc * (c * u + d * v) - du * (-b2 + u2 + v2) -
1181 dw * (b * d + u * w) + c * (dc * u + dd * v) -
1182 u * (-db2 + du2 + dv2) - w * (db * d + du * w) +
1183 c * (c * du + d * dv) - w * (b * dd + u * dw));
1184 dF(2, 1) += -dC1 * u - C1 * du +
1185 dC3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
1186 v * (c * d - u * v)) +
1187 C3 * (-db * (b * u - d * w) + du * (-c2 + u2 + w2) -
1188 dv * (c * d - u * v) - b * (db * u - dd * w) +
1189 u * (-dc2 + du2 + dw2) - v * (dc * d - du * v) -
1190 b * (b * du - d * dw) - v * (c * dd - u * dv));
1191
1192 dF(1, 3) = dF(3, 1) =
1193 dC2 * (b * d + u * w) + C2 * (db * d + b * dd + du * w + u * dw);
1194
1195 dF(1, 3) += dC1 * v + C1 * dv +
1196 dC3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
1197 w * (b * c - v * w)) +
1198 C3 * (dd * (c * u + d * v) - dv * (-b2 + u2 + v2) +
1199 dw * (b * c - v * w) + d * (dc * u + dd * v) -
1200 v * (-db2 + du2 + dv2) + w * (db * c - dv * w) +
1201 d * (c * du + d * dv) + w * (b * dc - v * dw));
1202 dF(3, 1) += -dC1 * v - C1 * dv +
1203 dC3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
1204 v * (-d2 + v2 + w2)) +
1205 C3 * (-db * (b * v + c * w) - du * (c * d - u * v) +
1206 dv * (-d2 + v2 + w2) - b * (db * v + dc * w) -
1207 u * (dc * d - du * v) + v * (-dd2 + dv2 + dw2) -
1208 b * (b * dv + c * dw) - u * (c * dd - u * dv));
1209
1210 dF(2, 3) = dF(3, 2) =
1211 dC2 * (c * d - u * v) + C2 * (dc * d + c * dd - du * v - u * dv);
1212
1213 dF(2, 3) += dC1 * w + C1 * dw +
1214 dC3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
1215 w * (-c2 + u2 + w2)) +
1216 C3 * (-dd * (b * u - d * w) + dv * (b * c - v * w) -
1217 dw * (-c2 + u2 + w2) - d * (db * u - dd * w) +
1218 v * (db * c - dv * w) - w * (-dc2 + du2 + dw2) -
1219 d * (b * du - d * dw) + v * (b * dc - v * dw));
1220 dF(3, 2) += -dC1 * w - C1 * dw +
1221 dC3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
1222 w * (-d2 + v2 + w2)) +
1223 C3 * (-dc * (b * v + c * w) + du * (b * d + u * w) +
1224 dw * (-d2 + v2 + w2) - c * (db * v + dc * w) +
1225 u * (db * d + du * w) + w * (-dd2 + dv2 + dw2) -
1226 c * (b * dv + c * dw) + u * (b * dd + u * dw));
1227
1228 dF *= exp_a;
1229
1230 // Finalize derivation by the chain rule
1231 dF(0, 0) += F(0, 0) * da;
1232 dF(0, 1) += F(0, 1) * da;
1233 dF(0, 2) += F(0, 2) * da;
1234 dF(0, 3) += F(0, 3) * da;
1235 dF(1, 0) += F(1, 0) * da;
1236 dF(1, 1) += F(1, 1) * da;
1237 dF(1, 2) += F(1, 2) * da;
1238 dF(1, 3) += F(1, 3) * da;
1239 dF(2, 0) += F(2, 0) * da;
1240 dF(2, 1) += F(2, 1) * da;
1241 dF(2, 2) += F(2, 2) * da;
1242 dF(2, 3) += F(2, 3) * da;
1243 dF(3, 0) += F(3, 0) * da;
1244 dF(3, 1) += F(3, 1) * da;
1245 dF(3, 2) += F(3, 2) * da;
1246 dF(3, 3) += F(3, 3) * da;
1247 }
1248
1249 for (Index j = 0; j < nppd; j++) {
1250 if (not dlower_level_dx[j].NumberOfFrequencies()) continue;
1251
1252 MatrixView dF = dT_dx_lower_level(j, i, joker, joker);
1253
1254 const Numeric
1255 da = -0.5 * (r * dlower_level_dx[j].Kjj(iz, ia)[i] +
1256 ((j == it) ? dr_dTl * (upper_level.Kjj(iz, ia)[i] +
1257 lower_level.Kjj(iz, ia)[i])
1258 : 0.0)),
1259 db = -0.5 * (r * dlower_level_dx[j].K12(iz, ia)[i] +
1260 ((j == it) ? dr_dTl * (upper_level.K12(iz, ia)[i] +
1261 lower_level.K12(iz, ia)[i])
1262 : 0.0)),
1263 dc = -0.5 * (r * dlower_level_dx[j].K13(iz, ia)[i] +
1264 ((j == it) ? dr_dTl * (upper_level.K13(iz, ia)[i] +
1265 lower_level.K13(iz, ia)[i])
1266 : 0.0)),
1267 dd = -0.5 * (r * dlower_level_dx[j].K14(iz, ia)[i] +
1268 ((j == it) ? dr_dTl * (upper_level.K14(iz, ia)[i] +
1269 lower_level.K14(iz, ia)[i])
1270 : 0.0)),
1271 du = -0.5 * (r * dlower_level_dx[j].K23(iz, ia)[i] +
1272 ((j == it) ? dr_dTl * (upper_level.K23(iz, ia)[i] +
1273 lower_level.K23(iz, ia)[i])
1274 : 0.0)),
1275 dv = -0.5 * (r * dlower_level_dx[j].K24(iz, ia)[i] +
1276 ((j == it) ? dr_dTl * (upper_level.K24(iz, ia)[i] +
1277 lower_level.K24(iz, ia)[i])
1278 : 0.0)),
1279 dw = -0.5 * (r * dlower_level_dx[j].K34(iz, ia)[i] +
1280 ((j == it) ? dr_dTl * (upper_level.K34(iz, ia)[i] +
1281 lower_level.K34(iz, ia)[i])
1282 : 0.0));
1283
1284 const Numeric db2 = 2 * db * b, dc2 = 2 * dc * c, dd2 = 2 * dd * d,
1285 du2 = 2 * du * u, dv2 = 2 * dv * v, dw2 = 2 * dw * w;
1286
1287 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1288
1289 Numeric dConst1;
1290 if (Const1 > 0.) {
1291 dConst1 = db2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
1292 dConst1 += b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2);
1293
1294 dConst1 += dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
1295 dConst1 += c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2);
1296
1297 dConst1 += dd2 * (d2 * 0.5 + u2 - v2 - w2);
1298 dConst1 += d2 * (dd2 * 0.5 + du2 - dv2 - dw2);
1299
1300 dConst1 += du2 * (u2 * 0.5 + v2 + w2);
1301 dConst1 += u2 * (du2 * 0.5 + dv2 + dw2);
1302
1303 dConst1 += dv2 * (v2 * 0.5 + w2);
1304 dConst1 += v2 * (dv2 * 0.5 + dw2);
1305
1306 dConst1 += 4 * ((db * d * u * w - db * c * v * w - dc * d * u * v +
1307 b * dd * u * w - b * dc * v * w - c * dd * u * v +
1308 b * d * du * w - b * c * dv * w - c * d * du * v +
1309 b * d * u * dw - b * c * v * dw - c * d * u * dv));
1310 dConst1 += dw2 * w2;
1311 dConst1 /= Const1;
1312 } else
1313 dConst1 = 0.0;
1314
1315 Numeric dC0, dC1, dC2, dC3;
1316 if (x == 0.0) {
1317 const Numeric dy =
1318 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1319
1320 dC0 = 0.0;
1321 dC1 = 0.0;
1322 dC2 = -2 * y * dy * C2 * inv_x2y2 + dy * sin_y * inv_x2y2;
1323 dC3 = -2 * y * dy * C3 * inv_x2y2 +
1324 (dy * sin_y * inv_y2 - cos_y * dy * inv_y) * inv_x2y2;
1325 ;
1326 } else if (y == 0.0) {
1327 const Numeric dx =
1328 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1329
1330 dC0 = 0.0;
1331 dC1 = 0.0;
1332 dC2 = -2 * x * dx * C2 * inv_x2y2 + dx * sinh_x * inv_x2y2;
1333 dC3 = -2 * x * dx * C3 * inv_x2y2 +
1334 (cosh_x * dx * inv_x - dx * sinh_x * inv_x2) * inv_x2y2;
1335 } else {
1336 const Numeric dx =
1337 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1338 const Numeric dy =
1339 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1340 const Numeric dy2 = 2 * y * dy;
1341 const Numeric dx2 = 2 * x * dx;
1342 const Numeric dx2dy2 = dx2 + dy2;
1343
1344 dC0 = -dx2dy2 * C0 * inv_x2y2 +
1345 (2 * cos_y * dx * x + 2 * cosh_x * dy * y + dx * sinh_x * y2 -
1346 dy * sin_y * x2) *
1347 inv_x2y2;
1348
1349 dC1 = -dx2dy2 * C1 * inv_x2y2 +
1350 (cos_y * dy * x2 * inv_y + dx2 * sin_y * inv_y -
1351 dy * sin_y * x2 * inv_y2 - dx * sinh_x * y2 * inv_x2 +
1352 cosh_x * dx * y2 * inv_x + dy2 * sinh_x * inv_x) *
1353 inv_x2y2;
1354
1355 dC2 =
1356 -dx2dy2 * C2 * inv_x2y2 + (dx * sinh_x + dy * sin_y) * inv_x2y2;
1357
1358 dC3 = -dx2dy2 * C3 * inv_x2y2 +
1359 (dy * sin_y * inv_y2 - cos_y * dy * inv_y +
1360 cosh_x * dx * inv_x - dx * sinh_x * inv_x2) *
1361 inv_x2y2;
1362 }
1363
1364 dF(0, 0) = dF(1, 1) = dF(2, 2) = dF(3, 3) = dC0;
1365 dF(0, 0) += dC2 * (b2 + c2 + d2) + C2 * (db2 + dc2 + dd2);
1366 dF(1, 1) += dC2 * (b2 - u2 - v2) + C2 * (db2 - du2 - dv2);
1367 dF(2, 2) += dC2 * (c2 - u2 - w2) + C2 * (dc2 - du2 - dw2);
1368 dF(3, 3) += dC2 * (d2 - v2 - w2) + C2 * (dd2 - dv2 - dw2);
1369
1370 dF(0, 1) = dF(1, 0) = db * C1 + b * dC1;
1371
1372 dF(0, 1) += dC2 * (-c * u - d * v) +
1373 C2 * (-dc * u - dd * v - c * du - d * dv) +
1374 dC3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) -
1375 v * (b * v + c * w)) +
1376 C3 * (db * (b2 + c2 + d2) - du * (b * u - d * w) -
1377 dv * (b * v + c * w) + b * (db2 + dc2 + dd2) -
1378 u * (db * u - dd * w) - v * (db * v + dc * w) -
1379 u * (b * du - d * dw) - v * (b * dv + c * dw));
1380 dF(1, 0) += dC2 * (c * u + d * v) +
1381 C2 * (dc * u + dd * v + c * du + d * dv) +
1382 dC3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) +
1383 d * (b * d + u * w)) +
1384 C3 * (-db * (-b2 + u2 + v2) + dc * (b * c - v * w) +
1385 dd * (b * d + u * w) - b * (-db2 + du2 + dv2) +
1386 c * (db * c - dv * w) + d * (db * d + du * w) +
1387 c * (b * dc - v * dw) + d * (b * dd + u * dw));
1388
1389 dF(0, 2) = dF(2, 0) = dC1 * c + C1 * dc;
1390
1391 dF(0, 2) += dC2 * (b * u - d * w) +
1392 C2 * (db * u - dd * w + b * du - d * dw) +
1393 dC3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) -
1394 w * (b * v + c * w)) +
1395 C3 * (dc * (b2 + c2 + d2) - du * (c * u + d * v) -
1396 dw * (b * v + c * w) + c * (db2 + dc2 + dd2) -
1397 u * (dc * u + dd * v) - w * (db * v + dc * w) -
1398 u * (c * du + d * dv) - w * (b * dv + c * dw));
1399 dF(2, 0) += dC2 * (-b * u + d * w) +
1400 C2 * (-db * u + dd * w - b * du + d * dw) +
1401 dC3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) +
1402 d * (c * d - u * v)) +
1403 C3 * (db * (b * c - v * w) - dc * (-c2 + u2 + w2) +
1404 dd * (c * d - u * v) + b * (db * c - dv * w) -
1405 c * (-dc2 + du2 + dw2) + d * (dc * d - du * v) +
1406 b * (b * dc - v * dw) + d * (c * dd - u * dv));
1407
1408 dF(0, 3) = dF(3, 0) = dC1 * d + C1 * dd;
1409
1410 dF(0, 3) += dC2 * (b * v + c * w) +
1411 C2 * (db * v + dc * w + b * dv + c * dw) +
1412 dC3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) +
1413 w * (b * u - d * w)) +
1414 C3 * (dd * (b2 + c2 + d2) - dv * (c * u + d * v) +
1415 dw * (b * u - d * w) + d * (db2 + dc2 + dd2) -
1416 v * (dc * u + dd * v) + w * (db * u - dd * w) -
1417 v * (c * du + d * dv) + w * (b * du - d * dw));
1418 dF(3, 0) += dC2 * (-b * v - c * w) +
1419 C2 * (-db * v - dc * w - b * dv - c * dw) +
1420 dC3 * (b * (b * d + u * w) + c * (c * d - u * v) -
1421 d * (-d2 + v2 + w2)) +
1422 C3 * (db * (b * d + u * w) + dc * (c * d - u * v) -
1423 dd * (-d2 + v2 + w2) + b * (db * d + du * w) +
1424 c * (dc * d - du * v) - d * (-dd2 + dv2 + dw2) +
1425 b * (b * dd + u * dw) + c * (c * dd - u * dv));
1426
1427 dF(1, 2) = dF(2, 1) =
1428 dC2 * (b * c - v * w) + C2 * (db * c + b * dc - dv * w - v * dw);
1429
1430 dF(1, 2) += dC1 * u + C1 * du +
1431 dC3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
1432 w * (b * d + u * w)) +
1433 C3 * (dc * (c * u + d * v) - du * (-b2 + u2 + v2) -
1434 dw * (b * d + u * w) + c * (dc * u + dd * v) -
1435 u * (-db2 + du2 + dv2) - w * (db * d + du * w) +
1436 c * (c * du + d * dv) - w * (b * dd + u * dw));
1437 dF(2, 1) += -dC1 * u - C1 * du +
1438 dC3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
1439 v * (c * d - u * v)) +
1440 C3 * (-db * (b * u - d * w) + du * (-c2 + u2 + w2) -
1441 dv * (c * d - u * v) - b * (db * u - dd * w) +
1442 u * (-dc2 + du2 + dw2) - v * (dc * d - du * v) -
1443 b * (b * du - d * dw) - v * (c * dd - u * dv));
1444
1445 dF(1, 3) = dF(3, 1) =
1446 dC2 * (b * d + u * w) + C2 * (db * d + b * dd + du * w + u * dw);
1447
1448 dF(1, 3) += dC1 * v + C1 * dv +
1449 dC3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
1450 w * (b * c - v * w)) +
1451 C3 * (dd * (c * u + d * v) - dv * (-b2 + u2 + v2) +
1452 dw * (b * c - v * w) + d * (dc * u + dd * v) -
1453 v * (-db2 + du2 + dv2) + w * (db * c - dv * w) +
1454 d * (c * du + d * dv) + w * (b * dc - v * dw));
1455 dF(3, 1) += -dC1 * v - C1 * dv +
1456 dC3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
1457 v * (-d2 + v2 + w2)) +
1458 C3 * (-db * (b * v + c * w) - du * (c * d - u * v) +
1459 dv * (-d2 + v2 + w2) - b * (db * v + dc * w) -
1460 u * (dc * d - du * v) + v * (-dd2 + dv2 + dw2) -
1461 b * (b * dv + c * dw) - u * (c * dd - u * dv));
1462
1463 dF(2, 3) = dF(3, 2) =
1464 dC2 * (c * d - u * v) + C2 * (dc * d + c * dd - du * v - u * dv);
1465
1466 dF(2, 3) += dC1 * w + C1 * dw +
1467 dC3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
1468 w * (-c2 + u2 + w2)) +
1469 C3 * (-dd * (b * u - d * w) + dv * (b * c - v * w) -
1470 dw * (-c2 + u2 + w2) - d * (db * u - dd * w) +
1471 v * (db * c - dv * w) - w * (-dc2 + du2 + dw2) -
1472 d * (b * du - d * dw) + v * (b * dc - v * dw));
1473 dF(3, 2) += -dC1 * w - C1 * dw +
1474 dC3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
1475 w * (-d2 + v2 + w2)) +
1476 C3 * (-dc * (b * v + c * w) + du * (b * d + u * w) +
1477 dw * (-d2 + v2 + w2) - c * (db * v + dc * w) +
1478 u * (db * d + du * w) + w * (-dd2 + dv2 + dw2) -
1479 c * (b * dv + c * dw) + u * (b * dd + u * dw));
1480
1481 dF *= exp_a;
1482
1483 // Finalize derivation by the chain rule
1484 dF(0, 0) += F(0, 0) * da;
1485 dF(0, 1) += F(0, 1) * da;
1486 dF(0, 2) += F(0, 2) * da;
1487 dF(0, 3) += F(0, 3) * da;
1488 dF(1, 0) += F(1, 0) * da;
1489 dF(1, 1) += F(1, 1) * da;
1490 dF(1, 2) += F(1, 2) * da;
1491 dF(1, 3) += F(1, 3) * da;
1492 dF(2, 0) += F(2, 0) * da;
1493 dF(2, 1) += F(2, 1) * da;
1494 dF(2, 2) += F(2, 2) * da;
1495 dF(2, 3) += F(2, 3) * da;
1496 dF(3, 0) += F(3, 0) * da;
1497 dF(3, 1) += F(3, 1) * da;
1498 dF(3, 2) += F(3, 2) * da;
1499 dF(3, 3) += F(3, 3) * da;
1500 }
1501 }
1502 }
1503 }
1504}
1505
1506void PropagationMatrix::MatrixInverseAtPosition(MatrixView ret,
1507 const Index iv,
1508 const Index iz,
1509 const Index ia) const {
1510 switch (mstokes_dim) {
1511 case 1:
1512 ret(0, 0) = 1.0 / Kjj(iz, ia)[iv];
1513 break;
1514 case 2: {
1515 const Numeric a = Kjj(iz, ia)[iv], a2 = a * a, b = K12(iz, ia)[iv],
1516 b2 = b * b;
1517
1518 const Numeric f = a2 - b2;
1519
1520 const Numeric div = 1.0 / f;
1521
1522 ret(1, 1) = ret(0, 0) = Kjj(iz, ia)[iv] * div;
1523 ret(1, 0) = ret(0, 1) = -K12(iz, ia)[iv] * div;
1524 } break;
1525 case 3: {
1526 const Numeric a = Kjj(iz, ia)[iv], a2 = a * a, b = K12(iz, ia)[iv],
1527 b2 = b * b, c = K13(iz, ia)[iv], c2 = c * c,
1528 u = K23(iz, ia)[iv], u2 = u * u;
1529
1530 const Numeric f = a * (a2 - b2 - c2 + u2);
1531
1532 const Numeric div = 1.0 / f;
1533
1534 ret(0, 0) = (a2 + u2) * div;
1535 ret(0, 1) = -(a * b + c * u) * div;
1536 ret(0, 2) = (-a * c + b * u) * div;
1537
1538 ret(1, 0) = (-a * b + c * u) * div;
1539 ret(1, 1) = (a2 - c2) * div;
1540 ret(1, 2) = (-a * u + b * c) * div;
1541
1542 ret(2, 0) = -(a * c + b * u) * div;
1543 ret(2, 1) = (a * u + b * c) * div;
1544 ret(2, 2) = (a2 - b2) * div;
1545 } break;
1546 case 4: {
1547 const Numeric a = Kjj(iz, ia)[iv], a2 = a * a, b = K12(iz, ia)[iv],
1548 b2 = b * b, c = K13(iz, ia)[iv], c2 = c * c,
1549 u = K23(iz, ia)[iv], u2 = u * u, d = K14(iz, ia)[iv],
1550 d2 = d * d, v = K24(iz, ia)[iv], v2 = v * v,
1551 w = K34(iz, ia)[iv], w2 = w * w;
1552
1553 const Numeric f = a2 * a2 - a2 * b2 - a2 * c2 - a2 * d2 + a2 * u2 +
1554 a2 * v2 + a2 * w2 - b2 * w2 + 2 * b * c * v * w -
1555 2 * b * d * u * w - c2 * v2 + 2 * c * d * u * v -
1556 d2 * u2;
1557
1558 const Numeric div = 1.0 / f;
1559
1560 ret(0, 0) = a * (a2 + u2 + v2 + w2) * div;
1561 ret(0, 1) =
1562 (-a2 * b - a * c * u - a * d * v - b * w2 + c * v * w - d * u * w) *
1563 div;
1564 ret(0, 2) =
1565 (-a2 * c + a * b * u - a * d * w + b * v * w - c * v2 + d * u * v) *
1566 div;
1567 ret(0, 3) =
1568 (-a2 * d + a * b * v + a * c * w - b * u * w + c * u * v - d * u2) *
1569 div;
1570
1571 ret(1, 0) =
1572 (-a2 * b + a * c * u + a * d * v - b * w2 + c * v * w - d * u * w) *
1573 div;
1574 ret(1, 1) = a * (a2 - c2 - d2 + w2) * div;
1575 ret(1, 2) =
1576 (-a2 * u + a * b * c - a * v * w + b * d * w - c * d * v + d2 * u) *
1577 div;
1578 ret(1, 3) =
1579 (-a2 * v + a * b * d + a * u * w - b * c * w + c2 * v - c * d * u) *
1580 div;
1581
1582 ret(2, 0) =
1583 (-a2 * c - a * b * u + a * d * w + b * v * w - c * v2 + d * u * v) *
1584 div;
1585 ret(2, 1) =
1586 (a2 * u + a * b * c - a * v * w - b * d * w + c * d * v - d2 * u) *
1587 div;
1588 ret(2, 2) = a * (a2 - b2 - d2 + v2) * div;
1589 ret(2, 3) =
1590 (-a2 * w + a * c * d - a * u * v + b2 * w - b * c * v + b * d * u) *
1591 div;
1592
1593 ret(3, 0) =
1594 (-a2 * d - a * b * v - a * c * w - b * u * w + c * u * v - d * u2) *
1595 div;
1596 ret(3, 1) =
1597 (a2 * v + a * b * d + a * u * w + b * c * w - c2 * v + c * d * u) *
1598 div;
1599 ret(3, 2) =
1600 (a2 * w + a * c * d - a * u * v - b2 * w + b * c * v - b * d * u) *
1601 div;
1602 ret(3, 3) = a * (a2 - b2 - c2 + u2) * div;
1603 } break;
1604 default:
1605 ARTS_ASSERT(false, "Strange stokes dimensions");
1606 break;
1607 }
1608}
1609
1610void PropagationMatrix::AddAverageAtPosition(const ConstMatrixView& mat1,
1611 const ConstMatrixView& mat2,
1612 const Index iv,
1613 const Index iz,
1614 const Index ia) {
1615 switch (mstokes_dim) {
1616 case 4:
1617 mdata(ia, iz, iv, 3) += (mat1(3, 0) + mat2(3, 0)) * 0.5;
1618 mdata(ia, iz, iv, 5) += (mat1(1, 3) + mat2(1, 3)) * 0.5;
1619 mdata(ia, iz, iv, 6) += (mat1(2, 3) + mat2(2, 3)) * 0.5; /* FALLTHROUGH */
1620 case 3:
1621 mdata(ia, iz, iv, 2) += (mat1(2, 0) + mat2(2, 0)) * 0.5;
1622 mdata(ia, iz, iv, mstokes_dim) +=
1623 (mat1(1, 2) + mat2(1, 2)) * 0.5; /* FALLTHROUGH */
1624 case 2:
1625 mdata(ia, iz, iv, 1) += (mat1(1, 0) + mat2(1, 0)) * 0.5; /* FALLTHROUGH */
1626 case 1:
1627 mdata(ia, iz, iv, 0) += (mat1(0, 0) + mat2(0, 0)) * 0.5; /* FALLTHROUGH */
1628 }
1629}
1630
1631void PropagationMatrix::MultiplyAndAdd(const Numeric x,
1632 const PropagationMatrix& y) {
1633 for (Index i = 0; i < maa; i++)
1634 for (Index j = 0; j < mza; j++)
1635 for (Index k = 0; k < mfreqs; k++)
1636 for (Index l = 0; l < NumberOfNeededVectors(); l++)
1637 mdata(i, j, k, l) += x * y.mdata(i, j, k, l);
1638}
1639
1640bool PropagationMatrix::FittingShape(const ConstMatrixView& x) const {
1641 Index nelem = x.nrows();
1642 if (mstokes_dim == nelem) {
1643 if (x.ncols() == nelem) {
1644 switch (mstokes_dim) {
1645 case 1:
1646 return true;
1647 break;
1648 case 2:
1649 if (x(0, 0) == x(1, 1)) return true;
1650 break;
1651 case 3:
1652 if (x(0, 0) == x(1, 1) and x(0, 0) == x(2, 2) and
1653 x(0, 1) == x(1, 0) and x(2, 0) == x(0, 2) and x(1, 2) == -x(2, 1))
1654 return true;
1655 break;
1656 case 4:
1657 if (x(0, 0) == x(1, 1) and x(0, 0) == x(2, 2) and
1658 x(0, 0) == x(3, 3) and x(0, 1) == x(1, 0) and
1659 x(2, 0) == x(0, 2) and x(3, 0) == x(0, 3) and
1660 x(1, 2) == -x(2, 1) and x(1, 3) == -x(3, 1) and
1661 x(3, 2) == -x(2, 3))
1662 return true;
1663 break;
1664 default:
1665 ARTS_ASSERT(false, "Stokes dimension does not agree with accepted values");
1666 break;
1667 }
1668 }
1669 }
1670 return false;
1671}
1672
1673void PropagationMatrix::GetTensor3(Tensor3View tensor3, Index iz, Index ia) {
1674 switch (mstokes_dim) {
1675 case 4:
1676 tensor3(joker, 3, 3) = mdata(ia, iz, joker, 0);
1677 tensor3(joker, 3, 2) = mdata(ia, iz, joker, 6);
1678 tensor3(joker, 3, 1) = mdata(ia, iz, joker, 5);
1679 tensor3(joker, 0, 3) = mdata(ia, iz, joker, 3);
1680 tensor3(joker, 3, 0) = mdata(ia, iz, joker, 3);
1681 tensor3(joker, 2, 3) = mdata(ia, iz, joker, 6);
1682 tensor3(joker, 1, 3) = mdata(ia, iz, joker, 5);
1683 tensor3(joker, 3, 2) *= -1;
1684 tensor3(joker, 3, 1) *= -1; /* FALLTHROUGH */
1685 case 3:
1686 tensor3(joker, 2, 2) = mdata(ia, iz, joker, 0);
1687 tensor3(joker, 2, 1) = mdata(ia, iz, joker, mstokes_dim);
1688 tensor3(joker, 2, 0) = mdata(ia, iz, joker, 2);
1689 tensor3(joker, 0, 2) = mdata(ia, iz, joker, 2);
1690 tensor3(joker, 1, 2) = mdata(ia, iz, joker, mstokes_dim);
1691 tensor3(joker, 2, 1) *= -1; /* FALLTHROUGH */
1692 case 2:
1693 tensor3(joker, 1, 1) = mdata(ia, iz, joker, 0);
1694 tensor3(joker, 1, 0) = mdata(ia, iz, joker, 1);
1695 tensor3(joker, 0, 1) = mdata(ia, iz, joker, 1); /* FALLTHROUGH */
1696 case 1:
1697 tensor3(joker, 0, 0) = mdata(ia, iz, joker, 0);
1698 break;
1699 default:
1700 ARTS_ASSERT(false, "Stokes dimension does not agree with accepted values");
1701 break;
1702 }
1703}
1704
1705void PropagationMatrix::LeftMultiplyAtPosition(MatrixView out,
1706 const ConstMatrixView& in,
1707 const Index iv,
1708 const Index iz,
1709 const Index ia) const {
1710 switch (mstokes_dim) {
1711 case 1:
1712 out(0, 0) = Kjj(iz, ia)[iv] * in(0, 0);
1713 break;
1714 case 2: {
1715 const Numeric a = Kjj(iz, ia)[iv], b = K12(iz, ia)[iv], m11 = in(0, 0),
1716 m12 = in(0, 1), m21 = in(1, 0), m22 = in(1, 1);
1717 out(0, 0) = a * m11 + b * m21;
1718 out(0, 1) = a * m12 + b * m22;
1719 out(1, 0) = a * m21 + b * m11;
1720 out(1, 1) = a * m22 + b * m12;
1721 } break;
1722 case 3: {
1723 const Numeric a = Kjj(iz, ia)[iv], b = K12(iz, ia)[iv],
1724 c = K13(iz, ia)[iv], u = K23(iz, ia)[iv], m11 = in(0, 0),
1725 m12 = in(0, 1), m13 = in(0, 2), m21 = in(1, 0),
1726 m22 = in(1, 1), m23 = in(1, 2), m31 = in(2, 0),
1727 m32 = in(2, 1), m33 = in(2, 2);
1728 out(0, 0) = a * m11 + b * m21 + c * m31;
1729 out(0, 1) = a * m12 + b * m22 + c * m32;
1730 out(0, 2) = a * m13 + b * m23 + c * m33;
1731 out(1, 0) = a * m21 + b * m11 + m31 * u;
1732 out(1, 1) = a * m22 + b * m12 + m32 * u;
1733 out(1, 2) = a * m23 + b * m13 + m33 * u;
1734 out(2, 0) = a * m31 + c * m11 - m21 * u;
1735 out(2, 1) = a * m32 + c * m12 - m22 * u;
1736 out(2, 2) = a * m33 + c * m13 - m23 * u;
1737 } break;
1738 case 4: {
1739 const Numeric a = Kjj(iz, ia)[iv], b = K12(iz, ia)[iv],
1740 c = K13(iz, ia)[iv], u = K23(iz, ia)[iv],
1741 d = K14(iz, ia)[iv], v = K24(iz, ia)[iv],
1742 w = K34(iz, ia)[iv], m11 = in(0, 0), m12 = in(0, 1),
1743 m13 = in(0, 2), m14 = in(0, 3), m21 = in(1, 0),
1744 m22 = in(1, 1), m23 = in(1, 2), m24 = in(1, 3),
1745 m31 = in(2, 0), m32 = in(2, 1), m33 = in(2, 2),
1746 m34 = in(2, 3), m41 = in(3, 0), m42 = in(3, 1),
1747 m43 = in(3, 2), m44 = in(3, 3);
1748 out(0, 0) = a * m11 + b * m21 + c * m31 + d * m41;
1749 out(0, 1) = a * m12 + b * m22 + c * m32 + d * m42;
1750 out(0, 2) = a * m13 + b * m23 + c * m33 + d * m43;
1751 out(0, 3) = a * m14 + b * m24 + c * m34 + d * m44;
1752 out(1, 0) = a * m21 + b * m11 + m31 * u + m41 * v;
1753 out(1, 1) = a * m22 + b * m12 + m32 * u + m42 * v;
1754 out(1, 2) = a * m23 + b * m13 + m33 * u + m43 * v;
1755 out(1, 3) = a * m24 + b * m14 + m34 * u + m44 * v;
1756 out(2, 0) = a * m31 + c * m11 - m21 * u + m41 * w;
1757 out(2, 1) = a * m32 + c * m12 - m22 * u + m42 * w;
1758 out(2, 2) = a * m33 + c * m13 - m23 * u + m43 * w;
1759 out(2, 3) = a * m34 + c * m14 - m24 * u + m44 * w;
1760 out(3, 0) = a * m41 + d * m11 - m21 * v - m31 * w;
1761 out(3, 1) = a * m42 + d * m12 - m22 * v - m32 * w;
1762 out(3, 2) = a * m43 + d * m13 - m23 * v - m33 * w;
1763 out(3, 3) = a * m44 + d * m14 - m24 * v - m34 * w;
1764 }
1765 }
1766}
1767
1768void PropagationMatrix::RightMultiplyAtPosition(MatrixView out,
1769 const ConstMatrixView& in,
1770 const Index iv,
1771 const Index iz,
1772 const Index ia) const {
1773 switch (mstokes_dim) {
1774 case 1:
1775 out(0, 0) = in(0, 0) * Kjj(iz, ia)[iv];
1776 break;
1777 case 2: {
1778 const Numeric a = Kjj(iz, ia)[iv], b = K12(iz, ia)[iv], m11 = in(0, 0),
1779 m12 = in(0, 1), m21 = in(1, 0), m22 = in(1, 1);
1780 out(0, 0) = a * m11 + b * m12;
1781 out(0, 1) = a * m12 + b * m11;
1782 out(1, 0) = a * m21 + b * m22;
1783 out(1, 1) = a * m22 + b * m21;
1784 } break;
1785 case 3: {
1786 const Numeric a = Kjj(iz, ia)[iv], b = K12(iz, ia)[iv],
1787 c = K13(iz, ia)[iv], u = K23(iz, ia)[iv], m11 = in(0, 0),
1788 m12 = in(0, 1), m13 = in(0, 2), m21 = in(1, 0),
1789 m22 = in(1, 1), m23 = in(1, 2), m31 = in(2, 0),
1790 m32 = in(2, 1), m33 = in(2, 2);
1791 out(0, 0) = a * m11 + b * m12 + c * m13;
1792 out(0, 1) = a * m12 + b * m11 - m13 * u;
1793 out(0, 2) = a * m13 + c * m11 + m12 * u;
1794 out(1, 0) = a * m21 + b * m22 + c * m23;
1795 out(1, 1) = a * m22 + b * m21 - m23 * u;
1796 out(1, 2) = a * m23 + c * m21 + m22 * u;
1797 out(2, 0) = a * m31 + b * m32 + c * m33;
1798 out(2, 1) = a * m32 + b * m31 - m33 * u;
1799 out(2, 2) = a * m33 + c * m31 + m32 * u;
1800 } break;
1801 case 4: {
1802 const Numeric a = Kjj(iz, ia)[iv], b = K12(iz, ia)[iv],
1803 c = K13(iz, ia)[iv], u = K23(iz, ia)[iv],
1804 d = K14(iz, ia)[iv], v = K24(iz, ia)[iv],
1805 w = K34(iz, ia)[iv], m11 = in(0, 0), m12 = in(0, 1),
1806 m13 = in(0, 2), m14 = in(0, 3), m21 = in(1, 0),
1807 m22 = in(1, 1), m23 = in(1, 2), m24 = in(1, 3),
1808 m31 = in(2, 0), m32 = in(2, 1), m33 = in(2, 2),
1809 m34 = in(2, 3), m41 = in(3, 0), m42 = in(3, 1),
1810 m43 = in(3, 2), m44 = in(3, 3);
1811 out(0, 0) = a * m11 + b * m12 + c * m13 + d * m14;
1812 out(0, 1) = a * m12 + b * m11 - m13 * u - m14 * v;
1813 out(0, 2) = a * m13 + c * m11 + m12 * u - m14 * w;
1814 out(0, 3) = a * m14 + d * m11 + m12 * v + m13 * w;
1815 out(1, 0) = a * m21 + b * m22 + c * m23 + d * m24;
1816 out(1, 1) = a * m22 + b * m21 - m23 * u - m24 * v;
1817 out(1, 2) = a * m23 + c * m21 + m22 * u - m24 * w;
1818 out(1, 3) = a * m24 + d * m21 + m22 * v + m23 * w;
1819 out(2, 0) = a * m31 + b * m32 + c * m33 + d * m34;
1820 out(2, 1) = a * m32 + b * m31 - m33 * u - m34 * v;
1821 out(2, 2) = a * m33 + c * m31 + m32 * u - m34 * w;
1822 out(2, 3) = a * m34 + d * m31 + m32 * v + m33 * w;
1823 out(3, 0) = a * m41 + b * m42 + c * m43 + d * m44;
1824 out(3, 1) = a * m42 + b * m41 - m43 * u - m44 * v;
1825 out(3, 2) = a * m43 + c * m41 + m42 * u - m44 * w;
1826 out(3, 3) = a * m44 + d * m41 + m42 * v + m43 * w;
1827 }
1828 }
1829}
1830
1831void PropagationMatrix::RemoveAtPosition(const ConstMatrixView& x,
1832 const Index iv,
1833 const Index iz,
1834 const Index ia) {
1835 switch (mstokes_dim) {
1836 case 4:
1837 mdata(ia, iz, iv, 5) -= x(1, 3);
1838 mdata(ia, iz, iv, 6) -= x(2, 3);
1839 mdata(ia, iz, iv, 3) -= x(0, 3); /* FALLTHROUGH */
1840 case 3:
1841 mdata(ia, iz, iv, 2) -= x(0, 2);
1842 mdata(ia, iz, iv, mstokes_dim) -= x(1, 2); /* FALLTHROUGH */
1843 case 2:
1844 mdata(ia, iz, iv, 1) -= x(0, 1); /* FALLTHROUGH */
1845 case 1:
1846 mdata(ia, iz, iv, 0) -= x(0, 0); /* FALLTHROUGH */
1847 }
1848}
1849
1850void PropagationMatrix::AddAtPosition(const ConstMatrixView& x,
1851 const Index iv,
1852 const Index iz,
1853 const Index ia) {
1854 switch (mstokes_dim) {
1855 case 4:
1856 mdata(ia, iz, iv, 5) += x(1, 3);
1857 mdata(ia, iz, iv, 6) += x(2, 3);
1858 mdata(ia, iz, iv, 3) += x(0, 3); /* FALLTHROUGH */
1859 case 3:
1860 mdata(ia, iz, iv, 2) += x(0, 2);
1861 mdata(ia, iz, iv, mstokes_dim) += x(1, 2); /* FALLTHROUGH */
1862 case 2:
1863 mdata(ia, iz, iv, 1) += x(0, 1); /* FALLTHROUGH */
1864 case 1:
1865 mdata(ia, iz, iv, 0) += x(0, 0); /* FALLTHROUGH */
1866 }
1867}
1868
1869void PropagationMatrix::MultiplyAtPosition(const ConstMatrixView& x,
1870 const Index iv,
1871 const Index iz,
1872 const Index ia) {
1873 switch (mstokes_dim) {
1874 case 4:
1875 mdata(ia, iz, iv, 5) *= x(1, 3);
1876 mdata(ia, iz, iv, 6) *= x(2, 3);
1877 mdata(ia, iz, iv, 3) *= x(0, 3); /* FALLTHROUGH */
1878 case 3:
1879 mdata(ia, iz, iv, 2) *= x(0, 2);
1880 mdata(ia, iz, iv, mstokes_dim) *= x(1, 2); /* FALLTHROUGH */
1881 case 2:
1882 mdata(ia, iz, iv, 1) *= x(0, 1); /* FALLTHROUGH */
1883 case 1:
1884 mdata(ia, iz, iv, 0) *= x(0, 0); /* FALLTHROUGH */
1885 }
1886}
1887
1888void PropagationMatrix::DivideAtPosition(const ConstMatrixView& x,
1889 const Index iv,
1890 const Index iz,
1891 const Index ia) {
1892 switch (mstokes_dim) {
1893 case 4:
1894 mdata(ia, iz, iv, 5) /= x(1, 3);
1895 mdata(ia, iz, iv, 6) /= x(2, 3);
1896 mdata(ia, iz, iv, 3) /= x(0, 3); /* FALLTHROUGH */
1897 case 3:
1898 mdata(ia, iz, iv, 2) /= x(0, 2);
1899 mdata(ia, iz, iv, mstokes_dim) /= x(1, 2); /* FALLTHROUGH */
1900 case 2:
1901 mdata(ia, iz, iv, 1) /= x(0, 1); /* FALLTHROUGH */
1902 case 1:
1903 mdata(ia, iz, iv, 0) /= x(0, 0); /* FALLTHROUGH */
1904 }
1905}
1906
1907void PropagationMatrix::SetAtPosition(const ConstMatrixView& x,
1908 const Index iv,
1909 const Index iz,
1910 const Index ia) {
1911 switch (mstokes_dim) {
1912 case 4:
1913 mdata(ia, iz, iv, 5) = x(1, 3);
1914 mdata(ia, iz, iv, 6) = x(2, 3);
1915 mdata(ia, iz, iv, 3) = x(0, 3); /* FALLTHROUGH */
1916 case 3:
1917 mdata(ia, iz, iv, 2) = x(0, 2);
1918 mdata(ia, iz, iv, mstokes_dim) = x(1, 2); /* FALLTHROUGH */
1919 case 2:
1920 mdata(ia, iz, iv, 1) = x(0, 1); /* FALLTHROUGH */
1921 case 1:
1922 mdata(ia, iz, iv, 0) = x(0, 0); /* FALLTHROUGH */
1923 }
1924}
1925
1926void PropagationMatrix::MatrixAtPosition(MatrixView ret,
1927 const Index iv,
1928 const Index iz,
1929 const Index ia) const {
1930 switch (mstokes_dim) {
1931 case 4:
1932 ret(3, 3) = mdata(ia, iz, iv, 0);
1933 ret(3, 1) = -mdata(ia, iz, iv, 5);
1934 ret(1, 3) = mdata(ia, iz, iv, 5);
1935 ret(3, 2) = -mdata(ia, iz, iv, 6);
1936 ret(2, 3) = mdata(ia, iz, iv, 6);
1937 ret(0, 3) = ret(3, 0) = mdata(ia, iz, iv, 3); /* FALLTHROUGH */
1938 case 3:
1939 ret(2, 2) = mdata(ia, iz, iv, 0);
1940 ret(2, 1) = -mdata(ia, iz, iv, 3);
1941 ret(1, 2) = mdata(ia, iz, iv, 3);
1942 ret(2, 0) = ret(0, 2) = mdata(ia, iz, iv, 2); /* FALLTHROUGH */
1943 case 2:
1944 ret(1, 1) = mdata(ia, iz, iv, 0);
1945 ret(1, 0) = ret(0, 1) = mdata(ia, iz, iv, 1); /* FALLTHROUGH */
1946 case 1:
1947 ret(0, 0) = mdata(ia, iz, iv, 0); /* FALLTHROUGH */
1948 }
1949}
1950
1951#pragma GCC diagnostic push
1952#pragma GCC diagnostic ignored "-Wreturn-type"
1953Numeric PropagationMatrix::operator()(const Index iv,
1954 const Index is1,
1955 const Index is2,
1956 const Index iz,
1957 const Index ia) const {
1958 switch (is1) {
1959 case 0:
1960 switch (is2) {
1961 case 0:
1962 return mdata(ia, iz, iv, 0);
1963 break;
1964 case 1:
1965 return mdata(ia, iz, iv, 1);
1966 break;
1967 case 2:
1968 return mdata(ia, iz, iv, 2);
1969 break;
1970 case 3:
1971 return mdata(ia, iz, iv, 3);
1972 break;
1973 default:
1974 ARTS_ASSERT(false, "out of bounds");
1975 }
1976 break;
1977 case 1:
1978 switch (is2) {
1979 case 0:
1980 return mdata(ia, iz, iv, 1);
1981 break;
1982 case 1:
1983 return mdata(ia, iz, iv, 0);
1984 break;
1985 case 2:
1986 return mdata(ia, iz, iv, mstokes_dim);
1987 break;
1988 case 3:
1989 return mdata(ia, iz, iv, 5);
1990 break;
1991 default:
1992 ARTS_ASSERT(false, "out of bounds");
1993 }
1994 case 2:
1995 switch (is2) {
1996 case 0:
1997 return mdata(ia, iz, iv, 2);
1998 break;
1999 case 1:
2000 return -mdata(ia, iz, iv, mstokes_dim);
2001 break;
2002 case 2:
2003 return mdata(ia, iz, iv, 0);
2004 break;
2005 case 3:
2006 return mdata(ia, iz, iv, 6);
2007 break;
2008 default:
2009 ARTS_ASSERT(false, "out of bounds");
2010 }
2011 break;
2012 case 3:
2013 switch (is2) {
2014 case 0:
2015 return mdata(ia, iz, iv, 3);
2016 break;
2017 case 1:
2018 return -mdata(ia, iz, iv, 5);
2019 break;
2020 case 2:
2021 return -mdata(ia, iz, iv, 6);
2022 break;
2023 case 3:
2024 return mdata(ia, iz, iv, 0);
2025 break;
2026 default:
2027 ARTS_ASSERT(false, "out of bounds");
2028 }
2029 break;
2030 default:
2031 ARTS_ASSERT(false, "out of bounds");
2032 }
2033}
2034#pragma GCC diagnostic pop
2035
2036// Needs to be implemented in this file!!!
2037std::ostream& operator<<(std::ostream& os, const PropagationMatrix& pm) {
2038 os << pm.Data() << "\n";
2039 return os;
2040}
2041
2042std::ostream& operator<<(std::ostream& os,
2043 const ArrayOfPropagationMatrix& apm) {
2044 for (auto& pm : apm) os << pm;
2045 return os;
2046}
2047
2048std::ostream& operator<<(std::ostream& os,
2049 const ArrayOfArrayOfPropagationMatrix& aapm) {
2050 for (auto& apm : aapm) os << apm;
2051 return os;
2052}
2053
2054// Needs to be implemented in this file!!!
2055std::ostream& operator<<(std::ostream& os, const StokesVector& sv) {
2056 os << sv.Data() << "\n";
2057 return os;
2058}
2059
2060std::ostream& operator<<(std::ostream& os, const ArrayOfStokesVector& asv) {
2061 for (auto& sv : asv) os << sv;
2062 return os;
2063}
2064
2065std::ostream& operator<<(std::ostream& os,
2066 const ArrayOfArrayOfStokesVector& aasv) {
2067 for (auto& asv : aasv) os << asv;
2068 return os;
2069}
arts_omp.h
Header file for helper functions for OpenMP.
Array
This can be used to make arrays out of anything.
Definition: array.h:108
Array::nelem
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:197
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:1050
ConstMatrixView::nrows
Index nrows() const noexcept
Definition: matpackI.h:1061
ConstMatrixView::ncols
Index ncols() const noexcept
Definition: matpackI.h:1062
MatrixView
The MatrixView class.
Definition: matpackI.h:1169
PropagationMatrix::DivideAtPosition
void DivideAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Divide at position.
Definition: propagationmatrix.h:575
PropagationMatrix::operator()
Numeric operator()(const Index iv=0, const Index is1=0, const Index is2=0, const Index iz=0, const Index ia=0) const
access operator.
Definition: propagationmatrix.cc:1953
PropagationMatrix::RightMultiplyAtPosition
void RightMultiplyAtPosition(MatrixView out, const ConstMatrixView &in, const Index iv=0, const Index iz=0, const Index ia=0) const
Multiply the matrix input from the right of this at position.
Definition: propagationmatrix.cc:1768
PropagationMatrix::NumberOfFrequencies
Index NumberOfFrequencies() const
The number of frequencies of the propagation matrix.
Definition: propagationmatrix.h:339
PropagationMatrix::MultiplyAndAdd
void MultiplyAndAdd(const Numeric x, const PropagationMatrix &y)
Multiply input by scalar and add to this.
Definition: propagationmatrix.cc:1631
PropagationMatrix::K24
VectorView K24(const Index iz=0, const Index ia=0)
Vector view to K(1, 3) elements.
Definition: propagationmatrix.h:955
PropagationMatrix::Data
Tensor4 & Data()
Get full view to data.
Definition: propagationmatrix.h:1043
PropagationMatrix::MatrixAtPosition
void MatrixAtPosition(MatrixView ret, const Index iv=0, const Index iz=0, const Index ia=0) const
Sets the dense matrix.
Definition: propagationmatrix.cc:1926
PropagationMatrix::K13
VectorView K13(const Index iz=0, const Index ia=0)
Vector view to K(0, 2) elements.
Definition: propagationmatrix.h:925
PropagationMatrix::NumberOfNeededVectors
Index NumberOfNeededVectors() const
The number of required vectors to fill this PropagationMatrix.
Definition: propagationmatrix.h:387
PropagationMatrix::Kjj
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
Definition: propagationmatrix.h:905
PropagationMatrix::StokesDimensions
Index StokesDimensions() const
The stokes dimension of the propagation matrix.
Definition: propagationmatrix.h:336
PropagationMatrix::K34
VectorView K34(const Index iz=0, const Index ia=0)
Vector view to K(2, 3) elements.
Definition: propagationmatrix.h:965
PropagationMatrix::SetAtPosition
void SetAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Set the At Position object.
Definition: propagationmatrix.h:498
PropagationMatrix::AddAtPosition
void AddAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Addition at position operator.
Definition: propagationmatrix.h:729
PropagationMatrix::FittingShape
bool FittingShape(const ConstMatrixView &x) const
Tests of the input matrix fits Propagation Matrix style.
Definition: propagationmatrix.cc:1640
PropagationMatrix::LeftMultiplyAtPosition
void LeftMultiplyAtPosition(MatrixView out, const ConstMatrixView &in, const Index iv=0, const Index iz=0, const Index ia=0) const
Multiply the matrix input from the left of this at position.
Definition: propagationmatrix.cc:1705
PropagationMatrix::MultiplyAtPosition
void MultiplyAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Multiply operator at position.
Definition: propagationmatrix.h:647
PropagationMatrix::K12
VectorView K12(const Index iz=0, const Index ia=0)
Vector view to K(0, 1) elements.
Definition: propagationmatrix.h:915
PropagationMatrix::K23
VectorView K23(const Index iz=0, const Index ia=0)
Vector view to K(1, 2) elements.
Definition: propagationmatrix.h:945
PropagationMatrix::K14
VectorView K14(const Index iz=0, const Index ia=0)
Vector view to K(0, 3) elements.
Definition: propagationmatrix.h:935
PropagationMatrix::GetTensor3
void GetTensor3(Tensor3View tensor3, const Index iz=0, const Index ia=0)
Get a Tensor3 object from this.
Definition: propagationmatrix.cc:1673
PropagationMatrix::RemoveAtPosition
void RemoveAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Subtraction at position.
Definition: propagationmatrix.h:805
PropagationMatrix::AddAverageAtPosition
void AddAverageAtPosition(const ConstMatrixView &mat1, const ConstMatrixView &mat2, const Index iv=0, const Index iz=0, const Index ia=0)
Add the average of the two input at position.
Definition: propagationmatrix.cc:1610
PropagationMatrix::MatrixInverseAtPosition
void MatrixInverseAtPosition(MatrixView ret, const Index iv=0, const Index iz=0, const Index ia=0) const
Return the matrix inverse at the position.
Definition: propagationmatrix.cc:1506
StokesVector
Stokes vector is as Propagation matrix but only has 4 possible values.
Definition: propagationmatrix.h:1181
Tensor3View
The Tensor3View class.
Definition: matpackIII.h:244
Tensor4View
The Tensor4View class.
Definition: matpackIV.h:290
ARTS_ASSERT
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
dx
#define dx
Definition: legacy_continua.cc:22121
lin_alg.h
Linear algebra functions.
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
joker
const Joker joker
a2
#define a2
Definition: matpack_complex.h:73
Complex
std::complex< Numeric > Complex
Definition: matpack_complex.h:35
b2
#define b2
Definition: matpack_complex.h:74
Absorption::nelem
Index nelem(const Lines &l)
Number of lines.
Definition: absorptionlines.h:951
Interpolation::l
constexpr Numeric l(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Index j, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
Definition: interpolation_lagrange.h:214
compute_transmission_matrix_and_derivative
void compute_transmission_matrix_and_derivative(Tensor3View T, Tensor4View dT_dx_upper_level, Tensor4View dT_dx_lower_level, const Numeric &r, const PropagationMatrix &upper_level, const PropagationMatrix &lower_level, const ArrayOfPropagationMatrix &dupper_level_dx, const ArrayOfPropagationMatrix &dlower_level_dx, const Numeric &dr_dTu, const Numeric &dr_dTl, const Index it, const Index iz, const Index ia)
Definition: propagationmatrix.cc:478
compute_transmission_matrix_from_averaged_matrix_at_frequency
void compute_transmission_matrix_from_averaged_matrix_at_frequency(MatrixView T, const Numeric &r, const PropagationMatrix &averaged_propagation_matrix, const Index iv, const Index iz, const Index ia)
Compute the matrix exponent as the transmission matrix of this propagation matrix.
Definition: propagationmatrix.cc:272
compute_transmission_matrix
void compute_transmission_matrix(Tensor3View T, const Numeric &r, const PropagationMatrix &upper_level, const PropagationMatrix &lower_level, const Index iz, const Index ia)
Compute the matrix exponent as the transmission matrix of this propagation matrix.
Definition: propagationmatrix.cc:33
operator<<
std::ostream & operator<<(std::ostream &os, const PropagationMatrix &pm)
output operator
Definition: propagationmatrix.cc:2037
propagationmatrix.h
Stuff related to the propagation matrix.
u
#define u
d
#define d
v
#define v
w
#define w
a
#define a
c
#define c
b
#define b
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:739