ARTS 2.5.4 (git: 31ce4f0e)
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
45template <class base>
46class LazyScale {
47 public:
48
54 LazyScale(const base& t, const Numeric& x) : bas(t), scale(x) {}
55
57 const base& bas;
58
60 const Numeric& scale;
61};
62
63template<bool matrix>
64constexpr Index need2stokes(Index nstokes_needed) {
65 if constexpr (matrix) {
66 if (nstokes_needed == 7) return 4;
67 if (nstokes_needed == 4) return 3;
68 if (nstokes_needed == 2) return 2;
69 if (nstokes_needed == 1) return 1;
70 }
71
72 if (nstokes_needed < 5 and nstokes_needed > 0) return nstokes_needed;
73
74 ARTS_ASSERT (false, "Cannot understand the input Stokes dimensions");
76}
77
78template<bool matrix>
79constexpr Index stokes2need(Index nstokes) {
80 if constexpr (matrix) {
81 if (nstokes == 4) return 7;
82 if (nstokes == 3) return 4;
83 if (nstokes == 2) return 2;
84 if (nstokes == 1) return 1;
85 }
86
87 if (nstokes < 5 and nstokes > 0) return nstokes;
88
89 ARTS_ASSERT (false, "Cannot understand the input Stokes dimensions");
91}
92
118 public:
130 PropagationMatrix(const Index nr_frequencies = 0,
131 const Index stokes_dim = 1,
132 const Index nr_za = 1,
133 const Index nr_aa = 1,
134 const Numeric v = 0.0)
135 : mfreqs(nr_frequencies),
136 mstokes_dim(stokes_dim),
137 mza(nr_za),
138 maa(nr_aa) {
139 ARTS_ASSERT(mstokes_dim < 5 and mstokes_dim > 0);
141 }
142
148
154 : mfreqs(pm.mfreqs),
155 mstokes_dim(pm.mstokes_dim),
156 mza(pm.mza),
157 maa(pm.maa),
158 mdata(std::move(pm.mdata)),
159 mvectortype(pm.mvectortype) {}
160
166 : mfreqs(x.nrows()),
167 mza(x.npages()),
168 maa(x.nbooks()),
169 mdata(std::move(x)) {
170 switch (x.ncols()) {
171 case 7:
172 mstokes_dim = 4;
173 break;
174 case 4:
175 mstokes_dim = 3;
176 break;
177 case 2:
178 mstokes_dim = 2;
179 break;
180 case 1:
181 mstokes_dim = 1;
182 break;
183 default:
184 ARTS_ASSERT(false, "Tensor4 not representative of PropagationMatrix");
185 }
186 }
187
193 explicit PropagationMatrix(const ConstMatrixView& x, const bool& assume_fit = false)
194 : mfreqs(1), mstokes_dim(x.ncols()), mza(1), maa(1) {
195 ARTS_ASSERT(mstokes_dim < 5 and mstokes_dim > 0);
196 mvectortype = false;
197
198 if (not assume_fit) {
200 "Matrix not fit as propagation matrix");
201 }
202
204
205 switch (mstokes_dim) {
206 case 4:
207 mdata(0, 0, 0, 5) = x(1, 3);
208 mdata(0, 0, 0, 6) = x(2, 3);
209 mdata(0, 0, 0, 3) = x(0, 3); /* FALLTHROUGH */
210 case 3:
211 mdata(0, 0, 0, mstokes_dim) = x(1, 2);
212 mdata(0, 0, 0, 2) = x(0, 2); /* FALLTHROUGH */
213 case 2:
214 mdata(0, 0, 0, 1) = x(0, 1); /* FALLTHROUGH */
215 case 1:
216 mdata(0, 0, 0, 0) = x(0, 0); /* FALLTHROUGH */
217 }
218 };
219
221 [[nodiscard]] Index StokesDimensions() const { return mstokes_dim; };
222
224 [[nodiscard]] Index NumberOfFrequencies() const { return mfreqs; };
225
227 [[nodiscard]] Index NumberOfZenithAngles() const { return mza; };
228
230 [[nodiscard]] Index NumberOfAzimuthAngles() const { return maa; };
231
232 [[nodiscard]] bool OK() const {return mdata.ncols() == NumberOfNeededVectors() and mdata.nrows() == mfreqs and mdata.npages() == mza and mdata.nbooks() == maa;}
233
235 [[nodiscard]] bool IsEmpty() const { return not mfreqs or not mza or not maa; };
236
244 [[nodiscard]] bool IsZero(const Index iv = 0,
245 const Index iz = 0,
246 const Index ia = 0) const {
247 // FIXME: matpack does not do pointers in a clear manner
248 // return std::any_of(mdata(ia, iz, iv, joker).begin(), mdata(ia, iz, iv, joker).end(), [](auto& x){return x not_eq 0;});
249 for (auto& n : mdata(ia, iz, iv, joker))
250 if (n not_eq 0) return false;
251 return true;
252 };
253
261 [[nodiscard]] bool IsRotational(const Index iv = 0,
262 const Index iz = 0,
263 const Index ia = 0) const {
264 if (mdata(ia, iz, iv, 0) == 0.0)
265 return true;
266 return false;
267 };
268
270 #pragma GCC diagnostic push
271 #pragma GCC diagnostic ignored "-Wreturn-type"
272 [[nodiscard]] Index NumberOfNeededVectors() const {
273 if (not mvectortype) {
274 return stokes2need<true>(mstokes_dim);
275 }
276 return stokes2need<false>(mstokes_dim);
277 }
278 #pragma GCC diagnostic pop
279
281 Numeric operator()(const Index iv = 0,
282 const Index is1 = 0,
283 const Index is2 = 0,
284 const Index iz = 0,
285 const Index ia = 0) const;
286
296 void AddFaraday(const Numeric& rot,
297 const Index iv = 0,
298 const Index iz = 0,
299 const Index ia = 0) {
300 mdata(ia, iz, iv, mstokes_dim) += rot;
301 }
302
312 void SetFaraday(const Numeric& rot,
313 const Index iv = 0,
314 const Index iz = 0,
315 const Index ia = 0) {
316 mdata(ia, iz, iv, mstokes_dim) = rot;
317 }
318
327 const Index iv = 0,
328 const Index iz = 0,
329 const Index ia = 0) const;
330
337 if (this != &pm) {
338 mfreqs = pm.mfreqs;
339 mstokes_dim = pm.mstokes_dim;
340 mza = pm.mza;
341 maa = pm.maa;
342 mdata = std::move(pm.mdata);
343 mvectortype = pm.mvectortype;
344 }
345 return *this;
346 }
347
354 operator=(lpms.bas);
355 mdata *= lpms.scale;
356 return *this;
357 }
358
365
372 mdata = x;
373 return *this;
374 }
375
384 const Index iv = 0,
385 const Index iz = 0,
386 const Index ia = 0) {
387 mdata(ia, iz, iv, joker) = x.mdata(ia, iz, iv, joker);
388 }
389
397 void SetAtPosition(const ConstMatrixView& x,
398 const Index iv = 0,
399 const Index iz = 0,
400 const Index ia = 0);
401
402
410 void SetAtPosition(const Numeric& x,
411 const Index iv = 0,
412 const Index iz = 0,
413 const Index ia = 0) {
414 mdata(ia, iz, iv, joker) = x;
415 }
416
423 mdata /= other.mdata;
424 return *this;
425 }
426
433 for (Index i = 0; i < NumberOfNeededVectors(); i++) {
434 for (Index j = 0; j < mza; j++) {
435 for (Index k = 0; k < maa; k++) {
436 mdata(k, j, joker, i) /= x;
437 }
438 }
439 }
440 return *this;
441 }
442
449 mdata /= x;
450 return *this;
451 }
452
461 const Index iv = 0,
462 const Index iz = 0,
463 const Index ia = 0) {
464 mdata(ia, iz, iv, joker) /= x.mdata(ia, iz, iv, joker);
465 }
466
474 void DivideAtPosition(const ConstMatrixView& x,
475 const Index iv = 0,
476 const Index iz = 0,
477 const Index ia = 0);
478
487 const Index iv = 0,
488 const Index iz = 0,
489 const Index ia = 0) {
490 mdata(ia, iz, iv, joker) /= x;
491 }
492
499 mdata *= other.mdata;
500 return *this;
501 }
502
509 for (Index i = 0; i < NumberOfNeededVectors(); i++)
510 for (Index j = 0; j < mza; j++)
511 for (Index k = 0; k < maa; k++) mdata(k, j, joker, i) *= x;
512 return *this;
513 }
514
521 mdata *= x;
522 return *this;
523 }
524
533 const Index iv = 0,
534 const Index iz = 0,
535 const Index ia = 0) {
536 mdata(ia, iz, iv, joker) *= x.mdata(ia, iz, iv, joker);
537 }
538
547 const Index iv = 0,
548 const Index iz = 0,
549 const Index ia = 0);
550
559 const Index iv = 0,
560 const Index iz = 0,
561 const Index ia = 0) {
562 mdata(ia, iz, iv, joker) *= x;
563 }
564
571 mdata += other.mdata;
572 return *this;
573 }
574
581 MultiplyAndAdd(lpms.scale, lpms.bas);
582 return *this;
583 }
584
591 for (Index i = 0; i < NumberOfNeededVectors(); i++)
592 for (Index j = 0; j < mza; j++)
593 for (Index k = 0; k < maa; k++) mdata(k, j, joker, i) += x;
594 return *this;
595 }
596
603 mdata += x;
604 return *this;
605 }
606
615 const Index iv = 0,
616 const Index iz = 0,
617 const Index ia = 0) {
618 mdata(ia, iz, iv, joker) += x.mdata(ia, iz, iv, joker);
619 }
620
628 void AddAtPosition(const ConstMatrixView& x,
629 const Index iv = 0,
630 const Index iz = 0,
631 const Index ia = 0);
632
640 void AddAtPosition(const Numeric& x,
641 const Index iv = 0,
642 const Index iz = 0,
643 const Index ia = 0) {
644 mdata(ia, iz, iv, joker) += x;
645 }
646
653 mdata -= other.mdata;
654 return *this;
655 }
656
663 for (Index i = 0; i < NumberOfNeededVectors(); i++) {
664 for (Index j = 0; j < mza; j++) {
665 for (Index k = 0; k < maa; k++) {
666 mdata(k, j, joker, i) -= x;
667 }
668 }
669 }
670 return *this;
671 }
672
679 mdata -= x;
680 return *this;
681 }
682
691 const Index iv = 0,
692 const Index iz = 0,
693 const Index ia = 0) {
694 mdata(ia, iz, iv, joker) -= x.mdata(ia, iz, iv, joker);
695 }
696
704 void RemoveAtPosition(const ConstMatrixView& x,
705 const Index iv = 0,
706 const Index iz = 0,
707 const Index ia = 0);
708
717 const Index iv = 0,
718 const Index iz = 0,
719 const Index ia = 0) {
720 mdata(ia, iz, iv, joker) -= x;
721 }
722
731 const Index iv = 0,
732 const Index iz = 0,
733 const Index ia = 0) {
734 for (Index i = 0; i < mstokes_dim; i++) mdata(ia, iz, iv, i) += x[i];
735 }
736
745 void AddAverageAtPosition(const ConstMatrixView& mat1,
746 const ConstMatrixView& mat2,
747 const Index iv = 0,
748 const Index iz = 0,
749 const Index ia = 0);
750
756 void MultiplyAndAdd(const Numeric x, const PropagationMatrix& y);
757
766 const Index iv = 0,
767 const Index iz = 0,
768 const Index ia = 0) const;
769
774 [[nodiscard]] bool FittingShape(const ConstMatrixView& x) const;
775
782 void GetTensor3(Tensor3View tensor3, const Index iz = 0, const Index ia = 0);
783
790 VectorView Kjj(const Index iz = 0, const Index ia = 0) {
791 return mdata(ia, iz, joker, 0);
792 }
793
800 VectorView K12(const Index iz = 0, const Index ia = 0) {
801 return mdata(ia, iz, joker, 1);
802 }
803
810 VectorView K13(const Index iz = 0, const Index ia = 0) {
811 return mdata(ia, iz, joker, 2);
812 }
813
820 VectorView K14(const Index iz = 0, const Index ia = 0) {
821 return mdata(ia, iz, joker, 3);
822 }
823
830 VectorView K23(const Index iz = 0, const Index ia = 0) {
831 return mdata(ia, iz, joker, mstokes_dim);
832 }
833
840 VectorView K24(const Index iz = 0, const Index ia = 0) {
841 return mdata(ia, iz, joker, 5);
842 }
843
850 VectorView K34(const Index iz = 0, const Index ia = 0) {
851 return mdata(ia, iz, joker, 6);
852 }
853
860 [[nodiscard]] ConstVectorView Kjj(const Index iz = 0, const Index ia = 0) const {
861 return mdata(ia, iz, joker, 0);
862 }
863
870 [[nodiscard]] ConstVectorView K12(const Index iz = 0, const Index ia = 0) const {
871 return mdata(ia, iz, joker, 1);
872 }
873
880 [[nodiscard]] ConstVectorView K13(const Index iz = 0, const Index ia = 0) const {
881 return mdata(ia, iz, joker, 2);
882 }
883
890 [[nodiscard]] ConstVectorView K14(const Index iz = 0, const Index ia = 0) const {
891 return mdata(ia, iz, joker, 3);
892 }
893
900 [[nodiscard]] ConstVectorView K23(const Index iz = 0, const Index ia = 0) const {
901 return mdata(ia, iz, joker, mstokes_dim);
902 }
903
910 [[nodiscard]] ConstVectorView K24(const Index iz = 0, const Index ia = 0) const {
911 return mdata(ia, iz, joker, 5);
912 }
913
920 [[nodiscard]] ConstVectorView K34(const Index iz = 0, const Index ia = 0) const {
921 return mdata(ia, iz, joker, 6);
922 }
923
925 void SetZero() { mdata = 0.0; }
926
928 Tensor4& Data() { return mdata; }
929
931 [[nodiscard]] const Tensor4& Data() const { return mdata; }
932
942 const ConstMatrixView& in,
943 const Index iv = 0,
944 const Index iz = 0,
945 const Index ia = 0) const;
946
947
957 const ConstMatrixView& in,
958 const Index iv = 0,
959 const Index iz = 0,
960 const Index ia = 0) const;
961
962 friend std::ostream& operator<<(std::ostream& os, const PropagationMatrix& pm);
963
964 protected:
968 bool mvectortype{false};
969};
970
973
991 const Numeric& r,
992 const PropagationMatrix& upper_level,
993 const PropagationMatrix& lower_level,
994 const Index iz = 0,
995 const Index ia = 0);
996
997
1012 MatrixView T,
1013 const Numeric& r,
1014 const PropagationMatrix& averaged_propagation_matrix,
1015 const Index iv,
1016 const Index iz = 0,
1017 const Index ia = 0);
1018
1043 Tensor3View T,
1044 Tensor4View dT_upper_level,
1045 Tensor4View dT_lower_level,
1046 const Numeric& r,
1047 const PropagationMatrix& upper_level,
1048 const PropagationMatrix& lower_level,
1049 const Array<PropagationMatrix>& dprop_mat_upper_level,
1050 const Array<PropagationMatrix>& dprop_mat_lower_level,
1051 const Numeric& dr_dTu = 0.0,
1052 const Numeric& dr_dTl = 0.0,
1053 const Index it = -1,
1054 const Index iz = 0,
1055 const Index ia = 0);
1056
1058class StokesVector final : public PropagationMatrix {
1059 public:
1060
1072 StokesVector(const Index nr_frequencies = 0,
1073 const Index stokes_dim = 1,
1074 const Index nr_za = 1,
1075 const Index nr_aa = 1,
1076 const Numeric& v = 0.0) {
1077 mvectortype = true;
1078 mfreqs = nr_frequencies;
1079 mstokes_dim = stokes_dim;
1080 mza = nr_za;
1081 maa = nr_aa;
1082 ARTS_ASSERT(mstokes_dim < 5 and mstokes_dim > 0);
1084 };
1085
1090 explicit StokesVector(Tensor4 x) {
1091 mfreqs = x.nrows();
1092 mstokes_dim = x.ncols();
1093 mza = x.npages();
1094 maa = x.nbooks();
1095 mdata = std::move(x);
1096 mvectortype = true;
1097 ARTS_ASSERT (not (mstokes_dim > 4 or mstokes_dim < 1),
1098 "Tensor4 is bad for StokesVector");
1099 }
1100
1105 explicit StokesVector(const ConstVectorView& x) {
1106 mfreqs = 1;
1107 mstokes_dim = x.nelem();
1108 mza = 1;
1109 maa = 1;
1110 ARTS_ASSERT(mstokes_dim < 5 and mstokes_dim > 0);
1111 mvectortype = true;
1112 mdata.resize(1, 1, 1, mstokes_dim);
1113 for (Index i = 0; i < mstokes_dim; i++) mdata(0, 0, 0, i) = x[i];
1114 };
1115
1121 mfreqs = x.nrows();
1122 mstokes_dim = x.ncols();
1123 mza = 1;
1124 maa = 1;
1125 ARTS_ASSERT(mstokes_dim < 5 and mstokes_dim > 0);
1126 mvectortype = true;
1128 for (Index j = 0; j < mstokes_dim; j++) {
1129 for (Index i = 0; i < mstokes_dim; i++) {
1130 mdata(0, 0, j, i) = x(j, i);
1131 }
1132 }
1133 };
1134
1142 const StokesVector& b,
1143 const Numeric& scale = 0.5) {
1144 mfreqs = a.NumberOfFrequencies();
1145 mstokes_dim = a.StokesDimensions();
1146 mza = a.NumberOfZenithAngles();
1147 maa = a.NumberOfAzimuthAngles();
1149 mvectortype = true;
1150
1151 for (Index i = 0; i < maa; i++)
1152 for (Index j = 0; j < mza; j++)
1153 for (Index k = 0; k < mfreqs; k++) {
1154 switch (mstokes_dim) {
1155 case 4:
1156 mdata(i, j, k, 3) = (a.mdata(i, j, k, 3) + b.mdata(i, j, k, 3)) *
1157 scale; /* FALLTHROUGH */
1158 case 3:
1159 mdata(i, j, k, 2) = (a.mdata(i, j, k, 2) + b.mdata(i, j, k, 2)) *
1160 scale; /* FALLTHROUGH */
1161 case 2:
1162 mdata(i, j, k, 1) = (a.mdata(i, j, k, 1) + b.mdata(i, j, k, 1)) *
1163 scale; /* FALLTHROUGH */
1164 case 1:
1165 mdata(i, j, k, 0) = (a.mdata(i, j, k, 0) + b.mdata(i, j, k, 0)) *
1166 scale; /* FALLTHROUGH */
1167 }
1168 }
1169 };
1170
1172 [[nodiscard]] Index NumberOfNeededVectors() const { return mstokes_dim; }
1173
1180 mdata += x.Data()(joker, joker, joker, Range(0, mstokes_dim, 1));
1181 return *this;
1182 }
1183
1190 MultiplyAndAdd(lpms.scale, lpms.bas);
1191 return *this;
1192 }
1193
1206 return *this;
1207 }
1208
1215 operator=(lpms.bas);
1216 mdata *= lpms.scale;
1217 return *this;
1218 }
1219
1226 mdata = x;
1227 return *this;
1228 }
1229
1235 void MultiplyAndAdd(const Numeric x, const PropagationMatrix& y) {
1240
1241 const Tensor4& data = y.Data();
1242
1243 for (Index i = 0; i < maa; i++) {
1244 for (Index j = 0; j < mza; j++) {
1245 for (Index k = 0; k < mfreqs; k++) {
1246 switch (mstokes_dim) {
1247 case 4:
1248 mdata(i, j, k, 3) += x * data(i, j, k, 3); /* FALLTHROUGH */
1249 case 3:
1250 mdata(i, j, k, 2) += x * data(i, j, k, 2); /* FALLTHROUGH */
1251 case 2:
1252 mdata(i, j, k, 1) += x * data(i, j, k, 1); /* FALLTHROUGH */
1253 case 1:
1254 mdata(i, j, k, 0) += x * data(i, j, k, 0); /* FALLTHROUGH */
1255 }
1256 }
1257 }
1258 }
1259 }
1260
1269 const Index iz = 0,
1270 const Index ia = 0) {
1271 return mdata(ia, iz, iv, joker);
1272 }
1273
1281 [[nodiscard]] ConstVectorView VectorAtPosition(const Index iv = 0,
1282 const Index iz = 0,
1283 const Index ia = 0) const {
1284 return mdata(ia, iz, iv, joker);
1285 }
1286
1288 const Index iv = 0,
1289 const Index iz = 0,
1290 const Index ia = 0) {
1291 mdata(ia, iz, iv, joker) = x;
1292 }
1293
1303 const ConstVectorView& vec2,
1304 const Index iv = 0,
1305 const Index iz = 0,
1306 const Index ia = 0) {
1307 switch (mstokes_dim) {
1308 case 4:
1309 mdata(ia, iz, iv, 3) += (vec1[3] + vec2[3]) * 0.5; /* FALLTHROUGH */
1310 case 3:
1311 mdata(ia, iz, iv, 2) += (vec1[2] + vec2[2]) * 0.5; /* FALLTHROUGH */
1312 case 2:
1313 mdata(ia, iz, iv, 1) += (vec1[1] + vec2[1]) * 0.5; /* FALLTHROUGH */
1314 case 1:
1315 mdata(ia, iz, iv, 0) += (vec1[0] + vec2[0]) * 0.5; /* FALLTHROUGH */
1316 }
1317 }
1318
1327 [[nodiscard]] bool IsPolarized(const Index iv = 0,
1328 const Index iz = 0,
1329 const Index ia = 0) const {
1330 switch (mstokes_dim) {
1331 case 4:
1332 if (K14(iz, ia)[iv] not_eq 0.0) return true; /* FALLTHROUGH */
1333 case 3:
1334 if (K13(iz, ia)[iv] not_eq 0.0) return true; /* FALLTHROUGH */
1335 case 2:
1336 if (K12(iz, ia)[iv] not_eq 0.0) return true; /* FALLTHROUGH */
1337 }
1338 return false;
1339 }
1340
1349 [[nodiscard]] bool IsUnpolarized(const Index iv = 0,
1350 const Index iz = 0,
1351 const Index ia = 0) const {
1352 return not IsPolarized(iv, iz, ia);
1353 }
1354
1355 [[nodiscard]] bool allZeroes() const {
1356 for (Index i=0; i<maa; i++)
1357 for (Index j=0; j<mza; j++)
1358 for (Index k=0; k<mfreqs; k++)
1359 for (Index m=0; m<mstokes_dim; m++)
1360 if (mdata(i, j, k, m) not_eq 0)
1361 return false;
1362 return true;
1363 }
1364
1365 friend std::ostream& operator<<(std::ostream& os, const StokesVector& pm);
1366};
1367
1371
1379 const Numeric& x) {
1380 return LazyScale<PropagationMatrix>(pm, x);
1381}
1382
1390 const PropagationMatrix& pm) {
1391 return LazyScale<PropagationMatrix>(pm, x);
1392}
1393
1394
1396template <class T>
1398 const Vector& f_grid,
1399 const Index sd = 1) noexcept {
1400 const Index nf = f_grid.nelem();
1401 for (auto& var : main) {
1402 const bool bad_stokes = sd > var.StokesDimensions();
1403 const bool bad_freq = nf not_eq var.NumberOfFrequencies();
1404 if (bad_freq or bad_stokes) return true;
1405 }
1406
1407 return false;
1408}
1409
1410#endif //propagationmatrix_h
int main(int argc, char **argv)
This can be used to make arrays out of anything.
Definition: array.h:48
A constant view of a Matrix.
Definition: matpackI.h:1043
Index nrows() const noexcept
Definition: matpackI.h:1055
Index ncols() const noexcept
Definition: matpackI.h:1056
Index ncols() const noexcept
Definition: matpackIV.h:142
Index nrows() const noexcept
Definition: matpackIV.h:141
Index nbooks() const noexcept
Definition: matpackIV.h:139
Index npages() const noexcept
Definition: matpackIV.h:140
A constant view of a Vector.
Definition: matpackI.h:512
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
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:1164
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.
friend std::ostream & operator<<(std::ostream &os, const PropagationMatrix &pm)
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:159
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
friend std::ostream & operator<<(std::ostream &os, const StokesVector &pm)
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.
StokesVector(ConstMatrixView x)
Construct a new Stokes Vector object.
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:246
The Tensor4View class.
Definition: matpackIV.h:292
The Tensor4 class.
Definition: matpackIV.h:429
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1054
The VectorView class.
Definition: matpackI.h:663
The Vector class.
Definition: matpackI.h:899
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define max(a, b)
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
constexpr Numeric k
Boltzmann constant convenience name [J/K].
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)
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.
LazyScale< PropagationMatrix > operator*(const PropagationMatrix &pm, const Numeric &x)
Returns a lazy multiplier.
constexpr Index stokes2need(Index nstokes)
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)
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.
#define v
#define a
#define b