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