ARTS  2.4.0(git:4fb77825)
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 
33 void 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 
272 void 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 
478 void 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 
1506 void 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  throw std::runtime_error("Strange stokes dimensions");
1606  break;
1607  }
1608 }
1609 
1610 void PropagationMatrix::AddAverageAtPosition(ConstMatrixView mat1,
1611  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 
1631 void 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 
1640 bool PropagationMatrix::FittingShape(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  throw std::runtime_error(
1666  "Stokes dimension does not agree with accepted values");
1667  break;
1668  }
1669  }
1670  }
1671  return false;
1672 }
1673 
1674 void 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  throw std::runtime_error(
1702  "Stokes dimension does not agree with accepted values");
1703  break;
1704  }
1705 }
1706 
1707 void PropagationMatrix::LeftMultiplyAtPosition(MatrixView out,
1708  ConstMatrixView in,
1709  const Index iv,
1710  const Index iz,
1711  const Index ia) const {
1712  switch (mstokes_dim) {
1713  case 1:
1714  out(0, 0) = Kjj(iz, ia)[iv] * in(0, 0);
1715  break;
1716  case 2: {
1717  const Numeric a = Kjj(iz, ia)[iv], b = K12(iz, ia)[iv], m11 = in(0, 0),
1718  m12 = in(0, 1), m21 = in(1, 0), m22 = in(1, 1);
1719  out(0, 0) = a * m11 + b * m21;
1720  out(0, 1) = a * m12 + b * m22;
1721  out(1, 0) = a * m21 + b * m11;
1722  out(1, 1) = a * m22 + b * m12;
1723  } break;
1724  case 3: {
1725  const Numeric a = Kjj(iz, ia)[iv], b = K12(iz, ia)[iv],
1726  c = K13(iz, ia)[iv], u = K23(iz, ia)[iv], m11 = in(0, 0),
1727  m12 = in(0, 1), m13 = in(0, 2), m21 = in(1, 0),
1728  m22 = in(1, 1), m23 = in(1, 2), m31 = in(2, 0),
1729  m32 = in(2, 1), m33 = in(2, 2);
1730  out(0, 0) = a * m11 + b * m21 + c * m31;
1731  out(0, 1) = a * m12 + b * m22 + c * m32;
1732  out(0, 2) = a * m13 + b * m23 + c * m33;
1733  out(1, 0) = a * m21 + b * m11 + m31 * u;
1734  out(1, 1) = a * m22 + b * m12 + m32 * u;
1735  out(1, 2) = a * m23 + b * m13 + m33 * u;
1736  out(2, 0) = a * m31 + c * m11 - m21 * u;
1737  out(2, 1) = a * m32 + c * m12 - m22 * u;
1738  out(2, 2) = a * m33 + c * m13 - m23 * u;
1739  } break;
1740  case 4: {
1741  const Numeric a = Kjj(iz, ia)[iv], b = K12(iz, ia)[iv],
1742  c = K13(iz, ia)[iv], u = K23(iz, ia)[iv],
1743  d = K14(iz, ia)[iv], v = K24(iz, ia)[iv],
1744  w = K34(iz, ia)[iv], m11 = in(0, 0), m12 = in(0, 1),
1745  m13 = in(0, 2), m14 = in(0, 3), m21 = in(1, 0),
1746  m22 = in(1, 1), m23 = in(1, 2), m24 = in(1, 3),
1747  m31 = in(2, 0), m32 = in(2, 1), m33 = in(2, 2),
1748  m34 = in(2, 3), m41 = in(3, 0), m42 = in(3, 1),
1749  m43 = in(3, 2), m44 = in(3, 3);
1750  out(0, 0) = a * m11 + b * m21 + c * m31 + d * m41;
1751  out(0, 1) = a * m12 + b * m22 + c * m32 + d * m42;
1752  out(0, 2) = a * m13 + b * m23 + c * m33 + d * m43;
1753  out(0, 3) = a * m14 + b * m24 + c * m34 + d * m44;
1754  out(1, 0) = a * m21 + b * m11 + m31 * u + m41 * v;
1755  out(1, 1) = a * m22 + b * m12 + m32 * u + m42 * v;
1756  out(1, 2) = a * m23 + b * m13 + m33 * u + m43 * v;
1757  out(1, 3) = a * m24 + b * m14 + m34 * u + m44 * v;
1758  out(2, 0) = a * m31 + c * m11 - m21 * u + m41 * w;
1759  out(2, 1) = a * m32 + c * m12 - m22 * u + m42 * w;
1760  out(2, 2) = a * m33 + c * m13 - m23 * u + m43 * w;
1761  out(2, 3) = a * m34 + c * m14 - m24 * u + m44 * w;
1762  out(3, 0) = a * m41 + d * m11 - m21 * v - m31 * w;
1763  out(3, 1) = a * m42 + d * m12 - m22 * v - m32 * w;
1764  out(3, 2) = a * m43 + d * m13 - m23 * v - m33 * w;
1765  out(3, 3) = a * m44 + d * m14 - m24 * v - m34 * w;
1766  }
1767  }
1768 }
1769 
1770 void PropagationMatrix::RightMultiplyAtPosition(MatrixView out,
1771  ConstMatrixView in,
1772  const Index iv,
1773  const Index iz,
1774  const Index ia) const {
1775  switch (mstokes_dim) {
1776  case 1:
1777  out(0, 0) = in(0, 0) * Kjj(iz, ia)[iv];
1778  break;
1779  case 2: {
1780  const Numeric a = Kjj(iz, ia)[iv], b = K12(iz, ia)[iv], m11 = in(0, 0),
1781  m12 = in(0, 1), m21 = in(1, 0), m22 = in(1, 1);
1782  out(0, 0) = a * m11 + b * m12;
1783  out(0, 1) = a * m12 + b * m11;
1784  out(1, 0) = a * m21 + b * m22;
1785  out(1, 1) = a * m22 + b * m21;
1786  } break;
1787  case 3: {
1788  const Numeric a = Kjj(iz, ia)[iv], b = K12(iz, ia)[iv],
1789  c = K13(iz, ia)[iv], u = K23(iz, ia)[iv], m11 = in(0, 0),
1790  m12 = in(0, 1), m13 = in(0, 2), m21 = in(1, 0),
1791  m22 = in(1, 1), m23 = in(1, 2), m31 = in(2, 0),
1792  m32 = in(2, 1), m33 = in(2, 2);
1793  out(0, 0) = a * m11 + b * m12 + c * m13;
1794  out(0, 1) = a * m12 + b * m11 - m13 * u;
1795  out(0, 2) = a * m13 + c * m11 + m12 * u;
1796  out(1, 0) = a * m21 + b * m22 + c * m23;
1797  out(1, 1) = a * m22 + b * m21 - m23 * u;
1798  out(1, 2) = a * m23 + c * m21 + m22 * u;
1799  out(2, 0) = a * m31 + b * m32 + c * m33;
1800  out(2, 1) = a * m32 + b * m31 - m33 * u;
1801  out(2, 2) = a * m33 + c * m31 + m32 * u;
1802  } break;
1803  case 4: {
1804  const Numeric a = Kjj(iz, ia)[iv], b = K12(iz, ia)[iv],
1805  c = K13(iz, ia)[iv], u = K23(iz, ia)[iv],
1806  d = K14(iz, ia)[iv], v = K24(iz, ia)[iv],
1807  w = K34(iz, ia)[iv], m11 = in(0, 0), m12 = in(0, 1),
1808  m13 = in(0, 2), m14 = in(0, 3), m21 = in(1, 0),
1809  m22 = in(1, 1), m23 = in(1, 2), m24 = in(1, 3),
1810  m31 = in(2, 0), m32 = in(2, 1), m33 = in(2, 2),
1811  m34 = in(2, 3), m41 = in(3, 0), m42 = in(3, 1),
1812  m43 = in(3, 2), m44 = in(3, 3);
1813  out(0, 0) = a * m11 + b * m12 + c * m13 + d * m14;
1814  out(0, 1) = a * m12 + b * m11 - m13 * u - m14 * v;
1815  out(0, 2) = a * m13 + c * m11 + m12 * u - m14 * w;
1816  out(0, 3) = a * m14 + d * m11 + m12 * v + m13 * w;
1817  out(1, 0) = a * m21 + b * m22 + c * m23 + d * m24;
1818  out(1, 1) = a * m22 + b * m21 - m23 * u - m24 * v;
1819  out(1, 2) = a * m23 + c * m21 + m22 * u - m24 * w;
1820  out(1, 3) = a * m24 + d * m21 + m22 * v + m23 * w;
1821  out(2, 0) = a * m31 + b * m32 + c * m33 + d * m34;
1822  out(2, 1) = a * m32 + b * m31 - m33 * u - m34 * v;
1823  out(2, 2) = a * m33 + c * m31 + m32 * u - m34 * w;
1824  out(2, 3) = a * m34 + d * m31 + m32 * v + m33 * w;
1825  out(3, 0) = a * m41 + b * m42 + c * m43 + d * m44;
1826  out(3, 1) = a * m42 + b * m41 - m43 * u - m44 * v;
1827  out(3, 2) = a * m43 + c * m41 + m42 * u - m44 * w;
1828  out(3, 3) = a * m44 + d * m41 + m42 * v + m43 * w;
1829  }
1830  }
1831 }
1832 
1833 void PropagationMatrix::RemoveAtPosition(ConstMatrixView x,
1834  const Index iv,
1835  const Index iz,
1836  const Index ia) {
1837  switch (mstokes_dim) {
1838  case 4:
1839  mdata(ia, iz, iv, 5) -= x(1, 3);
1840  mdata(ia, iz, iv, 6) -= x(2, 3);
1841  mdata(ia, iz, iv, 3) -= x(0, 3); /* FALLTHROUGH */
1842  case 3:
1843  mdata(ia, iz, iv, 2) -= x(0, 2);
1844  mdata(ia, iz, iv, mstokes_dim) -= x(1, 2); /* FALLTHROUGH */
1845  case 2:
1846  mdata(ia, iz, iv, 1) -= x(0, 1); /* FALLTHROUGH */
1847  case 1:
1848  mdata(ia, iz, iv, 0) -= x(0, 0); /* FALLTHROUGH */
1849  }
1850 }
1851 
1852 void PropagationMatrix::AddAtPosition(ConstMatrixView x,
1853  const Index iv,
1854  const Index iz,
1855  const Index ia) {
1856  switch (mstokes_dim) {
1857  case 4:
1858  mdata(ia, iz, iv, 5) += x(1, 3);
1859  mdata(ia, iz, iv, 6) += x(2, 3);
1860  mdata(ia, iz, iv, 3) += x(0, 3); /* FALLTHROUGH */
1861  case 3:
1862  mdata(ia, iz, iv, 2) += x(0, 2);
1863  mdata(ia, iz, iv, mstokes_dim) += x(1, 2); /* FALLTHROUGH */
1864  case 2:
1865  mdata(ia, iz, iv, 1) += x(0, 1); /* FALLTHROUGH */
1866  case 1:
1867  mdata(ia, iz, iv, 0) += x(0, 0); /* FALLTHROUGH */
1868  }
1869 }
1870 
1871 void PropagationMatrix::MultiplyAtPosition(ConstMatrixView x,
1872  const Index iv,
1873  const Index iz,
1874  const Index ia) {
1875  switch (mstokes_dim) {
1876  case 4:
1877  mdata(ia, iz, iv, 5) *= x(1, 3);
1878  mdata(ia, iz, iv, 6) *= x(2, 3);
1879  mdata(ia, iz, iv, 3) *= x(0, 3); /* FALLTHROUGH */
1880  case 3:
1881  mdata(ia, iz, iv, 2) *= x(0, 2);
1882  mdata(ia, iz, iv, mstokes_dim) *= x(1, 2); /* FALLTHROUGH */
1883  case 2:
1884  mdata(ia, iz, iv, 1) *= x(0, 1); /* FALLTHROUGH */
1885  case 1:
1886  mdata(ia, iz, iv, 0) *= x(0, 0); /* FALLTHROUGH */
1887  }
1888 }
1889 
1890 void PropagationMatrix::DivideAtPosition(ConstMatrixView x,
1891  const Index iv,
1892  const Index iz,
1893  const Index ia) {
1894  switch (mstokes_dim) {
1895  case 4:
1896  mdata(ia, iz, iv, 5) /= x(1, 3);
1897  mdata(ia, iz, iv, 6) /= x(2, 3);
1898  mdata(ia, iz, iv, 3) /= x(0, 3); /* FALLTHROUGH */
1899  case 3:
1900  mdata(ia, iz, iv, 2) /= x(0, 2);
1901  mdata(ia, iz, iv, mstokes_dim) /= x(1, 2); /* FALLTHROUGH */
1902  case 2:
1903  mdata(ia, iz, iv, 1) /= x(0, 1); /* FALLTHROUGH */
1904  case 1:
1905  mdata(ia, iz, iv, 0) /= x(0, 0); /* FALLTHROUGH */
1906  }
1907 }
1908 
1909 void PropagationMatrix::SetAtPosition(ConstMatrixView x,
1910  const Index iv,
1911  const Index iz,
1912  const Index ia) {
1913  switch (mstokes_dim) {
1914  case 4:
1915  mdata(ia, iz, iv, 5) = x(1, 3);
1916  mdata(ia, iz, iv, 6) = x(2, 3);
1917  mdata(ia, iz, iv, 3) = x(0, 3); /* FALLTHROUGH */
1918  case 3:
1919  mdata(ia, iz, iv, 2) = x(0, 2);
1920  mdata(ia, iz, iv, mstokes_dim) = x(1, 2); /* FALLTHROUGH */
1921  case 2:
1922  mdata(ia, iz, iv, 1) = x(0, 1); /* FALLTHROUGH */
1923  case 1:
1924  mdata(ia, iz, iv, 0) = x(0, 0); /* FALLTHROUGH */
1925  }
1926 }
1927 
1928 void PropagationMatrix::MatrixAtPosition(MatrixView ret,
1929  const Index iv,
1930  const Index iz,
1931  const Index ia) const {
1932  switch (mstokes_dim) {
1933  case 4:
1934  ret(3, 3) = mdata(ia, iz, iv, 0);
1935  ret(3, 1) = -mdata(ia, iz, iv, 5);
1936  ret(1, 3) = mdata(ia, iz, iv, 5);
1937  ret(3, 2) = -mdata(ia, iz, iv, 6);
1938  ret(2, 3) = mdata(ia, iz, iv, 6);
1939  ret(0, 3) = ret(3, 0) = mdata(ia, iz, iv, 3); /* FALLTHROUGH */
1940  case 3:
1941  ret(2, 2) = mdata(ia, iz, iv, 0);
1942  ret(2, 1) = -mdata(ia, iz, iv, 3);
1943  ret(1, 2) = mdata(ia, iz, iv, 3);
1944  ret(2, 0) = ret(0, 2) = mdata(ia, iz, iv, 2); /* FALLTHROUGH */
1945  case 2:
1946  ret(1, 1) = mdata(ia, iz, iv, 0);
1947  ret(1, 0) = ret(0, 1) = mdata(ia, iz, iv, 1); /* FALLTHROUGH */
1948  case 1:
1949  ret(0, 0) = mdata(ia, iz, iv, 0); /* FALLTHROUGH */
1950  }
1951 }
1952 
1953 Numeric 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  throw std::runtime_error("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  throw std::runtime_error("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  throw std::runtime_error("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  throw std::runtime_error("out of bounds");
2028  }
2029  break;
2030  default:
2031  throw std::runtime_error("out of bounds");
2032  }
2033 }
2034 
2035 // Needs to be implemented in this file!!!
2036 std::ostream& operator<<(std::ostream& os, const PropagationMatrix& pm) {
2037  os << pm.Data() << "\n";
2038  return os;
2039 }
2040 
2041 std::ostream& operator<<(std::ostream& os,
2042  const ArrayOfPropagationMatrix& apm) {
2043  for (auto& pm : apm) os << pm;
2044  return os;
2045 }
2046 
2047 std::ostream& operator<<(std::ostream& os,
2048  const ArrayOfArrayOfPropagationMatrix& aapm) {
2049  for (auto& apm : aapm) os << apm;
2050  return os;
2051 }
2052 
2053 // Needs to be implemented in this file!!!
2054 std::ostream& operator<<(std::ostream& os, const StokesVector& sv) {
2055  os << sv.Data() << "\n";
2056  return os;
2057 }
2058 
2059 std::ostream& operator<<(std::ostream& os, const ArrayOfStokesVector& asv) {
2060  for (auto& sv : asv) os << sv;
2061  return os;
2062 }
2063 
2064 std::ostream& operator<<(std::ostream& os,
2065  const ArrayOfArrayOfStokesVector& aasv) {
2066  for (auto& asv : aasv) os << asv;
2067  return os;
2068 }
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:1928
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
StokesVector
Stokes vector is as Propagation matrix but only has 4 possible values.
Definition: propagationmatrix.h:1075
PropagationMatrix::NumberOfNeededVectors
Index NumberOfNeededVectors() const
The number of required vectors to fill this PropagationMatrix.
Definition: propagationmatrix.h:257
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
w
Complex w(Complex z) noexcept
The Faddeeva function.
Definition: linefunctions.cc:42
joker
const Joker joker
PropagationMatrix::maa
Index maa
Definition: propagationmatrix.h:973
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
PropagationMatrix::FittingShape
bool FittingShape(ConstMatrixView x) const
Tests of the input matrix fits Propagation Matrix style.
Definition: propagationmatrix.cc:1640
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
PropagationMatrix
Definition: propagationmatrix.h:87
PropagationMatrix::mfreqs
Index mfreqs
Definition: propagationmatrix.h:972
b2
#define b2
Definition: complex.h:59
PropagationMatrix::K13
VectorView K13(const Index iz=0, const Index ia=0)
Vector view to K(0, 2) elements.
Definition: propagationmatrix.h:819
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:699
Complex
std::complex< Numeric > Complex
Definition: complex.h:33
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
PropagationMatrix::RightMultiplyAtPosition
void RightMultiplyAtPosition(MatrixView out, 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:1770
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
PropagationMatrix::AddAverageAtPosition
void AddAverageAtPosition(ConstMatrixView mat1, 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::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:541
PropagationMatrix::MultiplyAndAdd
void MultiplyAndAdd(const Numeric x, const PropagationMatrix &y)
Multiply input by scalar and add to this.
Definition: propagationmatrix.cc:1631
Tensor3View
The Tensor3View class.
Definition: matpackIII.h:239
Array
This can be used to make arrays out of anything.
Definition: array.h:108
PropagationMatrix::K24
VectorView K24(const Index iz=0, const Index ia=0)
Vector view to K(1, 3) elements.
Definition: propagationmatrix.h:849
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::K23
VectorView K23(const Index iz=0, const Index ia=0)
Vector view to K(1, 2) elements.
Definition: propagationmatrix.h:839
Absorption::nelem
Index nelem(const Lines &l)
Number of lines.
Definition: absorptionlines.h:1820
dw
constexpr Complex dw(Complex z, Complex w) noexcept
The Faddeeva function partial derivative.
Definition: linefunctions.cc:45
PropagationMatrix::K12
VectorView K12(const Index iz=0, const Index ia=0)
Vector view to K(0, 1) elements.
Definition: propagationmatrix.h:809
PropagationMatrix::K34
VectorView K34(const Index iz=0, const Index ia=0)
Vector view to K(2, 3) elements.
Definition: propagationmatrix.h:859
PropagationMatrix::mdata
Tensor4 mdata
Definition: propagationmatrix.h:974
PropagationMatrix::NumberOfFrequencies
Index NumberOfFrequencies() const
The number of frequencies of the propagation matrix.
Definition: propagationmatrix.h:204
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
lin_alg.h
Linear algebra functions.
PropagationMatrix::Data
Tensor4 & Data()
Get full view to data.
Definition: propagationmatrix.h:937
propagationmatrix.h
Stuff related to the propagation matrix.
PropagationMatrix::LeftMultiplyAtPosition
void LeftMultiplyAtPosition(MatrixView out, 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:1707
PropagationMatrix::K14
VectorView K14(const Index iz=0, const Index ia=0)
Vector view to K(0, 3) elements.
Definition: propagationmatrix.h:829
PropagationMatrix::mstokes_dim
Index mstokes_dim
Definition: propagationmatrix.h:972
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:982
PropagationMatrix::StokesDimensions
Index StokesDimensions() const
The stokes dimension of the propagation matrix.
Definition: propagationmatrix.h:201
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:469
PropagationMatrix::Kjj
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
Definition: propagationmatrix.h:799
Tensor4View
The Tensor4View class.
Definition: matpackIV.h:284
dx
#define dx
Definition: legacy_continua.cc:22180
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
PropagationMatrix::GetTensor3
void GetTensor3(Tensor3View tensor3, const Index iz=0, const Index ia=0)
Get a Tensor3 object from this.
Definition: propagationmatrix.cc:1674
operator<<
std::ostream & operator<<(std::ostream &os, const PropagationMatrix &pm)
output operator
Definition: propagationmatrix.cc:2036
a2
#define a2
Definition: complex.h:58
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
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:623
PropagationMatrix::mza
Index mza
Definition: propagationmatrix.h:973
arts_omp.h
Header file for helper functions for OpenMP.
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
x2
#define x2
Definition: lineshapemodel.h:415
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:392