ARTS 2.5.4 (git: 4c0d3b4d)
transmissionmatrix.h
Go to the documentation of this file.
1/* Copyright (C) 2018
2 * Richard Larsson <ric.larsson@gmail.com>
3 *
4 * This program is free software; you can redistribute it and/or modify it
5 * under the terms of the GNU General Public License as published by the
6 * Free Software Foundation; either version 2, or (at your option) any
7 * later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program; if not, write to the Free Software
16 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17 * USA. */
18
29#ifndef transmissionmatrix_h
30#define transmissionmatrix_h
31
32#pragma GCC diagnostic push
33#pragma GCC diagnostic ignored "-Wshadow"
34#pragma GCC diagnostic ignored "-Wconversion"
35#pragma GCC diagnostic ignored "-Wfloat-conversion"
36#include <Eigen/Dense>
37#include <Eigen/StdVector>
38#pragma GCC diagnostic pop
39
40#include "jacobian.h"
41#include "propagationmatrix.h"
42
46 std::vector<Eigen::Matrix4d, Eigen::aligned_allocator<Eigen::Matrix4d>> T4;
47 std::vector<Eigen::Matrix3d, Eigen::aligned_allocator<Eigen::Matrix3d>> T3;
48 std::vector<Eigen::Matrix2d, Eigen::aligned_allocator<Eigen::Matrix2d>> T2;
49 std::vector<Eigen::Matrix<double, 1, 1>,
50 Eigen::aligned_allocator<Eigen::Matrix<double, 1, 1>>>
52
58 TransmissionMatrix(Index nf = 0, Index stokes = 1)
59 : stokes_dim(stokes),
60 T4(stokes_dim == 4 ? nf : 0, Eigen::Matrix4d::Identity()),
61 T3(stokes_dim == 3 ? nf : 0, Eigen::Matrix3d::Identity()),
62 T2(stokes_dim == 2 ? nf : 0, Eigen::Matrix2d::Identity()),
63 T1(stokes_dim == 1 ? nf : 0, Eigen::Matrix<double, 1, 1>::Identity()) {
64 ARTS_ASSERT(stokes_dim < 5 and stokes_dim > 0);
65 }
66
72 : stokes_dim(tm.stokes_dim),
73 T4(std::move(tm.T4)),
74 T3(std::move(tm.T3)),
75 T2(std::move(tm.T2)),
76 T1(std::move(tm.T1)) {}
77
83
90
97 stokes_dim = tm.stokes_dim;
98 T4 = std::move(tm.T4);
99 T3 = std::move(tm.T3);
100 T2 = std::move(tm.T2);
101 T1 = std::move(tm.T1);
102 return *this;
103 }
104
110 explicit TransmissionMatrix(const PropagationMatrix& pm, const Numeric& r = 1.0);
111
112 operator Tensor3() const {
114 for (size_t i = 0; i < T4.size(); i++)
115 for (size_t j = 0; j < 4; j++)
116 for (size_t k = 0; k < 4; k++) T(i, j, k) = T4[i](j, k);
117 for (size_t i = 0; i < T3.size(); i++)
118 for (size_t j = 0; j < 3; j++)
119 for (size_t k = 0; k < 3; k++) T(i, j, k) = T3[i](j, k);
120 for (size_t i = 0; i < T2.size(); i++)
121 for (size_t j = 0; j < 2; j++)
122 for (size_t k = 0; k < 2; k++) T(i, j, k) = T2[i](j, k);
123 for (size_t i = 0; i < T1.size(); i++) T(i, 0, 0) = T1[i](0, 0);
124 return T;
125 }
126
132 [[nodiscard]] const Eigen::Matrix4d& Mat4(size_t i) const { return T4[i]; }
133
139 [[nodiscard]] const Eigen::Matrix3d& Mat3(size_t i) const { return T3[i]; }
140
146 [[nodiscard]] const Eigen::Matrix2d& Mat2(size_t i) const { return T2[i]; }
147
153 [[nodiscard]] const Eigen::Matrix<double, 1, 1>& Mat1(size_t i) const { return T1[i]; }
154
160 #pragma GCC diagnostic push
161 #pragma GCC diagnostic ignored "-Wreturn-type"
162 [[nodiscard]] Eigen::MatrixXd Mat(size_t i) const {
163 switch (stokes_dim) {
164 case 1:
165 return Mat1(i);
166 case 2:
167 return Mat2(i);
168 case 3:
169 return Mat3(i);
170 case 4:
171 return Mat4(i);
172 }
173 }
174 #pragma GCC diagnostic pop
175
181 Eigen::Matrix4d& Mat4(size_t i) { return T4[i]; }
182
188 Eigen::Matrix3d& Mat3(size_t i) { return T3[i]; }
189
195 Eigen::Matrix2d& Mat2(size_t i) { return T2[i]; }
196
202 Eigen::Matrix<double, 1, 1>& Mat1(size_t i) { return T1[i]; }
203
205 void setIdentity() {
206 std::fill(T4.begin(), T4.end(), Eigen::Matrix4d::Identity());
207 std::fill(T3.begin(), T3.end(), Eigen::Matrix3d::Identity());
208 std::fill(T2.begin(), T2.end(), Eigen::Matrix2d::Identity());
209 std::fill(T1.begin(), T1.end(), Eigen::Matrix<double, 1, 1>::Identity());
210 }
211
213 void setZero() {
214 std::fill(T4.begin(), T4.end(), Eigen::Matrix4d::Zero());
215 std::fill(T3.begin(), T3.end(), Eigen::Matrix3d::Zero());
216 std::fill(T2.begin(), T2.end(), Eigen::Matrix2d::Zero());
217 std::fill(T1.begin(), T1.end(), Eigen::Matrix<double, 1, 1>::Zero());
218 }
219
227 void mul(const TransmissionMatrix& A, const TransmissionMatrix& B) {
228 for (size_t i = 0; i < T4.size(); i++) T4[i].noalias() = A.T4[i] * B.T4[i];
229 for (size_t i = 0; i < T3.size(); i++) T3[i].noalias() = A.T3[i] * B.T3[i];
230 for (size_t i = 0; i < T2.size(); i++) T2[i].noalias() = A.T2[i] * B.T2[i];
231 for (size_t i = 0; i < T1.size(); i++) T1[i].noalias() = A.T1[i] * B.T1[i];
232 }
233
234
243 for (size_t i = 0; i < T4.size(); i++) T4[i] = A.T4[i] * B.T4[i];
244 for (size_t i = 0; i < T3.size(); i++) T3[i] = A.T3[i] * B.T3[i];
245 for (size_t i = 0; i < T2.size(); i++) T2[i] = A.T2[i] * B.T2[i];
246 for (size_t i = 0; i < T1.size(); i++) T1[i] = A.T1[i] * B.T1[i];
247 }
248
256 Numeric operator()(const Index i, const Index j, const Index k) const {
257 switch (stokes_dim) {
258 case 4:
259 return T4[i](j, k);
260 case 3:
261 return T3[i](j, k);
262 case 2:
263 return T2[i](j, k);
264 default:
265 return T1[i](j, k);
266 }
267 }
268
270 [[nodiscard]] Index Frequencies() const {
271 switch (stokes_dim) {
272 case 4:
273 return Index(T4.size());
274 case 3:
275 return Index(T3.size());
276 case 2:
277 return Index(T2.size());
278 default:
279 return Index(T1.size());
280 }
281 }
282
289 for (size_t i = 0; i < T4.size(); i++)
290 T4[i].noalias() = lstm.scale * lstm.bas.Mat4(i);
291 for (size_t i = 0; i < T3.size(); i++)
292 T3[i].noalias() = lstm.scale * lstm.bas.Mat3(i);
293 for (size_t i = 0; i < T2.size(); i++)
294 T2[i].noalias() = lstm.scale * lstm.bas.Mat2(i);
295 for (size_t i = 0; i < T1.size(); i++)
296 T1[i].noalias() = lstm.scale * lstm.bas.Mat1(i);
297 return *this;
298 }
299
306 std::transform(T4.begin(), T4.end(), T4.begin(), [scale](auto& T){return T*scale;});
307 std::transform(T3.begin(), T3.end(), T3.begin(), [scale](auto& T){return T*scale;});
308 std::transform(T2.begin(), T2.end(), T2.begin(), [scale](auto& T){return T*scale;});
309 std::transform(T1.begin(), T1.end(), T1.begin(), [scale](auto& T){return T*scale;});
310 return *this;
311 }
312
319 operator=(lstm.bas);
320 operator*=(lstm.scale);
321 return *this;
322 }
323
325 friend std::ostream& operator<<(std::ostream& os,
326 const TransmissionMatrix& tm);
327
329 friend std::istream& operator>>(std::istream& data, TransmissionMatrix& tm);
330
332 template <int N> auto& TraMat(size_t i) noexcept {
333 static_assert (N > 0 and N < 5, "Bad size N");
334 if constexpr (N == 1) return T1[i];
335 else if constexpr (N == 2) return T2[i];
336 else if constexpr (N == 3) return T3[i];
337 else if constexpr (N == 4) return T4[i];
338 }
339
341 template <int N> [[nodiscard]] auto& TraMat(size_t i) const noexcept {
342 static_assert (N > 0 and N < 5, "Bad size N");
343 if constexpr (N == 1) return T1[i];
344 else if constexpr (N == 2) return T2[i];
345 else if constexpr (N == 3) return T3[i];
346 else if constexpr (N == 4) return T4[i];
347 }
348
350 template <int N> [[nodiscard]] Eigen::Matrix<Numeric, N, N> OptDepth(size_t i) const noexcept {
351 return (-TraMat<N>(i).diagonal().array().log().matrix()).asDiagonal();
352 }
353
373 template <int N>
374 [[nodiscard]] Eigen::Matrix<Numeric, N, 1> second_order_integration_source(const Eigen::Matrix<Numeric, N, N> T,
375 const Eigen::Matrix<Numeric, N, 1> far,
376 const Eigen::Matrix<Numeric, N, 1> close,
377 const Eigen::Matrix<Numeric, N, N> Kfar,
378 const Eigen::Matrix<Numeric, N, N> Kclose,
379 const Numeric r)
380 const noexcept {
381 static_assert (N>0 and N<5, "Bad size N");
382
384 if (T(0, 0) < 0.99) {
385 const auto od = 0.5 * r * (Kfar + Kclose);
386 Eigen::Matrix<Numeric, N, 1> second = od.inverse() * ((I - (I + od) * T) * far + (od - I + T) * close);
387 if ((far[0] < close[0] and Kfar(0, 0) > Kclose(0, 0)) or (far[0] > close[0] and Kfar(0, 0) < Kclose(0, 0))) {
388 Eigen::Matrix<Numeric, N, 1> second_limit = 0.5 * (I - T) * (far + close);
389
390 // FIXME: This is the equation given in the source material...
391 // const Eigen::Matrix<Numeric, N, 1> second_limit = 0.25 * r * (Kfar * far + Kclose * close);
392
393 if (second_limit[0] > second[0])
394 return second;
395 return second_limit;
396 }
397
398 return second;
399 }
400
401 return 0.5 * (I - T) * (far + close);
402 }
403
404
405 template <int N>
406 [[nodiscard]] Eigen::Matrix<Numeric, N, 1> second_order_integration_dsource(size_t i,
407 const TransmissionMatrix& dx,
408 const Eigen::Matrix<Numeric, N, 1> far,
409 const Eigen::Matrix<Numeric, N, 1> close,
410 const Eigen::Matrix<Numeric, N, 1> d,
411 bool isfar) const noexcept {
412 static_assert (N > 0 and N < 5, "Bad size N");
413
415 const auto T = TraMat<N>(i);
416 const auto& dTdx = dx.TraMat<N>(i);
417 const Eigen::Matrix<Numeric, N, 1> first = 0.5 * (I - T) * (far + close);
418 if (T(0, 0) < 0.99) {
419 const auto od = OptDepth<N>(i);
420 const auto doddx = dx.OptDepth<N>(i);
421 const Eigen::Matrix<Numeric, N, 1> second = od.inverse() * ((I - (I + od) * T) * far + (od - I + T) * close);
422 if (first[0] > second[0]) {
423 if (isfar) {
424 return od.inverse() * ((I - (I + od) * T) * d
425 - (I + od) * dTdx * far
426 + (doddx + T) * close
427 - doddx * second);
428 }
429
430 return od.inverse() * (
431 - (I + od) * dTdx * far
432 + (od - I + T) * d
433 + (doddx + T) * close
434 - doddx * second);
435 }
436
437 return 0.5 * ((I - T) * d - dTdx * (far + close));
438 }
439
440 return 0.5 * ((I - T) * d - dTdx * (far + close));
441 }
442};
443
451 const Numeric& x) {
452 return LazyScale<TransmissionMatrix>(tm, x);
453}
454
462 const TransmissionMatrix& tm) {
463 return LazyScale<TransmissionMatrix>(tm, x);
464}
465
469 std::vector<Eigen::Vector4d, Eigen::aligned_allocator<Eigen::Vector4d>> R4;
470 std::vector<Eigen::Vector3d, Eigen::aligned_allocator<Eigen::Vector3d>> R3;
471 std::vector<Eigen::Vector2d, Eigen::aligned_allocator<Eigen::Vector2d>> R2;
472 std::vector<Eigen::Matrix<double, 1, 1>,
473 Eigen::aligned_allocator<Eigen::Matrix<double, 1, 1>>>
475
481 RadiationVector(Index nf = 0, Index stokes = 1)
482 : stokes_dim(stokes),
483 R4(stokes_dim == 4 ? nf : 0, Eigen::Vector4d::Zero()),
484 R3(stokes_dim == 3 ? nf : 0, Eigen::Vector3d::Zero()),
485 R2(stokes_dim == 2 ? nf : 0, Eigen::Vector2d::Zero()),
486 R1(stokes_dim == 1 ? nf : 0, Eigen::Matrix<double, 1, 1>::Zero()) {
487 ARTS_ASSERT(stokes_dim < 5 and stokes_dim > 0);
488 }
489
495 : stokes_dim(rv.stokes_dim),
496 R4(std::move(rv.R4)),
497 R3(std::move(rv.R3)),
498 R2(std::move(rv.R2)),
499 R1(std::move(rv.R1)) {}
500
505 RadiationVector(const RadiationVector& rv) = default;
506
513
520 stokes_dim = rv.stokes_dim;
521 R4 = std::move(rv.R4);
522 R3 = std::move(rv.R3);
523 R2 = std::move(rv.R2);
524 R1 = std::move(rv.R1);
525 return *this;
526 }
527
533 for (size_t i = 0; i < R4.size(); i++) R4[i] = T.Mat4(i) * R4[i];
534 for (size_t i = 0; i < R3.size(); i++) R3[i] = T.Mat3(i) * R3[i];
535 for (size_t i = 0; i < R2.size(); i++) R2[i] = T.Mat2(i) * R2[i];
536 for (size_t i = 0; i < R1.size(); i++) R1[i] = T.Mat1(i) * R1[i];
537 }
538
543 void SetZero(size_t i) {
544 switch (stokes_dim) {
545 case 4:
546 R4[i].noalias() = Eigen::Vector4d::Zero();
547 break;
548 case 3:
549 R3[i].noalias() = Eigen::Vector3d::Zero();
550 break;
551 case 2:
552 R2[i].noalias() = Eigen::Vector2d::Zero();
553 break;
554 case 1:
555 R1[i][0] = 0;
556 break;
557 }
558 }
559
565 [[nodiscard]] const Eigen::Vector4d& Vec4(size_t i) const { return R4[i]; }
566
572 [[nodiscard]] const Eigen::Vector3d& Vec3(size_t i) const { return R3[i]; }
573
579 [[nodiscard]] const Eigen::Vector2d& Vec2(size_t i) const { return R2[i]; }
580
586 [[nodiscard]] const Eigen::Matrix<double, 1, 1>& Vec1(size_t i) const { return R1[i]; }
587
593 [[nodiscard]] Eigen::VectorXd Vec(size_t i) const {
594 switch (stokes_dim) {
595 case 4:
596 return Vec4(i);
597 case 3:
598 return Vec3(i);
599 case 2:
600 return Vec2(i);
601 default:
602 return Vec1(i);
603 }
604 }
605
606
612 Eigen::Vector4d& Vec4(size_t i) { return R4[i]; }
613
619 Eigen::Vector3d& Vec3(size_t i) { return R3[i]; }
620
626 Eigen::Vector2d& Vec2(size_t i) { return R2[i]; }
627
633 Eigen::Matrix<double, 1, 1>& Vec1(size_t i) { return R1[i]; }
634
640 void rem_avg(const RadiationVector& O1, const RadiationVector& O2) {
641 for (size_t i = 0; i < R4.size(); i++)
642 R4[i].noalias() -= 0.5 * (O1.R4[i] + O2.R4[i]);
643 for (size_t i = 0; i < R3.size(); i++)
644 R3[i].noalias() -= 0.5 * (O1.R3[i] + O2.R3[i]);
645 for (size_t i = 0; i < R2.size(); i++)
646 R2[i].noalias() -= 0.5 * (O1.R2[i] + O2.R2[i]);
647 for (size_t i = 0; i < R1.size(); i++)
648 R1[i].noalias() -= 0.5 * (O1.R1[i] + O2.R1[i]);
649 }
650
651
657 void add_avg(const RadiationVector& O1, const RadiationVector& O2) {
658 for (size_t i = 0; i < R4.size(); i++)
659 R4[i].noalias() += 0.5 * (O1.R4[i] + O2.R4[i]);
660 for (size_t i = 0; i < R3.size(); i++)
661 R3[i].noalias() += 0.5 * (O1.R3[i] + O2.R3[i]);
662 for (size_t i = 0; i < R2.size(); i++)
663 R2[i].noalias() += 0.5 * (O1.R2[i] + O2.R2[i]);
664 for (size_t i = 0; i < R1.size(); i++)
665 R1[i].noalias() += 0.5 * (O1.R1[i] + O2.R1[i]);
666 }
667
674 void add_weighted(const TransmissionMatrix& T, const RadiationVector& far, const RadiationVector& close,
675 const ConstMatrixView& Kfar, const ConstMatrixView& Kclose, const Numeric r) {
676 for (size_t i = 0; i < R4.size(); i++) {
677 R4[i].noalias() += T.second_order_integration_source<4>(T.TraMat<4>(i), far.R4[i], close.R4[i], prop_matrix<4>(Kfar(i, joker)), prop_matrix<4>(Kclose(i, joker)), r);
678 }
679 for (size_t i = 0; i < R3.size(); i++) {
680 R3[i].noalias() += T.second_order_integration_source<3>(T.TraMat<3>(i), far.R3[i], close.R3[i], prop_matrix<3>(Kfar(i, joker)), prop_matrix<3>(Kclose(i, joker)), r);
681 }
682 for (size_t i = 0; i < R2.size(); i++) {
683 R2[i].noalias() += T.second_order_integration_source<2>(T.TraMat<2>(i), far.R2[i], close.R2[i], prop_matrix<2>(Kfar(i, joker)), prop_matrix<2>(Kclose(i, joker)), r);
684 }
685 for (size_t i = 0; i < R1.size(); i++) {
686 R1[i].noalias() += T.second_order_integration_source<1>(T.TraMat<1>(i), far.R1[i], close.R1[i], prop_matrix<1>(Kfar(i, joker)), prop_matrix<1>(Kclose(i, joker)), r);
687 }
688 }
689
699 const TransmissionMatrix& dT,
700 const TransmissionMatrix& T,
701 const RadiationVector& ImJ,
702 const RadiationVector& dJ) {
703 for (size_t i = 0; i < R4.size(); i++)
704 R4[i].noalias() += PiT.Mat4(i) * (dT.Mat4(i) * ImJ.R4[i] + dJ.R4[i] -
705 T.Mat4(i) * dJ.R4[i]);
706 for (size_t i = 0; i < R3.size(); i++)
707 R3[i].noalias() += PiT.Mat3(i) * (dT.Mat3(i) * ImJ.R3[i] + dJ.R3[i] -
708 T.Mat3(i) * dJ.R3[i]);
709 for (size_t i = 0; i < R2.size(); i++)
710 R2[i].noalias() += PiT.Mat2(i) * (dT.Mat2(i) * ImJ.R2[i] + dJ.R2[i] -
711 T.Mat2(i) * dJ.R2[i]);
712 for (size_t i = 0; i < R1.size(); i++)
713 R1[i].noalias() += PiT.Mat1(i) * (dT.Mat1(i) * ImJ.R1[i] + dJ.R1[i] -
714 T.Mat1(i) * dJ.R1[i]);
715 }
716
726 const TransmissionMatrix& dT,
727 const TransmissionMatrix& T,
728 const RadiationVector& I,
729 const RadiationVector& far,
730 const RadiationVector& close,
731 const RadiationVector& d,
732 bool isfar) {
733 for (size_t i = 0; i < R4.size(); i++)
734 R4[i].noalias() += PiT.Mat4(i) * (dT.Mat4(i) * I.R4[i] + T.second_order_integration_dsource<4>(i, dT, far.R4[i], close.R4[i], d.R4[i], isfar));
735 for (size_t i = 0; i < R3.size(); i++)
736 R3[i].noalias() += PiT.Mat3(i) * (dT.Mat3(i) * I.R3[i] + T.second_order_integration_dsource<3>(i, dT, far.R3[i], close.R3[i], d.R3[i], isfar));
737 for (size_t i = 0; i < R2.size(); i++)
738 R2[i].noalias() += PiT.Mat2(i) * (dT.Mat2(i) * I.R2[i] + T.second_order_integration_dsource<2>(i, dT, far.R2[i], close.R2[i], d.R2[i], isfar));
739 for (size_t i = 0; i < R1.size(); i++)
740 R1[i].noalias() += PiT.Mat1(i) * (dT.Mat1(i) * I.R1[i] + T.second_order_integration_dsource<1>(i, dT, far.R1[i], close.R1[i], d.R1[i], isfar));
741 }
742
743
751 const TransmissionMatrix& dT,
752 const RadiationVector& I) {
753 for (size_t i = 0; i < R4.size(); i++)
754 R4[i].noalias() += PiT.Mat4(i) * dT.Mat4(i) * I.R4[i];
755 for (size_t i = 0; i < R3.size(); i++)
756 R3[i].noalias() += PiT.Mat3(i) * dT.Mat3(i) * I.R3[i];
757 for (size_t i = 0; i < R2.size(); i++)
758 R2[i].noalias() += PiT.Mat2(i) * dT.Mat2(i) * I.R2[i];
759 for (size_t i = 0; i < R1.size(); i++)
760 R1[i].noalias() += PiT.Mat1(i) * dT.Mat1(i) * I.R1[i];
761 }
762
763
774 const RadiationVector& x) {
775 for (size_t i = 0; i < R4.size(); i++)
776 R4[i].noalias() += A.Mat4(i) * x.R4[i];
777 for (size_t i = 0; i < R3.size(); i++)
778 R3[i].noalias() += A.Mat3(i) * x.R3[i];
779 for (size_t i = 0; i < R2.size(); i++)
780 R2[i].noalias() += A.Mat2(i) * x.R2[i];
781 for (size_t i = 0; i < R1.size(); i++)
782 R1[i].noalias() += A.Mat1(i) * x.R1[i];
783 }
784
793 const TransmissionMatrix& PiT,
794 const TransmissionMatrix& Z,
795 const TransmissionMatrix& dZ) {
796 for (size_t i = 0; i < R4.size(); i++)
797 R4[i] = PiT.Mat4(i) * (Z.Mat4(i) * R4[i] + dZ.Mat4(i) * I.R4[i]);
798 for (size_t i = 0; i < R3.size(); i++)
799 R3[i] = PiT.Mat3(i) * (Z.Mat3(i) * R3[i] + dZ.Mat3(i) * I.R3[i]);
800 for (size_t i = 0; i < R2.size(); i++)
801 R2[i] = PiT.Mat2(i) * (Z.Mat2(i) * R2[i] + dZ.Mat2(i) * I.R2[i]);
802 for (size_t i = 0; i < R1.size(); i++)
803 R4[i][0] = PiT.Mat1(i)[0] *
804 (Z.Mat1(i)[0] * R1[i][0] + dZ.Mat1(i)[0] * I.R1[i][0]);
805 }
806
815 const TransmissionMatrix& Tr,
816 const TransmissionMatrix& Tf,
817 const TransmissionMatrix& Z) {
818 for (size_t i = 0; i < R4.size(); i++)
819 R4[i].noalias() = Tr.Mat4(i) * Z.Mat4(i) * Tf.Mat4(i) * I0.R4[i];
820 for (size_t i = 0; i < R3.size(); i++)
821 R3[i].noalias() = Tr.Mat3(i) * Z.Mat3(i) * Tf.Mat3(i) * I0.R3[i];
822 for (size_t i = 0; i < R2.size(); i++)
823 R2[i].noalias() = Tr.Mat2(i) * Z.Mat2(i) * Tf.Mat2(i) * I0.R2[i];
824 for (size_t i = 0; i < R1.size(); i++)
825 R1[i].noalias() = Tr.Mat1(i) * Z.Mat1(i) * Tf.Mat1(i) * I0.R1[i];
826 }
827
836 const TransmissionMatrix& Tr,
837 const TransmissionMatrix& Tf,
838 const TransmissionMatrix& dZ) {
839 for (size_t i = 0; i < R4.size(); i++)
840 R4[i].noalias() += Tr.Mat4(i) * dZ.Mat4(i) * Tf.Mat4(i) * I0.R4[i];
841 for (size_t i = 0; i < R3.size(); i++)
842 R3[i].noalias() += Tr.Mat3(i) * dZ.Mat3(i) * Tf.Mat3(i) * I0.R3[i];
843 for (size_t i = 0; i < R2.size(); i++)
844 R2[i].noalias() += Tr.Mat2(i) * dZ.Mat2(i) * Tf.Mat2(i) * I0.R2[i];
845 for (size_t i = 0; i < R1.size(); i++)
846 R1[i].noalias() += Tr.Mat1(i) * dZ.Mat1(i) * Tf.Mat1(i) * I0.R1[i];
847 }
848
855 ARTS_ASSERT(M.ncols() == stokes_dim and M.nrows() == Frequencies());
856 for (size_t i = 0; i < R4.size(); i++) {
857 R4[i][0] = M(i, 0);
858 R4[i][1] = M(i, 1);
859 R4[i][2] = M(i, 2);
860 R4[i][3] = M(i, 3);
861 }
862 for (size_t i = 0; i < R3.size(); i++) {
863 R3[i][0] = M(i, 0);
864 R3[i][1] = M(i, 1);
865 R3[i][2] = M(i, 2);
866 }
867 for (size_t i = 0; i < R2.size(); i++) {
868 R2[i][0] = M(i, 0);
869 R2[i][1] = M(i, 1);
870 }
871 for (size_t i = 0; i < R1.size(); i++) {
872 R1[i][0] = M(i, 0);
873 }
874 return *this;
875 }
876
883 const Numeric& operator()(const Index i, const Index j) const {
884 switch (stokes_dim) {
885 case 4:
886 return R4[i][j];
887 case 3:
888 return R3[i][j];
889 case 2:
890 return R2[i][j];
891 default:
892 return R1[i][j];
893 }
894 }
895
900 operator Matrix() const {
902 for (size_t i = 0; i < R4.size(); i++)
903 for (size_t j = 0; j < 4; j++) M(i, j) = R4[i](j);
904 for (size_t i = 0; i < R3.size(); i++)
905 for (size_t j = 0; j < 3; j++) M(i, j) = R3[i](j);
906 for (size_t i = 0; i < R2.size(); i++)
907 for (size_t j = 0; j < 2; j++) M(i, j) = R2[i](j);
908 for (size_t i = 0; i < R1.size(); i++) M(i, 0) = R1[i](0);
909 return M;
910 }
911
920 const ConstVectorView& B,
921 const StokesVector& S,
922 Index i) {
923 ARTS_ASSERT(a.NumberOfAzimuthAngles() == 1);
924 ARTS_ASSERT(a.NumberOfZenithAngles() == 1);
927 switch (stokes_dim) {
928 case 4:
929 if (not S.IsEmpty())
930 R4[i].noalias() =
931 Eigen::Vector4d(a.Kjj()[i], a.K12()[i], a.K13()[i], a.K14()[i]) *
932 B[i] +
933 Eigen::Vector4d(S.Kjj()[i], S.K12()[i], S.K13()[i], S.K14()[i]);
934 else
935 R4[i].noalias() =
936 Eigen::Vector4d(a.Kjj()[i], a.K12()[i], a.K13()[i], a.K14()[i]) *
937 B[i];
938 break;
939 case 3:
940 if (not S.IsEmpty())
941 R3[i].noalias() =
942 Eigen::Vector3d(a.Kjj()[i], a.K12()[i], a.K13()[i]) * B[i] +
943 Eigen::Vector3d(S.Kjj()[i], S.K12()[i], S.K13()[i]);
944 else
945 R3[i].noalias() =
946 Eigen::Vector3d(a.Kjj()[i], a.K12()[i], a.K13()[i]) * B[i];
947 break;
948 case 2:
949 if (not S.IsEmpty())
950 R2[i].noalias() = Eigen::Vector2d(a.Kjj()[i], a.K12()[i]) * B[i] +
951 Eigen::Vector2d(S.Kjj()[i], S.K12()[i]);
952 else
953 R2[i].noalias() = Eigen::Vector2d(a.Kjj()[i], a.K12()[i]) * B[i];
954 break;
955 default:
956 if (not S.IsEmpty())
957 R1[i][0] = a.Kjj()[i] * B[i] + S.Kjj()[i];
958 else
959 R1[i][0] = a.Kjj()[i] * B[i];
960 }
961 }
962
964 [[nodiscard]] Index Frequencies() const {
965 switch (stokes_dim) {
966 case 4:
967 return Index(R4.size());
968 case 3:
969 return Index(R3.size());
970 case 2:
971 return Index(R2.size());
972 default:
973 return Index(R1.size());
974 }
975 }
976
978 friend std::ostream& operator<<(std::ostream& os, const RadiationVector& rv);
979
981 friend std::istream& operator>>(std::istream& data, RadiationVector& rv);
982};
983
990
992std::ostream& operator<<(std::ostream& os, const TransmissionMatrix& tm);
993
995std::ostream& operator<<(std::ostream& os,
996 const ArrayOfTransmissionMatrix& atm);
997
999std::ostream& operator<<(std::ostream& os,
1001
1003std::ostream& operator<<(std::ostream& os, const RadiationVector& rv);
1004
1006std::ostream& operator<<(std::ostream& os, const ArrayOfRadiationVector& arv);
1007
1009std::ostream& operator<<(std::ostream& os,
1010 const ArrayOfArrayOfRadiationVector& aarv);
1011
1013std::istream& operator>>(std::istream& is, TransmissionMatrix& tm);
1014
1016std::istream& operator>>(std::istream& is, RadiationVector& rv);
1017
1022};
1023
1026 Forward,
1027 Reverse,
1028};
1029
1032 Emission,
1035};
1036
1055 const RadiationVector& J1,
1056 const RadiationVector& J2,
1057 const ArrayOfRadiationVector& dJ1,
1058 const ArrayOfRadiationVector& dJ2,
1059 const TransmissionMatrix& T,
1060 const TransmissionMatrix& PiT,
1061 const ArrayOfTransmissionMatrix& dT1,
1062 const ArrayOfTransmissionMatrix& dT2,
1063 const PropagationMatrix& K1,
1064 const PropagationMatrix& K2,
1065 const ArrayOfPropagationMatrix& dK1,
1066 const ArrayOfPropagationMatrix& dK2,
1067 const Numeric r,
1068 const Vector& dr1,
1069 const Vector& dr2,
1070 const Index ia,
1071 const Index iz,
1072 const RadiativeTransferSolver solver);
1073
1091 const PropagationMatrix& K,
1092 const StokesVector& a,
1093 const StokesVector& S,
1094 const ArrayOfPropagationMatrix& dK,
1095 const ArrayOfStokesVector& da,
1096 const ArrayOfStokesVector& dS,
1097 const ConstVectorView& B,
1098 const ConstVectorView& dB_dT,
1099 const ArrayOfRetrievalQuantity& jacobian_quantities,
1100 const bool& jacobian_do);
1101
1119 const PropagationMatrix& K1,
1120 const PropagationMatrix& K2,
1121 const ArrayOfPropagationMatrix& dK1,
1122 const ArrayOfPropagationMatrix& dK2,
1123 const Numeric& r,
1124 const Numeric& dr_dtemp1,
1125 const Numeric& dr_dtemp2,
1126 const Index temp_deriv_pos);
1127
1136 const CumulativeTransmission type) /*[[expects: T.nelem()>0]]*/;
1137
1155 const RadiationVector & I_incoming,
1157 const ArrayOfTransmissionMatrix& PiTf,
1158 const ArrayOfTransmissionMatrix& PiTr,
1163 const BackscatterSolver solver);
1164
1177 const ConstMatrixView& pnd);
1178
1189 const ConstTensor5View& Pe, const ArrayOfMatrix& dpnd_dx);
1190
1191#endif // transmissionmatrix_h
This can be used to make arrays out of anything.
Definition: array.h:108
A constant view of a Matrix.
Definition: matpackI.h:1050
A constant view of a Tensor5.
Definition: matpackV.h:141
A constant view of a Vector.
Definition: matpackI.h:517
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.
The Matrix class.
Definition: matpackI.h:1270
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 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 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.
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3 class.
Definition: matpackIII.h:344
The Vector class.
Definition: matpackI.h:908
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
Routines for setting up the jacobian.
#define dx
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1484
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
Matrix matrix(const std::array< Numeric, 196 > data)
Create the square matrix from the static data.
Definition: igrf13.cc:206
invlib::MatrixIdentity< Matrix > Identity
Definition: oem.h:38
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:33
Stuff related to the propagation matrix.
#define d
#define a
#define N
Definition: rng.cc:164
#define M
Definition: rng.cc:165
Radiation Vector for Stokes dimension 1-4.
std::vector< Eigen::Matrix< double, 1, 1 >, Eigen::aligned_allocator< Eigen::Matrix< double, 1, 1 > > > R1
RadiationVector(RadiationVector &&rv) noexcept
Construct a new Radiation Vector object.
RadiationVector(Index nf=0, Index stokes=1)
Construct a new Radiation Vector object.
Eigen::Vector2d & Vec2(size_t i)
Return Vector at position.
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.
Eigen::Matrix< double, 1, 1 > & Vec1(size_t i)
Return Vector at position.
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::Vector2d, Eigen::aligned_allocator< Eigen::Vector2d > > R2
void add_avg(const RadiationVector &O1, const RadiationVector &O2)
Add the average of two other RadiationVector to *this.
RadiationVector & operator=(RadiationVector &&rv) noexcept
Assign old radiation vector to this.
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.
friend std::ostream & operator<<(std::ostream &os, const RadiationVector &rv)
Output operator.
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.
std::vector< Eigen::Vector4d, Eigen::aligned_allocator< Eigen::Vector4d > > R4
const Eigen::Vector2d & Vec2(size_t i) const
Return Vector at position.
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.
Eigen::Vector4d & Vec4(size_t i)
Return Vector at position.
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 ConstMatrixView &M)
Set *this from matrix.
RadiationVector & operator=(const RadiationVector &rv)=default
Assign old radiation vector to this.
std::vector< Eigen::Vector3d, Eigen::aligned_allocator< Eigen::Vector3d > > R3
Eigen::Vector3d & Vec3(size_t i)
Return Vector at position.
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.
friend std::istream & operator>>(std::istream &data, RadiationVector &rv)
Input operator.
void setBackscatterTransmissionDerivative(const RadiationVector &I0, const TransmissionMatrix &Tr, const TransmissionMatrix &Tf, const TransmissionMatrix &dZ)
Set this to backscatter transmission scatter derivative.
RadiationVector(const RadiationVector &rv)=default
Construct a new Radiation Vector object.
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
TransmissionMatrix & operator=(const LazyScale< TransmissionMatrix > &lstm)
Assign lazily.
TransmissionMatrix(TransmissionMatrix &&tm) noexcept
Construct a new Transmission Matrix object.
Eigen::Matrix3d & Mat3(size_t i)
Get Matrix at position.
std::vector< Eigen::Matrix3d, Eigen::aligned_allocator< Eigen::Matrix3d > > T3
void mul_aliased(const TransmissionMatrix &A, const TransmissionMatrix &B)
Set this to a multiple of A by B.
Eigen::Matrix< Numeric, N, N > OptDepth(size_t i) const noexcept
Simple template access for the optical depth.
TransmissionMatrix & operator*=(const Numeric &scale)
Scale self.
Eigen::Matrix4d & Mat4(size_t i)
Get Matrix at position.
const Eigen::Matrix2d & Mat2(size_t i) const
Get Matrix at position.
TransmissionMatrix(const TransmissionMatrix &tm)=default
Construct a new Transmission Matrix object.
std::vector< Eigen::Matrix< double, 1, 1 >, Eigen::aligned_allocator< Eigen::Matrix< double, 1, 1 > > > T1
TransmissionMatrix & operator=(TransmissionMatrix &&tm) noexcept
Move operator.
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
auto & TraMat(size_t i) const noexcept
Simple template access for the transmission.
void mul(const TransmissionMatrix &A, const TransmissionMatrix &B)
Set this to a multiple of A by B.
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.
friend std::istream & operator>>(std::istream &data, TransmissionMatrix &tm)
Input operator.
const Eigen::Matrix3d & Mat3(size_t i) const
Get Matrix at position.
friend std::ostream & operator<<(std::ostream &os, const TransmissionMatrix &tm)
Output operator.
Numeric operator()(const Index i, const Index j, const Index k) const
Access value in matrix.
std::vector< Eigen::Matrix2d, Eigen::aligned_allocator< Eigen::Matrix2d > > T2
TransmissionMatrix & operator=(const TransmissionMatrix &tm)=default
Assignment operator.
TransmissionMatrix & operator+=(const LazyScale< TransmissionMatrix > &lstm)
Assign to *this lazily.
void setZero()
Set to zero matrix.
Eigen::Matrix< double, 1, 1 > & Mat1(size_t i)
Get Matrix at position.
TransmissionMatrix(Index nf=0, Index stokes=1)
Construct a new Transmission Matrix object.
Eigen::Matrix2d & Mat2(size_t i)
Get Matrix at position.
std::vector< Eigen::Matrix4d, Eigen::aligned_allocator< Eigen::Matrix4d > > T4
void stepwise_source(RadiationVector &J, ArrayOfRadiationVector &dJ, const PropagationMatrix &K, const StokesVector &a, const StokesVector &S, const ArrayOfPropagationMatrix &dK, const ArrayOfStokesVector &da, const ArrayOfStokesVector &dS, const ConstVectorView &B, const ConstVectorView &dB_dT, const ArrayOfRetrievalQuantity &jacobian_quantities, const bool &jacobian_do)
Set the stepwise source.
RadiativeTransferSolver
Intended to hold various forward solvers.
ArrayOfTransmissionMatrix bulk_backscatter(const ConstTensor5View &Pe, const ConstMatrixView &pnd)
Bulk back-scattering.
BackscatterSolver
Intended to hold various backscatter solvers.
CumulativeTransmission
Intended to hold various ways to accumulate the transmission matrix.
ArrayOfArrayOfTransmissionMatrix bulk_backscatter_derivative(const ConstTensor5View &Pe, const ArrayOfMatrix &dpnd_dx)
Derivatives of bulk back-scattering
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.
std::istream & operator>>(std::istream &is, TransmissionMatrix &tm)
Input operator.
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.
std::ostream & operator<<(std::ostream &os, const TransmissionMatrix &tm)
Output operator.
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.