42 if (mstokes_dim == 1) {
43 for (
Index i = 0; i < mfreqs; i++)
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++) {
51 (upper_level.
Kjj(iz, ia)[i] +
52 lower_level.
Kjj(iz, ia)[i]),
54 (upper_level.
K12(iz, ia)[i] +
55 lower_level.
K12(iz, ia)[i]);
60 F(0, 1) = F(1, 0) = 0.;
61 F(0, 0) = F(1, 1) = exp_a;
68 F(0, 0) = F(1, 1) = C0 + C1 *
a;
69 F(0, 1) = F(1, 0) = C1 *
b;
73 }
else if (mstokes_dim == 3) {
74 for (
Index i = 0; i < mfreqs; i++) {
79 (upper_level.
Kjj(iz, ia)[i] + lower_level.
Kjj(iz, ia)[i]),
81 (upper_level.
K12(iz, ia)[i] + lower_level.
K12(iz, ia)[i]),
83 (upper_level.
K13(iz, ia)[i] + lower_level.
K13(iz, ia)[i]),
85 (upper_level.
K23(iz, ia)[i] + lower_level.
K23(iz, ia)[i]);
89 if (
b == 0. and
c == 0. and
u == 0.) {
91 F(0, 0) = F(1, 1) = F(2, 2) = exp_a;
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);
100 const Numeric C0 = (
a2 * (cosh_x - 1) -
a * x * sinh_x + x2) *
102 const Numeric C1 = (2 *
a * (1 - cosh_x) + x * sinh_x) *
105 (cosh_x - 1) * inv_x2;
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);
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);
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);
120 F(1, 2) = C1 *
u + C2 * (2 *
a *
u +
b *
c);
121 F(2, 1) = -C1 *
u - C2 * (2 *
a *
u -
b *
c);
125 }
else if (mstokes_dim == 4) {
127 for (
Index i = 0; i < mfreqs; i++) {
132 (upper_level.
Kjj(iz, ia)[i] + lower_level.
Kjj(iz, ia)[i]),
134 (upper_level.
K12(iz, ia)[i] + lower_level.
K12(iz, ia)[i]),
136 (upper_level.
K13(iz, ia)[i] + lower_level.
K13(iz, ia)[i]),
138 (upper_level.
K14(iz, ia)[i] + lower_level.
K14(iz, ia)[i]),
140 (upper_level.
K23(iz, ia)[i] + lower_level.
K23(iz, ia)[i]),
142 (upper_level.
K24(iz, ia)[i] + lower_level.
K24(iz, ia)[i]),
144 (upper_level.
K34(iz, ia)[i] + lower_level.
K34(iz, ia)[i]);
148 if (
b == 0. and
c == 0. and
d == 0. and
u == 0. and
v == 0. and
w == 0.) {
150 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
157 const Numeric Const2 =
b2 + c2 + d2 - u2 - v2 - w2;
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) +
169 Const1 =
sqrt(Const1);
175 const Numeric x = sqrt_BpA.real() * sqrt_05;
176 const Numeric y = sqrt_BmA.imag() * sqrt_05;
181 const Numeric cosh_x = cosh(x);
182 const Numeric sinh_x = sinh(x);
184 const Numeric inv_x2y2 = 1.0 / x2y2;
187 Numeric inv_y = 0.0, inv_x = 0.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) {
200 C2 = (cosh_x - 1.0) * inv_x2y2;
201 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
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;
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);
220 F(0, 1) = F(1, 0) = C1 *
b;
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));
229 F(0, 2) = F(2, 0) = C1 *
c;
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));
238 F(0, 3) = F(3, 0) = C1 *
d;
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));
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));
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));
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));
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];
291 T(0, 1) = T(1, 0) = 0.;
292 T(0, 0) = T(1, 1) = exp_a;
297 T(0, 0) = T(1, 1) = C0 + C1 *
a;
298 T(0, 1) = T(1, 0) = C1 *
b;
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];
310 if (
b == 0. and
c == 0. and
u == 0.) {
312 T(0, 0) = T(1, 1) = T(2, 2) = exp_a;
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);
319 const Numeric C0 = (
a2 * (cosh_x - 1) -
a * x * sinh_x + x2) *
321 const Numeric C1 = (2 *
a * (1 - cosh_x) + x * sinh_x) *
324 (cosh_x - 1) * inv_x2;
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);
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);
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);
339 T(1, 2) = C1 *
u + C2 * (2 *
a *
u +
b *
c);
340 T(2, 1) = -C1 *
u - C2 * (2 *
a *
u -
b *
c);
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];
355 if (
b == 0. and
c == 0. and
d == 0. and
u == 0. and
v == 0. and
w == 0.) {
357 T(0, 0) = T(1, 1) = T(2, 2) = T(3, 3) = exp_a;
362 const Numeric Const2 =
b2 + c2 + d2 - u2 - v2 - w2;
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);
371 Const1 += 8 * (
b *
d *
u *
w -
b *
c *
v *
w -
c *
d *
u *
v);
375 Const1 =
sqrt(Const1);
381 const Numeric x = sqrt_BpA.real() * sqrt_05;
382 const Numeric y = sqrt_BmA.imag() * sqrt_05;
387 const Numeric cosh_x = cosh(x);
388 const Numeric sinh_x = sinh(x);
390 const Numeric inv_x2y2 = 1.0 / x2y2;
393 Numeric inv_y = 0.0, inv_x = 0.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) {
406 C2 = (cosh_x - 1.0) * inv_x2y2;
407 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
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;
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);
426 T(0, 1) = T(1, 0) = C1 *
b;
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));
435 T(0, 2) = T(2, 0) = C1 *
c;
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));
444 T(0, 3) = T(3, 0) = C1 *
d;
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));
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));
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));
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));
496 if (mstokes_dim == 1) {
497 for (
Index i = 0; i < mfreqs; i++) {
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()) {
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]))
507 dT_dx_upper_level(j, i, 0, 0) = T(i, 0, 0) * da;
510 if (dlower_level_dx[j].NumberOfFrequencies()) {
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]))
516 dT_dx_lower_level(j, i, 0, 0) = T(i, 0, 0) * da;
520 }
else if (mstokes_dim == 2) {
521 for (
Index i = 0; i < mfreqs; i++) {
525 (upper_level.
Kjj(iz, ia)[i] +
526 lower_level.
Kjj(iz, ia)[i]),
528 (upper_level.
K12(iz, ia)[i] +
529 lower_level.
K12(iz, ia)[i]);
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()) {
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])
544 dT_dx_upper_level(j, i, 0, 0) = dT_dx_upper_level(j, i, 1, 1) *= da;
547 if (dlower_level_dx[j].NumberOfFrequencies()) {
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])
554 dT_dx_lower_level(j, i, 0, 0) = dT_dx_lower_level(j, i, 1, 1) *= da;
563 F(0, 0) = F(1, 1) = C0 + C1 *
a;
564 F(0, 1) = F(1, 0) = C1 *
b;
568 for (
Index j = 0; j < nppd; j++) {
569 if (not dlower_level_dx[j].NumberOfFrequencies())
continue;
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])
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])
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;
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;
591 for (
Index j = 0; j < nppd; j++) {
592 if (not dupper_level_dx[j].NumberOfFrequencies())
continue;
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])
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])
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;
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;
615 }
else if (mstokes_dim == 3) {
616 for (
Index i = 0; i < mfreqs; i++) {
621 (upper_level.
Kjj(iz, ia)[i] + lower_level.
Kjj(iz, ia)[i]),
623 (upper_level.
K12(iz, ia)[i] + lower_level.
K12(iz, ia)[i]),
625 (upper_level.
K13(iz, ia)[i] + lower_level.
K13(iz, ia)[i]),
627 (upper_level.
K23(iz, ia)[i] + lower_level.
K23(iz, ia)[i]);
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()) {
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])
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;
646 if (dlower_level_dx[j].NumberOfFrequencies()) {
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])
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;
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);
665 const Numeric C0 = (
a2 * (cosh_x - 1) -
a * x * sinh_x) * inv_x2 +
667 const Numeric C1 = (2 *
a * (1 - cosh_x) + x * sinh_x) *
670 (cosh_x - 1) * inv_x2;
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);
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);
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);
685 F(1, 2) = C1 *
u + C2 * (2 *
a *
u +
b *
c);
686 F(2, 1) = -C1 *
u - C2 * (2 *
a *
u -
b *
c);
690 for (
Index j = 0; j < nppd; j++) {
691 if (not dlower_level_dx[j].NumberOfFrequencies())
continue;
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])
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])
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])
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])
716 const Numeric da2 = 2 *
a * da, db2 = 2 *
b * db, dc2 = 2 *
c * dc,
719 const Numeric dx = (db2 + dc2 - du2) / x / 2;
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) *
727 -2 * C1 *
dx / x + (2 * da * (1 - cosh_x) - 2 *
a * sinh_x *
dx +
728 x * cosh_x *
dx + sinh_x *
dx) *
730 const Numeric dC2 = (sinh_x / x - 2 * C2) *
dx / x;
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);
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);
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);
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);
755 for (
int s1 = 0; s1 < 3; s1++)
756 for (
int s2 = 0; s2 < 3; s2++) dF(s1, s2) += F(s1, s2) * da;
759 for (
Index j = 0; j < nppd; j++) {
760 if (not dupper_level_dx[j].NumberOfFrequencies())
continue;
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])
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])
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])
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])
785 const Numeric da2 = 2 *
a * da, db2 = 2 *
b * db, dc2 = 2 *
c * dc,
788 const Numeric dx = (db2 + dc2 - du2) / x / 2;
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) *
796 -2 * C1 *
dx / x + (2 * da * (1 - cosh_x) - 2 *
a * sinh_x *
dx +
797 x * cosh_x *
dx + sinh_x *
dx) *
799 const Numeric dC2 = (sinh_x / x - 2 * C2) *
dx / x;
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);
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);
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);
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);
824 for (
int s1 = 0; s1 < 3; s1++)
825 for (
int s2 = 0; s2 < 3; s2++) dF(s1, s2) += F(s1, s2) * da;
828 }
else if (mstokes_dim == 4) {
830 for (
Index i = 0; i < mfreqs; i++) {
835 (upper_level.
Kjj(iz, ia)[i] + lower_level.
Kjj(iz, ia)[i]),
837 (upper_level.
K12(iz, ia)[i] + lower_level.
K12(iz, ia)[i]),
839 (upper_level.
K13(iz, ia)[i] + lower_level.
K13(iz, ia)[i]),
841 (upper_level.
K14(iz, ia)[i] + lower_level.
K14(iz, ia)[i]),
843 (upper_level.
K23(iz, ia)[i] + lower_level.
K23(iz, ia)[i]),
845 (upper_level.
K24(iz, ia)[i] + lower_level.
K24(iz, ia)[i]),
847 (upper_level.
K34(iz, ia)[i] + lower_level.
K34(iz, ia)[i]);
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()) {
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])
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) *=
868 if (dlower_level_dx[j].NumberOfFrequencies()) {
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])
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) *=
886 const Numeric Const2 =
b2 + c2 + d2 - u2 - v2 - w2;
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);
895 Const1 += 8 * (
b *
d *
u *
w -
b *
c *
v *
w -
c *
d *
u *
v);
899 Const1 =
sqrt(Const1);
905 const Numeric x = sqrt_BpA.real() * sqrt_05;
906 const Numeric y = sqrt_BmA.imag() * sqrt_05;
911 const Numeric cosh_x = cosh(x);
912 const Numeric sinh_x = sinh(x);
914 const Numeric inv_x2y2 = 1.0 / x2y2;
917 Numeric inv_y = 0.0, inv_x = 0.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) {
930 C2 = (cosh_x - 1.0) * inv_x2y2;
931 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
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;
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);
948 F(0, 1) = F(1, 0) = C1 *
b;
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));
956 F(0, 2) = F(2, 0) = C1 *
c;
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));
964 F(0, 3) = F(3, 0) = C1 *
d;
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));
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));
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));
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));
993 const Numeric inv_x2 = inv_x * inv_x;
994 const Numeric inv_y2 = inv_y * inv_y;
996 for (
Index j = 0; j < nppd; j++) {
997 if (not dupper_level_dx[j].NumberOfFrequencies())
continue;
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])
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])
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])
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])
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])
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])
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])
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;
1034 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1038 dConst1 = db2 * (
b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
1039 dConst1 +=
b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2);
1041 dConst1 += dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
1042 dConst1 += c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2);
1044 dConst1 += dd2 * (d2 * 0.5 + u2 - v2 - w2);
1045 dConst1 += d2 * (dd2 * 0.5 + du2 - dv2 - dw2);
1047 dConst1 += du2 * (u2 * 0.5 + v2 + w2);
1048 dConst1 += u2 * (du2 * 0.5 + dv2 + dw2);
1050 dConst1 += dv2 * (v2 * 0.5 + w2);
1051 dConst1 += v2 * (dv2 * 0.5 + dw2);
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;
1065 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
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;
1073 }
else if (y == 0.0) {
1075 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
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;
1084 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1086 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1087 const Numeric dy2 = 2 * y * dy;
1089 const Numeric dx2dy2 = dx2 + dy2;
1091 dC0 = -dx2dy2 * C0 * inv_x2y2 +
1092 (2 * cos_y *
dx * x + 2 * cosh_x * dy * y +
dx * sinh_x * y2 -
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) *
1103 -dx2dy2 * C2 * inv_x2y2 + (
dx * sinh_x + dy * sin_y) * inv_x2y2;
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) *
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);
1117 dF(0, 1) = dF(1, 0) = db * C1 +
b * dC1;
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));
1136 dF(0, 2) = dF(2, 0) = dC1 *
c + C1 * dc;
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));
1155 dF(0, 3) = dF(3, 0) = dC1 *
d + C1 * dd;
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));
1174 dF(1, 2) = dF(2, 1) =
1175 dC2 * (
b *
c -
v *
w) + C2 * (db *
c +
b * dc - dv *
w -
v * dw);
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));
1192 dF(1, 3) = dF(3, 1) =
1193 dC2 * (
b *
d +
u *
w) + C2 * (db *
d +
b * dd + du *
w +
u * dw);
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));
1210 dF(2, 3) = dF(3, 2) =
1211 dC2 * (
c *
d -
u *
v) + C2 * (dc *
d +
c * dd - du *
v -
u * dv);
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));
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;
1249 for (
Index j = 0; j < nppd; j++) {
1250 if (not dlower_level_dx[j].NumberOfFrequencies())
continue;
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])
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])
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])
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])
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])
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])
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])
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;
1287 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1291 dConst1 = db2 * (
b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
1292 dConst1 +=
b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2);
1294 dConst1 += dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
1295 dConst1 += c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2);
1297 dConst1 += dd2 * (d2 * 0.5 + u2 - v2 - w2);
1298 dConst1 += d2 * (dd2 * 0.5 + du2 - dv2 - dw2);
1300 dConst1 += du2 * (u2 * 0.5 + v2 + w2);
1301 dConst1 += u2 * (du2 * 0.5 + dv2 + dw2);
1303 dConst1 += dv2 * (v2 * 0.5 + w2);
1304 dConst1 += v2 * (dv2 * 0.5 + dw2);
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;
1318 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
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;
1326 }
else if (y == 0.0) {
1328 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
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;
1337 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1339 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1340 const Numeric dy2 = 2 * y * dy;
1342 const Numeric dx2dy2 = dx2 + dy2;
1344 dC0 = -dx2dy2 * C0 * inv_x2y2 +
1345 (2 * cos_y *
dx * x + 2 * cosh_x * dy * y +
dx * sinh_x * y2 -
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) *
1356 -dx2dy2 * C2 * inv_x2y2 + (
dx * sinh_x + dy * sin_y) * inv_x2y2;
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) *
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);
1370 dF(0, 1) = dF(1, 0) = db * C1 +
b * dC1;
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));
1389 dF(0, 2) = dF(2, 0) = dC1 *
c + C1 * dc;
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));
1408 dF(0, 3) = dF(3, 0) = dC1 *
d + C1 * dd;
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));
1427 dF(1, 2) = dF(2, 1) =
1428 dC2 * (
b *
c -
v *
w) + C2 * (db *
c +
b * dc - dv *
w -
v * dw);
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));
1445 dF(1, 3) = dF(3, 1) =
1446 dC2 * (
b *
d +
u *
w) + C2 * (db *
d +
b * dd + du *
w +
u * dw);
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));
1463 dF(2, 3) = dF(3, 2) =
1464 dC2 * (
c *
d -
u *
v) + C2 * (dc *
d +
c * dd - du *
v -
u * dv);
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));
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;
1509 const Index ia)
const {
1512 ret(0, 0) = 1.0 /
Kjj(iz, ia)[iv];
1522 ret(1, 1) = ret(0, 0) =
Kjj(iz, ia)[iv] * div;
1523 ret(1, 0) = ret(0, 1) = -
K12(iz, ia)[iv] * div;
1528 u =
K23(iz, ia)[iv], u2 =
u *
u;
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;
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;
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;
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;
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 -
1560 ret(0, 0) =
a * (
a2 + u2 + v2 + w2) * div;
1574 ret(1, 1) =
a * (
a2 - c2 - d2 + w2) * div;
1588 ret(2, 2) =
a * (
a2 -
b2 - d2 + v2) * div;
1602 ret(3, 3) =
a * (
a2 -
b2 - c2 + u2) * div;
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;
1621 mdata(ia, iz, iv, 2) += (mat1(2, 0) + mat2(2, 0)) * 0.5;
1623 (mat1(1, 2) + mat2(1, 2)) * 0.5;
1625 mdata(ia, iz, iv, 1) += (mat1(1, 0) + mat2(1, 0)) * 0.5;
1627 mdata(ia, iz, iv, 0) += (mat1(0, 0) + mat2(0, 0)) * 0.5;
1649 if (x(0, 0) == x(1, 1))
return true;
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))
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))
1665 ARTS_ASSERT(
false,
"Stokes dimension does not agree with accepted values");
1683 tensor3(
joker, 3, 2) *= -1;
1684 tensor3(
joker, 3, 1) *= -1;
1691 tensor3(
joker, 2, 1) *= -1;
1700 ARTS_ASSERT(
false,
"Stokes dimension does not agree with accepted values");
1709 const Index ia)
const {
1712 out(0, 0) =
Kjj(iz, ia)[iv] * in(0, 0);
1715 const Numeric a =
Kjj(iz, ia)[iv],
b =
K12(iz, ia)[iv], m11 = in(0, 0),
1716 m12 = in(0, 1), m21 = in(1, 0), m22 = in(1, 1);
1717 out(0, 0) =
a * m11 +
b * m21;
1718 out(0, 1) =
a * m12 +
b * m22;
1719 out(1, 0) =
a * m21 +
b * m11;
1720 out(1, 1) =
a * m22 +
b * m12;
1724 c =
K13(iz, ia)[iv],
u =
K23(iz, ia)[iv], m11 = in(0, 0),
1725 m12 = in(0, 1), m13 = in(0, 2), m21 = in(1, 0),
1726 m22 = in(1, 1), m23 = in(1, 2), m31 = in(2, 0),
1727 m32 = in(2, 1), m33 = in(2, 2);
1728 out(0, 0) =
a * m11 +
b * m21 +
c * m31;
1729 out(0, 1) =
a * m12 +
b * m22 +
c * m32;
1730 out(0, 2) =
a * m13 +
b * m23 +
c * m33;
1731 out(1, 0) =
a * m21 +
b * m11 + m31 *
u;
1732 out(1, 1) =
a * m22 +
b * m12 + m32 *
u;
1733 out(1, 2) =
a * m23 +
b * m13 + m33 *
u;
1734 out(2, 0) =
a * m31 +
c * m11 - m21 *
u;
1735 out(2, 1) =
a * m32 +
c * m12 - m22 *
u;
1736 out(2, 2) =
a * m33 +
c * m13 - m23 *
u;
1740 c =
K13(iz, ia)[iv],
u =
K23(iz, ia)[iv],
1741 d =
K14(iz, ia)[iv],
v =
K24(iz, ia)[iv],
1742 w =
K34(iz, ia)[iv], m11 = in(0, 0), m12 = in(0, 1),
1743 m13 = in(0, 2), m14 = in(0, 3), m21 = in(1, 0),
1744 m22 = in(1, 1), m23 = in(1, 2), m24 = in(1, 3),
1745 m31 = in(2, 0), m32 = in(2, 1), m33 = in(2, 2),
1746 m34 = in(2, 3), m41 = in(3, 0), m42 = in(3, 1),
1747 m43 = in(3, 2), m44 = in(3, 3);
1748 out(0, 0) =
a * m11 +
b * m21 +
c * m31 +
d * m41;
1749 out(0, 1) =
a * m12 +
b * m22 +
c * m32 +
d * m42;
1750 out(0, 2) =
a * m13 +
b * m23 +
c * m33 +
d * m43;
1751 out(0, 3) =
a * m14 +
b * m24 +
c * m34 +
d * m44;
1752 out(1, 0) =
a * m21 +
b * m11 + m31 *
u + m41 *
v;
1753 out(1, 1) =
a * m22 +
b * m12 + m32 *
u + m42 *
v;
1754 out(1, 2) =
a * m23 +
b * m13 + m33 *
u + m43 *
v;
1755 out(1, 3) =
a * m24 +
b * m14 + m34 *
u + m44 *
v;
1756 out(2, 0) =
a * m31 +
c * m11 - m21 *
u + m41 *
w;
1757 out(2, 1) =
a * m32 +
c * m12 - m22 *
u + m42 *
w;
1758 out(2, 2) =
a * m33 +
c * m13 - m23 *
u + m43 *
w;
1759 out(2, 3) =
a * m34 +
c * m14 - m24 *
u + m44 *
w;
1760 out(3, 0) =
a * m41 +
d * m11 - m21 *
v - m31 *
w;
1761 out(3, 1) =
a * m42 +
d * m12 - m22 *
v - m32 *
w;
1762 out(3, 2) =
a * m43 +
d * m13 - m23 *
v - m33 *
w;
1763 out(3, 3) =
a * m44 +
d * m14 - m24 *
v - m34 *
w;
1772 const Index ia)
const {
1775 out(0, 0) = in(0, 0) *
Kjj(iz, ia)[iv];
1778 const Numeric a =
Kjj(iz, ia)[iv],
b =
K12(iz, ia)[iv], m11 = in(0, 0),
1779 m12 = in(0, 1), m21 = in(1, 0), m22 = in(1, 1);
1780 out(0, 0) =
a * m11 +
b * m12;
1781 out(0, 1) =
a * m12 +
b * m11;
1782 out(1, 0) =
a * m21 +
b * m22;
1783 out(1, 1) =
a * m22 +
b * m21;
1787 c =
K13(iz, ia)[iv],
u =
K23(iz, ia)[iv], m11 = in(0, 0),
1788 m12 = in(0, 1), m13 = in(0, 2), m21 = in(1, 0),
1789 m22 = in(1, 1), m23 = in(1, 2), m31 = in(2, 0),
1790 m32 = in(2, 1), m33 = in(2, 2);
1791 out(0, 0) =
a * m11 +
b * m12 +
c * m13;
1792 out(0, 1) =
a * m12 +
b * m11 - m13 *
u;
1793 out(0, 2) =
a * m13 +
c * m11 + m12 *
u;
1794 out(1, 0) =
a * m21 +
b * m22 +
c * m23;
1795 out(1, 1) =
a * m22 +
b * m21 - m23 *
u;
1796 out(1, 2) =
a * m23 +
c * m21 + m22 *
u;
1797 out(2, 0) =
a * m31 +
b * m32 +
c * m33;
1798 out(2, 1) =
a * m32 +
b * m31 - m33 *
u;
1799 out(2, 2) =
a * m33 +
c * m31 + m32 *
u;
1803 c =
K13(iz, ia)[iv],
u =
K23(iz, ia)[iv],
1804 d =
K14(iz, ia)[iv],
v =
K24(iz, ia)[iv],
1805 w =
K34(iz, ia)[iv], m11 = in(0, 0), m12 = in(0, 1),
1806 m13 = in(0, 2), m14 = in(0, 3), m21 = in(1, 0),
1807 m22 = in(1, 1), m23 = in(1, 2), m24 = in(1, 3),
1808 m31 = in(2, 0), m32 = in(2, 1), m33 = in(2, 2),
1809 m34 = in(2, 3), m41 = in(3, 0), m42 = in(3, 1),
1810 m43 = in(3, 2), m44 = in(3, 3);
1811 out(0, 0) =
a * m11 +
b * m12 +
c * m13 +
d * m14;
1812 out(0, 1) =
a * m12 +
b * m11 - m13 *
u - m14 *
v;
1813 out(0, 2) =
a * m13 +
c * m11 + m12 *
u - m14 *
w;
1814 out(0, 3) =
a * m14 +
d * m11 + m12 *
v + m13 *
w;
1815 out(1, 0) =
a * m21 +
b * m22 +
c * m23 +
d * m24;
1816 out(1, 1) =
a * m22 +
b * m21 - m23 *
u - m24 *
v;
1817 out(1, 2) =
a * m23 +
c * m21 + m22 *
u - m24 *
w;
1818 out(1, 3) =
a * m24 +
d * m21 + m22 *
v + m23 *
w;
1819 out(2, 0) =
a * m31 +
b * m32 +
c * m33 +
d * m34;
1820 out(2, 1) =
a * m32 +
b * m31 - m33 *
u - m34 *
v;
1821 out(2, 2) =
a * m33 +
c * m31 + m32 *
u - m34 *
w;
1822 out(2, 3) =
a * m34 +
d * m31 + m32 *
v + m33 *
w;
1823 out(3, 0) =
a * m41 +
b * m42 +
c * m43 +
d * m44;
1824 out(3, 1) =
a * m42 +
b * m41 - m43 *
u - m44 *
v;
1825 out(3, 2) =
a * m43 +
c * m41 + m42 *
u - m44 *
w;
1826 out(3, 3) =
a * m44 +
d * m41 + m42 *
v + m43 *
w;
1837 mdata(ia, iz, iv, 5) -= x(1, 3);
1838 mdata(ia, iz, iv, 6) -= x(2, 3);
1839 mdata(ia, iz, iv, 3) -= x(0, 3);
1841 mdata(ia, iz, iv, 2) -= x(0, 2);
1844 mdata(ia, iz, iv, 1) -= x(0, 1);
1846 mdata(ia, iz, iv, 0) -= x(0, 0);
1856 mdata(ia, iz, iv, 5) += x(1, 3);
1857 mdata(ia, iz, iv, 6) += x(2, 3);
1858 mdata(ia, iz, iv, 3) += x(0, 3);
1860 mdata(ia, iz, iv, 2) += x(0, 2);
1863 mdata(ia, iz, iv, 1) += x(0, 1);
1865 mdata(ia, iz, iv, 0) += x(0, 0);
1875 mdata(ia, iz, iv, 5) *= x(1, 3);
1876 mdata(ia, iz, iv, 6) *= x(2, 3);
1877 mdata(ia, iz, iv, 3) *= x(0, 3);
1879 mdata(ia, iz, iv, 2) *= x(0, 2);
1882 mdata(ia, iz, iv, 1) *= x(0, 1);
1884 mdata(ia, iz, iv, 0) *= x(0, 0);
1894 mdata(ia, iz, iv, 5) /= x(1, 3);
1895 mdata(ia, iz, iv, 6) /= x(2, 3);
1896 mdata(ia, iz, iv, 3) /= x(0, 3);
1898 mdata(ia, iz, iv, 2) /= x(0, 2);
1901 mdata(ia, iz, iv, 1) /= x(0, 1);
1903 mdata(ia, iz, iv, 0) /= x(0, 0);
1913 mdata(ia, iz, iv, 5) = x(1, 3);
1914 mdata(ia, iz, iv, 6) = x(2, 3);
1915 mdata(ia, iz, iv, 3) = x(0, 3);
1917 mdata(ia, iz, iv, 2) = x(0, 2);
1920 mdata(ia, iz, iv, 1) = x(0, 1);
1922 mdata(ia, iz, iv, 0) = x(0, 0);
1929 const Index ia)
const {
1932 ret(3, 3) =
mdata(ia, iz, iv, 0);
1933 ret(3, 1) = -
mdata(ia, iz, iv, 5);
1934 ret(1, 3) =
mdata(ia, iz, iv, 5);
1935 ret(3, 2) = -
mdata(ia, iz, iv, 6);
1936 ret(2, 3) =
mdata(ia, iz, iv, 6);
1937 ret(0, 3) = ret(3, 0) =
mdata(ia, iz, iv, 3);
1939 ret(2, 2) =
mdata(ia, iz, iv, 0);
1940 ret(2, 1) = -
mdata(ia, iz, iv, 3);
1941 ret(1, 2) =
mdata(ia, iz, iv, 3);
1942 ret(2, 0) = ret(0, 2) =
mdata(ia, iz, iv, 2);
1944 ret(1, 1) =
mdata(ia, iz, iv, 0);
1945 ret(1, 0) = ret(0, 1) =
mdata(ia, iz, iv, 1);
1947 ret(0, 0) =
mdata(ia, iz, iv, 0);
1951#pragma GCC diagnostic push
1952#pragma GCC diagnostic ignored "-Wreturn-type"
1957 const Index ia)
const {
1962 return mdata(ia, iz, iv, 0);
1965 return mdata(ia, iz, iv, 1);
1968 return mdata(ia, iz, iv, 2);
1971 return mdata(ia, iz, iv, 3);
1980 return mdata(ia, iz, iv, 1);
1983 return mdata(ia, iz, iv, 0);
1989 return mdata(ia, iz, iv, 5);
1997 return mdata(ia, iz, iv, 2);
2003 return mdata(ia, iz, iv, 0);
2006 return mdata(ia, iz, iv, 6);
2015 return mdata(ia, iz, iv, 3);
2018 return -
mdata(ia, iz, iv, 5);
2021 return -
mdata(ia, iz, iv, 6);
2024 return mdata(ia, iz, iv, 0);
2034#pragma GCC diagnostic pop
2038 os << pm.
Data() <<
"\n";
2044 for (
auto& pm : apm) os << pm;
2050 for (
auto& apm : aapm) os << apm;
2056 os << sv.
Data() <<
"\n";
2061 for (
auto& sv : asv) os << sv;
2067 for (
auto& asv : aasv) os << asv;
Header file for helper functions for OpenMP.
This can be used to make arrays out of anything.
Index nelem() const ARTS_NOEXCEPT
Number of elements.
A constant view of a Matrix.
Index nrows() const noexcept
Index ncols() const noexcept
void DivideAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Divide at position.
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.
void RightMultiplyAtPosition(MatrixView out, const ConstMatrixView &in, const Index iv=0, const Index iz=0, const Index ia=0) const
Multiply the matrix input from the right of this at position.
Index NumberOfFrequencies() const
The number of frequencies of the propagation matrix.
void MultiplyAndAdd(const Numeric x, const PropagationMatrix &y)
Multiply input by scalar and add to this.
VectorView K24(const Index iz=0, const Index ia=0)
Vector view to K(1, 3) elements.
Tensor4 & Data()
Get full view to data.
void MatrixAtPosition(MatrixView ret, const Index iv=0, const Index iz=0, const Index ia=0) const
Sets the dense matrix.
VectorView K13(const Index iz=0, const Index ia=0)
Vector view to K(0, 2) elements.
Index NumberOfNeededVectors() const
The number of required vectors to fill this PropagationMatrix.
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
Index StokesDimensions() const
The stokes dimension of the propagation matrix.
VectorView K34(const Index iz=0, const Index ia=0)
Vector view to K(2, 3) elements.
void SetAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Set the At Position object.
void AddAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Addition at position operator.
bool FittingShape(const ConstMatrixView &x) const
Tests of the input matrix fits Propagation Matrix style.
void LeftMultiplyAtPosition(MatrixView out, const ConstMatrixView &in, const Index iv=0, const Index iz=0, const Index ia=0) const
Multiply the matrix input from the left of this at position.
void MultiplyAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Multiply operator at position.
VectorView K12(const Index iz=0, const Index ia=0)
Vector view to K(0, 1) elements.
VectorView K23(const Index iz=0, const Index ia=0)
Vector view to K(1, 2) elements.
VectorView K14(const Index iz=0, const Index ia=0)
Vector view to K(0, 3) elements.
void GetTensor3(Tensor3View tensor3, const Index iz=0, const Index ia=0)
Get a Tensor3 object from this.
void RemoveAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Subtraction at position.
void AddAverageAtPosition(const ConstMatrixView &mat1, const ConstMatrixView &mat2, const Index iv=0, const Index iz=0, const Index ia=0)
Add the average of the two input at position.
void MatrixInverseAtPosition(MatrixView ret, const Index iv=0, const Index iz=0, const Index ia=0) const
Return the matrix inverse at the position.
Stokes vector is as Propagation matrix but only has 4 possible values.
#define ARTS_ASSERT(condition,...)
Linear algebra functions.
NUMERIC Numeric
The type to use for all floating point numbers.
INDEX Index
The type to use for all integer numbers and indices.
std::complex< Numeric > Complex
Index nelem(const Lines &l)
Number of lines.
constexpr Numeric l(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Index j, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
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)
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.
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.
std::ostream & operator<<(std::ostream &os, const PropagationMatrix &pm)
output operator
Stuff related to the propagation matrix.
Numeric sqrt(const Rational r)
Square root.