18 T4(stokes_dim == 4 ? nf : 0, Eigen::Matrix4d::Identity()),
19 T3(stokes_dim == 3 ? nf : 0, Eigen::Matrix3d::Identity()),
20 T2(stokes_dim == 2 ? nf : 0, Eigen::Matrix2d::Identity()),
21 T1(stokes_dim == 1 ? nf : 0, Eigen::Matrix<double, 1, 1>::Identity()) {
36 const LazyScale<TransmissionMatrix>& lstm) {
42TransmissionMatrix::operator Tensor3()
const {
43 Tensor3 T(Frequencies(), stokes_dim, stokes_dim);
44 for (
size_t i = 0; i < T4.size(); i++)
45 for (
size_t j = 0; j < 4; j++)
46 for (
size_t k = 0; k < 4; k++) T(i, j, k) = T4[i](j, k);
47 for (
size_t i = 0; i < T3.size(); i++)
48 for (
size_t j = 0; j < 3; j++)
49 for (
size_t k = 0; k < 3; k++) T(i, j, k) = T3[i](j, k);
50 for (
size_t i = 0; i < T2.size(); i++)
51 for (
size_t j = 0; j < 2; j++)
52 for (
size_t k = 0; k < 2; k++) T(i, j, k) = T2[i](j, k);
53 for (
size_t i = 0; i < T1.size(); i++) T(i, 0, 0) = T1[i](0, 0);
57#pragma GCC diagnostic push
58#pragma GCC diagnostic ignored "-Wreturn-type"
71#pragma GCC diagnostic pop
74 std::fill(
T4.begin(),
T4.end(), Eigen::Matrix4d::Identity());
75 std::fill(
T3.begin(),
T3.end(), Eigen::Matrix3d::Identity());
76 std::fill(
T2.begin(),
T2.end(), Eigen::Matrix2d::Identity());
77 std::fill(
T1.begin(),
T1.end(), Eigen::Matrix<double, 1, 1>::Identity());
81 std::fill(
T4.begin(),
T4.end(), Eigen::Matrix4d::Zero());
82 std::fill(
T3.begin(),
T3.end(), Eigen::Matrix3d::Zero());
83 std::fill(
T2.begin(),
T2.end(), Eigen::Matrix2d::Zero());
84 std::fill(
T1.begin(),
T1.end(), Eigen::Matrix<double, 1, 1>::Zero());
89 for (
size_t i = 0; i <
T4.size(); i++)
T4[i].noalias() = A.
T4[i] * B.
T4[i];
90 for (
size_t i = 0; i <
T3.size(); i++)
T3[i].noalias() = A.
T3[i] * B.
T3[i];
91 for (
size_t i = 0; i <
T2.size(); i++)
T2[i].noalias() = A.
T2[i] * B.
T2[i];
92 for (
size_t i = 0; i <
T1.size(); i++)
T1[i].noalias() = A.
T1[i] * B.
T1[i];
97 for (
size_t i = 0; i <
T4.size(); i++)
T4[i] = A.
T4[i] * B.
T4[i];
98 for (
size_t i = 0; i <
T3.size(); i++)
T3[i] = A.
T3[i] * B.
T3[i];
99 for (
size_t i = 0; i <
T2.size(); i++)
T2[i] = A.
T2[i] * B.
T2[i];
100 for (
size_t i = 0; i <
T1.size(); i++)
T1[i] = A.
T1[i] * B.
T1[i];
105 const Index k)
const {
134 return Index(
T4.size());
136 return Index(
T3.size());
138 return Index(
T2.size());
140 return Index(
T1.size());
148 operator()(0, i, j) = mat(i,j);
152 const LazyScale<TransmissionMatrix>& lstm) {
153 for (
size_t i = 0; i <
T4.size(); i++)
154 T4[i].noalias() = lstm.scale * lstm.bas.Mat4(i);
155 for (
size_t i = 0; i <
T3.size(); i++)
156 T3[i].noalias() = lstm.scale * lstm.bas.Mat3(i);
157 for (
size_t i = 0; i <
T2.size(); i++)
158 T2[i].noalias() = lstm.scale * lstm.bas.Mat2(i);
159 for (
size_t i = 0; i <
T1.size(); i++)
160 T1[i].noalias() = lstm.scale * lstm.bas.Mat1(i);
166 T4.begin(),
T4.end(),
T4.begin(), [scale](
auto& T) { return T * scale; });
168 T3.begin(),
T3.end(),
T3.begin(), [scale](
auto& T) { return T * scale; });
170 T2.begin(),
T2.end(),
T2.begin(), [scale](
auto& T) { return T * scale; });
172 T1.begin(),
T1.end(),
T1.begin(), [scale](
auto& T) { return T * scale; });
181LazyScale<TransmissionMatrix>
operator*(
const Numeric& x,
187 : stokes_dim(stokes),
188 R4(stokes_dim == 4 ? nf : 0, Eigen::Vector4d::Zero()),
189 R3(stokes_dim == 3 ? nf : 0, Eigen::Vector3d::Zero()),
190 R2(stokes_dim == 2 ? nf : 0, Eigen::Vector2d::Zero()),
191 R1(stokes_dim == 1 ? nf : 0, Eigen::Matrix<double, 1, 1>::Zero()) {
221 for (
size_t i = 0; i <
R4.size(); i++)
R4[i] = T.
Mat4(i) *
R4[i];
222 for (
size_t i = 0; i <
R3.size(); i++)
R3[i] = T.
Mat3(i) *
R3[i];
223 for (
size_t i = 0; i <
R2.size(); i++)
R2[i] = T.
Mat2(i) *
R2[i];
224 for (
size_t i = 0; i <
R1.size(); i++)
R1[i] = T.
Mat1(i) *
R1[i];
230 R4[i].noalias() = Eigen::Vector4d::Zero();
233 R3[i].noalias() = Eigen::Vector3d::Zero();
236 R2[i].noalias() = Eigen::Vector2d::Zero();
263 for (
size_t i = 0; i <
R4.size(); i++)
264 R4[i].noalias() -= 0.5 * (O1.
R4[i] + O2.
R4[i]);
265 for (
size_t i = 0; i <
R3.size(); i++)
266 R3[i].noalias() -= 0.5 * (O1.
R3[i] + O2.
R3[i]);
267 for (
size_t i = 0; i <
R2.size(); i++)
268 R2[i].noalias() -= 0.5 * (O1.
R2[i] + O2.
R2[i]);
269 for (
size_t i = 0; i <
R1.size(); i++)
270 R1[i].noalias() -= 0.5 * (O1.
R1[i] + O2.
R1[i]);
275 for (
size_t i = 0; i <
R4.size(); i++)
276 R4[i].noalias() += 0.5 * (O1.
R4[i] + O2.
R4[i]);
277 for (
size_t i = 0; i <
R3.size(); i++)
278 R3[i].noalias() += 0.5 * (O1.
R3[i] + O2.
R3[i]);
279 for (
size_t i = 0; i <
R2.size(); i++)
280 R2[i].noalias() += 0.5 * (O1.
R2[i] + O2.
R2[i]);
281 for (
size_t i = 0; i <
R1.size(); i++)
282 R1[i].noalias() += 0.5 * (O1.
R1[i] + O2.
R1[i]);
288 const ConstMatrixView& Kfar,
289 const ConstMatrixView& Kclose,
291 for (
size_t i = 0; i <
R4.size(); i++) {
296 prop_matrix<4>(Kfar(i, joker)),
297 prop_matrix<4>(Kclose(i, joker)),
300 for (
size_t i = 0; i <
R3.size(); i++) {
305 prop_matrix<3>(Kfar(i, joker)),
306 prop_matrix<3>(Kclose(i, joker)),
309 for (
size_t i = 0; i <
R2.size(); i++) {
314 prop_matrix<2>(Kfar(i, joker)),
315 prop_matrix<2>(Kclose(i, joker)),
318 for (
size_t i = 0; i <
R1.size(); i++) {
323 prop_matrix<1>(Kfar(i, joker)),
324 prop_matrix<1>(Kclose(i, joker)),
334 for (
size_t i = 0; i <
R4.size(); i++)
335 R4[i].noalias() += PiT.
Mat4(i) * (dT.
Mat4(i) * ImJ.
R4[i] + dJ.
R4[i] -
337 for (
size_t i = 0; i <
R3.size(); i++)
338 R3[i].noalias() += PiT.
Mat3(i) * (dT.
Mat3(i) * ImJ.
R3[i] + dJ.
R3[i] -
340 for (
size_t i = 0; i <
R2.size(); i++)
341 R2[i].noalias() += PiT.
Mat2(i) * (dT.
Mat2(i) * ImJ.
R2[i] + dJ.
R2[i] -
343 for (
size_t i = 0; i <
R1.size(); i++)
344 R1[i].noalias() += PiT.
Mat1(i) * (dT.
Mat1(i) * ImJ.
R1[i] + dJ.
R1[i] -
356 for (
size_t i = 0; i <
R4.size(); i++)
360 i, dT, far.
R4[i], close.
R4[i],
d.R4[i], isfar));
361 for (
size_t i = 0; i <
R3.size(); i++)
365 i, dT, far.
R3[i], close.
R3[i],
d.R3[i], isfar));
366 for (
size_t i = 0; i <
R2.size(); i++)
370 i, dT, far.
R2[i], close.
R2[i],
d.R2[i], isfar));
371 for (
size_t i = 0; i <
R1.size(); i++)
375 i, dT, far.
R1[i], close.
R1[i],
d.R1[i], isfar));
381 for (
size_t i = 0; i <
R4.size(); i++)
383 for (
size_t i = 0; i <
R3.size(); i++)
385 for (
size_t i = 0; i <
R2.size(); i++)
387 for (
size_t i = 0; i <
R1.size(); i++)
393 for (
size_t i = 0; i <
R4.size(); i++)
R4[i].noalias() += A.
Mat4(i) * x.
R4[i];
394 for (
size_t i = 0; i <
R3.size(); i++)
R3[i].noalias() += A.
Mat3(i) * x.
R3[i];
395 for (
size_t i = 0; i <
R2.size(); i++)
R2[i].noalias() += A.
Mat2(i) * x.
R2[i];
396 for (
size_t i = 0; i <
R1.size(); i++)
R1[i].noalias() += A.
Mat1(i) * x.
R1[i];
403 for (
size_t i = 0; i <
R4.size(); i++)
405 for (
size_t i = 0; i <
R3.size(); i++)
407 for (
size_t i = 0; i <
R2.size(); i++)
409 for (
size_t i = 0; i <
R1.size(); i++)
418 for (
size_t i = 0; i <
R4.size(); i++)
420 for (
size_t i = 0; i <
R3.size(); i++)
422 for (
size_t i = 0; i <
R2.size(); i++)
424 for (
size_t i = 0; i <
R1.size(); i++)
433 for (
size_t i = 0; i <
R4.size(); i++)
435 for (
size_t i = 0; i <
R3.size(); i++)
437 for (
size_t i = 0; i <
R2.size(); i++)
439 for (
size_t i = 0; i <
R1.size(); i++)
445 for (
size_t i = 0; i <
R4.size(); i++) {
451 for (
size_t i = 0; i <
R3.size(); i++) {
456 for (
size_t i = 0; i <
R2.size(); i++) {
460 for (
size_t i = 0; i <
R1.size(); i++) {
467 for (
size_t i = 0; i <
R4.size(); i++)
468 R4[i].noalias() += rv.
R4[i];
469 for (
size_t i = 0; i <
R3.size(); i++)
470 R3[i].noalias() += rv.
R3[i];
471 for (
size_t i = 0; i <
R2.size(); i++)
472 R2[i].noalias() += rv.
R2[i];
473 for (
size_t i = 0; i <
R1.size(); i++)
474 R1[i].noalias() += rv.
R1[i];
493RadiationVector::operator Matrix()
const {
494 Matrix M(Frequencies(), stokes_dim);
495 for (
size_t i = 0; i < R4.size(); i++)
496 for (
size_t j = 0; j < 4; j++) M(i, j) = R4[i](j);
497 for (
size_t i = 0; i < R3.size(); i++)
498 for (
size_t j = 0; j < 3; j++) M(i, j) = R3[i](j);
499 for (
size_t i = 0; i < R2.size(); i++)
500 for (
size_t j = 0; j < 2; j++) M(i, j) = R2[i](j);
501 for (
size_t i = 0; i < R1.size(); i++) M(i, 0) = R1[i](0);
506 const ConstVectorView& B,
507 const StokesVector& S,
517 Eigen::Vector4d(
a.Kjj()[i],
a.K12()[i],
a.K13()[i],
a.K14()[i]) *
519 Eigen::Vector4d(S.Kjj()[i], S.K12()[i], S.K13()[i], S.K14()[i]);
522 Eigen::Vector4d(
a.Kjj()[i],
a.K12()[i],
a.K13()[i],
a.K14()[i]) *
528 Eigen::Vector3d(
a.Kjj()[i],
a.K12()[i],
a.K13()[i]) * B[i] +
529 Eigen::Vector3d(S.Kjj()[i], S.K12()[i], S.K13()[i]);
532 Eigen::Vector3d(
a.Kjj()[i],
a.K12()[i],
a.K13()[i]) * B[i];
536 R2[i].noalias() = Eigen::Vector2d(
a.Kjj()[i],
a.K12()[i]) * B[i] +
537 Eigen::Vector2d(S.Kjj()[i], S.K12()[i]);
539 R2[i].noalias() = Eigen::Vector2d(
a.Kjj()[i],
a.K12()[i]) * B[i];
543 R1[i][0] =
a.Kjj()[i] * B[i] + S.Kjj()[i];
545 R1[i][0] =
a.Kjj()[i] * B[i];
552 return Index(
R4.size());
554 return Index(
R3.size());
556 return Index(
R2.size());
558 return Index(
R1.size());
566 const StokesVector&
a,
567 const ConstVectorView& B,
568 const StokesVector& da,
569 const ConstVectorView& dB_dT,
570 const StokesVector& dS,
573 static_assert(N > 0 and N < 5,
"Bad stokes dimensions");
575 if constexpr (N == 1) {
577 return Eigen::Matrix<Numeric, 1, 1>(dS.Kjj()[i] + da.Kjj()[i] * B[i] +
578 a.Kjj()[i] * dB_dT[i]);
579 return Eigen::Matrix<Numeric, 1, 1>(dS.Kjj()[i] + da.Kjj()[i] * B[i]);
582 if constexpr (N == 2) {
584 return Eigen::Vector2d(
585 dS.Kjj()[i] + da.Kjj()[i] * B[i] +
a.Kjj()[i] * dB_dT[i],
586 dS.K12()[i] + da.K12()[i] * B[i] +
a.K12()[i] * dB_dT[i]);
587 return Eigen::Vector2d(dS.Kjj()[i] + da.Kjj()[i] * B[i],
588 dS.K12()[i] + da.K12()[i] * B[i]);
591 if constexpr (N == 3) {
593 return Eigen::Vector3d(
594 dS.Kjj()[i] + da.Kjj()[i] * B[i] +
a.Kjj()[i] * dB_dT[i],
595 dS.K12()[i] + da.K12()[i] * B[i] +
a.K12()[i] * dB_dT[i],
596 dS.K13()[i] + da.K13()[i] * B[i] +
a.K13()[i] * dB_dT[i]);
597 return Eigen::Vector3d(dS.Kjj()[i] + da.Kjj()[i] * B[i],
598 dS.K12()[i] + da.K12()[i] * B[i],
599 dS.K13()[i] + da.K13()[i] * B[i]);
602 if constexpr (N == 4) {
604 return Eigen::Vector4d(
605 dS.Kjj()[i] + da.Kjj()[i] * B[i] +
a.Kjj()[i] * dB_dT[i],
606 dS.K12()[i] + da.K12()[i] * B[i] +
a.K12()[i] * dB_dT[i],
607 dS.K13()[i] + da.K13()[i] * B[i] +
a.K13()[i] * dB_dT[i],
608 dS.K14()[i] + da.K14()[i] * B[i] +
a.K14()[i] * dB_dT[i]);
609 return Eigen::Vector4d(dS.Kjj()[i] + da.Kjj()[i] * B[i],
610 dS.K12()[i] + da.K12()[i] * B[i],
611 dS.K13()[i] + da.K13()[i] * B[i],
612 dS.K14()[i] + da.K14()[i] * B[i]);
618 const StokesVector&
a,
619 const ConstVectorView& B,
620 const StokesVector& da,
621 const ConstVectorView& dB_dT,
624 static_assert(N > 0 and N < 5,
"Bad stokes dimensions");
626 if constexpr (N == 1) {
628 return Eigen::Matrix<Numeric, 1, 1>(da.Kjj()[i] * B[i] +
629 a.Kjj()[i] * dB_dT[i]);
630 return Eigen::Matrix<Numeric, 1, 1>(da.Kjj()[i] * B[i]);
633 if constexpr (N == 2) {
635 return Eigen::Vector2d(da.Kjj()[i] * B[i] +
a.Kjj()[i] * dB_dT[i],
636 da.K12()[i] * B[i] +
a.K12()[i] * dB_dT[i]);
637 return Eigen::Vector2d(da.Kjj()[i] * B[i], da.K12()[i] * B[i]);
640 if constexpr (N == 3) {
642 return Eigen::Vector3d(da.Kjj()[i] * B[i] +
a.Kjj()[i] * dB_dT[i],
643 da.K12()[i] * B[i] +
a.K12()[i] * dB_dT[i],
644 da.K13()[i] * B[i] +
a.K13()[i] * dB_dT[i]);
645 return Eigen::Vector3d(
646 da.Kjj()[i] * B[i], da.K12()[i] * B[i], da.K13()[i] * B[i]);
649 if constexpr (N == 4) {
651 return Eigen::Vector4d(da.Kjj()[i] * B[i] +
a.Kjj()[i] * dB_dT[i],
652 da.K12()[i] * B[i] +
a.K12()[i] * dB_dT[i],
653 da.K13()[i] * B[i] +
a.K13()[i] * dB_dT[i],
654 da.K14()[i] * B[i] +
a.K14()[i] * dB_dT[i]);
655 return Eigen::Vector4d(da.Kjj()[i] * B[i],
663 const PropagationMatrix& K1,
664 const PropagationMatrix& K2,
667 const Index ia = 0) noexcept {
668 for (Index i = 0; i < K1.NumberOfFrequencies(); i++)
670 std::exp(-0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]));
674 const PropagationMatrix& K1,
675 const PropagationMatrix& K2,
678 const Index ia = 0) noexcept {
679 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
680 const Numeric
a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
681 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]);
682 const Numeric exp_a = std::exp(
a);
683 const Numeric cb = std::cosh(
b), sb = std::sinh(
b);
684 T.
Mat2(i).noalias() =
685 (Eigen::Matrix2d() << cb, sb, sb, cb).finished() * exp_a;
690 const PropagationMatrix& K1,
691 const PropagationMatrix& K2,
694 const Index ia = 0) noexcept {
695 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
696 const Numeric
a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
697 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]),
698 c = -0.5 * r * (K1.K13(iz, ia)[i] + K2.K13(iz, ia)[i]),
699 u = -0.5 * r * (K1.K23(iz, ia)[i] + K2.K23(iz, ia)[i]);
700 const Numeric exp_a = std::exp(
a);
702 if (
b == 0. and
c == 0. and
u == 0.) {
703 T.
Mat3(i).noalias() = Eigen::Matrix3d::Identity() * exp_a;
705 const Numeric a2 =
a *
a, b2 =
b *
b, c2 =
c *
c, u2 =
u *
u;
706 const Numeric Const = b2 + c2 - u2;
708 const bool real = Const > 0;
709 const bool imag = Const < 0;
710 const bool either = real or imag;
713 std::sqrt(imag ? -Const : Const);
715 (real ? 1 : -1) * x * x;
716 const Numeric inv_x2 =
721 real ? std::sinh(x) : std::sin(-x);
723 real ? std::cosh(x) : std::cos(+x);
732 either ? a2 * (cx - 1.0) -
a * x * sx + x2 : 1.0 + 0.5 * a2 -
a;
733 const Numeric C1 = either ? 2.0 *
a * (1.0 - cx) + x * sx : 1.0 -
a;
734 const Numeric C2 = either ? cx - 1.0 : 0.5;
736 T.
Mat3(i).noalias() =
738 (Eigen::Matrix3d() << C0 + C1 *
a + C2 * (a2 + b2 + c2),
739 C1 *
b + C2 * (2 *
a *
b -
c *
u),
740 C1 *
c + C2 * (2 *
a *
c +
b *
u),
741 C1 *
b + C2 * (2 *
a *
b +
c *
u),
742 C0 + C1 *
a + C2 * (a2 + b2 - u2),
743 C1 *
u + C2 * (2 *
a *
u +
b *
c),
744 C1 *
c + C2 * (2 *
a *
c -
b *
u),
745 -C1 *
u - C2 * (2 *
a *
u -
b *
c),
746 C0 + C1 *
a + C2 * (a2 + c2 - u2))
753 const PropagationMatrix& K1,
754 const PropagationMatrix& K2,
757 const Index ia = 0) noexcept {
759 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
760 const Numeric
a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
761 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]),
762 c = -0.5 * r * (K1.K13(iz, ia)[i] + K2.K13(iz, ia)[i]),
763 d = -0.5 * r * (K1.K14(iz, ia)[i] + K2.K14(iz, ia)[i]),
764 u = -0.5 * r * (K1.K23(iz, ia)[i] + K2.K23(iz, ia)[i]),
765 v = -0.5 * r * (K1.K24(iz, ia)[i] + K2.K24(iz, ia)[i]),
766 w = -0.5 * r * (K1.K34(iz, ia)[i] + K2.K34(iz, ia)[i]);
767 const Numeric exp_a = std::exp(
a);
769 if (
b == 0. and
c == 0. and
d == 0. and
u == 0. and
v == 0. and
w == 0.)
770 T.
Mat4(i).noalias() = Eigen::Matrix4d::Identity() * exp_a;
772 const Numeric b2 =
b *
b, c2 =
c *
c, d2 =
d *
d, u2 =
u *
u, v2 =
v *
v,
776 w2 * w2 + 2 * (b2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
777 c2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
778 d2 * (d2 * 0.5 + u2 - v2 - w2) +
779 u2 * (u2 * 0.5 + v2 + w2) + v2 * (v2 * 0.5 + w2) +
781 const Complex Const1 = std::sqrt(Complex(tmp, 0));
782 const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
784 const Complex x = std::sqrt(Const2 + Const1) * sqrt_05;
785 const Complex y = std::sqrt(Const2 - Const1) * sqrt_05 * Complex(0, 1);
786 const Complex x2 = x * x;
787 const Complex y2 = y * y;
788 const Complex cy = std::cos(y);
789 const Complex sy = std::sin(y);
790 const Complex cx = std::cosh(x);
791 const Complex sx = std::sinh(x);
795 const bool both_zero = y_zero and x_zero;
796 const bool either_zero = y_zero or x_zero;
806 const Complex ix = x_zero ? 0.0 : 1.0 / x;
807 const Complex iy = y_zero ? 0.0 : 1.0 / y;
808 const Complex inv_x2y2 =
815 either_zero ? 1.0 : ((cy * x2 + cx * y2) * inv_x2y2).real();
817 either_zero ? 1.0 : ((sy * x2 * iy + sx * y2 * ix) * inv_x2y2).real();
818 const Numeric C2 = both_zero ? 0.5 : ((cx - cy) * inv_x2y2).real();
819 const Numeric C3 = both_zero ? 1.0 / 6.0
820 : ((x_zero ? 1.0 - sy * iy
821 : y_zero ? sx * ix - 1.0
822 : sx * ix - sy * iy) *
825 T.
Mat4(i).noalias() =
826 exp_a * (Eigen::Matrix4d() << C0 + C2 * (b2 + c2 + d2),
827 C1 *
b + C2 * (-
c *
u -
d *
v) +
828 C3 * (
b * (b2 + c2 + d2) -
u * (
b *
u -
d *
w) -
829 v * (
b *
v +
c *
w)),
830 C1 *
c + C2 * (
b *
u -
d *
w) +
831 C3 * (
c * (b2 + c2 + d2) -
u * (
c *
u +
d *
v) -
832 w * (
b *
v +
c *
w)),
833 C1 *
d + C2 * (
b *
v +
c *
w) +
834 C3 * (
d * (b2 + c2 + d2) -
v * (
c *
u +
d *
v) +
835 w * (
b *
u -
d *
w)),
837 C1 *
b + C2 * (
c *
u +
d *
v) +
838 C3 * (-
b * (-b2 + u2 + v2) +
c * (
b *
c -
v *
w) +
839 d * (
b *
d +
u *
w)),
840 C0 + C2 * (b2 - u2 - v2),
841 C2 * (
b *
c -
v *
w) + C1 *
u +
842 C3 * (
c * (
c *
u +
d *
v) -
u * (-b2 + u2 + v2) -
843 w * (
b *
d +
u *
w)),
844 C2 * (
b *
d +
u *
w) + C1 *
v +
845 C3 * (
d * (
c *
u +
d *
v) -
v * (-b2 + u2 + v2) +
846 w * (
b *
c -
v *
w)),
848 C1 *
c + C2 * (-
b *
u +
d *
w) +
849 C3 * (
b * (
b *
c -
v *
w) -
c * (-c2 + u2 + w2) +
850 d * (
c *
d -
u *
v)),
851 C2 * (
b *
c -
v *
w) - C1 *
u +
852 C3 * (-
b * (
b *
u -
d *
w) +
u * (-c2 + u2 + w2) -
853 v * (
c *
d -
u *
v)),
854 C0 + C2 * (c2 - u2 - w2),
855 C2 * (
c *
d -
u *
v) + C1 *
w +
856 C3 * (-
d * (
b *
u -
d *
w) +
v * (
b *
c -
v *
w) -
857 w * (-c2 + u2 + w2)),
859 C1 *
d + C2 * (-
b *
v -
c *
w) +
860 C3 * (
b * (
b *
d +
u *
w) +
c * (
c *
d -
u *
v) -
861 d * (-d2 + v2 + w2)),
862 C2 * (
b *
d +
u *
w) - C1 *
v +
863 C3 * (-
b * (
b *
v +
c *
w) -
u * (
c *
d -
u *
v) +
864 v * (-d2 + v2 + w2)),
865 C2 * (
c *
d -
u *
v) - C1 *
w +
866 C3 * (-
c * (
b *
v +
c *
w) +
u * (
b *
d +
u *
w) +
867 w * (-d2 + v2 + w2)),
868 C0 + C2 * (d2 - v2 - w2))
877 const PropagationMatrix& K1,
878 const PropagationMatrix& K2,
879 const ArrayOfPropagationMatrix& dK1,
880 const ArrayOfPropagationMatrix& dK2,
882 const Numeric& dr_dT1,
883 const Numeric& dr_dT2,
886 const Index ia)
noexcept {
887 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
889 std::exp(-0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]));
890 for (Index j = 0; j < dT1.nelem(); j++) {
891 if (dK1[j].NumberOfFrequencies())
892 dT1[j].Mat1(i)(0, 0) =
895 (r * dK1[j].Kjj(iz, ia)[i] +
896 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
898 if (dK2[j].NumberOfFrequencies())
899 dT2[j].Mat1(i)(0, 0) =
902 (r * dK2[j].Kjj(iz, ia)[i] +
903 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
912 const PropagationMatrix& K1,
913 const PropagationMatrix& K2,
914 const ArrayOfPropagationMatrix& dK1,
915 const ArrayOfPropagationMatrix& dK2,
917 const Numeric& dr_dT1,
918 const Numeric& dr_dT2,
921 const Index ia)
noexcept {
922 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
923 const Numeric
a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
924 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]);
925 const Numeric exp_a = std::exp(
a);
926 const Numeric cb = std::cosh(
b), sb = std::sinh(
b);
927 T.Mat2(i).noalias() =
928 (Eigen::Matrix2d() << cb, sb, sb, cb).finished() * exp_a;
929 for (Index j = 0; j < dT1.nelem(); j++) {
930 if (dK1[j].NumberOfFrequencies()) {
931 const Numeric da = -0.5 * (r * dK1[j].Kjj(iz, ia)[i] +
932 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] +
935 db = -0.5 * (r * dK1[j].K12(iz, ia)[i] +
936 ((j == it) ? dr_dT1 * (K1.K12(iz, ia)[i] +
939 dT1[j].Mat2(i).noalias() =
941 (Eigen::Matrix2d() << sb, cb, cb, sb).finished() * exp_a * db;
943 if (dK2[j].NumberOfFrequencies()) {
944 const Numeric da = -0.5 * (r * dK2[j].Kjj(iz, ia)[i] +
945 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] +
948 db = -0.5 * (r * dK2[j].K12(iz, ia)[i] +
949 ((j == it) ? dr_dT2 * (K1.K12(iz, ia)[i] +
952 dT2[j].Mat2(i).noalias() =
954 (Eigen::Matrix2d() << sb, cb, cb, sb).finished() * exp_a * db;
963 const PropagationMatrix& K1,
964 const PropagationMatrix& K2,
965 const ArrayOfPropagationMatrix& dK1,
966 const ArrayOfPropagationMatrix& dK2,
968 const Numeric& dr_dT1,
969 const Numeric& dr_dT2,
972 const Index ia)
noexcept {
973 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
974 const Numeric
a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
975 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]),
976 c = -0.5 * r * (K1.K13(iz, ia)[i] + K2.K13(iz, ia)[i]),
977 u = -0.5 * r * (K1.K23(iz, ia)[i] + K2.K23(iz, ia)[i]);
978 const Numeric exp_a = std::exp(
a);
980 if (
b == 0. and
c == 0. and
u == 0.) {
981 T.Mat3(i).noalias() = Eigen::Matrix3d::Identity() * exp_a;
982 for (Index j = 0; j < dT1.nelem(); j++) {
983 if (dK1[j].NumberOfFrequencies())
984 dT1[j].Mat3(i).noalias() =
987 (r * dK1[j].Kjj(iz, ia)[i] +
988 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
990 if (dK2[j].NumberOfFrequencies())
991 dT2[j].Mat3(i).noalias() =
994 (r * dK2[j].Kjj(iz, ia)[i] +
995 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
999 const Numeric a2 =
a *
a, b2 =
b *
b, c2 =
c *
c, u2 =
u *
u;
1000 const Numeric Const = b2 + c2 - u2;
1002 const bool real = Const > 0;
1003 const bool imag = Const < 0;
1004 const bool either = real or imag;
1007 std::sqrt(imag ? -Const : Const);
1009 (real ? 1 : -1) * x * x;
1010 const Numeric inv_x =
1013 const Numeric inv_x2 = inv_x * inv_x;
1016 real ? std::sinh(x) : std::sin(-x);
1018 real ? std::cosh(x) : std::cos(+x);
1027 either ? a2 * (cx - 1.0) -
a * x * sx + x2 : 1.0 + 0.5 * a2 -
a;
1028 const Numeric C1 = either ? 2.0 *
a * (1.0 - cx) + x * sx : 1.0 -
a;
1029 const Numeric C2 = either ? cx - 1.0 : 0.5;
1031 T.Mat3(i).noalias() =
1033 (Eigen::Matrix3d() << C0 + C1 *
a + C2 * (a2 + b2 + c2),
1034 C1 *
b + C2 * (2 *
a *
b -
c *
u),
1035 C1 *
c + C2 * (2 *
a *
c +
b *
u),
1036 C1 *
b + C2 * (2 *
a *
b +
c *
u),
1037 C0 + C1 *
a + C2 * (a2 + b2 - u2),
1038 C1 *
u + C2 * (2 *
a *
u +
b *
c),
1039 C1 *
c + C2 * (2 *
a *
c -
b *
u),
1040 -C1 *
u - C2 * (2 *
a *
u -
b *
c),
1041 C0 + C1 *
a + C2 * (a2 + c2 - u2))
1044 for (Index j = 0; j < dT1.nelem(); j++) {
1045 if (dK1[j].NumberOfFrequencies()) {
1046 const Numeric da = -0.5 * (r * dK1[j].Kjj(iz, ia)[i] +
1047 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] +
1050 db = -0.5 * (r * dK1[j].K12(iz, ia)[i] +
1051 ((j == it) ? dr_dT1 * (K1.K12(iz, ia)[i] +
1054 dc = -0.5 * (r * dK1[j].K13(iz, ia)[i] +
1055 ((j == it) ? dr_dT1 * (K1.K13(iz, ia)[i] +
1058 du = -0.5 * (r * dK1[j].K23(iz, ia)[i] +
1059 ((j == it) ? dr_dT1 * (K1.K23(iz, ia)[i] +
1062 const Numeric da2 = 2 *
a * da, db2 = 2 *
b * db, dc2 = 2 *
c * dc,
1064 const Numeric dx = either ? ((db2 + dc2 - du2) * inv_x * 0.5) : 0,
1066 const Numeric dsx = (real ? 1 : -1) * cx * dx;
1067 const Numeric dcx = sx * dx;
1069 const Numeric dC0 = either ? da2 * (cx - 1.0) + da2 * (dcx - 1.0) -
1070 da * x * sx -
a * dx * sx -
1074 either ? 2.0 * (da * (1.0 - cx) -
a * dcx) + dx * sx + x * dsx
1076 const Numeric dC2 = either ? dcx : 0;
1078 dT1[j].Mat3(i).noalias() =
1079 T.Mat3(i) * (da + dx2 * inv_x2) +
1081 (Eigen::Matrix3d() << dC0 + dC1 *
a + C1 * da +
1082 dC2 * (a2 + b2 + c2) +
1083 C2 * (da2 + db2 + dc2),
1084 dC1 *
b + C1 * db + dC2 * (2 *
a *
b -
c *
u) +
1085 C2 * (2 * da *
b - dc *
u + 2 *
a * db -
c * du),
1086 dC1 *
c + C1 * dc + dC2 * (2 *
a *
c +
b *
u) +
1087 C2 * (2 * da *
c + db *
u + 2 *
a * dc +
b * du),
1088 dC1 *
b + C1 * db + dC2 * (2 *
a *
b +
c *
u) +
1089 C2 * (2 * da *
b + dc *
u + 2 *
a * db +
c * du),
1090 dC0 + dC1 *
a + C1 * da + dC2 * (a2 + b2 - u2) +
1091 C2 * (da2 + db2 - du2),
1092 dC1 *
u + C1 * du + dC2 * (2 *
a *
u +
b *
c) +
1093 C2 * (2 * da *
u + db *
c + 2 *
a * du +
b * dc),
1094 dC1 *
c + C1 * dc + dC2 * (2 *
a *
c -
b *
u) +
1095 C2 * (2 * da *
c - db *
u + 2 *
a * dc -
b * du),
1096 -dC1 *
u - C1 * du - dC2 * (2 *
a *
u -
b *
c) -
1097 C2 * (2 * da *
u - db *
c + 2 *
a * du -
b * dc),
1098 dC0 + dC1 *
a + C1 * da + dC2 * (a2 + c2 - u2) +
1099 C2 * (da2 + dc2 - du2))
1102 if (dK2[j].NumberOfFrequencies()) {
1103 const Numeric da = -0.5 * (r * dK2[j].Kjj(iz, ia)[i] +
1104 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] +
1107 db = -0.5 * (r * dK2[j].K12(iz, ia)[i] +
1108 ((j == it) ? dr_dT2 * (K1.K12(iz, ia)[i] +
1111 dc = -0.5 * (r * dK2[j].K13(iz, ia)[i] +
1112 ((j == it) ? dr_dT2 * (K1.K13(iz, ia)[i] +
1115 du = -0.5 * (r * dK2[j].K23(iz, ia)[i] +
1116 ((j == it) ? dr_dT2 * (K1.K23(iz, ia)[i] +
1119 const Numeric da2 = 2 *
a * da, db2 = 2 *
b * db, dc2 = 2 *
c * dc,
1121 const Numeric dx = either ? ((db2 + dc2 - du2) * inv_x * 0.5) : 0,
1123 const Numeric dsx = (real ? 1 : -1) * cx * dx;
1124 const Numeric dcx = sx * dx;
1126 const Numeric dC0 = either ? da2 * (cx - 1.0) + da2 * (dcx - 1.0) -
1127 da * x * sx -
a * dx * sx -
1131 either ? 2.0 * (da * (1.0 - cx) -
a * dcx) + dx * sx + x * dsx
1133 const Numeric dC2 = either ? dcx : 0;
1135 dT2[j].Mat3(i).noalias() =
1136 T.Mat3(i) * (da + dx2 * inv_x2) +
1138 (Eigen::Matrix3d() << dC0 + dC1 *
a + C1 * da +
1139 dC2 * (a2 + b2 + c2) +
1140 C2 * (da2 + db2 + dc2),
1141 dC1 *
b + C1 * db + dC2 * (2 *
a *
b -
c *
u) +
1142 C2 * (2 * da *
b - dc *
u + 2 *
a * db -
c * du),
1143 dC1 *
c + C1 * dc + dC2 * (2 *
a *
c +
b *
u) +
1144 C2 * (2 * da *
c + db *
u + 2 *
a * dc +
b * du),
1145 dC1 *
b + C1 * db + dC2 * (2 *
a *
b +
c *
u) +
1146 C2 * (2 * da *
b + dc *
u + 2 *
a * db +
c * du),
1147 dC0 + dC1 *
a + C1 * da + dC2 * (a2 + b2 - u2) +
1148 C2 * (da2 + db2 - du2),
1149 dC1 *
u + C1 * du + dC2 * (2 *
a *
u +
b *
c) +
1150 C2 * (2 * da *
u + db *
c + 2 *
a * du +
b * dc),
1151 dC1 *
c + C1 * dc + dC2 * (2 *
a *
c -
b *
u) +
1152 C2 * (2 * da *
c - db *
u + 2 *
a * dc -
b * du),
1153 -dC1 *
u - C1 * du - dC2 * (2 *
a *
u -
b *
c) -
1154 C2 * (2 * da *
u - db *
c + 2 *
a * du -
b * dc),
1155 dC0 + dC1 *
a + C1 * da + dC2 * (a2 + c2 - u2) +
1156 C2 * (da2 + dc2 - du2))
1167 const PropagationMatrix& K1,
1168 const PropagationMatrix& K2,
1169 const ArrayOfPropagationMatrix& dK1,
1170 const ArrayOfPropagationMatrix& dK2,
1172 const Numeric& dr_dT1,
1173 const Numeric& dr_dT2,
1176 const Index ia)
noexcept {
1178 for (Index i = 0; i < K1.NumberOfFrequencies(); i++) {
1179 const Numeric
a = -0.5 * r * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i]),
1180 b = -0.5 * r * (K1.K12(iz, ia)[i] + K2.K12(iz, ia)[i]),
1181 c = -0.5 * r * (K1.K13(iz, ia)[i] + K2.K13(iz, ia)[i]),
1182 d = -0.5 * r * (K1.K14(iz, ia)[i] + K2.K14(iz, ia)[i]),
1183 u = -0.5 * r * (K1.K23(iz, ia)[i] + K2.K23(iz, ia)[i]),
1184 v = -0.5 * r * (K1.K24(iz, ia)[i] + K2.K24(iz, ia)[i]),
1185 w = -0.5 * r * (K1.K34(iz, ia)[i] + K2.K34(iz, ia)[i]);
1186 const Numeric exp_a = std::exp(
a);
1188 if (
b == 0. and
c == 0. and
d == 0. and
u == 0. and
v == 0. and
w == 0.) {
1189 T.Mat4(i).noalias() = Eigen::Matrix4d::Identity() * exp_a;
1190 for (Index j = 0; j < dK1.nelem(); j++) {
1191 if (dK1[j].NumberOfFrequencies())
1192 dT1[j].Mat4(i).noalias() =
1195 (r * dK1[j].Kjj(iz, ia)[i] +
1196 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
1198 if (dK2[j].NumberOfFrequencies())
1199 dT2[j].Mat4(i).noalias() =
1202 (r * dK2[j].Kjj(iz, ia)[i] +
1203 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] + K2.Kjj(iz, ia)[i])
1207 const Numeric b2 =
b *
b, c2 =
c *
c, d2 =
d *
d, u2 =
u *
u, v2 =
v *
v,
1210 w2 * w2 + 2 * (b2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
1211 c2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
1212 d2 * (d2 * 0.5 + u2 - v2 - w2) +
1213 u2 * (u2 * 0.5 + v2 + w2) + v2 * (v2 * 0.5 + w2) +
1215 const Complex Const1 = std::sqrt(Complex(tmp, 0));
1216 const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
1217 const Complex tmp_x_sqrt = std::sqrt(Const2 + Const1);
1218 const Complex tmp_y_sqrt = std::sqrt(Const2 - Const1);
1219 const Complex x = tmp_x_sqrt * sqrt_05;
1220 const Complex y = tmp_y_sqrt * sqrt_05 * Complex(0, 1);
1221 const Complex x2 = x * x;
1222 const Complex y2 = y * y;
1223 const Complex cy = std::cos(y);
1224 const Complex sy = std::sin(y);
1225 const Complex cx = std::cosh(x);
1226 const Complex sx = std::sinh(x);
1230 const bool both_zero = y_zero and x_zero;
1231 const bool either_zero = y_zero or x_zero;
1241 const Complex ix = x_zero ? 0.0 : 1.0 / x;
1242 const Complex iy = y_zero ? 0.0 : 1.0 / y;
1243 const Complex inv_x2y2 =
1248 const Complex C0c = either_zero ? 1.0 : (cy * x2 + cx * y2) * inv_x2y2;
1250 either_zero ? 1.0 : (sy * x2 * iy + sx * y2 * ix) * inv_x2y2;
1251 const Complex C2c = both_zero ? 0.5 : (cx - cy) * inv_x2y2;
1252 const Complex C3c = both_zero ? 1.0 / 6.0
1253 : (x_zero ? 1.0 - sy * iy
1254 : y_zero ? sx * ix - 1.0
1255 : sx * ix - sy * iy) *
1258 const Numeric& C0 = real_val(C0c);
1259 const Numeric& C1 = real_val(C1c);
1260 const Numeric& C2 = real_val(C2c);
1261 const Numeric& C3 = real_val(C3c);
1262 T.Mat4(i).noalias() =
1263 exp_a * (Eigen::Matrix4d() << C0 + C2 * (b2 + c2 + d2),
1264 C1 *
b + C2 * (-
c *
u -
d *
v) +
1265 C3 * (
b * (b2 + c2 + d2) -
u * (
b *
u -
d *
w) -
1266 v * (
b *
v +
c *
w)),
1267 C1 *
c + C2 * (
b *
u -
d *
w) +
1268 C3 * (
c * (b2 + c2 + d2) -
u * (
c *
u +
d *
v) -
1269 w * (
b *
v +
c *
w)),
1270 C1 *
d + C2 * (
b *
v +
c *
w) +
1271 C3 * (
d * (b2 + c2 + d2) -
v * (
c *
u +
d *
v) +
1272 w * (
b *
u -
d *
w)),
1274 C1 *
b + C2 * (
c *
u +
d *
v) +
1275 C3 * (-
b * (-b2 + u2 + v2) +
c * (
b *
c -
v *
w) +
1276 d * (
b *
d +
u *
w)),
1277 C0 + C2 * (b2 - u2 - v2),
1278 C2 * (
b *
c -
v *
w) + C1 *
u +
1279 C3 * (
c * (
c *
u +
d *
v) -
u * (-b2 + u2 + v2) -
1280 w * (
b *
d +
u *
w)),
1281 C2 * (
b *
d +
u *
w) + C1 *
v +
1282 C3 * (
d * (
c *
u +
d *
v) -
v * (-b2 + u2 + v2) +
1283 w * (
b *
c -
v *
w)),
1285 C1 *
c + C2 * (-
b *
u +
d *
w) +
1286 C3 * (
b * (
b *
c -
v *
w) -
c * (-c2 + u2 + w2) +
1287 d * (
c *
d -
u *
v)),
1288 C2 * (
b *
c -
v *
w) - C1 *
u +
1289 C3 * (-
b * (
b *
u -
d *
w) +
u * (-c2 + u2 + w2) -
1290 v * (
c *
d -
u *
v)),
1291 C0 + C2 * (c2 - u2 - w2),
1292 C2 * (
c *
d -
u *
v) + C1 *
w +
1293 C3 * (-
d * (
b *
u -
d *
w) +
v * (
b *
c -
v *
w) -
1294 w * (-c2 + u2 + w2)),
1296 C1 *
d + C2 * (-
b *
v -
c *
w) +
1297 C3 * (
b * (
b *
d +
u *
w) +
c * (
c *
d -
u *
v) -
1298 d * (-d2 + v2 + w2)),
1299 C2 * (
b *
d +
u *
w) - C1 *
v +
1300 C3 * (-
b * (
b *
v +
c *
w) -
u * (
c *
d -
u *
v) +
1301 v * (-d2 + v2 + w2)),
1302 C2 * (
c *
d -
u *
v) - C1 *
w +
1303 C3 * (-
c * (
b *
v +
c *
w) +
u * (
b *
d +
u *
w) +
1304 w * (-d2 + v2 + w2)),
1305 C0 + C2 * (d2 - v2 - w2))
1308 for (Index j = 0; j < dK1.nelem(); j++) {
1309 if (dK1[j].NumberOfFrequencies()) {
1310 const Numeric da = -0.5 * (r * dK1[j].Kjj(iz, ia)[i] +
1311 ((j == it) ? dr_dT1 * (K1.Kjj(iz, ia)[i] +
1314 db = -0.5 * (r * dK1[j].K12(iz, ia)[i] +
1315 ((j == it) ? dr_dT1 * (K1.K12(iz, ia)[i] +
1318 dc = -0.5 * (r * dK1[j].K13(iz, ia)[i] +
1319 ((j == it) ? dr_dT1 * (K1.K13(iz, ia)[i] +
1322 dd = -0.5 * (r * dK1[j].K14(iz, ia)[i] +
1323 ((j == it) ? dr_dT1 * (K1.K14(iz, ia)[i] +
1326 du = -0.5 * (r * dK1[j].K23(iz, ia)[i] +
1327 ((j == it) ? dr_dT1 * (K1.K23(iz, ia)[i] +
1330 dv = -0.5 * (r * dK1[j].K24(iz, ia)[i] +
1331 ((j == it) ? dr_dT1 * (K1.K24(iz, ia)[i] +
1334 dw = -0.5 * (r * dK1[j].K34(iz, ia)[i] +
1335 ((j == it) ? dr_dT1 * (K1.K34(iz, ia)[i] +
1338 const Numeric db2 = 2 * db *
b, dc2 = 2 * dc *
c, dd2 = 2 * dd *
d,
1339 du2 = 2 * du *
u, dv2 = 2 * dv *
v, dw2 = 2 * dw *
w;
1340 const Numeric dtmp =
1342 2 * (db2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
1343 b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2) +
1344 dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
1345 c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2) +
1346 dd2 * (d2 * 0.5 + u2 - v2 - w2) +
1347 d2 * (dd2 * 0.5 + du2 - dv2 - dw2) +
1348 du2 * (u2 * 0.5 + v2 + w2) + u2 * (du2 * 0.5 + dv2 + dw2) +
1349 dv2 * (v2 * 0.5 + w2) + v2 * (dv2 * 0.5 + dw2) +
1350 4 * (db *
d *
u *
w - db *
c *
v *
w - dc *
d *
u *
v +
1351 b * dd *
u *
w -
b * dc *
v *
w -
c * dd *
u *
v +
1352 b *
d * du *
w -
b *
c * dv *
w -
c *
d * du *
v +
1353 b *
d *
u * dw -
b *
c *
v * dw -
c *
d *
u * dv));
1354 const Complex dConst1 = 0.5 * dtmp / Const1;
1355 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1357 x_zero ? 0 : (0.5 * (dConst2 + dConst1) / tmp_x_sqrt) * sqrt_05;
1358 const Complex dy = y_zero ? 0
1359 : (0.5 * (dConst2 - dConst1) / tmp_y_sqrt) *
1360 sqrt_05 * Complex(0, 1);
1361 const Complex dx2 = 2 * x * dx;
1362 const Complex dy2 = 2 * y * dy;
1363 const Complex dcy = -sy * dy;
1364 const Complex dsy = cy * dy;
1365 const Complex dcx = sx * dx;
1366 const Complex dsx = cx * dx;
1367 const Complex dix = -dx * ix * ix;
1368 const Complex diy = -dy * iy * iy;
1369 const Complex dx2dy2 = dx2 + dy2;
1370 const Complex dC0c =
1373 : (dcy * x2 + cy * dx2 + dcx * y2 + cx * dy2 - C0c * dx2dy2) *
1375 const Complex dC1c =
1377 : (dsy * x2 * iy + sy * dx2 * iy + sy * x2 * diy +
1378 dsx * y2 * ix + sx * dy2 * ix + sx * y2 * dix -
1381 const Complex dC2c =
1382 both_zero ? 0.0 : (dcx - dcy - C2c * dx2dy2) * inv_x2y2;
1383 const Complex dC3c =
1386 : ((x_zero ? -dsy * iy - sy * diy
1387 : y_zero ? dsx * ix + sx * dix
1388 : dsx * ix + sx * dix - dsy * iy - sy * diy) -
1392 const Numeric& dC0 = real_val(dC0c);
1393 const Numeric& dC1 = real_val(dC1c);
1394 const Numeric& dC2 = real_val(dC2c);
1395 const Numeric& dC3 = real_val(dC3c);
1396 dT1[j].Mat4(i).noalias() =
1400 << dC0 + dC2 * (b2 + c2 + d2) + C2 * (db2 + dc2 + dd2),
1401 db * C1 +
b * dC1 + dC2 * (-
c *
u -
d *
v) +
1402 C2 * (-dc *
u - dd *
v -
c * du -
d * dv) +
1403 dC3 * (
b * (b2 + c2 + d2) -
u * (
b *
u -
d *
w) -
1404 v * (
b *
v +
c *
w)) +
1405 C3 * (db * (b2 + c2 + d2) - du * (
b *
u -
d *
w) -
1406 dv * (
b *
v +
c *
w) +
b * (db2 + dc2 + dd2) -
1407 u * (db *
u - dd *
w) -
v * (db *
v + dc *
w) -
1408 u * (
b * du -
d * dw) -
v * (
b * dv +
c * dw)),
1409 dC1 *
c + C1 * dc + dC2 * (
b *
u -
d *
w) +
1410 C2 * (db *
u - dd *
w +
b * du -
d * dw) +
1411 dC3 * (
c * (b2 + c2 + d2) -
u * (
c *
u +
d *
v) -
1412 w * (
b *
v +
c *
w)) +
1413 C3 * (dc * (b2 + c2 + d2) - du * (
c *
u +
d *
v) -
1414 dw * (
b *
v +
c *
w) +
c * (db2 + dc2 + dd2) -
1415 u * (dc *
u + dd *
v) -
w * (db *
v + dc *
w) -
1416 u * (
c * du +
d * dv) -
w * (
b * dv +
c * dw)),
1417 dC1 *
d + C1 * dd + dC2 * (
b *
v +
c *
w) +
1418 C2 * (db *
v + dc *
w +
b * dv +
c * dw) +
1419 dC3 * (
d * (b2 + c2 + d2) -
v * (
c *
u +
d *
v) +
1420 w * (
b *
u -
d *
w)) +
1421 C3 * (dd * (b2 + c2 + d2) - dv * (
c *
u +
d *
v) +
1422 dw * (
b *
u -
d *
w) +
d * (db2 + dc2 + dd2) -
1423 v * (dc *
u + dd *
v) +
w * (db *
u - dd *
w) -
1424 v * (
c * du +
d * dv) +
w * (
b * du -
d * dw)),
1426 db * C1 +
b * dC1 + dC2 * (
c *
u +
d *
v) +
1427 C2 * (dc *
u + dd *
v +
c * du +
d * dv) +
1428 dC3 * (-
b * (-b2 + u2 + v2) +
c * (
b *
c -
v *
w) +
1429 d * (
b *
d +
u *
w)) +
1430 C3 * (-db * (-b2 + u2 + v2) + dc * (
b *
c -
v *
w) +
1431 dd * (
b *
d +
u *
w) -
b * (-db2 + du2 + dv2) +
1432 c * (db *
c - dv *
w) +
d * (db *
d + du *
w) +
1433 c * (
b * dc -
v * dw) +
d * (
b * dd +
u * dw)),
1434 dC0 + dC2 * (b2 - u2 - v2) + C2 * (db2 - du2 - dv2),
1435 dC2 * (
b *
c -
v *
w) +
1436 C2 * (db *
c +
b * dc - dv *
w -
v * dw) + dC1 *
u +
1438 dC3 * (
c * (
c *
u +
d *
v) -
u * (-b2 + u2 + v2) -
1439 w * (
b *
d +
u *
w)) +
1440 C3 * (dc * (
c *
u +
d *
v) - du * (-b2 + u2 + v2) -
1441 dw * (
b *
d +
u *
w) +
c * (dc *
u + dd *
v) -
1442 u * (-db2 + du2 + dv2) -
w * (db *
d + du *
w) +
1443 c * (
c * du +
d * dv) -
w * (
b * dd +
u * dw)),
1444 dC2 * (
b *
d +
u *
w) +
1445 C2 * (db *
d +
b * dd + du *
w +
u * dw) + dC1 *
v +
1447 dC3 * (
d * (
c *
u +
d *
v) -
v * (-b2 + u2 + v2) +
1448 w * (
b *
c -
v *
w)) +
1449 C3 * (dd * (
c *
u +
d *
v) - dv * (-b2 + u2 + v2) +
1450 dw * (
b *
c -
v *
w) +
d * (dc *
u + dd *
v) -
1451 v * (-db2 + du2 + dv2) +
w * (db *
c - dv *
w) +
1452 d * (
c * du +
d * dv) +
w * (
b * dc -
v * dw)),
1454 dC1 *
c + C1 * dc + dC2 * (-
b *
u +
d *
w) +
1455 C2 * (-db *
u + dd *
w -
b * du +
d * dw) +
1456 dC3 * (
b * (
b *
c -
v *
w) -
c * (-c2 + u2 + w2) +
1457 d * (
c *
d -
u *
v)) +
1458 C3 * (db * (
b *
c -
v *
w) - dc * (-c2 + u2 + w2) +
1459 dd * (
c *
d -
u *
v) +
b * (db *
c - dv *
w) -
1460 c * (-dc2 + du2 + dw2) +
d * (dc *
d - du *
v) +
1461 b * (
b * dc -
v * dw) +
d * (
c * dd -
u * dv)),
1462 dC2 * (
b *
c -
v *
w) +
1463 C2 * (db *
c +
b * dc - dv *
w -
v * dw) - dC1 *
u -
1465 dC3 * (-
b * (
b *
u -
d *
w) +
u * (-c2 + u2 + w2) -
1466 v * (
c *
d -
u *
v)) +
1467 C3 * (-db * (
b *
u -
d *
w) + du * (-c2 + u2 + w2) -
1468 dv * (
c *
d -
u *
v) -
b * (db *
u - dd *
w) +
1469 u * (-dc2 + du2 + dw2) -
v * (dc *
d - du *
v) -
1470 b * (
b * du -
d * dw) -
v * (
c * dd -
u * dv)),
1471 dC0 + dC2 * (c2 - u2 - w2) + C2 * (dc2 - du2 - dw2),
1472 dC2 * (
c *
d -
u *
v) +
1473 C2 * (dc *
d +
c * dd - du *
v -
u * dv) + dC1 *
w +
1475 dC3 * (-
d * (
b *
u -
d *
w) +
v * (
b *
c -
v *
w) -
1476 w * (-c2 + u2 + w2)) +
1477 C3 * (-dd * (
b *
u -
d *
w) + dv * (
b *
c -
v *
w) -
1478 dw * (-c2 + u2 + w2) -
d * (db *
u - dd *
w) +
1479 v * (db *
c - dv *
w) -
w * (-dc2 + du2 + dw2) -
1480 d * (
b * du -
d * dw) +
v * (
b * dc -
v * dw)),
1482 dC1 *
d + C1 * dd + dC2 * (-
b *
v -
c *
w) +
1483 C2 * (-db *
v - dc *
w -
b * dv -
c * dw) +
1484 dC3 * (
b * (
b *
d +
u *
w) +
c * (
c *
d -
u *
v) -
1485 d * (-d2 + v2 + w2)) +
1486 C3 * (db * (
b *
d +
u *
w) + dc * (
c *
d -
u *
v) -
1487 dd * (-d2 + v2 + w2) +
b * (db *
d + du *
w) +
1488 c * (dc *
d - du *
v) -
d * (-dd2 + dv2 + dw2) +
1489 b * (
b * dd +
u * dw) +
c * (
c * dd -
u * dv)),
1490 dC2 * (
b *
d +
u *
w) +
1491 C2 * (db *
d +
b * dd + du *
w +
u * dw) - dC1 *
v -
1493 dC3 * (-
b * (
b *
v +
c *
w) -
u * (
c *
d -
u *
v) +
1494 v * (-d2 + v2 + w2)) +
1495 C3 * (-db * (
b *
v +
c *
w) - du * (
c *
d -
u *
v) +
1496 dv * (-d2 + v2 + w2) -
b * (db *
v + dc *
w) -
1497 u * (dc *
d - du *
v) +
v * (-dd2 + dv2 + dw2) -
1498 b * (
b * dv +
c * dw) -
u * (
c * dd -
u * dv)),
1499 dC2 * (
c *
d -
u *
v) +
1500 C2 * (dc *
d +
c * dd - du *
v -
u * dv) - dC1 *
w -
1502 dC3 * (-
c * (
b *
v +
c *
w) +
u * (
b *
d +
u *
w) +
1503 w * (-d2 + v2 + w2)) +
1504 C3 * (-dc * (
b *
v +
c *
w) + du * (
b *
d +
u *
w) +
1505 dw * (-d2 + v2 + w2) -
c * (db *
v + dc *
w) +
1506 u * (db *
d + du *
w) +
w * (-dd2 + dv2 + dw2) -
1507 c * (
b * dv +
c * dw) +
u * (
b * dd +
u * dw)),
1508 dC0 + dC2 * (d2 - v2 - w2) + C2 * (dd2 - dv2 - dw2))
1511 if (dK2[j].NumberOfFrequencies()) {
1512 const Numeric da = -0.5 * (r * dK2[j].Kjj(iz, ia)[i] +
1513 ((j == it) ? dr_dT2 * (K1.Kjj(iz, ia)[i] +
1516 db = -0.5 * (r * dK2[j].K12(iz, ia)[i] +
1517 ((j == it) ? dr_dT2 * (K1.K12(iz, ia)[i] +
1520 dc = -0.5 * (r * dK2[j].K13(iz, ia)[i] +
1521 ((j == it) ? dr_dT2 * (K1.K13(iz, ia)[i] +
1524 dd = -0.5 * (r * dK2[j].K14(iz, ia)[i] +
1525 ((j == it) ? dr_dT2 * (K1.K14(iz, ia)[i] +
1528 du = -0.5 * (r * dK2[j].K23(iz, ia)[i] +
1529 ((j == it) ? dr_dT2 * (K1.K23(iz, ia)[i] +
1532 dv = -0.5 * (r * dK2[j].K24(iz, ia)[i] +
1533 ((j == it) ? dr_dT2 * (K1.K24(iz, ia)[i] +
1536 dw = -0.5 * (r * dK2[j].K34(iz, ia)[i] +
1537 ((j == it) ? dr_dT2 * (K1.K34(iz, ia)[i] +
1540 const Numeric db2 = 2 * db *
b, dc2 = 2 * dc *
c, dd2 = 2 * dd *
d,
1541 du2 = 2 * du *
u, dv2 = 2 * dv *
v, dw2 = 2 * dw *
w;
1542 const Numeric dtmp =
1544 2 * (db2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2) +
1545 b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2) +
1546 dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2) +
1547 c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2) +
1548 dd2 * (d2 * 0.5 + u2 - v2 - w2) +
1549 d2 * (dd2 * 0.5 + du2 - dv2 - dw2) +
1550 du2 * (u2 * 0.5 + v2 + w2) + u2 * (du2 * 0.5 + dv2 + dw2) +
1551 dv2 * (v2 * 0.5 + w2) + v2 * (dv2 * 0.5 + dw2) +
1552 4 * (db *
d *
u *
w - db *
c *
v *
w - dc *
d *
u *
v +
1553 b * dd *
u *
w -
b * dc *
v *
w -
c * dd *
u *
v +
1554 b *
d * du *
w -
b *
c * dv *
w -
c *
d * du *
v +
1555 b *
d *
u * dw -
b *
c *
v * dw -
c *
d *
u * dv));
1556 const Complex dConst1 = 0.5 * dtmp / Const1;
1557 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1559 x_zero ? 0 : (0.5 * (dConst2 + dConst1) / tmp_x_sqrt) * sqrt_05;
1560 const Complex dy = y_zero ? 0
1561 : (0.5 * (dConst2 - dConst1) / tmp_y_sqrt) *
1562 sqrt_05 * Complex(0, 1);
1563 const Complex dx2 = 2 * x * dx;
1564 const Complex dy2 = 2 * y * dy;
1565 const Complex dcy = -sy * dy;
1566 const Complex dsy = cy * dy;
1567 const Complex dcx = sx * dx;
1568 const Complex dsx = cx * dx;
1569 const Complex dix = -dx * ix * ix;
1570 const Complex diy = -dy * iy * iy;
1571 const Complex dx2dy2 = dx2 + dy2;
1572 const Complex dC0c =
1575 : (dcy * x2 + cy * dx2 + dcx * y2 + cx * dy2 - C0c * dx2dy2) *
1577 const Complex dC1c =
1579 : (dsy * x2 * iy + sy * dx2 * iy + sy * x2 * diy +
1580 dsx * y2 * ix + sx * dy2 * ix + sx * y2 * dix -
1583 const Complex dC2c =
1584 both_zero ? 0.0 : (dcx - dcy - C2c * dx2dy2) * inv_x2y2;
1585 const Complex dC3c =
1588 : ((x_zero ? -dsy * iy - sy * diy
1589 : y_zero ? dsx * ix + sx * dix
1590 : dsx * ix + sx * dix - dsy * iy - sy * diy) -
1594 const Numeric& dC0 = real_val(dC0c);
1595 const Numeric& dC1 = real_val(dC1c);
1596 const Numeric& dC2 = real_val(dC2c);
1597 const Numeric& dC3 = real_val(dC3c);
1598 dT2[j].Mat4(i).noalias() =
1602 << dC0 + dC2 * (b2 + c2 + d2) + C2 * (db2 + dc2 + dd2),
1603 db * C1 +
b * dC1 + dC2 * (-
c *
u -
d *
v) +
1604 C2 * (-dc *
u - dd *
v -
c * du -
d * dv) +
1605 dC3 * (
b * (b2 + c2 + d2) -
u * (
b *
u -
d *
w) -
1606 v * (
b *
v +
c *
w)) +
1607 C3 * (db * (b2 + c2 + d2) - du * (
b *
u -
d *
w) -
1608 dv * (
b *
v +
c *
w) +
b * (db2 + dc2 + dd2) -
1609 u * (db *
u - dd *
w) -
v * (db *
v + dc *
w) -
1610 u * (
b * du -
d * dw) -
v * (
b * dv +
c * dw)),
1611 dC1 *
c + C1 * dc + dC2 * (
b *
u -
d *
w) +
1612 C2 * (db *
u - dd *
w +
b * du -
d * dw) +
1613 dC3 * (
c * (b2 + c2 + d2) -
u * (
c *
u +
d *
v) -
1614 w * (
b *
v +
c *
w)) +
1615 C3 * (dc * (b2 + c2 + d2) - du * (
c *
u +
d *
v) -
1616 dw * (
b *
v +
c *
w) +
c * (db2 + dc2 + dd2) -
1617 u * (dc *
u + dd *
v) -
w * (db *
v + dc *
w) -
1618 u * (
c * du +
d * dv) -
w * (
b * dv +
c * dw)),
1619 dC1 *
d + C1 * dd + dC2 * (
b *
v +
c *
w) +
1620 C2 * (db *
v + dc *
w +
b * dv +
c * dw) +
1621 dC3 * (
d * (b2 + c2 + d2) -
v * (
c *
u +
d *
v) +
1622 w * (
b *
u -
d *
w)) +
1623 C3 * (dd * (b2 + c2 + d2) - dv * (
c *
u +
d *
v) +
1624 dw * (
b *
u -
d *
w) +
d * (db2 + dc2 + dd2) -
1625 v * (dc *
u + dd *
v) +
w * (db *
u - dd *
w) -
1626 v * (
c * du +
d * dv) +
w * (
b * du -
d * dw)),
1628 db * C1 +
b * dC1 + dC2 * (
c *
u +
d *
v) +
1629 C2 * (dc *
u + dd *
v +
c * du +
d * dv) +
1630 dC3 * (-
b * (-b2 + u2 + v2) +
c * (
b *
c -
v *
w) +
1631 d * (
b *
d +
u *
w)) +
1632 C3 * (-db * (-b2 + u2 + v2) + dc * (
b *
c -
v *
w) +
1633 dd * (
b *
d +
u *
w) -
b * (-db2 + du2 + dv2) +
1634 c * (db *
c - dv *
w) +
d * (db *
d + du *
w) +
1635 c * (
b * dc -
v * dw) +
d * (
b * dd +
u * dw)),
1636 dC0 + dC2 * (b2 - u2 - v2) + C2 * (db2 - du2 - dv2),
1637 dC2 * (
b *
c -
v *
w) +
1638 C2 * (db *
c +
b * dc - dv *
w -
v * dw) + dC1 *
u +
1640 dC3 * (
c * (
c *
u +
d *
v) -
u * (-b2 + u2 + v2) -
1641 w * (
b *
d +
u *
w)) +
1642 C3 * (dc * (
c *
u +
d *
v) - du * (-b2 + u2 + v2) -
1643 dw * (
b *
d +
u *
w) +
c * (dc *
u + dd *
v) -
1644 u * (-db2 + du2 + dv2) -
w * (db *
d + du *
w) +
1645 c * (
c * du +
d * dv) -
w * (
b * dd +
u * dw)),
1646 dC2 * (
b *
d +
u *
w) +
1647 C2 * (db *
d +
b * dd + du *
w +
u * dw) + dC1 *
v +
1649 dC3 * (
d * (
c *
u +
d *
v) -
v * (-b2 + u2 + v2) +
1650 w * (
b *
c -
v *
w)) +
1651 C3 * (dd * (
c *
u +
d *
v) - dv * (-b2 + u2 + v2) +
1652 dw * (
b *
c -
v *
w) +
d * (dc *
u + dd *
v) -
1653 v * (-db2 + du2 + dv2) +
w * (db *
c - dv *
w) +
1654 d * (
c * du +
d * dv) +
w * (
b * dc -
v * dw)),
1656 dC1 *
c + C1 * dc + dC2 * (-
b *
u +
d *
w) +
1657 C2 * (-db *
u + dd *
w -
b * du +
d * dw) +
1658 dC3 * (
b * (
b *
c -
v *
w) -
c * (-c2 + u2 + w2) +
1659 d * (
c *
d -
u *
v)) +
1660 C3 * (db * (
b *
c -
v *
w) - dc * (-c2 + u2 + w2) +
1661 dd * (
c *
d -
u *
v) +
b * (db *
c - dv *
w) -
1662 c * (-dc2 + du2 + dw2) +
d * (dc *
d - du *
v) +
1663 b * (
b * dc -
v * dw) +
d * (
c * dd -
u * dv)),
1664 dC2 * (
b *
c -
v *
w) +
1665 C2 * (db *
c +
b * dc - dv *
w -
v * dw) - dC1 *
u -
1667 dC3 * (-
b * (
b *
u -
d *
w) +
u * (-c2 + u2 + w2) -
1668 v * (
c *
d -
u *
v)) +
1669 C3 * (-db * (
b *
u -
d *
w) + du * (-c2 + u2 + w2) -
1670 dv * (
c *
d -
u *
v) -
b * (db *
u - dd *
w) +
1671 u * (-dc2 + du2 + dw2) -
v * (dc *
d - du *
v) -
1672 b * (
b * du -
d * dw) -
v * (
c * dd -
u * dv)),
1673 dC0 + dC2 * (c2 - u2 - w2) + C2 * (dc2 - du2 - dw2),
1674 dC2 * (
c *
d -
u *
v) +
1675 C2 * (dc *
d +
c * dd - du *
v -
u * dv) + dC1 *
w +
1677 dC3 * (-
d * (
b *
u -
d *
w) +
v * (
b *
c -
v *
w) -
1678 w * (-c2 + u2 + w2)) +
1679 C3 * (-dd * (
b *
u -
d *
w) + dv * (
b *
c -
v *
w) -
1680 dw * (-c2 + u2 + w2) -
d * (db *
u - dd *
w) +
1681 v * (db *
c - dv *
w) -
w * (-dc2 + du2 + dw2) -
1682 d * (
b * du -
d * dw) +
v * (
b * dc -
v * dw)),
1684 dC1 *
d + C1 * dd + dC2 * (-
b *
v -
c *
w) +
1685 C2 * (-db *
v - dc *
w -
b * dv -
c * dw) +
1686 dC3 * (
b * (
b *
d +
u *
w) +
c * (
c *
d -
u *
v) -
1687 d * (-d2 + v2 + w2)) +
1688 C3 * (db * (
b *
d +
u *
w) + dc * (
c *
d -
u *
v) -
1689 dd * (-d2 + v2 + w2) +
b * (db *
d + du *
w) +
1690 c * (dc *
d - du *
v) -
d * (-dd2 + dv2 + dw2) +
1691 b * (
b * dd +
u * dw) +
c * (
c * dd -
u * dv)),
1692 dC2 * (
b *
d +
u *
w) +
1693 C2 * (db *
d +
b * dd + du *
w +
u * dw) - dC1 *
v -
1695 dC3 * (-
b * (
b *
v +
c *
w) -
u * (
c *
d -
u *
v) +
1696 v * (-d2 + v2 + w2)) +
1697 C3 * (-db * (
b *
v +
c *
w) - du * (
c *
d -
u *
v) +
1698 dv * (-d2 + v2 + w2) -
b * (db *
v + dc *
w) -
1699 u * (dc *
d - du *
v) +
v * (-dd2 + dv2 + dw2) -
1700 b * (
b * dv +
c * dw) -
u * (
c * dd -
u * dv)),
1701 dC2 * (
c *
d -
u *
v) +
1702 C2 * (dc *
d +
c * dd - du *
v -
u * dv) - dC1 *
w -
1704 dC3 * (-
c * (
b *
v +
c *
w) +
u * (
b *
d +
u *
w) +
1705 w * (-d2 + v2 + w2)) +
1706 C3 * (-dc * (
b *
v +
c *
w) + du * (
b *
d +
u *
w) +
1707 dw * (-d2 + v2 + w2) -
c * (db *
v + dc *
w) +
1708 u * (db *
d + du *
w) +
w * (-dd2 + dv2 + dw2) -
1709 c * (
b * dv +
c * dw) +
u * (
b * dd +
u * dw)),
1710 dC0 + dC2 * (d2 - v2 - w2) + C2 * (dd2 - dv2 - dw2))
1719 const PropagationMatrix& K1,
1720 const PropagationMatrix& K2,
1721 const Numeric& r)
noexcept {
1722 switch (K1.StokesDimensions()) {
1741 const PropagationMatrix& K1,
1742 const PropagationMatrix& K2,
1743 const ArrayOfPropagationMatrix& dK1,
1744 const ArrayOfPropagationMatrix& dK2,
1746 const Numeric& dr_dT1 = 0,
1747 const Numeric& dr_dT2 = 0,
1748 const Index it = -1,
1750 const Index ia = 0) noexcept {
1751 switch (K1.StokesDimensions()) {
1753 dtransmat4(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1756 dtransmat3(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1759 dtransmat2(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1762 dtransmat1(T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dT1, dr_dT2, it, iz, ia);
1770 const PropagationMatrix& K1,
1771 const PropagationMatrix& K2,
1772 const ArrayOfPropagationMatrix& dK1,
1773 const ArrayOfPropagationMatrix& dK2,
1775 const Numeric& dr_dtemp1,
1776 const Numeric& dr_dtemp2,
1777 const Index temp_deriv_pos) {
1778 if (not dT1.
nelem())
1782 T, dT1, dT2, K1, K2, dK1, dK2, r, dr_dtemp1, dr_dtemp2, temp_deriv_pos);
1788 const PropagationMatrix& K,
1789 const StokesVector&
a,
1790 const StokesVector& S,
1791 const ArrayOfPropagationMatrix& dK,
1792 const ArrayOfStokesVector& da,
1793 const ArrayOfStokesVector& dS,
1794 const ConstVectorView& B,
1795 const ConstVectorView& dB_dT,
1797 const bool& jacobian_do) {
1798 for (Index i = 0; i < K.NumberOfFrequencies(); i++) {
1799 if (K.IsRotational(i)) {
1802 for (Index j = 0; j < jacobian_quantities.
nelem(); j++)
1803 if (dJ[j].Frequencies()) dJ[j].SetZero(i);
1810 const auto invK = inv_prop_matrix<4>(K.Data()(0, 0, i, joker));
1813 for (Index j = 0; j < jacobian_quantities.
nelem(); j++) {
1814 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and
1815 dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1816 dJ[j].Vec4(i).noalias() =
1824 jacobian_quantities[j] == Jacobian::Atm::Temperature,
1826 prop_matrix<4>(dK[j].Data()(0, 0, i, joker)) * J.
Vec4(i));
1831 J_add.
Vec4(i) = invK * J_add.
Vec4(i);
1836 const auto invK = inv_prop_matrix<3>(K.Data()(0, 0, i, joker));
1839 for (Index j = 0; j < jacobian_quantities.
nelem(); j++) {
1841 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and
1842 dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1843 dJ[j].Vec3(i).noalias() =
1851 jacobian_quantities[j] == Jacobian::Atm::Temperature,
1853 prop_matrix<3>(dK[j].Data()(0, 0, i, joker)) * J.
Vec3(i));
1858 J_add.
Vec3(i) = invK * J_add.
Vec3(i);
1863 const auto invK = inv_prop_matrix<2>(K.Data()(0, 0, i, joker));
1866 for (Index j = 0; j < jacobian_quantities.
nelem(); j++) {
1868 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and
1869 dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1870 dJ[j].Vec2(i).noalias() =
1878 jacobian_quantities[j] == Jacobian::Atm::Temperature,
1880 prop_matrix<2>(dK[j].Data()(0, 0, i, joker)) * J.
Vec2(i));
1885 J_add.
Vec2(i) = invK * J_add.
Vec2(i);
1890 const auto invK = inv_prop_matrix<1>(K.Data()(0, 0, i, joker));
1893 for (Index j = 0; j < jacobian_quantities.
nelem(); j++) {
1895 if (dJ[j].Frequencies() == da[j].NumberOfFrequencies() and
1896 dJ[j].Frequencies() == dS[j].NumberOfFrequencies()) {
1897 dJ[j].Vec1(i).noalias() =
1905 jacobian_quantities[j] == Jacobian::Atm::Temperature,
1907 prop_matrix<1>(dK[j].Data()(0, 0, i, joker)) * J.
Vec1(i));
1912 J_add.
Vec1(i) = invK * J_add.
Vec1(i);
1937 [[maybe_unused]]
const PropagationMatrix& K1,
1938 [[maybe_unused]]
const PropagationMatrix& K2,
1939 [[maybe_unused]]
const ArrayOfPropagationMatrix& dK1,
1940 [[maybe_unused]]
const ArrayOfPropagationMatrix& dK2,
1941 [[maybe_unused]]
const Numeric r,
1942 [[maybe_unused]]
const Vector& dr1,
1943 [[maybe_unused]]
const Vector& dr2,
1944 [[maybe_unused]]
const Index ia,
1945 [[maybe_unused]]
const Index iz,
1950 for (
size_t i = 0; i < dI1.size(); i++) {
1951 dI1[i].addDerivEmission(PiT, dT1[i], T, I, dJ1[i]);
1952 dI2[i].addDerivEmission(PiT, dT2[i], T, I, dJ2[i]);
1959 for (
size_t i = 0; i < dI1.size(); i++) {
1960 dI1[i].addDerivTransmission(PiT, dT1[i], I);
1961 dI2[i].addDerivTransmission(PiT, dT2[i], I);
1969 "Cannot support derivatives with current integration method\n");
1975 K1.Data()(ia, iz, joker, joker),
1976 K2.Data()(ia, iz, joker, joker),
1987 const Index n = T.
nelem();
1988 const Index nf = n ? T[0].Frequencies() : 1;
1989 const Index ns = n ? T[0].stokes_dim : 1;
1993 for (Index i = 1; i < n; i++) PiT[i].
mul(PiT[i - 1], T[i]);
1996 for (Index i = 1; i < n; i++) PiT[i].
mul(T[i], PiT[i - 1]);
2016 const Index np = I.
nelem();
2017 const Index nv = np ? I[0].Frequencies() : 0;
2018 const Index ns = np ? I[0].stokes_dim : 1;
2019 const Index nq = np ? dI[0][0].
nelem() : 0;
2022 for (Index ip = 0; ip < np; ip++)
2023 I[ip].setBackscatterTransmission(I_incoming, PiTr[ip], PiTf[ip], Z[ip]);
2025 for (Index ip = 0; ip < np; ip++) {
2026 for (Index iq = 0; iq < nq; iq++) {
2027 dI[ip][ip][iq].setBackscatterTransmissionDerivative(
2028 I_incoming, PiTr[ip], PiTf[ip], dZ[ip][iq]);
2037 BackscatterSolverCommutativeTransmissionStokesDimOne:
2038 for (Index ip = 0; ip < np; ip++) {
2039 for (Index j = ip; j < np; j++) {
2040 for (Index iq = 0; iq < nq; iq++) {
2041 for (Index iv = 0; iv < nv; iv++) {
2042 dI[ip][j][iq].Vec1(iv).noalias() +=
2043 T[ip].Mat1(iv).inverse() *
2044 (dT1[ip][iq].Mat1(iv) + dT2[ip][iq].Mat1(iv)) *
2047 if (j < np - 1 and j > ip)
2048 dI[ip][j][iq].Vec1(iv).noalias() +=
2049 T[ip + 1].Mat1(iv).inverse() *
2050 (dT1[ip][iq].Mat1(iv) + dT2[ip][iq].Mat1(iv)) *
2058 for (Index ip = 0; ip < np; ip++) {
2059 for (Index j = ip; j < np; j++) {
2060 for (Index iq = 0; iq < nq; iq++) {
2061 for (Index iv = 0; iv < nv; iv++) {
2062 dI[ip][j][iq].Vec2(iv).noalias() +=
2063 T[ip].Mat2(iv).inverse() *
2064 (dT1[ip][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
2067 if (j < np - 1 and j > ip)
2068 dI[ip][j][iq].Vec2(iv).noalias() +=
2069 T[ip + 1].Mat2(iv).inverse() *
2070 (dT1[ip][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
2078 for (Index ip = 0; ip < np; ip++) {
2079 for (Index j = ip; j < np; j++) {
2080 for (Index iq = 0; iq < nq; iq++) {
2081 for (Index iv = 0; iv < nv; iv++) {
2082 dI[ip][j][iq].Vec3(iv).noalias() +=
2083 T[ip].Mat3(iv).inverse() *
2084 (dT1[ip][iq].Mat3(iv) + dT2[ip][iq].Mat3(iv)) *
2087 if (j < np - 1 and j > ip)
2088 dI[ip][j][iq].Vec3(iv).noalias() +=
2089 T[ip + 1].Mat3(iv).inverse() *
2090 (dT1[ip][iq].Mat3(iv) + dT2[ip][iq].Mat3(iv)) *
2098 for (Index ip = 0; ip < np; ip++) {
2099 for (Index j = ip; j < np; j++) {
2100 for (Index iq = 0; iq < nq; iq++) {
2101 for (Index iv = 0; iv < nv; iv++) {
2102 dI[ip][j][iq].Vec4(iv).noalias() +=
2103 T[ip].Mat4(iv).inverse() *
2104 (dT1[ip][iq].Mat4(iv) + dT2[ip][iq].Mat4(iv)) *
2107 if (j < np - 1 and j > ip)
2108 dI[ip][j][iq].Vec4(iv).noalias() +=
2109 T[ip + 1].Mat4(iv).inverse() *
2110 (dT1[ip][iq].Mat4(iv) + dT2[ip][iq].Mat4(iv)) *
2120 std::runtime_error(
"Do not activate this code. It is not ready yet");
2124 goto BackscatterSolverCommutativeTransmissionStokesDimOne;
2127 for (Index ip = 0; ip < np; ip++) {
2128 for (Index j = ip; j < np; j++) {
2129 for (Index iq = 0; iq < nq; iq++) {
2130 for (Index iv = 0; iv < nv; iv++) {
2132 dI[ip][j][iq].Vec2(iv).noalias() +=
2133 PiTr[ip - 2].Mat2(iv) *
2134 (dT1[ip - 1][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
2135 PiTr[ip - 1].Mat2(iv).inverse() * I[j].Vec2(iv);
2138 dI[ip][j][iq].Vec2(iv).noalias() +=
2139 PiTr[ip].Mat2(iv) * Z[ip].Mat2(iv) *
2140 PiTf[ip - 2].Mat2(iv) *
2141 (dT1[ip][iq].Mat2(iv) + dT2[ip + 1][iq].Mat2(iv)) *
2142 PiTf[ip - 1].Mat2(iv).inverse() * PiTf[ip].Mat2(iv) *
2145 dI[ip][j][iq].Vec2(iv).noalias() +=
2146 (dT1[ip - 1][iq].Mat2(iv) + dT2[ip][iq].Mat2(iv)) *
2147 PiTr[ip - 1].Mat2(iv).inverse() * I[j].Vec2(iv);
2150 dI[ip][j][iq].Vec2(iv).noalias() +=
2151 PiTr[ip].Mat2(iv) * Z[ip].Mat2(iv) *
2152 (dT1[ip][iq].Mat2(iv) + dT2[ip + 1][iq].Mat2(iv)) *
2153 PiTf[ip - 1].Mat2(iv).inverse() * PiTf[ip].Mat2(iv) *
2171 const ConstMatrixView& pnd) {
2172 const Index ns = Pe.ncols();
2173 const Index nv = Pe.npages();
2174 const Index np = Pe.nbooks();
2175 const Index ne = Pe.nshelves();
2179 for (Index ip = 0; ip < np; ip++) {
2184 for (Index iv = 0; iv < nv; iv++)
2185 for (Index ie = 0; ie < ne; ie++)
2186 aotm[ip].Mat4(iv).noalias() +=
2187 pnd(ie, ip) * prop_matrix<4>(Pe(ie, ip, iv, joker, joker));
2190 for (Index iv = 0; iv < nv; iv++)
2191 for (Index ie = 0; ie < ne; ie++)
2192 aotm[ip].Mat3(iv).noalias() +=
2193 pnd(ie, ip) * prop_matrix<3>(Pe(ie, ip, iv, joker, joker));
2196 for (Index iv = 0; iv < nv; iv++)
2197 for (Index ie = 0; ie < ne; ie++)
2198 aotm[ip].Mat2(iv).noalias() +=
2199 pnd(ie, ip) * prop_matrix<2>(Pe(ie, ip, iv, joker, joker));
2202 for (Index iv = 0; iv < nv; iv++)
2203 for (Index ie = 0; ie < ne; ie++)
2204 aotm[ip].Mat1(iv).noalias() +=
2205 pnd(ie, ip) * prop_matrix<1>(Pe(ie, ip, iv, joker, joker));
2213 const ConstTensor5View& Pe,
const ArrayOfMatrix& dpnd_dx) {
2214 const Index ns = Pe.ncols();
2215 const Index nv = Pe.npages();
2216 const Index np = Pe.nbooks();
2217 const Index ne = Pe.nshelves();
2218 const Index nq = dpnd_dx.nelem();
2223 for (Index ip = 0; ip < np; ip++) {
2224 for (Index iq = 0; iq < nq; iq++) {
2225 aoaotm[ip][iq].setZero();
2227 if (not dpnd_dx[iq].empty()) {
2230 for (Index iv = 0; iv < nv; iv++)
2231 for (Index ie = 0; ie < ne; ie++)
2232 aoaotm[ip][iq].Mat4(iv).noalias() +=
2233 dpnd_dx[iq](ie, ip) *
2234 prop_matrix<4>(Pe(ie, ip, iv, joker, joker));
2237 for (Index iv = 0; iv < nv; iv++)
2238 for (Index ie = 0; ie < ne; ie++)
2239 aoaotm[ip][iq].Mat3(iv).noalias() +=
2240 dpnd_dx[iq](ie, ip) *
2241 prop_matrix<3>(Pe(ie, ip, iv, joker, joker));
2244 for (Index iv = 0; iv < nv; iv++)
2245 for (Index ie = 0; ie < ne; ie++)
2246 aoaotm[ip][iq].Mat2(iv).noalias() +=
2247 dpnd_dx[iq](ie, ip) *
2248 prop_matrix<2>(Pe(ie, ip, iv, joker, joker));
2251 for (Index iv = 0; iv < nv; iv++)
2252 for (Index ie = 0; ie < ne; ie++)
2253 aoaotm[ip][iq].Mat1(iv).noalias() +=
2254 dpnd_dx[iq](ie, ip) *
2255 prop_matrix<1>(Pe(ie, ip, iv, joker, joker));
2267 for (
const auto& T : tm.
T4) os << T <<
'\n';
2268 for (
const auto& T : tm.
T3) os << T <<
'\n';
2269 for (
const auto& T : tm.
T2) os << T <<
'\n';
2270 for (
const auto& T : tm.
T1) os << T <<
'\n';
2276 for (
const auto& R : rv.
R4) os << R.transpose() <<
'\n';
2277 for (
const auto& R : rv.
R3) os << R.transpose() <<
'\n';
2278 for (
const auto& R : rv.
R2) os << R.transpose() <<
'\n';
2279 for (
const auto& R : rv.
R1) os << R.transpose() <<
'\n';
2284 for (
auto& T : tm.
T4)
2285 is >>
double_imanip() >> T(0, 0) >> T(0, 1) >> T(0, 2) >> T(0, 3) >>
2286 T(1, 0) >> T(1, 1) >> T(1, 2) >> T(1, 3) >> T(2, 0) >> T(2, 1) >>
2287 T(2, 2) >> T(2, 3) >> T(3, 0) >> T(3, 1) >> T(3, 2) >> T(3, 3);
2288 for (
auto& T : tm.
T3)
2289 is >>
double_imanip() >> T(0, 0) >> T(0, 1) >> T(0, 2) >> T(1, 0) >>
2290 T(1, 1) >> T(1, 2) >> T(2, 0) >> T(2, 1) >> T(2, 2);
2291 for (
auto& T : tm.
T2)
2292 is >>
double_imanip() >> T(0, 0) >> T(0, 1) >> T(1, 0) >> T(1, 1);
2298 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
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
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