ARTS 2.5.0 (git: 9ee3ac6c)
propagationmatrix.h
Go to the documentation of this file.
1/* Copyright (C) 2017
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
34#ifndef propagationmatrix_h
35#define propagationmatrix_h
36
37#include "matpack_complex.h"
38#include "matpackIV.h"
39
40template <int N>
41Eigen::Matrix<Numeric, N, N> prop_matrix(const ConstVectorView& vec) {
42#define a vec[0]
43#define b vec[1]
44#define c vec[2]
45#define d vec[3]
46#define u vec[N]
47#define v vec[5]
48#define w vec[6]
49 static_assert (N>0 and N<5, "Bad size N");
50 if constexpr (N == 1) {
51 return Eigen::Matrix<Numeric, 1, 1>(a);
52 } else if constexpr (N==2) {
53 return (Eigen::Matrix2d() << a, b, b, a).finished();
54 } else if constexpr (N==3) {
55 return (Eigen::Matrix3d() << a, b, c, b, a, u, c, -u, a).finished();
56 } else if constexpr (N==4) {
57 return (Eigen::Matrix4d() << a, b, c, d, b, a, u, v, c, -u, a, w, d, -v, -w, a).finished();
58 }
59#undef a
60#undef b
61#undef c
62#undef d
63#undef u
64#undef v
65#undef w
66}
67
68template <int N>
69Eigen::Matrix<Numeric, N, N> prop_matrix(const ConstMatrixView& m) {
70 static_assert (N>0 and N<5, "Bad size N");
71 if constexpr (N == 1) {
72 return Eigen::Matrix<Numeric, 1, 1>(m(0, 0));
73 } else if constexpr (N == 2) {
74 return (Eigen::Matrix2d() << m(0, 0), m(0, 1),
75 m(1, 0), m(1, 1)).finished();
76 } else if constexpr (N == 3) {
77 return (Eigen::Matrix3d() << m(0, 0), m(0, 1), m(0, 2),
78 m(1, 0), m(1, 1), m(1, 2),
79 m(2, 0), m(2, 1), m(2, 2)) .finished();
80 } else if constexpr (N == 4) {
81 return (Eigen::Matrix4d() << m(0, 0), m(0, 1), m(0, 2), m(0, 3),
82 m(1, 0), m(1, 1), m(1, 2), m(1, 3),
83 m(2, 0), m(2, 1), m(2, 2), m(2, 3),
84 m(3, 0), m(3, 1), m(3, 2), m(3, 3)) .finished();
85 }
86}
87
88template <int N>
89Eigen::Matrix<Numeric, N, N> inv_prop_matrix(const ConstVectorView& vec) {
90#define a vec[0]
91#define b vec[1]
92#define c vec[2]
93#define d vec[3]
94#define u vec[N]
95#define v vec[5]
96#define w vec[6]
97 static_assert (N>0 and N<5, "Bad size N");
98 if constexpr (N == 1) {
99 return Eigen::Matrix<double, 1, 1>(1 / a);
100 } else if constexpr (N==2) {
101 return (Eigen::Matrix2d() << a, -b, -b, a).finished() / (a * a - b * b);
102 } else if constexpr (N==3) {
103 return (Eigen::Matrix3d() << a * a + u * u,
104 -a * b - c * u,
105 -a * c + b * u,
106 -a * b + c * u,
107 a * a - c * c,
108 -a * u + b * c,
109 -a * c - b * u,
110 a * u + b * c,
111 a * a - b * b)
112 .finished() /
113 (a * a * a - a * b * b - a * c * c + a * u * u);
114 } else if constexpr (N==4) {
115 return (Eigen::Matrix4d() << a * a * a + a * u * u + a * v * v + a * w * w,
116 -a * a * b - a * c * u - a * d * v - b * w * w + c * v * w -
117 d * u * w,
118 -a * a * c + a * b * u - a * d * w + b * v * w - c * v * v +
119 d * u * v,
120 -a * a * d + a * b * v + a * c * w - b * u * w + c * u * v -
121 d * u * u,
122 -a * a * b + a * c * u + a * d * v - b * w * w + c * v * w -
123 d * u * w,
124 a * a * a - a * c * c - a * d * d + a * w * w,
125 -a * a * u + a * b * c - a * v * w + b * d * w - c * d * v +
126 d * d * u,
127 -a * a * v + a * b * d + a * u * w - b * c * w + c * c * v -
128 c * d * u,
129 -a * a * c - a * b * u + a * d * w + b * v * w - c * v * v +
130 d * u * v,
131 a * a * u + a * b * c - a * v * w - b * d * w + c * d * v - d * d * u,
132 a * a * a - a * b * b - a * d * d + a * v * v,
133 -a * a * w + a * c * d - a * u * v + b * b * w - b * c * v +
134 b * d * u,
135 -a * a * d - a * b * v - a * c * w - b * u * w + c * u * v -
136 d * u * u,
137 a * a * v + a * b * d + a * u * w + b * c * w - c * c * v + c * d * u,
138 a * a * w + a * c * d - a * u * v - b * b * w + b * c * v - b * d * u,
139 a * a * a - a * b * b - a * c * c + a * u * u)
140 .finished() /
141 (a * a * a * a - a * a * b * b - a * a * c * c - a * a * d * d +
142 a * a * u * u + a * a * v * v + a * a * w * w - b * b * w * w +
143 2 * b * c * v * w - 2 * b * d * u * w - c * c * v * v +
144 2 * c * d * u * v - d * d * u * u);
145 }
146#undef a
147#undef b
148#undef c
149#undef d
150#undef u
151#undef v
152#undef w
153}
154
160template <class base>
162 public:
163
169 LazyScale(const base& t, const Numeric& x) : bas(t), scale(x) {}
170
172 const base& bas;
173
176};
177
178template<bool matrix>
179constexpr Index need2stokes(Index nstokes_needed) {
180 if constexpr (matrix) {
181 if (nstokes_needed == 7) return 4;
182 if (nstokes_needed == 4) return 3;
183 if (nstokes_needed == 2) return 2;
184 if (nstokes_needed == 1) return 1;
185 }
186
187 if (nstokes_needed < 5 and nstokes_needed > 0) return nstokes_needed;
188
189 ARTS_ASSERT (false, "Cannot understand the input Stokes dimensions");
191}
192
193template<bool matrix>
194constexpr Index stokes2need(Index nstokes) {
195 if constexpr (matrix) {
196 if (nstokes == 4) return 7;
197 if (nstokes == 3) return 4;
198 if (nstokes == 2) return 2;
199 if (nstokes == 1) return 1;
200 }
201
202 if (nstokes < 5 and nstokes > 0) return nstokes;
203
204 ARTS_ASSERT (false, "Cannot understand the input Stokes dimensions");
206}
207
233 public:
245 PropagationMatrix(const Index nr_frequencies = 0,
246 const Index stokes_dim = 1,
247 const Index nr_za = 1,
248 const Index nr_aa = 1,
249 const Numeric v = 0.0)
250 : mfreqs(nr_frequencies),
251 mstokes_dim(stokes_dim),
252 mza(nr_za),
253 maa(nr_aa) {
254 ARTS_ASSERT(mstokes_dim < 5 and mstokes_dim > 0);
256 }
257
263
269 : mfreqs(pm.mfreqs),
270 mstokes_dim(pm.mstokes_dim),
271 mza(pm.mza),
272 maa(pm.maa),
273 mdata(std::move(pm.mdata)),
274 mvectortype(pm.mvectortype) {}
275
281 : mfreqs(x.nrows()),
282 mza(x.npages()),
283 maa(x.nbooks()),
284 mdata(std::move(x)) {
285 switch (x.ncols()) {
286 case 7:
287 mstokes_dim = 4;
288 break;
289 case 4:
290 mstokes_dim = 3;
291 break;
292 case 2:
293 mstokes_dim = 2;
294 break;
295 case 1:
296 mstokes_dim = 1;
297 break;
298 default:
299 ARTS_ASSERT(false, "Tensor4 not representative of PropagationMatrix");
300 }
301 }
302
308 explicit PropagationMatrix(const ConstMatrixView& x, const bool& assume_fit = false)
309 : mfreqs(1), mstokes_dim(x.ncols()), mza(1), maa(1) {
310 ARTS_ASSERT(mstokes_dim < 5 and mstokes_dim > 0);
311 mvectortype = false;
312
313 if (not assume_fit) {
315 "Matrix not fit as propagation matrix");
316 }
317
319
320 switch (mstokes_dim) {
321 case 4:
322 mdata(0, 0, 0, 5) = x(1, 3);
323 mdata(0, 0, 0, 6) = x(2, 3);
324 mdata(0, 0, 0, 3) = x(0, 3); /* FALLTHROUGH */
325 case 3:
326 mdata(0, 0, 0, mstokes_dim) = x(1, 2);
327 mdata(0, 0, 0, 2) = x(0, 2); /* FALLTHROUGH */
328 case 2:
329 mdata(0, 0, 0, 1) = x(0, 1); /* FALLTHROUGH */
330 case 1:
331 mdata(0, 0, 0, 0) = x(0, 0); /* FALLTHROUGH */
332 }
333 };
334
336 [[nodiscard]] Index StokesDimensions() const { return mstokes_dim; };
337
339 [[nodiscard]] Index NumberOfFrequencies() const { return mfreqs; };
340
342 [[nodiscard]] Index NumberOfZenithAngles() const { return mza; };
343
345 [[nodiscard]] Index NumberOfAzimuthAngles() const { return maa; };
346
347 [[nodiscard]] bool OK() const {return mdata.ncols() == NumberOfNeededVectors() and mdata.nrows() == mfreqs and mdata.npages() == mza and mdata.nbooks() == maa;}
348
350 [[nodiscard]] bool IsEmpty() const { return not mfreqs or not mza or not maa; };
351
359 [[nodiscard]] bool IsZero(const Index iv = 0,
360 const Index iz = 0,
361 const Index ia = 0) const {
362 // FIXME: matpack does not do pointers in a clear manner
363 // return std::any_of(mdata(ia, iz, iv, joker).begin(), mdata(ia, iz, iv, joker).end(), [](auto& x){return x not_eq 0;});
364 for (auto& n : mdata(ia, iz, iv, joker))
365 if (n not_eq 0) return false;
366 return true;
367 };
368
376 [[nodiscard]] bool IsRotational(const Index iv = 0,
377 const Index iz = 0,
378 const Index ia = 0) const {
379 if (mdata(ia, iz, iv, 0) == 0.0)
380 return true;
381 return false;
382 };
383
385 #pragma GCC diagnostic push
386 #pragma GCC diagnostic ignored "-Wreturn-type"
387 [[nodiscard]] Index NumberOfNeededVectors() const {
388 if (not mvectortype) {
389 return stokes2need<true>(mstokes_dim);
390 }
391 return stokes2need<false>(mstokes_dim);
392 }
393 #pragma GCC diagnostic pop
394
396 Numeric operator()(const Index iv = 0,
397 const Index is1 = 0,
398 const Index is2 = 0,
399 const Index iz = 0,
400 const Index ia = 0) const;
401
411 void AddFaraday(const Numeric& rot,
412 const Index iv = 0,
413 const Index iz = 0,
414 const Index ia = 0) {
415 mdata(ia, iz, iv, mstokes_dim) += rot;
416 }
417
427 void SetFaraday(const Numeric& rot,
428 const Index iv = 0,
429 const Index iz = 0,
430 const Index ia = 0) {
431 mdata(ia, iz, iv, mstokes_dim) = rot;
432 }
433
442 const Index iv = 0,
443 const Index iz = 0,
444 const Index ia = 0) const;
445
452 if (this != &pm) {
453 mfreqs = pm.mfreqs;
454 mstokes_dim = pm.mstokes_dim;
455 mza = pm.mza;
456 maa = pm.maa;
457 mdata = std::move(pm.mdata);
458 mvectortype = pm.mvectortype;
459 }
460 return *this;
461 }
462
469 operator=(lpms.bas);
470 mdata *= lpms.scale;
471 return *this;
472 }
473
480
487 mdata = x;
488 return *this;
489 }
490
499 const Index iv = 0,
500 const Index iz = 0,
501 const Index ia = 0) {
502 mdata(ia, iz, iv, joker) = x.mdata(ia, iz, iv, joker);
503 }
504
512 void SetAtPosition(const ConstMatrixView& x,
513 const Index iv = 0,
514 const Index iz = 0,
515 const Index ia = 0);
516
517
525 void SetAtPosition(const Numeric& x,
526 const Index iv = 0,
527 const Index iz = 0,
528 const Index ia = 0) {
529 mdata(ia, iz, iv, joker) = x;
530 }
531
538 mdata /= other.mdata;
539 return *this;
540 }
541
548 for (Index i = 0; i < NumberOfNeededVectors(); i++) {
549 for (Index j = 0; j < mza; j++) {
550 for (Index k = 0; k < maa; k++) {
551 mdata(k, j, joker, i) /= x;
552 }
553 }
554 }
555 return *this;
556 }
557
564 mdata /= x;
565 return *this;
566 }
567
576 const Index iv = 0,
577 const Index iz = 0,
578 const Index ia = 0) {
579 mdata(ia, iz, iv, joker) /= x.mdata(ia, iz, iv, joker);
580 }
581
589 void DivideAtPosition(const ConstMatrixView& x,
590 const Index iv = 0,
591 const Index iz = 0,
592 const Index ia = 0);
593
602 const Index iv = 0,
603 const Index iz = 0,
604 const Index ia = 0) {
605 mdata(ia, iz, iv, joker) /= x;
606 }
607
614 mdata *= other.mdata;
615 return *this;
616 }
617
624 for (Index i = 0; i < NumberOfNeededVectors(); i++)
625 for (Index j = 0; j < mza; j++)
626 for (Index k = 0; k < maa; k++) mdata(k, j, joker, i) *= x;
627 return *this;
628 }
629
636 mdata *= x;
637 return *this;
638 }
639
648 const Index iv = 0,
649 const Index iz = 0,
650 const Index ia = 0) {
651 mdata(ia, iz, iv, joker) *= x.mdata(ia, iz, iv, joker);
652 }
653
662 const Index iv = 0,
663 const Index iz = 0,
664 const Index ia = 0);
665
674 const Index iv = 0,
675 const Index iz = 0,
676 const Index ia = 0) {
677 mdata(ia, iz, iv, joker) *= x;
678 }
679
686 mdata += other.mdata;
687 return *this;
688 }
689
696 MultiplyAndAdd(lpms.scale, lpms.bas);
697 return *this;
698 }
699
706 for (Index i = 0; i < NumberOfNeededVectors(); i++)
707 for (Index j = 0; j < mza; j++)
708 for (Index k = 0; k < maa; k++) mdata(k, j, joker, i) += x;
709 return *this;
710 }
711
718 mdata += x;
719 return *this;
720 }
721
730 const Index iv = 0,
731 const Index iz = 0,
732 const Index ia = 0) {
733 mdata(ia, iz, iv, joker) += x.mdata(ia, iz, iv, joker);
734 }
735
743 void AddAtPosition(const ConstMatrixView& x,
744 const Index iv = 0,
745 const Index iz = 0,
746 const Index ia = 0);
747
755 void AddAtPosition(const Numeric& x,
756 const Index iv = 0,
757 const Index iz = 0,
758 const Index ia = 0) {
759 mdata(ia, iz, iv, joker) += x;
760 }
761
768 mdata -= other.mdata;
769 return *this;
770 }
771
778 for (Index i = 0; i < NumberOfNeededVectors(); i++) {
779 for (Index j = 0; j < mza; j++) {
780 for (Index k = 0; k < maa; k++) {
781 mdata(k, j, joker, i) -= x;
782 }
783 }
784 }
785 return *this;
786 }
787
794 mdata -= x;
795 return *this;
796 }
797
806 const Index iv = 0,
807 const Index iz = 0,
808 const Index ia = 0) {
809 mdata(ia, iz, iv, joker) -= x.mdata(ia, iz, iv, joker);
810 }
811
819 void RemoveAtPosition(const ConstMatrixView& x,
820 const Index iv = 0,
821 const Index iz = 0,
822 const Index ia = 0);
823
832 const Index iv = 0,
833 const Index iz = 0,
834 const Index ia = 0) {
835 mdata(ia, iz, iv, joker) -= x;
836 }
837
846 const Index iv = 0,
847 const Index iz = 0,
848 const Index ia = 0) {
849 for (Index i = 0; i < mstokes_dim; i++) mdata(ia, iz, iv, i) += x[i];
850 }
851
860 void AddAverageAtPosition(const ConstMatrixView& mat1,
861 const ConstMatrixView& mat2,
862 const Index iv = 0,
863 const Index iz = 0,
864 const Index ia = 0);
865
871 void MultiplyAndAdd(const Numeric x, const PropagationMatrix& y);
872
881 const Index iv = 0,
882 const Index iz = 0,
883 const Index ia = 0) const;
884
889 [[nodiscard]] bool FittingShape(const ConstMatrixView& x) const;
890
897 void GetTensor3(Tensor3View tensor3, const Index iz = 0, const Index ia = 0);
898
905 VectorView Kjj(const Index iz = 0, const Index ia = 0) {
906 return mdata(ia, iz, joker, 0);
907 }
908
915 VectorView K12(const Index iz = 0, const Index ia = 0) {
916 return mdata(ia, iz, joker, 1);
917 }
918
925 VectorView K13(const Index iz = 0, const Index ia = 0) {
926 return mdata(ia, iz, joker, 2);
927 }
928
935 VectorView K14(const Index iz = 0, const Index ia = 0) {
936 return mdata(ia, iz, joker, 3);
937 }
938
945 VectorView K23(const Index iz = 0, const Index ia = 0) {
946 return mdata(ia, iz, joker, mstokes_dim);
947 }
948
955 VectorView K24(const Index iz = 0, const Index ia = 0) {
956 return mdata(ia, iz, joker, 5);
957 }
958
965 VectorView K34(const Index iz = 0, const Index ia = 0) {
966 return mdata(ia, iz, joker, 6);
967 }
968
975 [[nodiscard]] ConstVectorView Kjj(const Index iz = 0, const Index ia = 0) const {
976 return mdata(ia, iz, joker, 0);
977 }
978
985 [[nodiscard]] ConstVectorView K12(const Index iz = 0, const Index ia = 0) const {
986 return mdata(ia, iz, joker, 1);
987 }
988
995 [[nodiscard]] ConstVectorView K13(const Index iz = 0, const Index ia = 0) const {
996 return mdata(ia, iz, joker, 2);
997 }
998
1005 [[nodiscard]] ConstVectorView K14(const Index iz = 0, const Index ia = 0) const {
1006 return mdata(ia, iz, joker, 3);
1007 }
1008
1015 [[nodiscard]] ConstVectorView K23(const Index iz = 0, const Index ia = 0) const {
1016 return mdata(ia, iz, joker, mstokes_dim);
1017 }
1018
1025 [[nodiscard]] ConstVectorView K24(const Index iz = 0, const Index ia = 0) const {
1026 return mdata(ia, iz, joker, 5);
1027 }
1028
1035 [[nodiscard]] ConstVectorView K34(const Index iz = 0, const Index ia = 0) const {
1036 return mdata(ia, iz, joker, 6);
1037 }
1038
1040 void SetZero() { mdata = 0.0; }
1041
1043 Tensor4& Data() { return mdata; }
1044
1046 [[nodiscard]] const Tensor4& Data() const { return mdata; }
1047
1057 const ConstMatrixView& in,
1058 const Index iv = 0,
1059 const Index iz = 0,
1060 const Index ia = 0) const;
1061
1062
1072 const ConstMatrixView& in,
1073 const Index iv = 0,
1074 const Index iz = 0,
1075 const Index ia = 0) const;
1076
1077 protected:
1081 bool mvectortype{false};
1082};
1083
1086
1104 const Numeric& r,
1105 const PropagationMatrix& upper_level,
1106 const PropagationMatrix& lower_level,
1107 const Index iz = 0,
1108 const Index ia = 0);
1109
1110
1125 MatrixView T,
1126 const Numeric& r,
1127 const PropagationMatrix& averaged_propagation_matrix,
1128 const Index iv,
1129 const Index iz = 0,
1130 const Index ia = 0);
1131
1156 Tensor3View T,
1157 Tensor4View dT_upper_level,
1158 Tensor4View dT_lower_level,
1159 const Numeric& r,
1160 const PropagationMatrix& upper_level,
1161 const PropagationMatrix& lower_level,
1162 const Array<PropagationMatrix>& dprop_mat_upper_level,
1163 const Array<PropagationMatrix>& dprop_mat_lower_level,
1164 const Numeric& dr_dTu = 0.0,
1165 const Numeric& dr_dTl = 0.0,
1166 const Index it = -1,
1167 const Index iz = 0,
1168 const Index ia = 0);
1169
1171std::ostream& operator<<(std::ostream& os, const PropagationMatrix& pm);
1172
1174std::ostream& operator<<(std::ostream& os, const ArrayOfPropagationMatrix& apm);
1175
1177std::ostream& operator<<(std::ostream& os,
1179
1181class StokesVector final : public PropagationMatrix {
1182 public:
1183
1195 StokesVector(const Index nr_frequencies = 0,
1196 const Index stokes_dim = 1,
1197 const Index nr_za = 1,
1198 const Index nr_aa = 1,
1199 const Numeric& v = 0.0) {
1200 mvectortype = true;
1201 mfreqs = nr_frequencies;
1202 mstokes_dim = stokes_dim;
1203 mza = nr_za;
1204 maa = nr_aa;
1205 ARTS_ASSERT(mstokes_dim < 5 and mstokes_dim > 0);
1207 };
1208
1213 explicit StokesVector(Tensor4 x) {
1214 mfreqs = x.nrows();
1215 mstokes_dim = x.ncols();
1216 mza = x.npages();
1217 maa = x.nbooks();
1218 mdata = std::move(x);
1219 mvectortype = true;
1220 ARTS_ASSERT (not (mstokes_dim > 4 or mstokes_dim < 1),
1221 "Tensor4 is bad for StokesVector");
1222 }
1223
1228 explicit StokesVector(const ConstVectorView& x) {
1229 mfreqs = 1;
1230 mstokes_dim = x.nelem();
1231 mza = 1;
1232 maa = 1;
1233 ARTS_ASSERT(mstokes_dim < 5 and mstokes_dim > 0);
1234 mvectortype = true;
1235 mdata.resize(1, 1, 1, mstokes_dim);
1236 for (Index i = 0; i < mstokes_dim; i++) mdata(0, 0, 0, i) = x[i];
1237 };
1238
1246 const StokesVector& b,
1247 const Numeric& scale = 0.5) {
1248 mfreqs = a.NumberOfFrequencies();
1249 mstokes_dim = a.StokesDimensions();
1250 mza = a.NumberOfZenithAngles();
1251 maa = a.NumberOfAzimuthAngles();
1253 mvectortype = true;
1254
1255 for (Index i = 0; i < maa; i++)
1256 for (Index j = 0; j < mza; j++)
1257 for (Index k = 0; k < mfreqs; k++) {
1258 switch (mstokes_dim) {
1259 case 4:
1260 mdata(i, j, k, 3) = (a.mdata(i, j, k, 3) + b.mdata(i, j, k, 3)) *
1261 scale; /* FALLTHROUGH */
1262 case 3:
1263 mdata(i, j, k, 2) = (a.mdata(i, j, k, 2) + b.mdata(i, j, k, 2)) *
1264 scale; /* FALLTHROUGH */
1265 case 2:
1266 mdata(i, j, k, 1) = (a.mdata(i, j, k, 1) + b.mdata(i, j, k, 1)) *
1267 scale; /* FALLTHROUGH */
1268 case 1:
1269 mdata(i, j, k, 0) = (a.mdata(i, j, k, 0) + b.mdata(i, j, k, 0)) *
1270 scale; /* FALLTHROUGH */
1271 }
1272 }
1273 };
1274
1276 [[nodiscard]] Index NumberOfNeededVectors() const { return mstokes_dim; }
1277
1284 mdata += x.Data()(joker, joker, joker, Range(0, mstokes_dim, 1));
1285 return *this;
1286 }
1287
1294 MultiplyAndAdd(lpms.scale, lpms.bas);
1295 return *this;
1296 }
1297
1310 return *this;
1311 }
1312
1319 operator=(lpms.bas);
1320 mdata *= lpms.scale;
1321 return *this;
1322 }
1323
1330 mdata = x;
1331 return *this;
1332 }
1333
1339 void MultiplyAndAdd(const Numeric x, const PropagationMatrix& y) {
1344
1345 const Tensor4& data = y.Data();
1346
1347 for (Index i = 0; i < maa; i++) {
1348 for (Index j = 0; j < mza; j++) {
1349 for (Index k = 0; k < mfreqs; k++) {
1350 switch (mstokes_dim) {
1351 case 4:
1352 mdata(i, j, k, 3) += x * data(i, j, k, 3); /* FALLTHROUGH */
1353 case 3:
1354 mdata(i, j, k, 2) += x * data(i, j, k, 2); /* FALLTHROUGH */
1355 case 2:
1356 mdata(i, j, k, 1) += x * data(i, j, k, 1); /* FALLTHROUGH */
1357 case 1:
1358 mdata(i, j, k, 0) += x * data(i, j, k, 0); /* FALLTHROUGH */
1359 }
1360 }
1361 }
1362 }
1363 }
1364
1373 const Index iz = 0,
1374 const Index ia = 0) {
1375 return mdata(ia, iz, iv, joker);
1376 }
1377
1385 [[nodiscard]] ConstVectorView VectorAtPosition(const Index iv = 0,
1386 const Index iz = 0,
1387 const Index ia = 0) const {
1388 return mdata(ia, iz, iv, joker);
1389 }
1390
1392 const Index iv = 0,
1393 const Index iz = 0,
1394 const Index ia = 0) {
1395 mdata(ia, iz, iv, joker) = x;
1396 }
1397
1407 const ConstVectorView& vec2,
1408 const Index iv = 0,
1409 const Index iz = 0,
1410 const Index ia = 0) {
1411 switch (mstokes_dim) {
1412 case 4:
1413 mdata(ia, iz, iv, 3) += (vec1[3] + vec2[3]) * 0.5; /* FALLTHROUGH */
1414 case 3:
1415 mdata(ia, iz, iv, 2) += (vec1[2] + vec2[2]) * 0.5; /* FALLTHROUGH */
1416 case 2:
1417 mdata(ia, iz, iv, 1) += (vec1[1] + vec2[1]) * 0.5; /* FALLTHROUGH */
1418 case 1:
1419 mdata(ia, iz, iv, 0) += (vec1[0] + vec2[0]) * 0.5; /* FALLTHROUGH */
1420 }
1421 }
1422
1431 [[nodiscard]] bool IsPolarized(const Index iv = 0,
1432 const Index iz = 0,
1433 const Index ia = 0) const {
1434 switch (mstokes_dim) {
1435 case 4:
1436 if (K14(iz, ia)[iv] not_eq 0.0) return true; /* FALLTHROUGH */
1437 case 3:
1438 if (K13(iz, ia)[iv] not_eq 0.0) return true; /* FALLTHROUGH */
1439 case 2:
1440 if (K12(iz, ia)[iv] not_eq 0.0) return true; /* FALLTHROUGH */
1441 }
1442 return false;
1443 }
1444
1453 [[nodiscard]] bool IsUnpolarized(const Index iv = 0,
1454 const Index iz = 0,
1455 const Index ia = 0) const {
1456 return not IsPolarized(iv, iz, ia);
1457 }
1458
1459 [[nodiscard]] bool allZeroes() const {
1460 for (Index i=0; i<maa; i++)
1461 for (Index j=0; j<mza; j++)
1462 for (Index k=0; k<mfreqs; k++)
1463 for (Index m=0; m<mstokes_dim; m++)
1464 if (mdata(i, j, k, m) not_eq 0)
1465 return false;
1466 return true;
1467 }
1468};
1469
1473
1474std::ostream& operator<<(std::ostream& os, const StokesVector& pm);
1475std::ostream& operator<<(std::ostream& os, const ArrayOfStokesVector& apm);
1476std::ostream& operator<<(std::ostream& os,
1477 const ArrayOfArrayOfStokesVector& aapm);
1478
1486 const Numeric& x) {
1487 return LazyScale<PropagationMatrix>(pm, x);
1488}
1489
1497 const PropagationMatrix& pm) {
1498 return LazyScale<PropagationMatrix>(pm, x);
1499}
1500
1501
1503template <class T>
1505 const Vector& f_grid,
1506 const Index sd = 1) noexcept {
1507 const Index nf = f_grid.nelem();
1508 for (auto& var : main) {
1509 const bool bad_stokes = sd > var.StokesDimensions();
1510 const bool bad_freq = nf not_eq var.NumberOfFrequencies();
1511 if (bad_freq or bad_stokes) return true;
1512 }
1513
1514 return false;
1515}
1516
1517#endif //propagationmatrix_h
Index nrows
Index nbooks
Index npages
void * data
Index ncols
This can be used to make arrays out of anything.
Definition: array.h:107
A constant view of a Matrix.
Definition: matpackI.h:1014
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:58
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:55
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:64
Index nrows() const
Returns the number of rows.
Definition: matpackIV.cc:61
A constant view of a Vector.
Definition: matpackI.h:489
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
Class to help with hidden temporary variables for operations of type Numeric times Class.
const base & bas
A const reference to a value.
LazyScale(const base &t, const Numeric &x)
Construct a new Lazy Scale object.
const Numeric & scale
A const reference to a Numeric.
The MatrixView class.
Definition: matpackI.h:1125
void DivideAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Divide at position.
PropagationMatrix & operator-=(const Numeric &x)
Subtraction operator.
Numeric operator()(const Index iv=0, const Index is1=0, const Index is2=0, const Index iz=0, const Index ia=0) const
access operator.
const Tensor4 & Data() const
Get full const view to data.
ConstVectorView K13(const Index iz=0, const Index ia=0) const
Vector view to K(0, 2) elements.
void RightMultiplyAtPosition(MatrixView out, const ConstMatrixView &in, const Index iv=0, const Index iz=0, const Index ia=0) const
Multiply the matrix input from the right of this at position.
void AddAbsorptionVectorAtPosition(const ConstVectorView &x, const Index iv=0, const Index iz=0, const Index ia=0)
Adds as a Stokes vector at position.
bool IsZero(const Index iv=0, const Index iz=0, const Index ia=0) const
False if any non-zeroes in internal Matrix representation.
Index NumberOfFrequencies() const
The number of frequencies of the propagation matrix.
PropagationMatrix(Tensor4 x)
Construct a new Propagation Matrix object.
void MultiplyAndAdd(const Numeric x, const PropagationMatrix &y)
Multiply input by scalar and add to this.
PropagationMatrix & operator-=(const ConstVectorView &x)
Subtraction operator.
ConstVectorView K24(const Index iz=0, const Index ia=0) const
Vector view to K(1, 3) elements.
VectorView K24(const Index iz=0, const Index ia=0)
Vector view to K(1, 3) elements.
Tensor4 & Data()
Get full view to data.
void MatrixAtPosition(MatrixView ret, const Index iv=0, const Index iz=0, const Index ia=0) const
Sets the dense matrix.
ConstVectorView K12(const Index iz=0, const Index ia=0) const
Vector view to K(0, 1) elements.
VectorView K13(const Index iz=0, const Index ia=0)
Vector view to K(0, 2) elements.
PropagationMatrix(const PropagationMatrix &pm)=default
Construct a new Propagation Matrix object.
PropagationMatrix & operator/=(const ConstVectorView &x)
Divide operator.
void SetFaraday(const Numeric &rot, const Index iv=0, const Index iz=0, const Index ia=0)
Sets the Faraday rotation to the PropagationMatrix at required position.
Index NumberOfNeededVectors() const
The number of required vectors to fill this PropagationMatrix.
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
ConstVectorView K23(const Index iz=0, const Index ia=0) const
Vector view to K(1, 3) elements.
void SetZero()
Sets all data to zero.
PropagationMatrix & operator=(const Numeric &x)
Sets all data to constant.
PropagationMatrix & operator+=(const LazyScale< PropagationMatrix > &lpms)
Addition operator.
Index StokesDimensions() const
The stokes dimension of the propagation matrix.
PropagationMatrix & operator+=(const ConstVectorView &x)
Addition operator.
PropagationMatrix & operator=(const PropagationMatrix &other)=default
Copy operator.
void AddFaraday(const Numeric &rot, const Index iv=0, const Index iz=0, const Index ia=0)
Adds the Faraday rotation to the PropagationMatrix at required position.
VectorView K34(const Index iz=0, const Index ia=0)
Vector view to K(2, 3) elements.
void SetAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Set the At Position object.
PropagationMatrix & operator*=(const PropagationMatrix &other)
Multiply operator.
Index NumberOfAzimuthAngles() const
The number of azimuth angles of the propagation matrix.
void AddAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Addition at position operator.
ConstVectorView Kjj(const Index iz=0, const Index ia=0) const
Vector view to diagonal elements.
void RemoveAtPosition(const Numeric &x, const Index iv=0, const Index iz=0, const Index ia=0)
Subtraction at position.
bool FittingShape(const ConstMatrixView &x) const
Tests of the input matrix fits Propagation Matrix style.
PropagationMatrix & operator=(const LazyScale< PropagationMatrix > &lpms)
Laze equal to opeartor.
ConstVectorView K14(const Index iz=0, const Index ia=0) const
Vector view to K(0, 3) elements.
PropagationMatrix & operator*=(const Numeric &x)
Multiply operator.
PropagationMatrix(PropagationMatrix &&pm) noexcept
Construct a new Propagation Matrix object.
void LeftMultiplyAtPosition(MatrixView out, const ConstMatrixView &in, const Index iv=0, const Index iz=0, const Index ia=0) const
Multiply the matrix input from the left of this at position.
void MultiplyAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Multiply operator at position.
PropagationMatrix & operator=(PropagationMatrix &&pm) noexcept
Move operator.
PropagationMatrix & operator/=(const Numeric &x)
Divide operator.
PropagationMatrix & operator+=(const Numeric &x)
Addition operator.
VectorView K12(const Index iz=0, const Index ia=0)
Vector view to K(0, 1) elements.
PropagationMatrix & operator+=(const PropagationMatrix &other)
Addition operator.
PropagationMatrix & operator/=(const PropagationMatrix &other)
Divide operator.
VectorView K23(const Index iz=0, const Index ia=0)
Vector view to K(1, 2) elements.
PropagationMatrix & operator*=(const ConstVectorView &x)
Multiply operator.
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.
void GetTensor3(Tensor3View tensor3, const Index iz=0, const Index ia=0)
Get a Tensor3 object from this.
PropagationMatrix & operator-=(const PropagationMatrix &other)
Subtraction operator.
ConstVectorView K34(const Index iz=0, const Index ia=0) const
Vector view to diagonal elements.
PropagationMatrix(const Index nr_frequencies=0, const Index stokes_dim=1, const Index nr_za=1, const Index nr_aa=1, const Numeric v=0.0)
Initialize variable sizes.
void RemoveAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Subtraction at position.
void MultiplyAtPosition(const Numeric &x, const Index iv=0, const Index iz=0, const Index ia=0)
Multiply operator at position.
void SetAtPosition(const Numeric &x, const Index iv=0, const Index iz=0, const Index ia=0)
Set the At Position object.
PropagationMatrix(const ConstMatrixView &x, const bool &assume_fit=false)
Initialize from single stokes_dim-by-stokes_dim matrix.
Index NumberOfZenithAngles() const
The number of zenith angles of the propagation matrix.
void AddAverageAtPosition(const ConstMatrixView &mat1, const ConstMatrixView &mat2, const Index iv=0, const Index iz=0, const Index ia=0)
Add the average of the two input at position.
void AddAtPosition(const Numeric &x, const Index iv=0, const Index iz=0, const Index ia=0)
Addition at position operator.
void DivideAtPosition(const Numeric &x, const Index iv=0, const Index iz=0, const Index ia=0)
Divide at position.
void MatrixInverseAtPosition(MatrixView ret, const Index iv=0, const Index iz=0, const Index ia=0) const
Return the matrix inverse at the position.
bool IsRotational(const Index iv=0, const Index iz=0, const Index ia=0) const
False if diagonal element is non-zero in internal Matrix representation.
The range class.
Definition: matpackI.h:165
Stokes vector is as Propagation matrix but only has 4 possible values.
bool IsPolarized(const Index iv=0, const Index iz=0, const Index ia=0) const
Checks if vector is polarized.
StokesVector & operator=(const LazyScale< PropagationMatrix > &lpms)
Set this lazily.
bool allZeroes() const
StokesVector & operator+=(const PropagationMatrix &x)
Addition operator.
StokesVector & operator+=(const LazyScale< PropagationMatrix > &lpms)
Addition operator.
void SetAtPosition(const ConstVectorView &x, const Index iv=0, const Index iz=0, const Index ia=0)
StokesVector & operator=(const PropagationMatrix &x)
Set *this from a Propagation matrix.
StokesVector(const ConstVectorView &x)
Construct a new Stokes Vector object.
VectorView VectorAtPosition(const Index iv=0, const Index iz=0, const Index ia=0)
Get a vectorview to the position.
ConstVectorView VectorAtPosition(const Index iv=0, const Index iz=0, const Index ia=0) const
Get a vectorview to the position.
StokesVector(const StokesVector &a, const StokesVector &b, const Numeric &scale=0.5)
Construct a new Stokes Vector as a scale between two others.
StokesVector(Tensor4 x)
Construct a new Propagation Matrix object.
bool IsUnpolarized(const Index iv=0, const Index iz=0, const Index ia=0) const
Checks if vector is polarized.
StokesVector(const Index nr_frequencies=0, const Index stokes_dim=1, const Index nr_za=1, const Index nr_aa=1, const Numeric &v=0.0)
Initialize variable sizes.
void MultiplyAndAdd(const Numeric x, const PropagationMatrix &y)
Add a scaled version of the input.
Index NumberOfNeededVectors() const
The number of required vectors to fill this StokesVector.
void AddAverageAtPosition(const ConstVectorView &vec1, const ConstVectorView &vec2, const Index iv=0, const Index iz=0, const Index ia=0)
Add the average of both inputs at position.
StokesVector & operator=(const Numeric &x)
Set this to constant value.
The Tensor3View class.
Definition: matpackIII.h:239
The Tensor4View class.
Definition: matpackIV.h:284
The Tensor4 class.
Definition: matpackIV.h:421
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1058
The VectorView class.
Definition: matpackI.h:626
The Vector class.
Definition: matpackI.h:876
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
int main(int argc, char **argv)
#define max(a, b)
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
Matrix matrix(const std::array< Numeric, 196 > data)
Create the square matrix from the static data.
Definition: igrf13.cc:206
VectorView var(VectorView var, const Vector &y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the variance of the ranged ys.
Definition: raw.cc:159
VectorView std(VectorView std, const Vector &y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the standard deviation of the ranged ys.
Definition: raw.cc:205
constexpr Index need2stokes(Index nstokes_needed)
#define u
void compute_transmission_matrix(Tensor3View T, const Numeric &r, const PropagationMatrix &upper_level, const PropagationMatrix &lower_level, const Index iz=0, const Index ia=0)
Compute the matrix exponent as the transmission matrix of this propagation matrix.
#define d
LazyScale< PropagationMatrix > operator*(const PropagationMatrix &pm, const Numeric &x)
Returns a lazy multiplier.
#define v
#define w
#define a
constexpr Index stokes2need(Index nstokes)
std::ostream & operator<<(std::ostream &os, const PropagationMatrix &pm)
output operator
void compute_transmission_matrix_and_derivative(Tensor3View T, Tensor4View dT_upper_level, Tensor4View dT_lower_level, const Numeric &r, const PropagationMatrix &upper_level, const PropagationMatrix &lower_level, const Array< PropagationMatrix > &dprop_mat_upper_level, const Array< PropagationMatrix > &dprop_mat_lower_level, const Numeric &dr_dTu=0.0, const Numeric &dr_dTl=0.0, const Index it=-1, const Index iz=0, const Index ia=0)
Eigen::Matrix< Numeric, N, N > inv_prop_matrix(const ConstVectorView &vec)
#define c
#define b
bool bad_propmat(const Array< T > &main, const Vector &f_grid, const Index sd=1) noexcept
Checks if a Propagation Matrix or something similar has good grids.
void compute_transmission_matrix_from_averaged_matrix_at_frequency(MatrixView T, const Numeric &r, const PropagationMatrix &averaged_propagation_matrix, const Index iv, const Index iz=0, const Index ia=0)
Compute the matrix exponent as the transmission matrix of this propagation matrix.
Eigen::Matrix< Numeric, N, N > prop_matrix(const ConstVectorView &vec)
#define N
Definition: rng.cc:164