ARTS 2.5.11 (git: 725533f0)
zeeman.cc
Go to the documentation of this file.
1
13#include "zeeman.h"
14
15#include "arts_conversions.h"
16#include "linescaling.h"
17#include "lineshape.h"
18#include "species_info.h"
19
21 PropagationMatrix& propmat_clearsky,
22 StokesVector& nlte_source,
23 ArrayOfPropagationMatrix& dpropmat_clearsky_dx,
24 ArrayOfStokesVector& dnlte_source_dx,
25 const ArrayOfArrayOfSpeciesTag& abs_species,
26 const ArrayOfSpeciesTag& select_abs_species,
27 const ArrayOfRetrievalQuantity& jacobian_quantities,
28 const ArrayOfArrayOfAbsorptionLines& abs_lines_per_species,
29 const SpeciesIsotopologueRatios& isotopologue_ratios,
30 const Vector& f_grid,
31 const Vector& rtp_vmr,
32 const EnergyLevelMap& rtp_nlte,
33 const Vector& rtp_mag,
34 const Vector& rtp_los,
35 const Numeric& rtp_pressure,
36 const Numeric& rtp_temperature,
37 const Index& nlte_do,
38 const Index& manual_tag,
39 const Numeric& H0,
40 const Numeric& theta0,
41 const Numeric& eta0) {
42
43 // Size of problem
44 const Index nf = f_grid.nelem();
45 const Index nq = jacobian_quantities.nelem();
46 const Index ns = abs_species.nelem();
47
48 // Possible things that can go wrong in this code (excluding line parameters)
49 check_abs_species(abs_species);
50 ARTS_USER_ERROR_IF((rtp_mag.nelem() not_eq 3) and (not manual_tag),
51 "Only for 3D *rtp_mag* or a manual magnetic field")
52 ARTS_USER_ERROR_IF(rtp_vmr.nelem() not_eq abs_species.nelem(),
53 "*rtp_vmr* must match *abs_species*")
54 ARTS_USER_ERROR_IF(propmat_clearsky.NumberOfFrequencies() not_eq nf,
55 "*f_grid* must match *propmat_clearsky*")
56 ARTS_USER_ERROR_IF(propmat_clearsky.StokesDimensions() not_eq 4,
57 "*propmat_clearsky* must have *stokes_dim* 4")
58 ARTS_USER_ERROR_IF(nlte_source.NumberOfFrequencies() not_eq nf,
59 "*f_grid* must match *nlte_source*")
60 ARTS_USER_ERROR_IF(nlte_source.StokesDimensions() not_eq 4,
61 "*nlte_source* must have *stokes_dim* 4")
62 ARTS_USER_ERROR_IF(not nq and (nq not_eq dpropmat_clearsky_dx.nelem()),
63 "*dpropmat_clearsky_dx* must match derived form of *jacobian_quantities*")
64 ARTS_USER_ERROR_IF(not nq and bad_propmat(dpropmat_clearsky_dx, f_grid, 4),
65 "*dpropmat_clearsky_dx* must have Stokes dim 4 and frequency dim same as *f_grid*")
66 ARTS_USER_ERROR_IF(nlte_do and (nq not_eq dnlte_source_dx.nelem()),
67 "*dnlte_source_dx* must match derived form of *jacobian_quantities* when non-LTE is on")
68 ARTS_USER_ERROR_IF(nlte_do and bad_propmat(dnlte_source_dx, f_grid, 4),
69 "*dnlte_source_dx* must have Stokes dim 4 and frequency dim same as *f_grid* when non-LTE is on")
70 ARTS_USER_ERROR_IF(any_negative(f_grid), "Negative frequency (at least one value).")
71 ARTS_USER_ERROR_IF(any_negative(rtp_vmr), "Negative VMR (at least one value).")
72 ARTS_USER_ERROR_IF(any_negative(rtp_nlte.value), "Negative NLTE (at least one value).")
73 ARTS_USER_ERROR_IF(rtp_temperature <= 0, "Non-positive temperature")
74 ARTS_USER_ERROR_IF(rtp_pressure <= 0, "Non-positive pressure")
75 ARTS_USER_ERROR_IF(manual_tag and H0 < 0, "Negative manual magnetic field strength")
76
77 // Magnetic field internals and derivatives...
78 const auto X =
79 manual_tag
82 : Zeeman::FromGrids(rtp_mag[0],
83 rtp_mag[1],
84 rtp_mag[2],
85 Conversion::deg2rad(rtp_los[0]),
86 Conversion::deg2rad(rtp_los[1]));
87
88 // Polarization
89 const auto polarization_scale_data = Zeeman::AllPolarization(X.theta, X.eta);
90 const auto polarization_scale_dtheta_data =
91 Zeeman::AllPolarization_dtheta(X.theta, X.eta);
92 const auto polarization_scale_deta_data =
93 Zeeman::AllPolarization_deta(X.theta, X.eta);
94
95 // Deal with sparse computational grid
96 const Vector f_grid_sparse(0);
97 const Numeric sparse_limit = 0;
98
99 for (auto polar : {Zeeman::Polarization::SigmaMinus,
102
103 // Calculations data
104 LineShape::ComputeData com(f_grid, jacobian_quantities, nlte_do);
105 LineShape::ComputeData sparse_com(f_grid_sparse, jacobian_quantities, nlte_do);
106
107 auto& pol = Zeeman::SelectPolarization(polarization_scale_data, polar);
108 auto& dpol_dtheta =
109 Zeeman::SelectPolarization(polarization_scale_dtheta_data, polar);
110 auto& dpol_deta =
111 Zeeman::SelectPolarization(polarization_scale_deta_data, polar);
112
113 for (Index ispecies = 0; ispecies < ns; ispecies++) {
114 // Skip it if there are no species or there is no Zeeman
115 if (not abs_species[ispecies].nelem() or
116 not abs_species[ispecies].Zeeman() or
117 not abs_lines_per_species[ispecies].nelem())
118 continue;
119 if (select_abs_species.nelem() and
120 select_abs_species not_eq abs_species[ispecies])
121 continue;
122
123 for (auto& band : abs_lines_per_species[ispecies]) {
124 LineShape::compute(com, sparse_com,
125 band, jacobian_quantities, rtp_nlte,
126 band.BroadeningSpeciesVMR(rtp_vmr, abs_species), abs_species[ispecies], rtp_vmr[ispecies],
127 isotopologue_ratios[band.Isotopologue()], rtp_pressure, rtp_temperature, X.H, sparse_limit,
128 polar, Options::LblSpeedup::None, false);
129
130 }
131 }
132
133 // Sum up the propagation matrix
134 Zeeman::sum(propmat_clearsky, com.F, pol);
135
136 // Sum up the Jacobian
137 for (Index j=0; j<nq; j++) {
138 auto& deriv = jacobian_quantities[j];
139
140 if (not deriv.propmattype()) continue;
141
142 if (deriv == Jacobian::Atm::MagneticU) {
143 Zeeman::dsum(dpropmat_clearsky_dx[j], com.F, com.dF(joker, j),
144 pol, dpol_dtheta, dpol_deta,
145 X.dH_du, X.dtheta_du, X.deta_du);
146 } else if (deriv == Jacobian::Atm::MagneticV) {
147 Zeeman::dsum(dpropmat_clearsky_dx[j], com.F, com.dF(joker, j),
148 pol, dpol_dtheta, dpol_deta,
149 X.dH_dv, X.dtheta_dv, X.deta_dv);
150 } else if (deriv == Jacobian::Atm::MagneticW) {
151 Zeeman::dsum(dpropmat_clearsky_dx[j], com.F, com.dF(joker, j),
152 pol, dpol_dtheta, dpol_deta,
153 X.dH_dw, X.dtheta_dw, X.deta_dw);
154 } else {
155 Zeeman::sum(dpropmat_clearsky_dx[j], com.dF(joker, j), pol);
156 }
157 }
158
159 if (nlte_do) {
160 // Sum up the source vector
161 Zeeman::sum(nlte_source, com.N, pol, false);
162
163 // Sum up the Jacobian
164 for (Index j=0; j<nq; j++) {
165 auto& deriv = jacobian_quantities[j];
166
167 if (not deriv.propmattype()) continue;
168
169 if (deriv == Jacobian::Atm::MagneticU) {
170 Zeeman::dsum(dnlte_source_dx[j], com.N, com.dN(joker, j),
171 pol, dpol_dtheta, dpol_deta,
172 X.dH_du, X.dtheta_du, X.deta_du, false);
173 } else if (deriv == Jacobian::Atm::MagneticV) {
174 Zeeman::dsum(dnlte_source_dx[j], com.N, com.dN(joker, j),
175 pol, dpol_dtheta, dpol_deta,
176 X.dH_dv, X.dtheta_dv, X.deta_dv, false);
177 } else if (deriv == Jacobian::Atm::MagneticW) {
178 Zeeman::dsum(dnlte_source_dx[j], com.N, com.dN(joker, j),
179 pol, dpol_dtheta, dpol_deta,
180 X.dH_dw, X.dtheta_dw, X.deta_dw, false);
181 } else {
182 Zeeman::sum(dnlte_source_dx[j], com.dN(joker, j), pol, false);
183 }
184 }
185 }
186 }
187}
Common ARTS conversions.
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
Constains various line scaling functions.
constexpr bool any_negative(const MatpackType &var) noexcept
Checks for negative values.
Definition math_funcs.h:132
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
void compute(ComputeData &com, ComputeData &sparse_com, const AbsorptionLines &band, const ArrayOfRetrievalQuantity &jacobian_quantities, const EnergyLevelMap &nlte, const Vector &vmrs, const ArrayOfSpeciesTag &self_tag, const Numeric &self_vmr, const Numeric &isot_ratio, const Numeric &P, const Numeric &T, const Numeric &H, const Numeric &sparse_lim, const Zeeman::Polarization zeeman_polarization, const Options::LblSpeedup speedup_type, const bool robust) ARTS_NOEXCEPT
Compute the absorption of an absorption band.
Implements Zeeman modeling.
const PolarizationVector & SelectPolarization(const AllPolarizationVectors &data, Polarization type) noexcept
Selects the polarization vector depending on polarization type.
constexpr Derived FromPreDerived(Numeric H, Numeric theta, Numeric eta) noexcept
Sets Derived from predefined Derived parameters.
Definition zeemandata.h:551
Derived FromGrids(Numeric u, Numeric v, Numeric w, Numeric z, Numeric a) noexcept
Computes the derived plane from ARTS grids.
AllPolarizationVectors AllPolarization_deta(Numeric theta, Numeric eta) noexcept
The derivative of AllPolarization wrt eta.
AllPolarizationVectors AllPolarization_dtheta(Numeric theta, const Numeric eta) noexcept
The derivative of AllPolarization wrt theta.
void dsum(PropagationMatrix &pm, const ComplexVectorView &abs, const ComplexVectorView &dabs, const PolarizationVector &polvec, const PolarizationVector &dpolvec_dtheta, const PolarizationVector &dpolvec_deta, const Numeric dH, const Numeric dt, const Numeric de, const bool do_phase) ARTS_NOEXCEPT
Sums the Zeeman components derivatives into a propagation matrix.
void sum(PropagationMatrix &pm, const ComplexVectorView &abs, const PolarizationVector &polvec, const bool do_phase) ARTS_NOEXCEPT
Sums the Zeeman components into a propagation matrix.
AllPolarizationVectors AllPolarization(Numeric theta, Numeric eta) noexcept
Computes the polarization of each polarization type.
Some molecular constants.
void check_abs_species(const ArrayOfArrayOfSpeciesTag &abs_species)
Main computational data for the line shape and strength calculations.
Definition lineshape.h:840
void zeeman_on_the_fly(PropagationMatrix &propmat_clearsky, StokesVector &nlte_source, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, ArrayOfStokesVector &dnlte_source_dx, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfSpeciesTag &select_abs_species, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfAbsorptionLines &abs_lines_per_species, const SpeciesIsotopologueRatios &isotopologue_ratios, const Vector &f_grid, const Vector &rtp_vmr, const EnergyLevelMap &rtp_nlte, const Vector &rtp_mag, const Vector &rtp_los, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const Index &nlte_do, const Index &manual_tag, const Numeric &H0, const Numeric &theta0, const Numeric &eta0)
Main and only way to compute Zeeman effect.
Definition zeeman.cc:20