34#ifndef propagationmatrix_h
35#define propagationmatrix_h
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;
72 if (nstokes_needed < 5 and nstokes_needed > 0)
return nstokes_needed;
74 ARTS_ASSERT (
false,
"Cannot understand the input Stokes dimensions");
75 return std::numeric_limits<Index>::max();
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;
87 if (nstokes < 5 and nstokes > 0)
return nstokes;
89 ARTS_ASSERT (
false,
"Cannot understand the input Stokes dimensions");
90 return std::numeric_limits<Index>::max();
131 const Index stokes_dim = 1,
132 const Index nr_za = 1,
133 const Index nr_aa = 1,
158 mdata(std::move(pm.mdata)),
169 mdata(std::move(x)) {
184 ARTS_ASSERT(
false,
"Tensor4 not representative of PropagationMatrix");
198 if (not assume_fit) {
200 "Matrix not fit as propagation matrix");
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);
212 mdata(0, 0, 0, 2) = x(0, 2);
214 mdata(0, 0, 0, 1) = x(0, 1);
216 mdata(0, 0, 0, 0) = x(0, 0);
246 const Index ia = 0)
const {
250 if (n not_eq 0)
return false;
263 const Index ia = 0)
const {
264 if (
mdata(ia, iz, iv, 0) == 0.0)
270 #pragma GCC diagnostic push
271 #pragma GCC diagnostic ignored "-Wreturn-type"
278 #pragma GCC diagnostic pop
285 const Index ia = 0)
const;
299 const Index ia = 0) {
315 const Index ia = 0) {
329 const Index ia = 0)
const;
342 mdata = std::move(pm.mdata);
386 const Index ia = 0) {
413 const Index ia = 0) {
463 const Index ia = 0) {
489 const Index ia = 0) {
535 const Index ia = 0) {
561 const Index ia = 0) {
617 const Index ia = 0) {
643 const Index ia = 0) {
693 const Index ia = 0) {
719 const Index ia = 0) {
733 const Index ia = 0) {
768 const Index ia = 0)
const;
945 const Index ia = 0)
const;
960 const Index ia = 0)
const;
1020 const Index ia = 0);
1056 const Index it = -1,
1058 const Index ia = 0);
1076 const Index stokes_dim = 1,
1077 const Index nr_za = 1,
1078 const Index nr_aa = 1,
1098 mdata = std::move(x);
1101 "Tensor4 is bad for StokesVector");
1133 mdata(0, 0, j, i) = x(j, i);
1147 mfreqs =
a.NumberOfFrequencies();
1149 mza =
a.NumberOfZenithAngles();
1150 maa =
a.NumberOfAzimuthAngles();
1159 mdata(i, j, k, 3) = (
a.mdata(i, j, k, 3) +
b.mdata(i, j, k, 3)) *
1162 mdata(i, j, k, 2) = (
a.mdata(i, j, k, 2) +
b.mdata(i, j, k, 2)) *
1165 mdata(i, j, k, 1) = (
a.mdata(i, j, k, 1) +
b.mdata(i, j, k, 1)) *
1168 mdata(i, j, k, 0) = (
a.mdata(i, j, k, 0) +
b.mdata(i, j, k, 0)) *
1208 mdata(
joker,
joker,
joker,
Range(0,
mstokes_dim, 1)) = x.
Data()(
joker,
joker,
joker,
Range(0,
mstokes_dim, 1));
1251 mdata(i, j, k, 3) += x * data(i, j, k, 3);
1253 mdata(i, j, k, 2) += x * data(i, j, k, 2);
1255 mdata(i, j, k, 1) += x * data(i, j, k, 1);
1257 mdata(i, j, k, 0) += x * data(i, j, k, 0);
1273 const Index ia = 0) {
1286 const Index ia = 0)
const {
1293 const Index ia = 0) {
1309 const Index ia = 0) {
1312 mdata(ia, iz, iv, 3) += (vec1[3] + vec2[3]) * 0.5;
1314 mdata(ia, iz, iv, 2) += (vec1[2] + vec2[2]) * 0.5;
1316 mdata(ia, iz, iv, 1) += (vec1[1] + vec2[1]) * 0.5;
1318 mdata(ia, iz, iv, 0) += (vec1[0] + vec2[0]) * 0.5;
1332 const Index ia = 0)
const {
1335 if (
K14(iz, ia)[iv] not_eq 0.0)
return true;
1337 if (
K13(iz, ia)[iv] not_eq 0.0)
return true;
1339 if (
K12(iz, ia)[iv] not_eq 0.0)
return true;
1354 const Index ia = 0)
const {
1363 if (
mdata(i, j, k, m) not_eq 0)
1397 const Index sd = 1) noexcept {
1399 for (
auto& var :
main) {
1400 const bool bad_stokes = sd > var.StokesDimensions();
1401 const bool bad_freq = nf not_eq var.NumberOfFrequencies();
1402 if (bad_freq or bad_stokes)
return true;
This can be used to make arrays out of anything.
A constant view of a Matrix.
Index nrows() const noexcept
Index ncols() const noexcept
Index ncols() const noexcept
Index nrows() const noexcept
Index nbooks() const noexcept
Index npages() const noexcept
A constant view of a Vector.
Index nelem() const noexcept
Returns the number of elements.
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.
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.
friend void swap(PropagationMatrix &, PropagationMatrix &) noexcept
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.
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.
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.
void resize(Index b, Index p, Index r, Index c)
Resize function.
#define ARTS_ASSERT(condition,...)
NUMERIC Numeric
The type to use for all floating point numbers.
INDEX Index
The type to use for all integer numbers and indices.
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.