36 T4(stokes_dim == 4 ? nf : 0, Eigen::Matrix4d::Identity()),
37 T3(stokes_dim == 3 ? nf : 0, Eigen::Matrix3d::Identity()),
38 T2(stokes_dim == 2 ? nf : 0, Eigen::Matrix2d::Identity()),
39 T1(stokes_dim == 1 ? nf : 0, Eigen::
Matrix<double, 1, 1>::Identity()) {
60TransmissionMatrix::operator
Tensor3()
const {
61 Tensor3 T(Frequencies(), stokes_dim, stokes_dim);
62 for (
size_t i = 0; i < T4.size(); i++)
63 for (
size_t j = 0; j < 4; j++)
64 for (
size_t k = 0; k < 4; k++) T(i, j, k) = T4[i](j, k);
65 for (
size_t i = 0; i < T3.size(); i++)
66 for (
size_t j = 0; j < 3; j++)
67 for (
size_t k = 0; k < 3; k++) T(i, j, k) = T3[i](j, k);
68 for (
size_t i = 0; i < T2.size(); i++)
69 for (
size_t j = 0; j < 2; j++)
70 for (
size_t k = 0; k < 2; k++) T(i, j, k) = T2[i](j, k);
71 for (
size_t i = 0; i < T1.size(); i++) T(i, 0, 0) = T1[i](0, 0);
75#pragma GCC diagnostic push
76#pragma GCC diagnostic ignored "-Wreturn-type"
89#pragma GCC diagnostic pop
92 std::fill(
T4.begin(),
T4.end(), Eigen::Matrix4d::Identity());
93 std::fill(
T3.begin(),
T3.end(), Eigen::Matrix3d::Identity());
94 std::fill(
T2.begin(),
T2.end(), Eigen::Matrix2d::Identity());
95 std::fill(
T1.begin(),
T1.end(), Eigen::Matrix<double, 1, 1>::Identity());
99 std::fill(
T4.begin(),
T4.end(), Eigen::Matrix4d::Zero());
100 std::fill(
T3.begin(),
T3.end(), Eigen::Matrix3d::Zero());
101 std::fill(
T2.begin(),
T2.end(), Eigen::Matrix2d::Zero());
102 std::fill(
T1.begin(),
T1.end(), Eigen::Matrix<double, 1, 1>::Zero());
107 for (
size_t i = 0; i <
T4.size(); i++)
T4[i].noalias() = A.
T4[i] * B.
T4[i];
108 for (
size_t i = 0; i <
T3.size(); i++)
T3[i].noalias() = A.
T3[i] * B.
T3[i];
109 for (
size_t i = 0; i <
T2.size(); i++)
T2[i].noalias() = A.
T2[i] * B.
T2[i];
110 for (
size_t i = 0; i <
T1.size(); i++)
T1[i].noalias() = A.
T1[i] * B.
T1[i];
115 for (
size_t i = 0; i <
T4.size(); i++)
T4[i] = A.
T4[i] * B.
T4[i];
116 for (
size_t i = 0; i <
T3.size(); i++)
T3[i] = A.
T3[i] * B.
T3[i];
117 for (
size_t i = 0; i <
T2.size(); i++)
T2[i] = A.
T2[i] * B.
T2[i];
118 for (
size_t i = 0; i <
T1.size(); i++)
T1[i] = A.
T1[i] * B.
T1[i];
123 const Index k)
const {
166 operator()(0, i, j) = mat(i,j);
171 for (
size_t i = 0; i <
T4.size(); i++)
172 T4[i].noalias() = lstm.
scale * lstm.
bas.Mat4(i);
173 for (
size_t i = 0; i <
T3.size(); i++)
174 T3[i].noalias() = lstm.
scale * lstm.
bas.Mat3(i);
175 for (
size_t i = 0; i <
T2.size(); i++)
176 T2[i].noalias() = lstm.
scale * lstm.
bas.Mat2(i);
177 for (
size_t i = 0; i <
T1.size(); i++)
178 T1[i].noalias() = lstm.
scale * lstm.
bas.Mat1(i);
184 T4.begin(),
T4.end(),
T4.begin(), [scale](
auto& T) { return T * scale; });
186 T3.begin(),
T3.end(),
T3.begin(), [scale](
auto& T) { return T * scale; });
188 T2.begin(),
T2.end(),
T2.begin(), [scale](
auto& T) { return T * scale; });
190 T1.begin(),
T1.end(),
T1.begin(), [scale](
auto& T) { return T * scale; });
205 : stokes_dim(stokes),
206 R4(stokes_dim == 4 ? nf : 0, Eigen::Vector4d::Zero()),
207 R3(stokes_dim == 3 ? nf : 0, Eigen::Vector3d::Zero()),
208 R2(stokes_dim == 2 ? nf : 0, Eigen::Vector2d::Zero()),
209 R1(stokes_dim == 1 ? nf : 0, Eigen::
Matrix<double, 1, 1>::Zero()) {
239 for (
size_t i = 0; i <
R4.size(); i++)
R4[i] = T.
Mat4(i) *
R4[i];
240 for (
size_t i = 0; i <
R3.size(); i++)
R3[i] = T.
Mat3(i) *
R3[i];
241 for (
size_t i = 0; i <
R2.size(); i++)
R2[i] = T.
Mat2(i) *
R2[i];
242 for (
size_t i = 0; i <
R1.size(); i++)
R1[i] = T.
Mat1(i) *
R1[i];
248 R4[i].noalias() = Eigen::Vector4d::Zero();
251 R3[i].noalias() = Eigen::Vector3d::Zero();
254 R2[i].noalias() = Eigen::Vector2d::Zero();
281 for (
size_t i = 0; i <
R4.size(); i++)
282 R4[i].noalias() -= 0.5 * (O1.
R4[i] + O2.
R4[i]);
283 for (
size_t i = 0; i <
R3.size(); i++)
284 R3[i].noalias() -= 0.5 * (O1.
R3[i] + O2.
R3[i]);
285 for (
size_t i = 0; i <
R2.size(); i++)
286 R2[i].noalias() -= 0.5 * (O1.
R2[i] + O2.
R2[i]);
287 for (
size_t i = 0; i <
R1.size(); i++)
288 R1[i].noalias() -= 0.5 * (O1.
R1[i] + O2.
R1[i]);
293 for (
size_t i = 0; i <
R4.size(); i++)
294 R4[i].noalias() += 0.5 * (O1.
R4[i] + O2.
R4[i]);
295 for (
size_t i = 0; i <
R3.size(); i++)
296 R3[i].noalias() += 0.5 * (O1.
R3[i] + O2.
R3[i]);
297 for (
size_t i = 0; i <
R2.size(); i++)
298 R2[i].noalias() += 0.5 * (O1.
R2[i] + O2.
R2[i]);
299 for (
size_t i = 0; i <
R1.size(); i++)
300 R1[i].noalias() += 0.5 * (O1.
R1[i] + O2.
R1[i]);
309 for (
size_t i = 0; i <
R4.size(); i++) {
314 prop_matrix<4>(Kfar(i,
joker)),
315 prop_matrix<4>(Kclose(i,
joker)),
318 for (
size_t i = 0; i <
R3.size(); i++) {
323 prop_matrix<3>(Kfar(i,
joker)),
324 prop_matrix<3>(Kclose(i,
joker)),
327 for (
size_t i = 0; i <
R2.size(); i++) {
332 prop_matrix<2>(Kfar(i,
joker)),
333 prop_matrix<2>(Kclose(i,
joker)),
336 for (
size_t i = 0; i <
R1.size(); i++) {
341 prop_matrix<1>(Kfar(i,
joker)),
342 prop_matrix<1>(Kclose(i,
joker)),
352 for (
size_t i = 0; i <
R4.size(); i++)
353 R4[i].noalias() += PiT.
Mat4(i) * (dT.
Mat4(i) * ImJ.
R4[i] + dJ.
R4[i] -
355 for (
size_t i = 0; i <
R3.size(); i++)
356 R3[i].noalias() += PiT.
Mat3(i) * (dT.
Mat3(i) * ImJ.
R3[i] + dJ.
R3[i] -
358 for (
size_t i = 0; i <
R2.size(); i++)
359 R2[i].noalias() += PiT.
Mat2(i) * (dT.
Mat2(i) * ImJ.
R2[i] + dJ.
R2[i] -
361 for (
size_t i = 0; i <
R1.size(); i++)
362 R1[i].noalias() += PiT.
Mat1(i) * (dT.
Mat1(i) * ImJ.
R1[i] + dJ.
R1[i] -
374 for (
size_t i = 0; i <
R4.size(); i++)
378 i, dT, far.
R4[i], close.
R4[i],
d.R4[i], isfar));
379 for (
size_t i = 0; i <
R3.size(); i++)
383 i, dT, far.
R3[i], close.
R3[i],
d.R3[i], isfar));
384 for (
size_t i = 0; i <
R2.size(); i++)
388 i, dT, far.
R2[i], close.
R2[i],
d.R2[i], isfar));
389 for (
size_t i = 0; i <
R1.size(); i++)
393 i, dT, far.
R1[i], close.
R1[i],
d.R1[i], isfar));
399 for (
size_t i = 0; i <
R4.size(); i++)
401 for (
size_t i = 0; i <
R3.size(); i++)
403 for (
size_t i = 0; i <
R2.size(); i++)
405 for (
size_t i = 0; i <
R1.size(); i++)
411 for (
size_t i = 0; i <
R4.size(); i++)
R4[i].noalias() += A.
Mat4(i) * x.
R4[i];
412 for (
size_t i = 0; i <
R3.size(); i++)
R3[i].noalias() += A.
Mat3(i) * x.
R3[i];
413 for (
size_t i = 0; i <
R2.size(); i++)
R2[i].noalias() += A.
Mat2(i) * x.
R2[i];
414 for (
size_t i = 0; i <
R1.size(); i++)
R1[i].noalias() += A.
Mat1(i) * x.
R1[i];
421 for (
size_t i = 0; i <
R4.size(); i++)
423 for (
size_t i = 0; i <
R3.size(); i++)
425 for (
size_t i = 0; i <
R2.size(); i++)
427 for (
size_t i = 0; i <
R1.size(); i++)
436 for (
size_t i = 0; i <
R4.size(); i++)
438 for (
size_t i = 0; i <
R3.size(); i++)
440 for (
size_t i = 0; i <
R2.size(); i++)
442 for (
size_t i = 0; i <
R1.size(); i++)
451 for (
size_t i = 0; i <
R4.size(); i++)
453 for (
size_t i = 0; i <
R3.size(); i++)
455 for (
size_t i = 0; i <
R2.size(); i++)
457 for (
size_t i = 0; i <
R1.size(); i++)
463 for (
size_t i = 0; i <
R4.size(); i++) {
469 for (
size_t i = 0; i <
R3.size(); i++) {
474 for (
size_t i = 0; i <
R2.size(); i++) {
478 for (
size_t i = 0; i <
R1.size(); i++) {
485 for (
size_t i = 0; i <
R4.size(); i++)
486 R4[i].noalias() += rv.
R4[i];
487 for (
size_t i = 0; i <
R3.size(); i++)
488 R3[i].noalias() += rv.
R3[i];
489 for (
size_t i = 0; i <
R2.size(); i++)
490 R2[i].noalias() += rv.
R2[i];
491 for (
size_t i = 0; i <
R1.size(); i++)
492 R1[i].noalias() += rv.
R1[i];
511RadiationVector::operator
Matrix()
const {
512 Matrix M(Frequencies(), stokes_dim);
513 for (
size_t i = 0; i < R4.size(); i++)
514 for (
size_t j = 0; j < 4; j++)
M(i, j) = R4[i](j);
515 for (
size_t i = 0; i < R3.size(); i++)
516 for (
size_t j = 0; j < 3; j++)
M(i, j) = R3[i](j);
517 for (
size_t i = 0; i < R2.size(); i++)
518 for (
size_t j = 0; j < 2; j++)
M(i, j) = R2[i](j);
519 for (
size_t i = 0; i < R1.size(); i++)
M(i, 0) = R1[i](0);
535 Eigen::Vector4d(
a.Kjj()[i],
a.K12()[i],
a.K13()[i],
a.K14()[i]) *
537 Eigen::Vector4d(S.
Kjj()[i], S.
K12()[i], S.
K13()[i], S.
K14()[i]);
540 Eigen::Vector4d(
a.Kjj()[i],
a.K12()[i],
a.K13()[i],
a.K14()[i]) *
546 Eigen::Vector3d(
a.Kjj()[i],
a.K12()[i],
a.K13()[i]) * B[i] +
547 Eigen::Vector3d(S.
Kjj()[i], S.
K12()[i], S.
K13()[i]);
550 Eigen::Vector3d(
a.Kjj()[i],
a.K12()[i],
a.K13()[i]) * B[i];
554 R2[i].noalias() = Eigen::Vector2d(
a.Kjj()[i],
a.K12()[i]) * B[i] +
555 Eigen::Vector2d(S.
Kjj()[i], S.
K12()[i]);
557 R2[i].noalias() = Eigen::Vector2d(
a.Kjj()[i],
a.K12()[i]) * B[i];
561 R1[i][0] =
a.Kjj()[i] * B[i] + S.
Kjj()[i];
563 R1[i][0] =
a.Kjj()[i] * B[i];
591 static_assert(
N > 0 and
N < 5,
"Bad stokes dimensions");
593 if constexpr (
N == 1) {
595 return Eigen::Matrix<Numeric, 1, 1>(dS.
Kjj()[i] + da.
Kjj()[i] * B[i] +
596 a.Kjj()[i] * dB_dT[i]);
597 return Eigen::Matrix<Numeric, 1, 1>(dS.
Kjj()[i] + da.
Kjj()[i] * B[i]);
600 if constexpr (
N == 2) {
602 return Eigen::Vector2d(
603 dS.
Kjj()[i] + da.
Kjj()[i] * B[i] +
a.Kjj()[i] * dB_dT[i],
604 dS.
K12()[i] + da.
K12()[i] * B[i] +
a.K12()[i] * dB_dT[i]);
605 return Eigen::Vector2d(dS.
Kjj()[i] + da.
Kjj()[i] * B[i],
606 dS.
K12()[i] + da.
K12()[i] * B[i]);
609 if constexpr (
N == 3) {
611 return Eigen::Vector3d(
612 dS.
Kjj()[i] + da.
Kjj()[i] * B[i] +
a.Kjj()[i] * dB_dT[i],
613 dS.
K12()[i] + da.
K12()[i] * B[i] +
a.K12()[i] * dB_dT[i],
614 dS.
K13()[i] + da.
K13()[i] * B[i] +
a.K13()[i] * dB_dT[i]);
615 return Eigen::Vector3d(dS.
Kjj()[i] + da.
Kjj()[i] * B[i],
616 dS.
K12()[i] + da.
K12()[i] * B[i],
617 dS.
K13()[i] + da.
K13()[i] * B[i]);
620 if constexpr (
N == 4) {
622 return Eigen::Vector4d(
623 dS.
Kjj()[i] + da.
Kjj()[i] * B[i] +
a.Kjj()[i] * dB_dT[i],
624 dS.
K12()[i] + da.
K12()[i] * B[i] +
a.K12()[i] * dB_dT[i],
625 dS.
K13()[i] + da.
K13()[i] * B[i] +
a.K13()[i] * dB_dT[i],
626 dS.
K14()[i] + da.
K14()[i] * B[i] +
a.K14()[i] * dB_dT[i]);
627 return Eigen::Vector4d(dS.
Kjj()[i] + da.
Kjj()[i] * B[i],
628 dS.
K12()[i] + da.
K12()[i] * B[i],
629 dS.
K13()[i] + da.
K13()[i] * B[i],
630 dS.
K14()[i] + da.
K14()[i] * B[i]);
642 static_assert(
N > 0 and
N < 5,
"Bad stokes dimensions");
644 if constexpr (
N == 1) {
646 return Eigen::Matrix<Numeric, 1, 1>(da.
Kjj()[i] * B[i] +
647 a.Kjj()[i] * dB_dT[i]);
648 return Eigen::Matrix<Numeric, 1, 1>(da.
Kjj()[i] * B[i]);
651 if constexpr (
N == 2) {
653 return Eigen::Vector2d(da.
Kjj()[i] * B[i] +
a.Kjj()[i] * dB_dT[i],
654 da.
K12()[i] * B[i] +
a.K12()[i] * dB_dT[i]);
655 return Eigen::Vector2d(da.
Kjj()[i] * B[i], da.
K12()[i] * B[i]);
658 if constexpr (
N == 3) {
660 return Eigen::Vector3d(da.
Kjj()[i] * B[i] +
a.Kjj()[i] * dB_dT[i],
661 da.
K12()[i] * B[i] +
a.K12()[i] * dB_dT[i],
662 da.
K13()[i] * B[i] +
a.K13()[i] * dB_dT[i]);
663 return Eigen::Vector3d(
664 da.
Kjj()[i] * B[i], da.
K12()[i] * B[i], da.
K13()[i] * B[i]);
667 if constexpr (
N == 4) {
669 return Eigen::Vector4d(da.
Kjj()[i] * B[i] +
a.Kjj()[i] * dB_dT[i],
670 da.
K12()[i] * B[i] +
a.K12()[i] * dB_dT[i],
671 da.
K13()[i] * B[i] +
a.K13()[i] * dB_dT[i],
672 da.
K14()[i] * B[i] +
a.K14()[i] * dB_dT[i]);
673 return Eigen::Vector4d(da.
Kjj()[i] * B[i],
685 const Index ia = 0) noexcept {
688 std::exp(-0.5 * r * (K1.
Kjj(iz, ia)[i] + K2.
Kjj(iz, ia)[i]));
696 const Index ia = 0) noexcept {
698 const Numeric a = -0.5 * r * (K1.
Kjj(iz, ia)[i] + K2.
Kjj(iz, ia)[i]),
699 b = -0.5 * r * (K1.
K12(iz, ia)[i] + K2.
K12(iz, ia)[i]);
701 const Numeric cb = std::cosh(
b), sb = std::sinh(
b);
702 T.
Mat2(i).noalias() =
703 (Eigen::Matrix2d() << cb, sb, sb, cb).finished() * exp_a;
712 const Index ia = 0) noexcept {
714 const Numeric a = -0.5 * r * (K1.
Kjj(iz, ia)[i] + K2.
Kjj(iz, ia)[i]),
715 b = -0.5 * r * (K1.
K12(iz, ia)[i] + K2.
K12(iz, ia)[i]),
716 c = -0.5 * r * (K1.
K13(iz, ia)[i] + K2.
K13(iz, ia)[i]),
717 u = -0.5 * r * (K1.
K23(iz, ia)[i] + K2.
K23(iz, ia)[i]);
720 if (
b == 0. and
c == 0. and
u == 0.) {
721 T.
Mat3(i).noalias() = Eigen::Matrix3d::Identity() * exp_a;
726 const bool real = Const > 0;
727 const bool imag = Const < 0;
731 std::sqrt(
imag ? -Const : Const);
733 (
real ? 1 : -1) * x * x;
739 real ? std::sinh(x) : std::sin(-x);
741 real ? std::cosh(x) : std::cos(+x);
750 either ?
a2 * (cx - 1.0) -
a * x * sx + x2 : 1.0 + 0.5 *
a2 -
a;
751 const Numeric C1 = either ? 2.0 *
a * (1.0 - cx) + x * sx : 1.0 -
a;
752 const Numeric C2 = either ? cx - 1.0 : 0.5;
754 T.
Mat3(i).noalias() =
756 (Eigen::Matrix3d() << C0 + C1 *
a + C2 * (
a2 +
b2 + c2),
757 C1 *
b + C2 * (2 *
a *
b -
c *
u),
758 C1 *
c + C2 * (2 *
a *
c +
b *
u),
759 C1 *
b + C2 * (2 *
a *
b +
c *
u),
760 C0 + C1 *
a + C2 * (
a2 +
b2 - u2),
761 C1 *
u + C2 * (2 *
a *
u +
b *
c),
762 C1 *
c + C2 * (2 *
a *
c -
b *
u),
763 -C1 *
u - C2 * (2 *
a *
u -
b *
c),
764 C0 + C1 *
a + C2 * (
a2 + c2 - u2))
775 const Index ia = 0) noexcept {
778 const Numeric a = -0.5 * r * (K1.
Kjj(iz, ia)[i] + K2.
Kjj(iz, ia)[i]),
779 b = -0.5 * r * (K1.
K12(iz, ia)[i] + K2.
K12(iz, ia)[i]),
780 c = -0.5 * r * (K1.
K13(iz, ia)[i] + K2.
K13(iz, ia)[i]),
781 d = -0.5 * r * (K1.
K14(iz, ia)[i] + K2.
K14(iz, ia)[i]),
782 u = -0.5 * r * (K1.
K23(iz, ia)[i] + K2.
K23(iz, ia)[i]),
783 v = -0.5 * r * (K1.
K24(iz, ia)[i] + K2.
K24(iz, ia)[i]),
784 w = -0.5 * r * (K1.
K34(iz, ia)[i] + K2.
K34(iz, ia)[i]);
787 if (
b == 0. and
c == 0. and
d == 0. and
u == 0. and
v == 0. and
w == 0.)
788 T.
Mat4(i).noalias() = Eigen::Matrix4d::Identity() * exp_a;
794 w2 * w2 + 2 * (
b2 * (
b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
795 c2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
796 d2 * (d2 * 0.5 + u2 - v2 - w2) +
797 u2 * (u2 * 0.5 + v2 + w2) + v2 * (v2 * 0.5 + w2) +
800 const Numeric Const2 =
b2 + c2 + d2 - u2 - v2 - w2;
802 const Complex x = std::sqrt(Const2 + Const1) * sqrt_05;
803 const Complex y = std::sqrt(Const2 - Const1) * sqrt_05 *
Complex(0, 1);
806 const Complex cy = std::cos(y);
807 const Complex sy = std::sin(y);
808 const Complex cx = std::cosh(x);
809 const Complex sx = std::sinh(x);
813 const bool both_zero = y_zero and x_zero;
814 const bool either_zero = y_zero or x_zero;
824 const Complex ix = x_zero ? 0.0 : 1.0 / x;
825 const Complex iy = y_zero ? 0.0 : 1.0 / y;
833 either_zero ? 1.0 : ((cy * x2 + cx * y2) * inv_x2y2).real();
835 either_zero ? 1.0 : ((sy * x2 * iy + sx * y2 * ix) * inv_x2y2).real();
836 const Numeric C2 = both_zero ? 0.5 : ((cx - cy) * inv_x2y2).real();
837 const Numeric C3 = both_zero ? 1.0 / 6.0
838 : ((x_zero ? 1.0 - sy * iy
839 : y_zero ? sx * ix - 1.0
840 : sx * ix - sy * iy) *
843 T.
Mat4(i).noalias() =
844 exp_a * (Eigen::Matrix4d() << C0 + C2 * (
b2 + c2 + d2),
845 C1 *
b + C2 * (-
c *
u -
d *
v) +
846 C3 * (
b * (
b2 + c2 + d2) -
u * (
b *
u -
d *
w) -
847 v * (
b *
v +
c *
w)),
848 C1 *
c + C2 * (
b *
u -
d *
w) +
849 C3 * (
c * (
b2 + c2 + d2) -
u * (
c *
u +
d *
v) -
850 w * (
b *
v +
c *
w)),
851 C1 *
d + C2 * (
b *
v +
c *
w) +
852 C3 * (
d * (
b2 + c2 + d2) -
v * (
c *
u +
d *
v) +
853 w * (
b *
u -
d *
w)),
855 C1 *
b + C2 * (
c *
u +
d *
v) +
856 C3 * (-
b * (-
b2 + u2 + v2) +
c * (
b *
c -
v *
w) +
857 d * (
b *
d +
u *
w)),
858 C0 + C2 * (
b2 - u2 - v2),
859 C2 * (
b *
c -
v *
w) + C1 *
u +
860 C3 * (
c * (
c *
u +
d *
v) -
u * (-
b2 + u2 + v2) -
861 w * (
b *
d +
u *
w)),
862 C2 * (
b *
d +
u *
w) + C1 *
v +
863 C3 * (
d * (
c *
u +
d *
v) -
v * (-
b2 + u2 + v2) +
864 w * (
b *
c -
v *
w)),
866 C1 *
c + C2 * (-
b *
u +
d *
w) +
867 C3 * (
b * (
b *
c -
v *
w) -
c * (-c2 + u2 + w2) +
868 d * (
c *
d -
u *
v)),
869 C2 * (
b *
c -
v *
w) - C1 *
u +
870 C3 * (-
b * (
b *
u -
d *
w) +
u * (-c2 + u2 + w2) -
871 v * (
c *
d -
u *
v)),
872 C0 + C2 * (c2 - u2 - w2),
873 C2 * (
c *
d -
u *
v) + C1 *
w +
874 C3 * (-
d * (
b *
u -
d *
w) +
v * (
b *
c -
v *
w) -
875 w * (-c2 + u2 + w2)),
877 C1 *
d + C2 * (-
b *
v -
c *
w) +
878 C3 * (
b * (
b *
d +
u *
w) +
c * (
c *
d -
u *
v) -
879 d * (-d2 + v2 + w2)),
880 C2 * (
b *
d +
u *
w) - C1 *
v +
881 C3 * (-
b * (
b *
v +
c *
w) -
u * (
c *
d -
u *
v) +
882 v * (-d2 + v2 + w2)),
883 C2 * (
c *
d -
u *
v) - C1 *
w +
884 C3 * (-
c * (
b *
v +
c *
w) +
u * (
b *
d +
u *
w) +
885 w * (-d2 + v2 + w2)),
886 C0 + C2 * (d2 - v2 - w2))
904 const Index ia)
noexcept {
905 for (
Index i = 0; i < K1.NumberOfFrequencies(); i++) {
907 std::exp(-0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]));
908 for (
Index j = 0; j < dT1.nelem(); j++) {
909 if (dK1[j].NumberOfFrequencies())
910 dT1[j].Mat1(i)(0, 0) =
913 (r * dK1[j].Kjj(iz, ia)[i] +
914 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
916 if (dK2[j].NumberOfFrequencies())
917 dT2[j].Mat1(i)(0, 0) =
920 (r * dK2[j].Kjj(iz, ia)[i] +
921 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
939 const Index ia)
noexcept {
940 for (
Index i = 0; i < K1.NumberOfFrequencies(); i++) {
941 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
942 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]);
944 const Numeric cb = std::cosh(
b), sb = std::sinh(
b);
945 T.Mat2(i).noalias() =
946 (Eigen::Matrix2d() << cb, sb, sb, cb).finished() * exp_a;
947 for (
Index j = 0; j < dT1.nelem(); j++) {
948 if (dK1[j].NumberOfFrequencies()) {
949 const Numeric da = -0.5 * (r * dK1[j].Kjj(iz, ia)[i] +
950 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] +
953 db = -0.5 * (r * dK1[j].K12(iz, ia)[i] +
954 ((j == it) ? dr_dT1 * (K1.K12(iz, ia)[i] +
957 dT1[j].Mat2(i).noalias() =
959 (Eigen::Matrix2d() << sb, cb, cb, sb).finished() * exp_a * db;
961 if (dK2[j].NumberOfFrequencies()) {
962 const Numeric da = -0.5 * (r * dK2[j].Kjj(iz, ia)[i] +
963 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] +
966 db = -0.5 * (r * dK2[j].K12(iz, ia)[i] +
967 ((j == it) ? dr_dT2 * (K1.K12(iz, ia)[i] +
970 dT2[j].Mat2(i).noalias() =
972 (Eigen::Matrix2d() << sb, cb, cb, sb).finished() * exp_a * db;
990 const Index ia)
noexcept {
991 for (
Index i = 0; i < K1.NumberOfFrequencies(); i++) {
992 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
993 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]),
994 c = -0.5 * r * (K1.K13(iz, ia)[i] + K2.K13(iz, ia)[i]),
995 u = -0.5 * r * (K1.K23(iz, ia)[i] + K2.K23(iz, ia)[i]);
998 if (
b == 0. and
c == 0. and
u == 0.) {
999 T.Mat3(i).noalias() = Eigen::Matrix3d::Identity() * exp_a;
1000 for (
Index j = 0; j < dT1.nelem(); j++) {
1001 if (dK1[j].NumberOfFrequencies())
1002 dT1[j].Mat3(i).noalias() =
1005 (r * dK1[j].Kjj(iz, ia)[i] +
1006 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
1008 if (dK2[j].NumberOfFrequencies())
1009 dT2[j].Mat3(i).noalias() =
1012 (r * dK2[j].Kjj(iz, ia)[i] +
1013 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
1020 const bool real = Const > 0;
1021 const bool imag = Const < 0;
1025 std::sqrt(
imag ? -Const : Const);
1027 (
real ? 1 : -1) * x * x;
1031 const Numeric inv_x2 = inv_x * inv_x;
1034 real ? std::sinh(x) : std::sin(-x);
1036 real ? std::cosh(x) : std::cos(+x);
1045 either ?
a2 * (cx - 1.0) -
a * x * sx + x2 : 1.0 + 0.5 *
a2 -
a;
1046 const Numeric C1 = either ? 2.0 *
a * (1.0 - cx) + x * sx : 1.0 -
a;
1047 const Numeric C2 = either ? cx - 1.0 : 0.5;
1049 T.Mat3(i).noalias() =
1051 (Eigen::Matrix3d() << C0 + C1 *
a + C2 * (
a2 +
b2 + c2),
1052 C1 *
b + C2 * (2 *
a *
b -
c *
u),
1053 C1 *
c + C2 * (2 *
a *
c +
b *
u),
1054 C1 *
b + C2 * (2 *
a *
b +
c *
u),
1055 C0 + C1 *
a + C2 * (
a2 +
b2 - u2),
1056 C1 *
u + C2 * (2 *
a *
u +
b *
c),
1057 C1 *
c + C2 * (2 *
a *
c -
b *
u),
1058 -C1 *
u - C2 * (2 *
a *
u -
b *
c),
1059 C0 + C1 *
a + C2 * (
a2 + c2 - u2))
1062 for (
Index j = 0; j < dT1.nelem(); j++) {
1063 if (dK1[j].NumberOfFrequencies()) {
1064 const Numeric da = -0.5 * (r * dK1[j].Kjj(iz, ia)[i] +
1065 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] +
1068 db = -0.5 * (r * dK1[j].K12(iz, ia)[i] +
1069 ((j == it) ? dr_dT1 * (K1.K12(iz, ia)[i] +
1072 dc = -0.5 * (r * dK1[j].K13(iz, ia)[i] +
1073 ((j == it) ? dr_dT1 * (K1.K13(iz, ia)[i] +
1076 du = -0.5 * (r * dK1[j].K23(iz, ia)[i] +
1077 ((j == it) ? dr_dT1 * (K1.K23(iz, ia)[i] +
1080 const Numeric da2 = 2 *
a * da, db2 = 2 *
b * db, dc2 = 2 *
c * dc,
1082 const Numeric dx = either ? ((db2 + dc2 - du2) * inv_x * 0.5) : 0,
1087 const Numeric dC0 = either ? da2 * (cx - 1.0) + da2 * (dcx - 1.0) -
1088 da * x * sx -
a * dx * sx -
1092 either ? 2.0 * (da * (1.0 - cx) -
a * dcx) + dx * sx + x * dsx
1094 const Numeric dC2 = either ? dcx : 0;
1096 dT1[j].Mat3(i).noalias() =
1097 T.Mat3(i) * (da + dx2 * inv_x2) +
1099 (Eigen::Matrix3d() << dC0 + dC1 *
a + C1 * da +
1100 dC2 * (
a2 +
b2 + c2) +
1101 C2 * (da2 + db2 + dc2),
1102 dC1 *
b + C1 * db + dC2 * (2 *
a *
b -
c *
u) +
1103 C2 * (2 * da *
b - dc *
u + 2 *
a * db -
c * du),
1104 dC1 *
c + C1 * dc + dC2 * (2 *
a *
c +
b *
u) +
1105 C2 * (2 * da *
c + db *
u + 2 *
a * dc +
b * du),
1106 dC1 *
b + C1 * db + dC2 * (2 *
a *
b +
c *
u) +
1107 C2 * (2 * da *
b + dc *
u + 2 *
a * db +
c * du),
1108 dC0 + dC1 *
a + C1 * da + dC2 * (
a2 +
b2 - u2) +
1109 C2 * (da2 + db2 - du2),
1110 dC1 *
u + C1 * du + dC2 * (2 *
a *
u +
b *
c) +
1111 C2 * (2 * da *
u + db *
c + 2 *
a * du +
b * dc),
1112 dC1 *
c + C1 * dc + dC2 * (2 *
a *
c -
b *
u) +
1113 C2 * (2 * da *
c - db *
u + 2 *
a * dc -
b * du),
1114 -dC1 *
u - C1 * du - dC2 * (2 *
a *
u -
b *
c) -
1115 C2 * (2 * da *
u - db *
c + 2 *
a * du -
b * dc),
1116 dC0 + dC1 *
a + C1 * da + dC2 * (
a2 + c2 - u2) +
1117 C2 * (da2 + dc2 - du2))
1120 if (dK2[j].NumberOfFrequencies()) {
1121 const Numeric da = -0.5 * (r * dK2[j].Kjj(iz, ia)[i] +
1122 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] +
1125 db = -0.5 * (r * dK2[j].K12(iz, ia)[i] +
1126 ((j == it) ? dr_dT2 * (K1.K12(iz, ia)[i] +
1129 dc = -0.5 * (r * dK2[j].K13(iz, ia)[i] +
1130 ((j == it) ? dr_dT2 * (K1.K13(iz, ia)[i] +
1133 du = -0.5 * (r * dK2[j].K23(iz, ia)[i] +
1134 ((j == it) ? dr_dT2 * (K1.K23(iz, ia)[i] +
1137 const Numeric da2 = 2 *
a * da, db2 = 2 *
b * db, dc2 = 2 *
c * dc,
1139 const Numeric dx = either ? ((db2 + dc2 - du2) * inv_x * 0.5) : 0,
1144 const Numeric dC0 = either ? da2 * (cx - 1.0) + da2 * (dcx - 1.0) -
1145 da * x * sx -
a * dx * sx -
1149 either ? 2.0 * (da * (1.0 - cx) -
a * dcx) + dx * sx + x * dsx
1151 const Numeric dC2 = either ? dcx : 0;
1153 dT2[j].Mat3(i).noalias() =
1154 T.Mat3(i) * (da + dx2 * inv_x2) +
1156 (Eigen::Matrix3d() << dC0 + dC1 *
a + C1 * da +
1157 dC2 * (
a2 +
b2 + c2) +
1158 C2 * (da2 + db2 + dc2),
1159 dC1 *
b + C1 * db + dC2 * (2 *
a *
b -
c *
u) +
1160 C2 * (2 * da *
b - dc *
u + 2 *
a * db -
c * du),
1161 dC1 *
c + C1 * dc + dC2 * (2 *
a *
c +
b *
u) +
1162 C2 * (2 * da *
c + db *
u + 2 *
a * dc +
b * du),
1163 dC1 *
b + C1 * db + dC2 * (2 *
a *
b +
c *
u) +
1164 C2 * (2 * da *
b + dc *
u + 2 *
a * db +
c * du),
1165 dC0 + dC1 *
a + C1 * da + dC2 * (
a2 +
b2 - u2) +
1166 C2 * (da2 + db2 - du2),
1167 dC1 *
u + C1 * du + dC2 * (2 *
a *
u +
b *
c) +
1168 C2 * (2 * da *
u + db *
c + 2 *
a * du +
b * dc),
1169 dC1 *
c + C1 * dc + dC2 * (2 *
a *
c -
b *
u) +
1170 C2 * (2 * da *
c - db *
u + 2 *
a * dc -
b * du),
1171 -dC1 *
u - C1 * du - dC2 * (2 *
a *
u -
b *
c) -
1172 C2 * (2 * da *
u - db *
c + 2 *
a * du -
b * dc),
1173 dC0 + dC1 *
a + C1 * da + dC2 * (
a2 + c2 - u2) +
1174 C2 * (da2 + dc2 - du2))
1194 const Index ia)
noexcept {
1196 for (
Index i = 0; i < K1.NumberOfFrequencies(); i++) {
1197 const Numeric a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
1198 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]),
1199 c = -0.5 * r * (K1.K13(iz, ia)[i] + K2.K13(iz, ia)[i]),
1200 d = -0.5 * r * (K1.K14(iz, ia)[i] + K2.K14(iz, ia)[i]),
1201 u = -0.5 * r * (K1.K23(iz, ia)[i] + K2.K23(iz, ia)[i]),
1202 v = -0.5 * r * (K1.K24(iz, ia)[i] + K2.K24(iz, ia)[i]),
1203 w = -0.5 * r * (K1.K34(iz, ia)[i] + K2.K34(iz, ia)[i]);
1206 if (
b == 0. and
c == 0. and
d == 0. and
u == 0. and
v == 0. and
w == 0.) {
1207 T.Mat4(i).noalias() = Eigen::Matrix4d::Identity() * exp_a;
1208 for (
Index j = 0; j < dK1.nelem(); j++) {
1209 if (dK1[j].NumberOfFrequencies())
1210 dT1[j].Mat4(i).noalias() =
1213 (r * dK1[j].Kjj(iz, ia)[i] +
1214 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
1216 if (dK2[j].NumberOfFrequencies())
1217 dT2[j].Mat4(i).noalias() =
1220 (r * dK2[j].Kjj(iz, ia)[i] +
1221 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
1228 w2 * w2 + 2 * (
b2 * (
b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
1229 c2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
1230 d2 * (d2 * 0.5 + u2 - v2 - w2) +
1231 u2 * (u2 * 0.5 + v2 + w2) + v2 * (v2 * 0.5 + w2) +
1234 const Numeric Const2 =
b2 + c2 + d2 - u2 - v2 - w2;
1235 const Complex tmp_x_sqrt = std::sqrt(Const2 + Const1);
1236 const Complex tmp_y_sqrt = std::sqrt(Const2 - Const1);
1237 const Complex x = tmp_x_sqrt * sqrt_05;
1241 const Complex cy = std::cos(y);
1242 const Complex sy = std::sin(y);
1243 const Complex cx = std::cosh(x);
1244 const Complex sx = std::sinh(x);
1248 const bool both_zero = y_zero and x_zero;
1249 const bool either_zero = y_zero or x_zero;
1259 const Complex ix = x_zero ? 0.0 : 1.0 / x;
1260 const Complex iy = y_zero ? 0.0 : 1.0 / y;
1266 const Complex C0c = either_zero ? 1.0 : (cy * x2 + cx * y2) * inv_x2y2;
1268 either_zero ? 1.0 : (sy * x2 * iy + sx * y2 * ix) * inv_x2y2;
1269 const Complex C2c = both_zero ? 0.5 : (cx - cy) * inv_x2y2;
1270 const Complex C3c = both_zero ? 1.0 / 6.0
1271 : (x_zero ? 1.0 - sy * iy
1272 : y_zero ? sx * ix - 1.0
1273 : sx * ix - sy * iy) *
1280 T.Mat4(i).noalias() =
1281 exp_a * (Eigen::Matrix4d() << C0 + C2 * (
b2 + c2 + d2),
1282 C1 *
b + C2 * (-
c *
u -
d *
v) +
1283 C3 * (
b * (
b2 + c2 + d2) -
u * (
b *
u -
d *
w) -
1284 v * (
b *
v +
c *
w)),
1285 C1 *
c + C2 * (
b *
u -
d *
w) +
1286 C3 * (
c * (
b2 + c2 + d2) -
u * (
c *
u +
d *
v) -
1287 w * (
b *
v +
c *
w)),
1288 C1 *
d + C2 * (
b *
v +
c *
w) +
1289 C3 * (
d * (
b2 + c2 + d2) -
v * (
c *
u +
d *
v) +
1290 w * (
b *
u -
d *
w)),
1292 C1 *
b + C2 * (
c *
u +
d *
v) +
1293 C3 * (-
b * (-
b2 + u2 + v2) +
c * (
b *
c -
v *
w) +
1294 d * (
b *
d +
u *
w)),
1295 C0 + C2 * (
b2 - u2 - v2),
1296 C2 * (
b *
c -
v *
w) + C1 *
u +
1297 C3 * (
c * (
c *
u +
d *
v) -
u * (-
b2 + u2 + v2) -
1298 w * (
b *
d +
u *
w)),
1299 C2 * (
b *
d +
u *
w) + C1 *
v +
1300 C3 * (
d * (
c *
u +
d *
v) -
v * (-
b2 + u2 + v2) +
1301 w * (
b *
c -
v *
w)),
1303 C1 *
c + C2 * (-
b *
u +
d *
w) +
1304 C3 * (
b * (
b *
c -
v *
w) -
c * (-c2 + u2 + w2) +
1305 d * (
c *
d -
u *
v)),
1306 C2 * (
b *
c -
v *
w) - C1 *
u +
1307 C3 * (-
b * (
b *
u -
d *
w) +
u * (-c2 + u2 + w2) -
1308 v * (
c *
d -
u *
v)),
1309 C0 + C2 * (c2 - u2 - w2),
1310 C2 * (
c *
d -
u *
v) + C1 *
w +
1311 C3 * (-
d * (
b *
u -
d *
w) +
v * (
b *
c -
v *
w) -
1312 w * (-c2 + u2 + w2)),
1314 C1 *
d + C2 * (-
b *
v -
c *
w) +
1315 C3 * (
b * (
b *
d +
u *
w) +
c * (
c *
d -
u *
v) -
1316 d * (-d2 + v2 + w2)),
1317 C2 * (
b *
d +
u *
w) - C1 *
v +
1318 C3 * (-
b * (
b *
v +
c *
w) -
u * (
c *
d -
u *
v) +
1319 v * (-d2 + v2 + w2)),
1320 C2 * (
c *
d -
u *
v) - C1 *
w +
1321 C3 * (-
c * (
b *
v +
c *
w) +
u * (
b *
d +
u *
w) +
1322 w * (-d2 + v2 + w2)),
1323 C0 + C2 * (d2 - v2 - w2))
1326 for (
Index j = 0; j < dK1.nelem(); j++) {
1327 if (dK1[j].NumberOfFrequencies()) {
1328 const Numeric da = -0.5 * (r * dK1[j].Kjj(iz, ia)[i] +
1329 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] +
1332 db = -0.5 * (r * dK1[j].K12(iz, ia)[i] +
1333 ((j == it) ? dr_dT1 * (K1.K12(iz, ia)[i] +
1336 dc = -0.5 * (r * dK1[j].K13(iz, ia)[i] +
1337 ((j == it) ? dr_dT1 * (K1.K13(iz, ia)[i] +
1340 dd = -0.5 * (r * dK1[j].K14(iz, ia)[i] +
1341 ((j == it) ? dr_dT1 * (K1.K14(iz, ia)[i] +
1344 du = -0.5 * (r * dK1[j].K23(iz, ia)[i] +
1345 ((j == it) ? dr_dT1 * (K1.K23(iz, ia)[i] +
1348 dv = -0.5 * (r * dK1[j].K24(iz, ia)[i] +
1349 ((j == it) ? dr_dT1 * (K1.K24(iz, ia)[i] +
1352 dw = -0.5 * (r * dK1[j].K34(iz, ia)[i] +
1353 ((j == it) ? dr_dT1 * (K1.K34(iz, ia)[i] +
1356 const Numeric db2 = 2 * db *
b, dc2 = 2 * dc *
c, dd2 = 2 * dd *
d,
1357 du2 = 2 * du *
u, dv2 = 2 * dv *
v, dw2 = 2 * dw *
w;
1360 2 * (db2 * (
b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
1361 b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2) +
1362 dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
1363 c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2) +
1364 dd2 * (d2 * 0.5 + u2 - v2 - w2) +
1365 d2 * (dd2 * 0.5 + du2 - dv2 - dw2) +
1366 du2 * (u2 * 0.5 + v2 + w2) + u2 * (du2 * 0.5 + dv2 + dw2) +
1367 dv2 * (v2 * 0.5 + w2) + v2 * (dv2 * 0.5 + dw2) +
1368 4 * (db *
d *
u *
w - db *
c *
v *
w - dc *
d *
u *
v +
1369 b * dd *
u *
w -
b * dc *
v *
w -
c * dd *
u *
v +
1370 b *
d * du *
w -
b *
c * dv *
w -
c *
d * du *
v +
1371 b *
d *
u * dw -
b *
c *
v * dw -
c *
d *
u * dv));
1372 const Complex dConst1 = 0.5 * dtmp / Const1;
1373 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1375 x_zero ? 0 : (0.5 * (dConst2 + dConst1) / tmp_x_sqrt) * sqrt_05;
1377 : (0.5 * (dConst2 - dConst1) / tmp_y_sqrt) *
1379 const Complex dx2 = 2 * x * dx;
1380 const Complex dy2 = 2 * y * dy;
1385 const Complex dix = -dx * ix * ix;
1386 const Complex diy = -dy * iy * iy;
1387 const Complex dx2dy2 = dx2 + dy2;
1391 : (dcy * x2 + cy * dx2 + dcx * y2 + cx * dy2 - C0c * dx2dy2) *
1395 : (dsy * x2 * iy + sy * dx2 * iy + sy * x2 * diy +
1396 dsx * y2 * ix + sx * dy2 * ix + sx * y2 * dix -
1400 both_zero ? 0.0 : (dcx - dcy - C2c * dx2dy2) * inv_x2y2;
1404 : ((x_zero ? -dsy * iy - sy * diy
1405 : y_zero ? dsx * ix + sx * dix
1406 : dsx * ix + sx * dix - dsy * iy - sy * diy) -
1414 dT1[j].Mat4(i).noalias() =
1418 << dC0 + dC2 * (
b2 + c2 + d2) + C2 * (db2 + dc2 + dd2),
1419 db * C1 +
b * dC1 + dC2 * (-
c *
u -
d *
v) +
1420 C2 * (-dc *
u - dd *
v -
c * du -
d * dv) +
1421 dC3 * (
b * (
b2 + c2 + d2) -
u * (
b *
u -
d *
w) -
1422 v * (
b *
v +
c *
w)) +
1423 C3 * (db * (
b2 + c2 + d2) - du * (
b *
u -
d *
w) -
1424 dv * (
b *
v +
c *
w) +
b * (db2 + dc2 + dd2) -
1425 u * (db *
u - dd *
w) -
v * (db *
v + dc *
w) -
1426 u * (
b * du -
d * dw) -
v * (
b * dv +
c * dw)),
1427 dC1 *
c + C1 * dc + dC2 * (
b *
u -
d *
w) +
1428 C2 * (db *
u - dd *
w +
b * du -
d * dw) +
1429 dC3 * (
c * (
b2 + c2 + d2) -
u * (
c *
u +
d *
v) -
1430 w * (
b *
v +
c *
w)) +
1431 C3 * (dc * (
b2 + c2 + d2) - du * (
c *
u +
d *
v) -
1432 dw * (
b *
v +
c *
w) +
c * (db2 + dc2 + dd2) -
1433 u * (dc *
u + dd *
v) -
w * (db *
v + dc *
w) -
1434 u * (
c * du +
d * dv) -
w * (
b * dv +
c * dw)),
1435 dC1 *
d + C1 * dd + dC2 * (
b *
v +
c *
w) +
1436 C2 * (db *
v + dc *
w +
b * dv +
c * dw) +
1437 dC3 * (
d * (
b2 + c2 + d2) -
v * (
c *
u +
d *
v) +
1438 w * (
b *
u -
d *
w)) +
1439 C3 * (dd * (
b2 + c2 + d2) - dv * (
c *
u +
d *
v) +
1440 dw * (
b *
u -
d *
w) +
d * (db2 + dc2 + dd2) -
1441 v * (dc *
u + dd *
v) +
w * (db *
u - dd *
w) -
1442 v * (
c * du +
d * dv) +
w * (
b * du -
d * dw)),
1444 db * C1 +
b * dC1 + dC2 * (
c *
u +
d *
v) +
1445 C2 * (dc *
u + dd *
v +
c * du +
d * dv) +
1446 dC3 * (-
b * (-
b2 + u2 + v2) +
c * (
b *
c -
v *
w) +
1447 d * (
b *
d +
u *
w)) +
1448 C3 * (-db * (-
b2 + u2 + v2) + dc * (
b *
c -
v *
w) +
1449 dd * (
b *
d +
u *
w) -
b * (-db2 + du2 + dv2) +
1450 c * (db *
c - dv *
w) +
d * (db *
d + du *
w) +
1451 c * (
b * dc -
v * dw) +
d * (
b * dd +
u * dw)),
1452 dC0 + dC2 * (
b2 - u2 - v2) + C2 * (db2 - du2 - dv2),
1453 dC2 * (
b *
c -
v *
w) +
1454 C2 * (db *
c +
b * dc - dv *
w -
v * dw) + dC1 *
u +
1456 dC3 * (
c * (
c *
u +
d *
v) -
u * (-
b2 + u2 + v2) -
1457 w * (
b *
d +
u *
w)) +
1458 C3 * (dc * (
c *
u +
d *
v) - du * (-
b2 + u2 + v2) -
1459 dw * (
b *
d +
u *
w) +
c * (dc *
u + dd *
v) -
1460 u * (-db2 + du2 + dv2) -
w * (db *
d + du *
w) +
1461 c * (
c * du +
d * dv) -
w * (
b * dd +
u * dw)),
1462 dC2 * (
b *
d +
u *
w) +
1463 C2 * (db *
d +
b * dd + du *
w +
u * dw) + dC1 *
v +
1465 dC3 * (
d * (
c *
u +
d *
v) -
v * (-
b2 + u2 + v2) +
1466 w * (
b *
c -
v *
w)) +
1467 C3 * (dd * (
c *
u +
d *
v) - dv * (-
b2 + u2 + v2) +
1468 dw * (
b *
c -
v *
w) +
d * (dc *
u + dd *
v) -
1469 v * (-db2 + du2 + dv2) +
w * (db *
c - dv *
w) +
1470 d * (
c * du +
d * dv) +
w * (
b * dc -
v * dw)),
1472 dC1 *
c + C1 * dc + dC2 * (-
b *
u +
d *
w) +
1473 C2 * (-db *
u + dd *
w -
b * du +
d * dw) +
1474 dC3 * (
b * (
b *
c -
v *
w) -
c * (-c2 + u2 + w2) +
1475 d * (
c *
d -
u *
v)) +
1476 C3 * (db * (
b *
c -
v *
w) - dc * (-c2 + u2 + w2) +
1477 dd * (
c *
d -
u *
v) +
b * (db *
c - dv *
w) -
1478 c * (-dc2 + du2 + dw2) +
d * (dc *
d - du *
v) +
1479 b * (
b * dc -
v * dw) +
d * (
c * dd -
u * dv)),
1480 dC2 * (
b *
c -
v *
w) +
1481 C2 * (db *
c +
b * dc - dv *
w -
v * dw) - dC1 *
u -
1483 dC3 * (-
b * (
b *
u -
d *
w) +
u * (-c2 + u2 + w2) -
1484 v * (
c *
d -
u *
v)) +
1485 C3 * (-db * (
b *
u -
d *
w) + du * (-c2 + u2 + w2) -
1486 dv * (
c *
d -
u *
v) -
b * (db *
u - dd *
w) +
1487 u * (-dc2 + du2 + dw2) -
v * (dc *
d - du *
v) -
1488 b * (
b * du -
d * dw) -
v * (
c * dd -
u * dv)),
1489 dC0 + dC2 * (c2 - u2 - w2) + C2 * (dc2 - du2 - dw2),
1490 dC2 * (
c *
d -
u *
v) +
1491 C2 * (dc *
d +
c * dd - du *
v -
u * dv) + dC1 *
w +
1493 dC3 * (-
d * (
b *
u -
d *
w) +
v * (
b *
c -
v *
w) -
1494 w * (-c2 + u2 + w2)) +
1495 C3 * (-dd * (
b *
u -
d *
w) + dv * (
b *
c -
v *
w) -
1496 dw * (-c2 + u2 + w2) -
d * (db *
u - dd *
w) +
1497 v * (db *
c - dv *
w) -
w * (-dc2 + du2 + dw2) -
1498 d * (
b * du -
d * dw) +
v * (
b * dc -
v * dw)),
1500 dC1 *
d + C1 * dd + dC2 * (-
b *
v -
c *
w) +
1501 C2 * (-db *
v - dc *
w -
b * dv -
c * dw) +
1502 dC3 * (
b * (
b *
d +
u *
w) +
c * (
c *
d -
u *
v) -
1503 d * (-d2 + v2 + w2)) +
1504 C3 * (db * (
b *
d +
u *
w) + dc * (
c *
d -
u *
v) -
1505 dd * (-d2 + v2 + w2) +
b * (db *
d + du *
w) +
1506 c * (dc *
d - du *
v) -
d * (-dd2 + dv2 + dw2) +
1507 b * (
b * dd +
u * dw) +
c * (
c * dd -
u * dv)),
1508 dC2 * (
b *
d +
u *
w) +
1509 C2 * (db *
d +
b * dd + du *
w +
u * dw) - dC1 *
v -
1511 dC3 * (-
b * (
b *
v +
c *
w) -
u * (
c *
d -
u *
v) +
1512 v * (-d2 + v2 + w2)) +
1513 C3 * (-db * (
b *
v +
c *
w) - du * (
c *
d -
u *
v) +
1514 dv * (-d2 + v2 + w2) -
b * (db *
v + dc *
w) -
1515 u * (dc *
d - du *
v) +
v * (-dd2 + dv2 + dw2) -
1516 b * (
b * dv +
c * dw) -
u * (
c * dd -
u * dv)),
1517 dC2 * (
c *
d -
u *
v) +
1518 C2 * (dc *
d +
c * dd - du *
v -
u * dv) - dC1 *
w -
1520 dC3 * (-
c * (
b *
v +
c *
w) +
u * (
b *
d +
u *
w) +
1521 w * (-d2 + v2 + w2)) +
1522 C3 * (-dc * (
b *
v +
c *
w) + du * (
b *
d +
u *
w) +
1523 dw * (-d2 + v2 + w2) -
c * (db *
v + dc *
w) +
1524 u * (db *
d + du *
w) +
w * (-dd2 + dv2 + dw2) -
1525 c * (
b * dv +
c * dw) +
u * (
b * dd +
u * dw)),
1526 dC0 + dC2 * (d2 - v2 - w2) + C2 * (dd2 - dv2 - dw2))
1529 if (dK2[j].NumberOfFrequencies()) {
1530 const Numeric da = -0.5 * (r * dK2[j].Kjj(iz, ia)[i] +
1531 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] +
1534 db = -0.5 * (r * dK2[j].K12(iz, ia)[i] +
1535 ((j == it) ? dr_dT2 * (K1.K12(iz, ia)[i] +
1538 dc = -0.5 * (r * dK2[j].K13(iz, ia)[i] +
1539 ((j == it) ? dr_dT2 * (K1.K13(iz, ia)[i] +
1542 dd = -0.5 * (r * dK2[j].K14(iz, ia)[i] +
1543 ((j == it) ? dr_dT2 * (K1.K14(iz, ia)[i] +
1546 du = -0.5 * (r * dK2[j].K23(iz, ia)[i] +
1547 ((j == it) ? dr_dT2 * (K1.K23(iz, ia)[i] +
1550 dv = -0.5 * (r * dK2[j].K24(iz, ia)[i] +
1551 ((j == it) ? dr_dT2 * (K1.K24(iz, ia)[i] +
1554 dw = -0.5 * (r * dK2[j].K34(iz, ia)[i] +
1555 ((j == it) ? dr_dT2 * (K1.K34(iz, ia)[i] +
1558 const Numeric db2 = 2 * db *
b, dc2 = 2 * dc *
c, dd2 = 2 * dd *
d,
1559 du2 = 2 * du *
u, dv2 = 2 * dv *
v, dw2 = 2 * dw *
w;
1562 2 * (db2 * (
b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
1563 b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2) +
1564 dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
1565 c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2) +
1566 dd2 * (d2 * 0.5 + u2 - v2 - w2) +
1567 d2 * (dd2 * 0.5 + du2 - dv2 - dw2) +
1568 du2 * (u2 * 0.5 + v2 + w2) + u2 * (du2 * 0.5 + dv2 + dw2) +
1569 dv2 * (v2 * 0.5 + w2) + v2 * (dv2 * 0.5 + dw2) +
1570 4 * (db *
d *
u *
w - db *
c *
v *
w - dc *
d *
u *
v +
1571 b * dd *
u *
w -
b * dc *
v *
w -
c * dd *
u *
v +
1572 b *
d * du *
w -
b *
c * dv *
w -
c *
d * du *
v +
1573 b *
d *
u * dw -
b *
c *
v * dw -
c *
d *
u * dv));
1574 const Complex dConst1 = 0.5 * dtmp / Const1;
1575 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1577 x_zero ? 0 : (0.5 * (dConst2 + dConst1) / tmp_x_sqrt) * sqrt_05;
1579 : (0.5 * (dConst2 - dConst1) / tmp_y_sqrt) *
1581 const Complex dx2 = 2 * x * dx;
1582 const Complex dy2 = 2 * y * dy;
1587 const Complex dix = -dx * ix * ix;
1588 const Complex diy = -dy * iy * iy;
1589 const Complex dx2dy2 = dx2 + dy2;
1593 : (dcy * x2 + cy * dx2 + dcx * y2 + cx * dy2 - C0c * dx2dy2) *
1597 : (dsy * x2 * iy + sy * dx2 * iy + sy * x2 * diy +
1598 dsx * y2 * ix + sx * dy2 * ix + sx * y2 * dix -
1602 both_zero ? 0.0 : (dcx - dcy - C2c * dx2dy2) * inv_x2y2;
1606 : ((x_zero ? -dsy * iy - sy * diy
1607 : y_zero ? dsx * ix + sx * dix
1608 : dsx * ix + sx * dix - dsy * iy - sy * diy) -
1616 dT2[j].Mat4(i).noalias() =
1620 << dC0 + dC2 * (
b2 + c2 + d2) + C2 * (db2 + dc2 + dd2),
1621 db * C1 +
b * dC1 + dC2 * (-
c *
u -
d *
v) +
1622 C2 * (-dc *
u - dd *
v -
c * du -
d * dv) +
1623 dC3 * (
b * (
b2 + c2 + d2) -
u * (
b *
u -
d *
w) -
1624 v * (
b *
v +
c *
w)) +
1625 C3 * (db * (
b2 + c2 + d2) - du * (
b *
u -
d *
w) -
1626 dv * (
b *
v +
c *
w) +
b * (db2 + dc2 + dd2) -
1627 u * (db *
u - dd *
w) -
v * (db *
v + dc *
w) -
1628 u * (
b * du -
d * dw) -
v * (
b * dv +
c * dw)),
1629 dC1 *
c + C1 * dc + dC2 * (
b *
u -
d *
w) +
1630 C2 * (db *
u - dd *
w +
b * du -
d * dw) +
1631 dC3 * (
c * (
b2 + c2 + d2) -
u * (
c *
u +
d *
v) -
1632 w * (
b *
v +
c *
w)) +
1633 C3 * (dc * (
b2 + c2 + d2) - du * (
c *
u +
d *
v) -
1634 dw * (
b *
v +
c *
w) +
c * (db2 + dc2 + dd2) -
1635 u * (dc *
u + dd *
v) -
w * (db *
v + dc *
w) -
1636 u * (
c * du +
d * dv) -
w * (
b * dv +
c * dw)),
1637 dC1 *
d + C1 * dd + dC2 * (
b *
v +
c *
w) +
1638 C2 * (db *
v + dc *
w +
b * dv +
c * dw) +
1639 dC3 * (
d * (
b2 + c2 + d2) -
v * (
c *
u +
d *
v) +
1640 w * (
b *
u -
d *
w)) +
1641 C3 * (dd * (
b2 + c2 + d2) - dv * (
c *
u +
d *
v) +
1642 dw * (
b *
u -
d *
w) +
d * (db2 + dc2 + dd2) -
1643 v * (dc *
u + dd *
v) +
w * (db *
u - dd *
w) -
1644 v * (
c * du +
d * dv) +
w * (
b * du -
d * dw)),
1646 db * C1 +
b * dC1 + dC2 * (
c *
u +
d *
v) +
1647 C2 * (dc *
u + dd *
v +
c * du +
d * dv) +
1648 dC3 * (-
b * (-
b2 + u2 + v2) +
c * (
b *
c -
v *
w) +
1649 d * (
b *
d +
u *
w)) +
1650 C3 * (-db * (-
b2 + u2 + v2) + dc * (
b *
c -
v *
w) +
1651 dd * (
b *
d +
u *
w) -
b * (-db2 + du2 + dv2) +
1652 c * (db *
c - dv *
w) +
d * (db *
d + du *
w) +
1653 c * (
b * dc -
v * dw) +
d * (
b * dd +
u * dw)),
1654 dC0 + dC2 * (
b2 - u2 - v2) + C2 * (db2 - du2 - dv2),
1655 dC2 * (
b *
c -
v *
w) +
1656 C2 * (db *
c +
b * dc - dv *
w -
v * dw) + dC1 *
u +
1658 dC3 * (
c * (
c *
u +
d *
v) -
u * (-
b2 + u2 + v2) -
1659 w * (
b *
d +
u *
w)) +
1660 C3 * (dc * (
c *
u +
d *
v) - du * (-
b2 + u2 + v2) -
1661 dw * (
b *
d +
u *
w) +
c * (dc *
u + dd *
v) -
1662 u * (-db2 + du2 + dv2) -
w * (db *
d + du *
w) +
1663 c * (
c * du +
d * dv) -
w * (
b * dd +
u * dw)),
1664 dC2 * (
b *
d +
u *
w) +
1665 C2 * (db *
d +
b * dd + du *
w +
u * dw) + dC1 *
v +
1667 dC3 * (
d * (
c *
u +
d *
v) -
v * (-
b2 + u2 + v2) +
1668 w * (
b *
c -
v *
w)) +
1669 C3 * (dd * (
c *
u +
d *
v) - dv * (-
b2 + u2 + v2) +
1670 dw * (
b *
c -
v *
w) +
d * (dc *
u + dd *
v) -
1671 v * (-db2 + du2 + dv2) +
w * (db *
c - dv *
w) +
1672 d * (
c * du +
d * dv) +
w * (
b * dc -
v * dw)),
1674 dC1 *
c + C1 * dc + dC2 * (-
b *
u +
d *
w) +
1675 C2 * (-db *
u + dd *
w -
b * du +
d * dw) +
1676 dC3 * (
b * (
b *
c -
v *
w) -
c * (-c2 + u2 + w2) +
1677 d * (
c *
d -
u *
v)) +
1678 C3 * (db * (
b *
c -
v *
w) - dc * (-c2 + u2 + w2) +
1679 dd * (
c *
d -
u *
v) +
b * (db *
c - dv *
w) -
1680 c * (-dc2 + du2 + dw2) +
d * (dc *
d - du *
v) +
1681 b * (
b * dc -
v * dw) +
d * (
c * dd -
u * dv)),
1682 dC2 * (
b *
c -
v *
w) +
1683 C2 * (db *
c +
b * dc - dv *
w -
v * dw) - dC1 *
u -
1685 dC3 * (-
b * (
b *
u -
d *
w) +
u * (-c2 + u2 + w2) -
1686 v * (
c *
d -
u *
v)) +
1687 C3 * (-db * (
b *
u -
d *
w) + du * (-c2 + u2 + w2) -
1688 dv * (
c *
d -
u *
v) -
b * (db *
u - dd *
w) +
1689 u * (-dc2 + du2 + dw2) -
v * (dc *
d - du *
v) -
1690 b * (
b * du -
d * dw) -
v * (
c * dd -
u * dv)),
1691 dC0 + dC2 * (c2 - u2 - w2) + C2 * (dc2 - du2 - dw2),
1692 dC2 * (
c *
d -
u *
v) +
1693 C2 * (dc *
d +
c * dd - du *
v -
u * dv) + dC1 *
w +
1695 dC3 * (-
d * (
b *
u -
d *
w) +
v * (
b *
c -
v *
w) -
1696 w * (-c2 + u2 + w2)) +
1697 C3 * (-dd * (
b *
u -
d *
w) + dv * (
b *
c -
v *
w) -
1698 dw * (-c2 + u2 + w2) -
d * (db *
u - dd *
w) +
1699 v * (db *
c - dv *
w) -
w * (-dc2 + du2 + dw2) -
1700 d * (
b * du -
d * dw) +
v * (
b * dc -
v * dw)),
1702 dC1 *
d + C1 * dd + dC2 * (-
b *
v -
c *
w) +
1703 C2 * (-db *
v - dc *
w -
b * dv -
c * dw) +
1704 dC3 * (
b * (
b *
d +
u *
w) +
c * (
c *
d -
u *
v) -
1705 d * (-d2 + v2 + w2)) +
1706 C3 * (db * (
b *
d +
u *
w) + dc * (
c *
d -
u *
v) -
1707 dd * (-d2 + v2 + w2) +
b * (db *
d + du *
w) +
1708 c * (dc *
d - du *
v) -
d * (-dd2 + dv2 + dw2) +
1709 b * (
b * dd +
u * dw) +
c * (
c * dd -
u * dv)),
1710 dC2 * (
b *
d +
u *
w) +
1711 C2 * (db *
d +
b * dd + du *
w +
u * dw) - dC1 *
v -
1713 dC3 * (-
b * (
b *
v +
c *
w) -
u * (
c *
d -
u *
v) +
1714 v * (-d2 + v2 + w2)) +
1715 C3 * (-db * (
b *
v +
c *
w) - du * (
c *
d -
u *
v) +
1716 dv * (-d2 + v2 + w2) -
b * (db *
v + dc *
w) -
1717 u * (dc *
d - du *
v) +
v * (-dd2 + dv2 + dw2) -
1718 b * (
b * dv +
c * dw) -
u * (
c * dd -
u * dv)),
1719 dC2 * (
c *
d -
u *
v) +
1720 C2 * (dc *
d +
c * dd - du *
v -
u * dv) - dC1 *
w -
1722 dC3 * (-
c * (
b *
v +
c *
w) +
u * (
b *
d +
u *
w) +
1723 w * (-d2 + v2 + w2)) +
1724 C3 * (-dc * (
b *
v +
c *
w) + du * (
b *
d +
u *
w) +
1725 dw * (-d2 + v2 + w2) -
c * (db *
v + dc *
w) +
1726 u * (db *
d + du *
w) +
w * (-dd2 + dv2 + dw2) -
1727 c * (
b * dv +
c * dw) +
u * (
b * dd +
u * dw)),
1728 dC0 + dC2 * (d2 - v2 - w2) + C2 * (dd2 - dv2 - dw2))
1740 switch (K1.StokesDimensions()) {
1766 const Index it = -1,
1768 const Index ia = 0) noexcept {
1771 dtransmat4(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1774 dtransmat3(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1777 dtransmat2(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1780 dtransmat1(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1795 const Index temp_deriv_pos) {
1796 if (not dT1.
nelem())
1800 T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dtemp1, dr_dtemp2, temp_deriv_pos);
1815 const bool& jacobian_do) {
1820 for (
Index j = 0; j < jacobian_quantities.
nelem(); j++)
1821 if (dJ[j].Frequencies()) dJ[j].SetZero(i);
1828 const auto invK = inv_prop_matrix<4>(K.
Data()(0, 0, i,
joker));
1831 for (
Index j = 0; j < jacobian_quantities.
nelem(); j++) {
1832 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and
1833 dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1834 dJ[j].Vec4(i).noalias() =
1842 jacobian_quantities[j] == Jacobian::Atm::Temperature,
1844 prop_matrix<4>(dK[j].Data()(0, 0, i,
joker)) * J.
Vec4(i));
1849 J_add.
Vec4(i) = invK * J_add.
Vec4(i);
1854 const auto invK = inv_prop_matrix<3>(K.
Data()(0, 0, i,
joker));
1857 for (
Index j = 0; j < jacobian_quantities.
nelem(); j++) {
1859 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and
1860 dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1861 dJ[j].Vec3(i).noalias() =
1869 jacobian_quantities[j] == Jacobian::Atm::Temperature,
1871 prop_matrix<3>(dK[j].Data()(0, 0, i,
joker)) * J.
Vec3(i));
1876 J_add.
Vec3(i) = invK * J_add.
Vec3(i);
1881 const auto invK = inv_prop_matrix<2>(K.
Data()(0, 0, i,
joker));
1884 for (
Index j = 0; j < jacobian_quantities.
nelem(); j++) {
1886 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and
1887 dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1888 dJ[j].Vec2(i).noalias() =
1896 jacobian_quantities[j] == Jacobian::Atm::Temperature,
1898 prop_matrix<2>(dK[j].Data()(0, 0, i,
joker)) * J.
Vec2(i));
1903 J_add.
Vec2(i) = invK * J_add.
Vec2(i);
1908 const auto invK = inv_prop_matrix<1>(K.
Data()(0, 0, i,
joker));
1911 for (
Index j = 0; j < jacobian_quantities.
nelem(); j++) {
1913 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and
1914 dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1915 dJ[j].Vec1(i).noalias() =
1923 jacobian_quantities[j] == Jacobian::Atm::Temperature,
1925 prop_matrix<1>(dK[j].Data()(0, 0, i,
joker)) * J.
Vec1(i));
1930 J_add.
Vec1(i) = invK * J_add.
Vec1(i);
1959 [[maybe_unused]]
const Numeric r,
1960 [[maybe_unused]]
const Vector& dr1,
1961 [[maybe_unused]]
const Vector& dr2,
1962 [[maybe_unused]]
const Index ia,
1963 [[maybe_unused]]
const Index iz,
1968 for (
size_t i = 0; i < dI1.size(); i++) {
1969 dI1[i].addDerivEmission(PiT, dT1[i], T, I, dJ1[i]);
1970 dI2[i].addDerivEmission(PiT, dT2[i], T, I, dJ2[i]);
1977 for (
size_t i = 0; i < dI1.size(); i++) {
1978 dI1[i].addDerivTransmission(PiT, dT1[i], I);
1979 dI2[i].addDerivTransmission(PiT, dT2[i], I);
1987 "Cannot support derivatives with current integration method\n");
2006 const Index nf = n ? T[0].Frequencies() : 1;
2007 const Index ns = n ? T[0].stokes_dim : 1;
2011 for (
Index i = 1; i < n; i++) PiT[i].
mul(PiT[i - 1], T[i]);
2014 for (
Index i = 1; i < n; i++) PiT[i].
mul(T[i], PiT[i - 1]);
2035 const Index nv = np ? I[0].Frequencies() : 0;
2036 const Index ns = np ? I[0].stokes_dim : 1;
2040 for (
Index ip = 0; ip < np; ip++)
2041 I[ip].setBackscatterTransmission(I_incoming, PiTr[ip], PiTf[ip], Z[ip]);
2043 for (
Index ip = 0; ip < np; ip++) {
2044 for (
Index iq = 0; iq < nq; iq++) {
2045 dI[ip][ip][iq].setBackscatterTransmissionDerivative(
2046 I_incoming, PiTr[ip], PiTf[ip], dZ[ip][iq]);
2055 BackscatterSolverCommutativeTransmissionStokesDimOne:
2056 for (
Index ip = 0; ip < np; ip++) {
2057 for (
Index j = ip; j < np; j++) {
2058 for (
Index iq = 0; iq < nq; iq++) {
2059 for (
Index iv = 0; iv < nv; iv++) {
2060 dI[ip][j][iq].Vec1(iv).noalias() +=
2061 T[ip].Mat1(iv).inverse() *
2062 (dT1[ip][iq].Mat1(iv) + dT2[ip][iq].Mat1(iv)) *
2065 if (j < np - 1 and j > ip)
2066 dI[ip][j][iq].Vec1(iv).noalias() +=
2067 T[ip + 1].Mat1(iv).inverse() *
2068 (dT1[ip][iq].Mat1(iv) + dT2[ip][iq].Mat1(iv)) *
2076 for (
Index ip = 0; ip < np; ip++) {
2077 for (
Index j = ip; j < np; j++) {
2078 for (
Index iq = 0; iq < nq; iq++) {
2079 for (
Index iv = 0; iv < nv; iv++) {
2080 dI[ip][j][iq].Vec2(iv).noalias() +=
2081 T[ip].Mat2(iv).inverse() *
2082 (dT1[ip][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
2085 if (j < np - 1 and j > ip)
2086 dI[ip][j][iq].Vec2(iv).noalias() +=
2087 T[ip + 1].Mat2(iv).inverse() *
2088 (dT1[ip][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
2096 for (
Index ip = 0; ip < np; ip++) {
2097 for (
Index j = ip; j < np; j++) {
2098 for (
Index iq = 0; iq < nq; iq++) {
2099 for (
Index iv = 0; iv < nv; iv++) {
2100 dI[ip][j][iq].Vec3(iv).noalias() +=
2101 T[ip].Mat3(iv).inverse() *
2102 (dT1[ip][iq].Mat3(iv) + dT2[ip][iq].Mat3(iv)) *
2105 if (j < np - 1 and j > ip)
2106 dI[ip][j][iq].Vec3(iv).noalias() +=
2107 T[ip + 1].Mat3(iv).inverse() *
2108 (dT1[ip][iq].Mat3(iv) + dT2[ip][iq].Mat3(iv)) *
2116 for (
Index ip = 0; ip < np; ip++) {
2117 for (
Index j = ip; j < np; j++) {
2118 for (
Index iq = 0; iq < nq; iq++) {
2119 for (
Index iv = 0; iv < nv; iv++) {
2120 dI[ip][j][iq].Vec4(iv).noalias() +=
2121 T[ip].Mat4(iv).inverse() *
2122 (dT1[ip][iq].Mat4(iv) + dT2[ip][iq].Mat4(iv)) *
2125 if (j < np - 1 and j > ip)
2126 dI[ip][j][iq].Vec4(iv).noalias() +=
2127 T[ip + 1].Mat4(iv).inverse() *
2128 (dT1[ip][iq].Mat4(iv) + dT2[ip][iq].Mat4(iv)) *
2138 std::runtime_error(
"Do not activate this code. It is not ready yet");
2142 goto BackscatterSolverCommutativeTransmissionStokesDimOne;
2145 for (
Index ip = 0; ip < np; ip++) {
2146 for (
Index j = ip; j < np; j++) {
2147 for (
Index iq = 0; iq < nq; iq++) {
2148 for (
Index iv = 0; iv < nv; iv++) {
2150 dI[ip][j][iq].Vec2(iv).noalias() +=
2151 PiTr[ip - 2].Mat2(iv) *
2152 (dT1[ip - 1][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
2153 PiTr[ip - 1].Mat2(iv).inverse() * I[j].Vec2(iv);
2156 dI[ip][j][iq].Vec2(iv).noalias() +=
2157 PiTr[ip].Mat2(iv) * Z[ip].Mat2(iv) *
2158 PiTf[ip - 2].Mat2(iv) *
2159 (dT1[ip][iq].Mat2(iv) + dT2[ip + 1][iq].Mat2(iv)) *
2160 PiTf[ip - 1].Mat2(iv).inverse() * PiTf[ip].Mat2(iv) *
2163 dI[ip][j][iq].Vec2(iv).noalias() +=
2164 (dT1[ip - 1][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
2165 PiTr[ip - 1].Mat2(iv).inverse() * I[j].Vec2(iv);
2168 dI[ip][j][iq].Vec2(iv).noalias() +=
2169 PiTr[ip].Mat2(iv) * Z[ip].Mat2(iv) *
2170 (dT1[ip][iq].Mat2(iv) + dT2[ip + 1][iq].Mat2(iv)) *
2171 PiTf[ip - 1].Mat2(iv).inverse() * PiTf[ip].Mat2(iv) *
2197 for (
Index ip = 0; ip < np; ip++) {
2202 for (
Index iv = 0; iv < nv; iv++)
2203 for (
Index ie = 0; ie < ne; ie++)
2204 aotm[ip].Mat4(iv).noalias() +=
2205 pnd(ie, ip) * prop_matrix<4>(Pe(ie, ip, iv,
joker,
joker));
2208 for (
Index iv = 0; iv < nv; iv++)
2209 for (
Index ie = 0; ie < ne; ie++)
2210 aotm[ip].Mat3(iv).noalias() +=
2211 pnd(ie, ip) * prop_matrix<3>(Pe(ie, ip, iv,
joker,
joker));
2214 for (
Index iv = 0; iv < nv; iv++)
2215 for (
Index ie = 0; ie < ne; ie++)
2216 aotm[ip].Mat2(iv).noalias() +=
2217 pnd(ie, ip) * prop_matrix<2>(Pe(ie, ip, iv,
joker,
joker));
2220 for (
Index iv = 0; iv < nv; iv++)
2221 for (
Index ie = 0; ie < ne; ie++)
2222 aotm[ip].Mat1(iv).noalias() +=
2223 pnd(ie, ip) * prop_matrix<1>(Pe(ie, ip, iv,
joker,
joker));
2241 for (
Index ip = 0; ip < np; ip++) {
2242 for (
Index iq = 0; iq < nq; iq++) {
2243 aoaotm[ip][iq].setZero();
2245 if (not dpnd_dx[iq].empty()) {
2248 for (
Index iv = 0; iv < nv; iv++)
2249 for (
Index ie = 0; ie < ne; ie++)
2250 aoaotm[ip][iq].Mat4(iv).noalias() +=
2251 dpnd_dx[iq](ie, ip) *
2252 prop_matrix<4>(Pe(ie, ip, iv,
joker,
joker));
2255 for (
Index iv = 0; iv < nv; iv++)
2256 for (
Index ie = 0; ie < ne; ie++)
2257 aoaotm[ip][iq].Mat3(iv).noalias() +=
2258 dpnd_dx[iq](ie, ip) *
2259 prop_matrix<3>(Pe(ie, ip, iv,
joker,
joker));
2262 for (
Index iv = 0; iv < nv; iv++)
2263 for (
Index ie = 0; ie < ne; ie++)
2264 aoaotm[ip][iq].Mat2(iv).noalias() +=
2265 dpnd_dx[iq](ie, ip) *
2266 prop_matrix<2>(Pe(ie, ip, iv,
joker,
joker));
2269 for (
Index iv = 0; iv < nv; iv++)
2270 for (
Index ie = 0; ie < ne; ie++)
2271 aoaotm[ip][iq].Mat1(iv).noalias() +=
2272 dpnd_dx[iq](ie, ip) *
2273 prop_matrix<1>(Pe(ie, ip, iv,
joker,
joker));
2285 for (
const auto& T : tm.
T4) os << T <<
'\n';
2286 for (
const auto& T : tm.
T3) os << T <<
'\n';
2287 for (
const auto& T : tm.
T2) os << T <<
'\n';
2288 for (
const auto& T : tm.
T1) os << T <<
'\n';
2294 for (
const auto& R : rv.
R4) os << R.transpose() <<
'\n';
2295 for (
const auto& R : rv.
R3) os << R.transpose() <<
'\n';
2296 for (
const auto& R : rv.
R2) os << R.transpose() <<
'\n';
2297 for (
const auto& R : rv.
R1) os << R.transpose() <<
'\n';
2302 for (
auto& T : tm.
T4)
2303 is >>
double_imanip() >> T(0, 0) >> T(0, 1) >> T(0, 2) >> T(0, 3) >>
2304 T(1, 0) >> T(1, 1) >> T(1, 2) >> T(1, 3) >> T(2, 0) >> T(2, 1) >>
2305 T(2, 2) >> T(2, 3) >> T(3, 0) >> T(3, 1) >> T(3, 2) >> T(3, 3);
2306 for (
auto& T : tm.
T3)
2307 is >>
double_imanip() >> T(0, 0) >> T(0, 1) >> T(0, 2) >> T(1, 0) >>
2308 T(1, 1) >> T(1, 2) >> T(2, 0) >> T(2, 1) >> T(2, 2);
2309 for (
auto& T : tm.
T2)
2310 is >>
double_imanip() >> T(0, 0) >> T(0, 1) >> T(1, 0) >> T(1, 1);
2316 for (
auto& R : rv.
R4) is >>
double_imanip() >> R[0] >> R[1] >> R[2] >> R[3];
This can be used to make arrays out of anything.
Index nelem() const ARTS_NOEXCEPT
A constant view of a Matrix.
A constant view of a Tensor5.
Index ncols() const noexcept
Index npages() const noexcept
Index nbooks() const noexcept
Index nshelves() const noexcept
A constant view of a Vector.
Class to help with hidden temporary variables for operations of type Numeric times Class.
const base & bas
A const reference to a value.
const Numeric & scale
A const reference to a Numeric.
Index NumberOfFrequencies() const
The number of frequencies of the propagation matrix.
VectorView K24(const Index iz=0, const Index ia=0)
Vector view to K(1, 3) elements.
Tensor4 & Data()
Get full view to data.
VectorView K13(const Index iz=0, const Index ia=0)
Vector view to K(0, 2) elements.
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.
Index NumberOfAzimuthAngles() const
The number of azimuth angles of the propagation matrix.
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.
bool IsEmpty() const
Asks if the class is empty.
Index NumberOfZenithAngles() const
The number of zenith angles of the propagation matrix.
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.
Stokes vector is as Propagation matrix but only has 4 possible values.
Input manipulator class for doubles to enable nan and inf parsing.
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR_IF(condition,...)
Fast double input stream with support for parsing nan and inf.
constexpr std::size_t mul(Inds... inds) noexcept
NUMERIC Numeric
The type to use for all floating point numbers.
INDEX Index
The type to use for all integer numbers and indices.
Numeric & real_val(Complex &c) noexcept
Return a reference to the real value of c.
constexpr Numeric imag(Complex c) noexcept
imag
std::complex< Numeric > Complex
constexpr Numeric real(Complex c) noexcept
real
constexpr Numeric inv_sqrt_2
Inverse of the square root of 2.
Radiation Vector for Stokes dimension 1-4.
std::vector< Eigen::Vector2d > R2
RadiationVector(Index nf=0, Index stokes=1)
Construct a new Radiation Vector object.
Index Frequencies() const
Get frequency count.
const Eigen::Matrix< double, 1, 1 > & Vec1(size_t i) const
Return Vector at position.
void setSource(const StokesVector &a, const ConstVectorView &B, const StokesVector &S, Index i)
Set this to source vector at position.
void addMultiplied(const TransmissionMatrix &A, const RadiationVector &x)
Add multiply.
const Eigen::Vector4d & Vec4(size_t i) const
Return Vector at position.
void leftMul(const TransmissionMatrix &T)
Multiply radiation vector from the left.
const Numeric & operator()(const Index i, const Index j) const
Access operator.
void SetZero(size_t i)
Set Radiation Vector to Zero at position.
void rem_avg(const RadiationVector &O1, const RadiationVector &O2)
Remove the average of two other RadiationVector from *this.
std::vector< Eigen::Matrix< double, 1, 1 > > R1
void add_avg(const RadiationVector &O1, const RadiationVector &O2)
Add the average of two other RadiationVector to *this.
std::vector< Eigen::Vector4d > R4
RadiationVector & operator+=(const RadiationVector &rv)
Addition operator.
Eigen::VectorXd Vec(size_t i) const
Return Vector at position by copy.
void addDerivTransmission(const TransmissionMatrix &PiT, const TransmissionMatrix &dT, const RadiationVector &I)
Add the transmission derivative to this.
void setDerivReflection(const RadiationVector &I, const TransmissionMatrix &PiT, const TransmissionMatrix &Z, const TransmissionMatrix &dZ)
Sets *this to the reflection derivative.
const Eigen::Vector3d & Vec3(size_t i) const
Return Vector at position.
const Eigen::Vector2d & Vec2(size_t i) const
Return Vector at position.
std::vector< Eigen::Vector3d > R3
void addWeightedDerivEmission(const TransmissionMatrix &PiT, const TransmissionMatrix &dT, const TransmissionMatrix &T, const RadiationVector &I, const RadiationVector &far, const RadiationVector &close, const RadiationVector &d, bool isfar)
Add the emission derivative to this.
void setBackscatterTransmission(const RadiationVector &I0, const TransmissionMatrix &Tr, const TransmissionMatrix &Tf, const TransmissionMatrix &Z)
Set this to backscatter transmission.
void addDerivEmission(const TransmissionMatrix &PiT, const TransmissionMatrix &dT, const TransmissionMatrix &T, const RadiationVector &ImJ, const RadiationVector &dJ)
Add the emission derivative to this.
RadiationVector & operator=(const RadiationVector &rv)=default
Assign old radiation vector to this.
void add_weighted(const TransmissionMatrix &T, const RadiationVector &far, const RadiationVector &close, const ConstMatrixView &Kfar, const ConstMatrixView &Kclose, const Numeric r)
Add the weighted source of two RadiationVector to *this.
void setBackscatterTransmissionDerivative(const RadiationVector &I0, const TransmissionMatrix &Tr, const TransmissionMatrix &Tf, const TransmissionMatrix &dZ)
Set this to backscatter transmission scatter derivative.
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
void mul_aliased(const TransmissionMatrix &A, const TransmissionMatrix &B)
Set this to a multiple of A by B.
std::vector< Eigen::Matrix< double, 1, 1 > > T1
std::vector< Eigen::Matrix4d > T4
TransmissionMatrix & operator*=(const Numeric &scale)
Scale self.
const Eigen::Matrix2d & Mat2(size_t i) const
Get Matrix at position.
Eigen::MatrixXd Mat(size_t i) const
Get Matrix at position by copy.
auto & TraMat(size_t i) noexcept
Simple template access for the transmission.
Index Frequencies() const
Number of frequencies.
Eigen::Matrix< Numeric, N, 1 > second_order_integration_dsource(size_t i, const TransmissionMatrix &dx, const Eigen::Matrix< Numeric, N, 1 > far, const Eigen::Matrix< Numeric, N, 1 > close, const Eigen::Matrix< Numeric, N, 1 > d, bool isfar) const noexcept
void mul(const TransmissionMatrix &A, const TransmissionMatrix &B)
Set this to a multiple of A by B.
std::vector< Eigen::Matrix3d > T3
const Eigen::Matrix4d & Mat4(size_t i) const
Get Matrix at position.
Eigen::Matrix< Numeric, N, 1 > second_order_integration_source(const Eigen::Matrix< Numeric, N, N > T, const Eigen::Matrix< Numeric, N, 1 > far, const Eigen::Matrix< Numeric, N, 1 > close, const Eigen::Matrix< Numeric, N, N > Kfar, const Eigen::Matrix< Numeric, N, N > Kclose, const Numeric r) const noexcept
void setIdentity()
Set to identity matrix.
const Eigen::Matrix< double, 1, 1 > & Mat1(size_t i) const
Get Matrix at position.
std::vector< Eigen::Matrix2d > T2
const Eigen::Matrix3d & Mat3(size_t i) const
Get Matrix at position.
Numeric operator()(const Index i, const Index j, const Index k) const
Access value in matrix.
TransmissionMatrix & operator=(const TransmissionMatrix &tm)=default
Assignment operator.
TransmissionMatrix & operator+=(const LazyScale< TransmissionMatrix > &lstm)
Assign to *this lazily.
void setZero()
Set to zero matrix.
TransmissionMatrix(Index nf=0, Index stokes=1)
Construct a new Transmission Matrix object.
void transmat3(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r, const Index iz=0, const Index ia=0) noexcept
constexpr Numeric lower_is_considered_zero_for_sinc_likes
ArrayOfTransmissionMatrix bulk_backscatter(const ConstTensor5View &Pe, const ConstMatrixView &pnd)
Bulk back-scattering.
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
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
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
void transmat4(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r, const Index iz=0, const Index ia=0) noexcept
void transmat2(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r, const Index iz=0, const Index ia=0) noexcept
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
ArrayOfArrayOfTransmissionMatrix bulk_backscatter_derivative(const ConstTensor5View &Pe, const ArrayOfMatrix &dpnd_dx)
Derivatives of bulk back-scattering
constexpr Eigen::Matrix< Numeric, N, 1 > source_vector(const StokesVector &a, const ConstVectorView &B, const StokesVector &da, const ConstVectorView &dB_dT, const StokesVector &dS, bool dT, Index i)
LazyScale< TransmissionMatrix > operator*(const TransmissionMatrix &tm, const Numeric &x)
Lazy scale of Transmission Matrix.
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
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.
void stepwise_source(RadiationVector &J, ArrayOfRadiationVector &dJ, RadiationVector &J_add, 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.
std::istream & operator>>(std::istream &is, TransmissionMatrix &tm)
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 PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric r, const Vector &dr1, const Vector &dr2, const Index ia, const Index iz, const RadiativeTransferSolver solver)
Update the Radiation Vector.
void transmat(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r) noexcept
std::ostream & operator<<(std::ostream &os, const TransmissionMatrix &tm)
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.
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
void transmat1(TransmissionMatrix &T, const PropagationMatrix &K1, const PropagationMatrix &K2, const Numeric &r, const Index iz=0, const Index ia=0) noexcept
Stuff related to the transmission matrix.
RadiativeTransferSolver
Intended to hold various forward solvers.
BackscatterSolver
Intended to hold various backscatter solvers.
@ CommutativeTransmission
CumulativeTransmission
Intended to hold various ways to accumulate the transmission matrix.
Array< TransmissionMatrix > ArrayOfTransmissionMatrix