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