ARTS 2.5.0 (git: 9ee3ac6c)
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
45 private:
47 std::vector<Eigen::Matrix4d, Eigen::aligned_allocator<Eigen::Matrix4d>> T4;
48 std::vector<Eigen::Matrix3d, Eigen::aligned_allocator<Eigen::Matrix3d>> T3;
49 std::vector<Eigen::Matrix2d, Eigen::aligned_allocator<Eigen::Matrix2d>> T2;
50 std::vector<Eigen::Matrix<double, 1, 1>,
51 Eigen::aligned_allocator<Eigen::Matrix<double, 1, 1>>>
53
54 public:
55
61 TransmissionMatrix(Index nf = 0, Index stokes = 1)
62 : stokes_dim(stokes),
63 T4(stokes_dim == 4 ? nf : 0, Eigen::Matrix4d::Identity()),
64 T3(stokes_dim == 3 ? nf : 0, Eigen::Matrix3d::Identity()),
65 T2(stokes_dim == 2 ? nf : 0, Eigen::Matrix2d::Identity()),
66 T1(stokes_dim == 1 ? nf : 0, Eigen::Matrix<double, 1, 1>::Identity()) {
67 ARTS_ASSERT(stokes_dim < 5 and stokes_dim > 0);
68 }
69
75 : stokes_dim(tm.stokes_dim),
76 T4(std::move(tm.T4)),
77 T3(std::move(tm.T3)),
78 T2(std::move(tm.T2)),
79 T1(std::move(tm.T1)) {}
80
86
93
100 stokes_dim = tm.stokes_dim;
101 T4 = std::move(tm.T4);
102 T3 = std::move(tm.T3);
103 T2 = std::move(tm.T2);
104 T1 = std::move(tm.T1);
105 return *this;
106 }
107
113 explicit TransmissionMatrix(const PropagationMatrix& pm, const Numeric& r = 1.0);
114
115 operator Tensor3() const {
117 for (size_t i = 0; i < T4.size(); i++)
118 for (size_t j = 0; j < 4; j++)
119 for (size_t k = 0; k < 4; k++) T(i, j, k) = T4[i](j, k);
120 for (size_t i = 0; i < T3.size(); i++)
121 for (size_t j = 0; j < 3; j++)
122 for (size_t k = 0; k < 3; k++) T(i, j, k) = T3[i](j, k);
123 for (size_t i = 0; i < T2.size(); i++)
124 for (size_t j = 0; j < 2; j++)
125 for (size_t k = 0; k < 2; k++) T(i, j, k) = T2[i](j, k);
126 for (size_t i = 0; i < T1.size(); i++) T(i, 0, 0) = T1[i](0, 0);
127 return T;
128 }
129
135 [[nodiscard]] const Eigen::Matrix4d& Mat4(size_t i) const { return T4[i]; }
136
142 [[nodiscard]] const Eigen::Matrix3d& Mat3(size_t i) const { return T3[i]; }
143
149 [[nodiscard]] const Eigen::Matrix2d& Mat2(size_t i) const { return T2[i]; }
150
156 [[nodiscard]] const Eigen::Matrix<double, 1, 1>& Mat1(size_t i) const { return T1[i]; }
157
163 #pragma GCC diagnostic push
164 #pragma GCC diagnostic ignored "-Wreturn-type"
165 [[nodiscard]] Eigen::MatrixXd Mat(size_t i) const {
166 switch (stokes_dim) {
167 case 1:
168 return Mat1(i);
169 case 2:
170 return Mat2(i);
171 case 3:
172 return Mat3(i);
173 case 4:
174 return Mat4(i);
175 }
176 }
177 #pragma GCC diagnostic pop
178
184 Eigen::Matrix4d& Mat4(size_t i) { return T4[i]; }
185
191 Eigen::Matrix3d& Mat3(size_t i) { return T3[i]; }
192
198 Eigen::Matrix2d& Mat2(size_t i) { return T2[i]; }
199
205 Eigen::Matrix<double, 1, 1>& Mat1(size_t i) { return T1[i]; }
206
208 void setIdentity() {
209 std::fill(T4.begin(), T4.end(), Eigen::Matrix4d::Identity());
210 std::fill(T3.begin(), T3.end(), Eigen::Matrix3d::Identity());
211 std::fill(T2.begin(), T2.end(), Eigen::Matrix2d::Identity());
212 std::fill(T1.begin(), T1.end(), Eigen::Matrix<double, 1, 1>::Identity());
213 }
214
216 void setZero() {
217 std::fill(T4.begin(), T4.end(), Eigen::Matrix4d::Zero());
218 std::fill(T3.begin(), T3.end(), Eigen::Matrix3d::Zero());
219 std::fill(T2.begin(), T2.end(), Eigen::Matrix2d::Zero());
220 std::fill(T1.begin(), T1.end(), Eigen::Matrix<double, 1, 1>::Zero());
221 }
222
230 void mul(const TransmissionMatrix& A, const TransmissionMatrix& B) {
231 for (size_t i = 0; i < T4.size(); i++) T4[i].noalias() = A.T4[i] * B.T4[i];
232 for (size_t i = 0; i < T3.size(); i++) T3[i].noalias() = A.T3[i] * B.T3[i];
233 for (size_t i = 0; i < T2.size(); i++) T2[i].noalias() = A.T2[i] * B.T2[i];
234 for (size_t i = 0; i < T1.size(); i++) T1[i].noalias() = A.T1[i] * B.T1[i];
235 }
236
237
246 for (size_t i = 0; i < T4.size(); i++) T4[i] = A.T4[i] * B.T4[i];
247 for (size_t i = 0; i < T3.size(); i++) T3[i] = A.T3[i] * B.T3[i];
248 for (size_t i = 0; i < T2.size(); i++) T2[i] = A.T2[i] * B.T2[i];
249 for (size_t i = 0; i < T1.size(); i++) T1[i] = A.T1[i] * B.T1[i];
250 }
251
259 Numeric operator()(const Index i, const Index j, const Index k) const {
260 switch (stokes_dim) {
261 case 4:
262 return T4[i](j, k);
263 case 3:
264 return T3[i](j, k);
265 case 2:
266 return T2[i](j, k);
267 default:
268 return T1[i](j, k);
269 }
270 }
271
273 [[nodiscard]] Index StokesDim() const { return stokes_dim; }
274
276 [[nodiscard]] Index Frequencies() const {
277 switch (stokes_dim) {
278 case 4:
279 return Index(T4.size());
280 case 3:
281 return Index(T3.size());
282 case 2:
283 return Index(T2.size());
284 default:
285 return Index(T1.size());
286 }
287 }
288
295 for (size_t i = 0; i < T4.size(); i++)
296 T4[i].noalias() = lstm.scale * lstm.bas.Mat4(i);
297 for (size_t i = 0; i < T3.size(); i++)
298 T3[i].noalias() = lstm.scale * lstm.bas.Mat3(i);
299 for (size_t i = 0; i < T2.size(); i++)
300 T2[i].noalias() = lstm.scale * lstm.bas.Mat2(i);
301 for (size_t i = 0; i < T1.size(); i++)
302 T1[i].noalias() = lstm.scale * lstm.bas.Mat1(i);
303 return *this;
304 }
305
312 std::transform(T4.begin(), T4.end(), T4.begin(), [scale](auto& T){return T*scale;});
313 std::transform(T3.begin(), T3.end(), T3.begin(), [scale](auto& T){return T*scale;});
314 std::transform(T2.begin(), T2.end(), T2.begin(), [scale](auto& T){return T*scale;});
315 std::transform(T1.begin(), T1.end(), T1.begin(), [scale](auto& T){return T*scale;});
316 return *this;
317 }
318
325 operator=(lstm.bas);
326 operator*=(lstm.scale);
327 return *this;
328 }
329
331 friend std::ostream& operator<<(std::ostream& os,
332 const TransmissionMatrix& tm);
333
335 friend std::istream& operator>>(std::istream& data, TransmissionMatrix& tm);
336
338 template <int N> auto& TraMat(size_t i) noexcept {
339 static_assert (N > 0 and N < 5, "Bad size N");
340 if constexpr (N == 1) return T1[i];
341 else if constexpr (N == 2) return T2[i];
342 else if constexpr (N == 3) return T3[i];
343 else if constexpr (N == 4) return T4[i];
344 }
345
347 template <int N> [[nodiscard]] auto& TraMat(size_t i) const noexcept {
348 static_assert (N > 0 and N < 5, "Bad size N");
349 if constexpr (N == 1) return T1[i];
350 else if constexpr (N == 2) return T2[i];
351 else if constexpr (N == 3) return T3[i];
352 else if constexpr (N == 4) return T4[i];
353 }
354
356 template <int N> [[nodiscard]] Eigen::Matrix<Numeric, N, N> OptDepth(size_t i) const noexcept {
357 return (-TraMat<N>(i).diagonal().array().log().matrix()).asDiagonal();
358 }
359
379 template <int N>
380 [[nodiscard]] Eigen::Matrix<Numeric, N, 1> second_order_integration_source(const Eigen::Matrix<Numeric, N, N> T,
381 const Eigen::Matrix<Numeric, N, 1> far,
382 const Eigen::Matrix<Numeric, N, 1> close,
383 const Eigen::Matrix<Numeric, N, N> Kfar,
384 const Eigen::Matrix<Numeric, N, N> Kclose,
385 const Numeric r)
386 const noexcept {
387 static_assert (N>0 and N<5, "Bad size N");
388
390 if (T(0, 0) < 0.99) {
391 const auto od = 0.5 * r * (Kfar + Kclose);
392 Eigen::Matrix<Numeric, N, 1> second = od.inverse() * ((I - (I + od) * T) * far + (od - I + T) * close);
393 if ((far[0] < close[0] and Kfar(0, 0) > Kclose(0, 0)) or (far[0] > close[0] and Kfar(0, 0) < Kclose(0, 0))) {
394 Eigen::Matrix<Numeric, N, 1> second_limit = 0.5 * (I - T) * (far + close);
395
396 // FIXME: This is the equation given in the source material...
397 // const Eigen::Matrix<Numeric, N, 1> second_limit = 0.25 * r * (Kfar * far + Kclose * close);
398
399 if (second_limit[0] > second[0])
400 return second;
401 return second_limit;
402 }
403
404 return second;
405 }
406
407 return 0.5 * (I - T) * (far + close);
408 }
409
410
411 template <int N>
412 [[nodiscard]] Eigen::Matrix<Numeric, N, 1> second_order_integration_dsource(size_t i,
413 const TransmissionMatrix& dx,
414 const Eigen::Matrix<Numeric, N, 1> far,
415 const Eigen::Matrix<Numeric, N, 1> close,
416 const Eigen::Matrix<Numeric, N, 1> d,
417 bool isfar) const noexcept {
418 static_assert (N > 0 and N < 5, "Bad size N");
419
421 const auto T = TraMat<N>(i);
422 const auto& dTdx = dx.TraMat<N>(i);
423 const Eigen::Matrix<Numeric, N, 1> first = 0.5 * (I - T) * (far + close);
424 if (T(0, 0) < 0.99) {
425 const auto od = OptDepth<N>(i);
426 const auto doddx = dx.OptDepth<N>(i);
427 const Eigen::Matrix<Numeric, N, 1> second = od.inverse() * ((I - (I + od) * T) * far + (od - I + T) * close);
428 if (first[0] > second[0]) {
429 if (isfar) {
430 return od.inverse() * ((I - (I + od) * T) * d
431 - (I + od) * dTdx * far
432 + (doddx + T) * close
433 - doddx * second);
434 }
435
436 return od.inverse() * (
437 - (I + od) * dTdx * far
438 + (od - I + T) * d
439 + (doddx + T) * close
440 - doddx * second);
441 }
442
443 return 0.5 * ((I - T) * d - dTdx * (far + close));
444 }
445
446 return 0.5 * ((I - T) * d - dTdx * (far + close));
447 }
448};
449
457 const Numeric& x) {
458 return LazyScale<TransmissionMatrix>(tm, x);
459}
460
468 const TransmissionMatrix& tm) {
469 return LazyScale<TransmissionMatrix>(tm, x);
470}
471
474 private:
476 std::vector<Eigen::Vector4d, Eigen::aligned_allocator<Eigen::Vector4d>> R4;
477 std::vector<Eigen::Vector3d, Eigen::aligned_allocator<Eigen::Vector3d>> R3;
478 std::vector<Eigen::Vector2d, Eigen::aligned_allocator<Eigen::Vector2d>> R2;
479 std::vector<Eigen::Matrix<double, 1, 1>,
480 Eigen::aligned_allocator<Eigen::Matrix<double, 1, 1>>>
482
483 public:
484
490 RadiationVector(Index nf = 0, Index stokes = 1)
491 : stokes_dim(stokes),
492 R4(stokes_dim == 4 ? nf : 0, Eigen::Vector4d::Zero()),
493 R3(stokes_dim == 3 ? nf : 0, Eigen::Vector3d::Zero()),
494 R2(stokes_dim == 2 ? nf : 0, Eigen::Vector2d::Zero()),
495 R1(stokes_dim == 1 ? nf : 0, Eigen::Matrix<double, 1, 1>::Zero()) {
496 ARTS_ASSERT(stokes_dim < 5 and stokes_dim > 0);
497 }
498
504 : stokes_dim(rv.stokes_dim),
505 R4(std::move(rv.R4)),
506 R3(std::move(rv.R3)),
507 R2(std::move(rv.R2)),
508 R1(std::move(rv.R1)) {}
509
514 RadiationVector(const RadiationVector& rv) = default;
515
522
529 stokes_dim = rv.stokes_dim;
530 R4 = std::move(rv.R4);
531 R3 = std::move(rv.R3);
532 R2 = std::move(rv.R2);
533 R1 = std::move(rv.R1);
534 return *this;
535 }
536
542 for (size_t i = 0; i < R4.size(); i++) R4[i] = T.Mat4(i) * R4[i];
543 for (size_t i = 0; i < R3.size(); i++) R3[i] = T.Mat3(i) * R3[i];
544 for (size_t i = 0; i < R2.size(); i++) R2[i] = T.Mat2(i) * R2[i];
545 for (size_t i = 0; i < R1.size(); i++) R1[i] = T.Mat1(i) * R1[i];
546 }
547
552 void SetZero(size_t i) {
553 switch (stokes_dim) {
554 case 4:
555 R4[i].noalias() = Eigen::Vector4d::Zero();
556 break;
557 case 3:
558 R3[i].noalias() = Eigen::Vector3d::Zero();
559 break;
560 case 2:
561 R2[i].noalias() = Eigen::Vector2d::Zero();
562 break;
563 case 1:
564 R1[i][0] = 0;
565 break;
566 }
567 }
568
574 [[nodiscard]] const Eigen::Vector4d& Vec4(size_t i) const { return R4[i]; }
575
581 [[nodiscard]] const Eigen::Vector3d& Vec3(size_t i) const { return R3[i]; }
582
588 [[nodiscard]] const Eigen::Vector2d& Vec2(size_t i) const { return R2[i]; }
589
595 [[nodiscard]] const Eigen::Matrix<double, 1, 1>& Vec1(size_t i) const { return R1[i]; }
596
602 [[nodiscard]] Eigen::VectorXd Vec(size_t i) const {
603 switch (stokes_dim) {
604 case 4:
605 return Vec4(i);
606 case 3:
607 return Vec3(i);
608 case 2:
609 return Vec2(i);
610 default:
611 return Vec1(i);
612 }
613 }
614
615
621 Eigen::Vector4d& Vec4(size_t i) { return R4[i]; }
622
628 Eigen::Vector3d& Vec3(size_t i) { return R3[i]; }
629
635 Eigen::Vector2d& Vec2(size_t i) { return R2[i]; }
636
642 Eigen::Matrix<double, 1, 1>& Vec1(size_t i) { return R1[i]; }
643
649 void rem_avg(const RadiationVector& O1, const RadiationVector& O2) {
650 for (size_t i = 0; i < R4.size(); i++)
651 R4[i].noalias() -= 0.5 * (O1.R4[i] + O2.R4[i]);
652 for (size_t i = 0; i < R3.size(); i++)
653 R3[i].noalias() -= 0.5 * (O1.R3[i] + O2.R3[i]);
654 for (size_t i = 0; i < R2.size(); i++)
655 R2[i].noalias() -= 0.5 * (O1.R2[i] + O2.R2[i]);
656 for (size_t i = 0; i < R1.size(); i++)
657 R1[i].noalias() -= 0.5 * (O1.R1[i] + O2.R1[i]);
658 }
659
660
666 void add_avg(const RadiationVector& O1, const RadiationVector& O2) {
667 for (size_t i = 0; i < R4.size(); i++)
668 R4[i].noalias() += 0.5 * (O1.R4[i] + O2.R4[i]);
669 for (size_t i = 0; i < R3.size(); i++)
670 R3[i].noalias() += 0.5 * (O1.R3[i] + O2.R3[i]);
671 for (size_t i = 0; i < R2.size(); i++)
672 R2[i].noalias() += 0.5 * (O1.R2[i] + O2.R2[i]);
673 for (size_t i = 0; i < R1.size(); i++)
674 R1[i].noalias() += 0.5 * (O1.R1[i] + O2.R1[i]);
675 }
676
683 void add_weighted(const TransmissionMatrix& T, const RadiationVector& far, const RadiationVector& close,
684 const ConstMatrixView& Kfar, const ConstMatrixView& Kclose, const Numeric r) {
685 for (size_t i = 0; i < R4.size(); i++) {
686 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);
687 }
688 for (size_t i = 0; i < R3.size(); i++) {
689 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);
690 }
691 for (size_t i = 0; i < R2.size(); i++) {
692 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);
693 }
694 for (size_t i = 0; i < R1.size(); i++) {
695 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);
696 }
697 }
698
708 const TransmissionMatrix& dT,
709 const TransmissionMatrix& T,
710 const RadiationVector& ImJ,
711 const RadiationVector& dJ) {
712 for (size_t i = 0; i < R4.size(); i++)
713 R4[i].noalias() += PiT.Mat4(i) * (dT.Mat4(i) * ImJ.R4[i] + dJ.R4[i] -
714 T.Mat4(i) * dJ.R4[i]);
715 for (size_t i = 0; i < R3.size(); i++)
716 R3[i].noalias() += PiT.Mat3(i) * (dT.Mat3(i) * ImJ.R3[i] + dJ.R3[i] -
717 T.Mat3(i) * dJ.R3[i]);
718 for (size_t i = 0; i < R2.size(); i++)
719 R2[i].noalias() += PiT.Mat2(i) * (dT.Mat2(i) * ImJ.R2[i] + dJ.R2[i] -
720 T.Mat2(i) * dJ.R2[i]);
721 for (size_t i = 0; i < R1.size(); i++)
722 R1[i].noalias() += PiT.Mat1(i) * (dT.Mat1(i) * ImJ.R1[i] + dJ.R1[i] -
723 T.Mat1(i) * dJ.R1[i]);
724 }
725
735 const TransmissionMatrix& dT,
736 const TransmissionMatrix& T,
737 const RadiationVector& I,
738 const RadiationVector& far,
739 const RadiationVector& close,
740 const RadiationVector& d,
741 bool isfar) {
742 for (size_t i = 0; i < R4.size(); i++)
743 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));
744 for (size_t i = 0; i < R3.size(); i++)
745 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));
746 for (size_t i = 0; i < R2.size(); i++)
747 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));
748 for (size_t i = 0; i < R1.size(); i++)
749 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));
750 }
751
752
760 const TransmissionMatrix& dT,
761 const RadiationVector& I) {
762 for (size_t i = 0; i < R4.size(); i++)
763 R4[i].noalias() += PiT.Mat4(i) * dT.Mat4(i) * I.R4[i];
764 for (size_t i = 0; i < R3.size(); i++)
765 R3[i].noalias() += PiT.Mat3(i) * dT.Mat3(i) * I.R3[i];
766 for (size_t i = 0; i < R2.size(); i++)
767 R2[i].noalias() += PiT.Mat2(i) * dT.Mat2(i) * I.R2[i];
768 for (size_t i = 0; i < R1.size(); i++)
769 R1[i].noalias() += PiT.Mat1(i) * dT.Mat1(i) * I.R1[i];
770 }
771
772
783 const RadiationVector& x) {
784 for (size_t i = 0; i < R4.size(); i++)
785 R4[i].noalias() += A.Mat4(i) * x.R4[i];
786 for (size_t i = 0; i < R3.size(); i++)
787 R3[i].noalias() += A.Mat3(i) * x.R3[i];
788 for (size_t i = 0; i < R2.size(); i++)
789 R2[i].noalias() += A.Mat2(i) * x.R2[i];
790 for (size_t i = 0; i < R1.size(); i++)
791 R1[i].noalias() += A.Mat1(i) * x.R1[i];
792 }
793
802 const TransmissionMatrix& PiT,
803 const TransmissionMatrix& Z,
804 const TransmissionMatrix& dZ) {
805 for (size_t i = 0; i < R4.size(); i++)
806 R4[i] = PiT.Mat4(i) * (Z.Mat4(i) * R4[i] + dZ.Mat4(i) * I.R4[i]);
807 for (size_t i = 0; i < R3.size(); i++)
808 R3[i] = PiT.Mat3(i) * (Z.Mat3(i) * R3[i] + dZ.Mat3(i) * I.R3[i]);
809 for (size_t i = 0; i < R2.size(); i++)
810 R2[i] = PiT.Mat2(i) * (Z.Mat2(i) * R2[i] + dZ.Mat2(i) * I.R2[i]);
811 for (size_t i = 0; i < R1.size(); i++)
812 R4[i][0] = PiT.Mat1(i)[0] *
813 (Z.Mat1(i)[0] * R1[i][0] + dZ.Mat1(i)[0] * I.R1[i][0]);
814 }
815
824 const TransmissionMatrix& Tr,
825 const TransmissionMatrix& Tf,
826 const TransmissionMatrix& Z) {
827 for (size_t i = 0; i < R4.size(); i++)
828 R4[i].noalias() = Tr.Mat4(i) * Z.Mat4(i) * Tf.Mat4(i) * I0.R4[i];
829 for (size_t i = 0; i < R3.size(); i++)
830 R3[i].noalias() = Tr.Mat3(i) * Z.Mat3(i) * Tf.Mat3(i) * I0.R3[i];
831 for (size_t i = 0; i < R2.size(); i++)
832 R2[i].noalias() = Tr.Mat2(i) * Z.Mat2(i) * Tf.Mat2(i) * I0.R2[i];
833 for (size_t i = 0; i < R1.size(); i++)
834 R1[i].noalias() = Tr.Mat1(i) * Z.Mat1(i) * Tf.Mat1(i) * I0.R1[i];
835 }
836
845 const TransmissionMatrix& Tr,
846 const TransmissionMatrix& Tf,
847 const TransmissionMatrix& dZ) {
848 for (size_t i = 0; i < R4.size(); i++)
849 R4[i].noalias() += Tr.Mat4(i) * dZ.Mat4(i) * Tf.Mat4(i) * I0.R4[i];
850 for (size_t i = 0; i < R3.size(); i++)
851 R3[i].noalias() += Tr.Mat3(i) * dZ.Mat3(i) * Tf.Mat3(i) * I0.R3[i];
852 for (size_t i = 0; i < R2.size(); i++)
853 R2[i].noalias() += Tr.Mat2(i) * dZ.Mat2(i) * Tf.Mat2(i) * I0.R2[i];
854 for (size_t i = 0; i < R1.size(); i++)
855 R1[i].noalias() += Tr.Mat1(i) * dZ.Mat1(i) * Tf.Mat1(i) * I0.R1[i];
856 }
857
864 ARTS_ASSERT(M.ncols() == stokes_dim and M.nrows() == Frequencies());
865 for (size_t i = 0; i < R4.size(); i++) {
866 R4[i][0] = M(i, 0);
867 R4[i][1] = M(i, 1);
868 R4[i][2] = M(i, 2);
869 R4[i][3] = M(i, 3);
870 }
871 for (size_t i = 0; i < R3.size(); i++) {
872 R3[i][0] = M(i, 0);
873 R3[i][1] = M(i, 1);
874 R3[i][2] = M(i, 2);
875 }
876 for (size_t i = 0; i < R2.size(); i++) {
877 R2[i][0] = M(i, 0);
878 R2[i][1] = M(i, 1);
879 }
880 for (size_t i = 0; i < R1.size(); i++) {
881 R1[i][0] = M(i, 0);
882 }
883 return *this;
884 }
885
892 const Numeric& operator()(const Index i, const Index j) const {
893 switch (stokes_dim) {
894 case 4:
895 return R4[i][j];
896 case 3:
897 return R3[i][j];
898 case 2:
899 return R2[i][j];
900 default:
901 return R1[i][j];
902 }
903 }
904
909 operator Matrix() const {
911 for (size_t i = 0; i < R4.size(); i++)
912 for (size_t j = 0; j < 4; j++) M(i, j) = R4[i](j);
913 for (size_t i = 0; i < R3.size(); i++)
914 for (size_t j = 0; j < 3; j++) M(i, j) = R3[i](j);
915 for (size_t i = 0; i < R2.size(); i++)
916 for (size_t j = 0; j < 2; j++) M(i, j) = R2[i](j);
917 for (size_t i = 0; i < R1.size(); i++) M(i, 0) = R1[i](0);
918 return M;
919 }
920
929 const ConstVectorView& B,
930 const StokesVector& S,
931 Index i) {
932 ARTS_ASSERT(a.NumberOfAzimuthAngles() == 1);
933 ARTS_ASSERT(a.NumberOfZenithAngles() == 1);
936 switch (stokes_dim) {
937 case 4:
938 if (not S.IsEmpty())
939 R4[i].noalias() =
940 Eigen::Vector4d(a.Kjj()[i], a.K12()[i], a.K13()[i], a.K14()[i]) *
941 B[i] +
942 Eigen::Vector4d(S.Kjj()[i], S.K12()[i], S.K13()[i], S.K14()[i]);
943 else
944 R4[i].noalias() =
945 Eigen::Vector4d(a.Kjj()[i], a.K12()[i], a.K13()[i], a.K14()[i]) *
946 B[i];
947 break;
948 case 3:
949 if (not S.IsEmpty())
950 R3[i].noalias() =
951 Eigen::Vector3d(a.Kjj()[i], a.K12()[i], a.K13()[i]) * B[i] +
952 Eigen::Vector3d(S.Kjj()[i], S.K12()[i], S.K13()[i]);
953 else
954 R3[i].noalias() =
955 Eigen::Vector3d(a.Kjj()[i], a.K12()[i], a.K13()[i]) * B[i];
956 break;
957 case 2:
958 if (not S.IsEmpty())
959 R2[i].noalias() = Eigen::Vector2d(a.Kjj()[i], a.K12()[i]) * B[i] +
960 Eigen::Vector2d(S.Kjj()[i], S.K12()[i]);
961 else
962 R2[i].noalias() = Eigen::Vector2d(a.Kjj()[i], a.K12()[i]) * B[i];
963 break;
964 default:
965 if (not S.IsEmpty())
966 R1[i][0] = a.Kjj()[i] * B[i] + S.Kjj()[i];
967 else
968 R1[i][0] = a.Kjj()[i] * B[i];
969 }
970 }
971
973 [[nodiscard]] Index StokesDim() const { return stokes_dim; }
974
976 [[nodiscard]] Index Frequencies() const {
977 switch (stokes_dim) {
978 case 4:
979 return Index(R4.size());
980 case 3:
981 return Index(R3.size());
982 case 2:
983 return Index(R2.size());
984 default:
985 return Index(R1.size());
986 }
987 }
988
990 friend std::ostream& operator<<(std::ostream& os, const RadiationVector& rv);
991
993 friend std::istream& operator>>(std::istream& data, RadiationVector& rv);
994};
995
1002
1004std::ostream& operator<<(std::ostream& os, const TransmissionMatrix& tm);
1005
1007std::ostream& operator<<(std::ostream& os,
1008 const ArrayOfTransmissionMatrix& atm);
1009
1011std::ostream& operator<<(std::ostream& os,
1013
1015std::ostream& operator<<(std::ostream& os, const RadiationVector& rv);
1016
1018std::ostream& operator<<(std::ostream& os, const ArrayOfRadiationVector& arv);
1019
1021std::ostream& operator<<(std::ostream& os,
1022 const ArrayOfArrayOfRadiationVector& aarv);
1023
1025std::istream& operator>>(std::istream& is, TransmissionMatrix& tm);
1026
1028std::istream& operator>>(std::istream& is, RadiationVector& rv);
1029
1034};
1035
1038 Forward,
1039 Reverse,
1040};
1041
1044 Emission,
1047};
1048
1067 const RadiationVector& J1,
1068 const RadiationVector& J2,
1069 const ArrayOfRadiationVector& dJ1,
1070 const ArrayOfRadiationVector& dJ2,
1071 const TransmissionMatrix& T,
1072 const TransmissionMatrix& PiT,
1073 const ArrayOfTransmissionMatrix& dT1,
1074 const ArrayOfTransmissionMatrix& dT2,
1075 const PropagationMatrix& K1,
1076 const PropagationMatrix& K2,
1077 const ArrayOfPropagationMatrix& dK1,
1078 const ArrayOfPropagationMatrix& dK2,
1079 const Numeric r,
1080 const Vector& dr1,
1081 const Vector& dr2,
1082 const Index ia,
1083 const Index iz,
1084 const RadiativeTransferSolver solver);
1085
1103 const PropagationMatrix& K,
1104 const StokesVector& a,
1105 const StokesVector& S,
1106 const ArrayOfPropagationMatrix& dK,
1107 const ArrayOfStokesVector& da,
1108 const ArrayOfStokesVector& dS,
1109 const ConstVectorView& B,
1110 const ConstVectorView& dB_dT,
1111 const ArrayOfRetrievalQuantity& jacobian_quantities,
1112 const bool& jacobian_do);
1113
1131 const PropagationMatrix& K1,
1132 const PropagationMatrix& K2,
1133 const ArrayOfPropagationMatrix& dK1,
1134 const ArrayOfPropagationMatrix& dK2,
1135 const Numeric& r,
1136 const Numeric& dr_dtemp1,
1137 const Numeric& dr_dtemp2,
1138 const Index temp_deriv_pos);
1139
1148 const CumulativeTransmission type) /*[[expects: T.nelem()>0]]*/;
1149
1167 const RadiationVector & I_incoming,
1169 const ArrayOfTransmissionMatrix& PiTf,
1170 const ArrayOfTransmissionMatrix& PiTr,
1175 const BackscatterSolver solver);
1176
1189 const ConstMatrixView& pnd);
1190
1201 const ConstTensor5View& Pe, const ArrayOfMatrix& dpnd_dx);
1202
1203#endif // transmissionmatrix_h
void * data
This can be used to make arrays out of anything.
Definition: array.h:107
A constant view of a Matrix.
Definition: matpackI.h:1014
A constant view of a Tensor5.
Definition: matpackV.h:143
A constant view of a Vector.
Definition: matpackI.h:489
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:1225
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.
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
Index StokesDim() const
Get Stokes dimension.
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.
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3 class.
Definition: matpackIII.h:339
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
TransmissionMatrix & operator=(const LazyScale< TransmissionMatrix > &lstm)
Assign lazily.
Index StokesDim() const
Stokes dimensionaility.
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
The Vector class.
Definition: matpackI.h:876
#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:1472
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
const Joker joker
char Type type
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
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.