ARTS 2.5.9 (git: 825fa5f2)
Go to the documentation of this file.
1/* Copyright (C) 2014
2 Richard Larsson <>
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.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 GNU General Public License for more details.
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. */
30#include "zeeman.h"
32#include "arts_conversions.h"
33#include "linescaling.h"
34#include "lineshape.h"
35#include "species_info.h"
38 PropagationMatrix& propmat_clearsky,
39 StokesVector& nlte_source,
40 ArrayOfPropagationMatrix& dpropmat_clearsky_dx,
41 ArrayOfStokesVector& dnlte_source_dx,
42 const ArrayOfArrayOfSpeciesTag& abs_species,
43 const ArrayOfSpeciesTag& select_abs_species,
44 const ArrayOfRetrievalQuantity& jacobian_quantities,
45 const ArrayOfArrayOfAbsorptionLines& abs_lines_per_species,
46 const SpeciesIsotopologueRatios& isotopologue_ratios,
47 const Vector& f_grid,
48 const Vector& rtp_vmr,
49 const EnergyLevelMap& rtp_nlte,
50 const Vector& rtp_mag,
51 const Vector& rtp_los,
52 const Numeric& rtp_pressure,
53 const Numeric& rtp_temperature,
54 const Index& nlte_do,
55 const Index& manual_tag,
56 const Numeric& H0,
57 const Numeric& theta0,
58 const Numeric& eta0) {
60 // Size of problem
61 const Index nf = f_grid.nelem();
62 const Index nq = jacobian_quantities.nelem();
63 const Index ns = abs_species.nelem();
65 // Possible things that can go wrong in this code (excluding line parameters)
66 check_abs_species(abs_species);
67 ARTS_USER_ERROR_IF((rtp_mag.nelem() not_eq 3) and (not manual_tag),
68 "Only for 3D *rtp_mag* or a manual magnetic field")
69 ARTS_USER_ERROR_IF(rtp_vmr.nelem() not_eq abs_species.nelem(),
70 "*rtp_vmr* must match *abs_species*")
71 ARTS_USER_ERROR_IF(propmat_clearsky.NumberOfFrequencies() not_eq nf,
72 "*f_grid* must match *propmat_clearsky*")
73 ARTS_USER_ERROR_IF(propmat_clearsky.StokesDimensions() not_eq 4,
74 "*propmat_clearsky* must have *stokes_dim* 4")
75 ARTS_USER_ERROR_IF(nlte_source.NumberOfFrequencies() not_eq nf,
76 "*f_grid* must match *nlte_source*")
77 ARTS_USER_ERROR_IF(nlte_source.StokesDimensions() not_eq 4,
78 "*nlte_source* must have *stokes_dim* 4")
79 ARTS_USER_ERROR_IF(not nq and (nq not_eq dpropmat_clearsky_dx.nelem()),
80 "*dpropmat_clearsky_dx* must match derived form of *jacobian_quantities*")
81 ARTS_USER_ERROR_IF(not nq and bad_propmat(dpropmat_clearsky_dx, f_grid, 4),
82 "*dpropmat_clearsky_dx* must have Stokes dim 4 and frequency dim same as *f_grid*")
83 ARTS_USER_ERROR_IF(nlte_do and (nq not_eq dnlte_source_dx.nelem()),
84 "*dnlte_source_dx* must match derived form of *jacobian_quantities* when non-LTE is on")
85 ARTS_USER_ERROR_IF(nlte_do and bad_propmat(dnlte_source_dx, f_grid, 4),
86 "*dnlte_source_dx* must have Stokes dim 4 and frequency dim same as *f_grid* when non-LTE is on")
87 ARTS_USER_ERROR_IF(any_negative(f_grid), "Negative frequency (at least one value).")
88 ARTS_USER_ERROR_IF(any_negative(rtp_vmr), "Negative VMR (at least one value).")
89 ARTS_USER_ERROR_IF(any_negative(rtp_nlte.value), "Negative NLTE (at least one value).")
90 ARTS_USER_ERROR_IF(rtp_temperature <= 0, "Non-positive temperature")
91 ARTS_USER_ERROR_IF(rtp_pressure <= 0, "Non-positive pressure")
92 ARTS_USER_ERROR_IF(manual_tag and H0 < 0, "Negative manual magnetic field strength")
94 // Magnetic field internals and derivatives...
95 const auto X =
96 manual_tag
99 : Zeeman::FromGrids(rtp_mag[0],
100 rtp_mag[1],
101 rtp_mag[2],
102 Conversion::deg2rad(rtp_los[0]),
103 Conversion::deg2rad(rtp_los[1]));
105 // Polarization
106 const auto polarization_scale_data = Zeeman::AllPolarization(X.theta, X.eta);
107 const auto polarization_scale_dtheta_data =
108 Zeeman::AllPolarization_dtheta(X.theta, X.eta);
109 const auto polarization_scale_deta_data =
110 Zeeman::AllPolarization_deta(X.theta, X.eta);
112 // Deal with sparse computational grid
113 const Vector f_grid_sparse(0);
114 const Numeric sparse_limit = 0;
116 for (auto polar : {Zeeman::Polarization::SigmaMinus,
120 // Calculations data
121 LineShape::ComputeData com(f_grid, jacobian_quantities, nlte_do);
122 LineShape::ComputeData sparse_com(f_grid_sparse, jacobian_quantities, nlte_do);
124 auto& pol = Zeeman::SelectPolarization(polarization_scale_data, polar);
125 auto& dpol_dtheta =
126 Zeeman::SelectPolarization(polarization_scale_dtheta_data, polar);
127 auto& dpol_deta =
128 Zeeman::SelectPolarization(polarization_scale_deta_data, polar);
130 for (Index ispecies = 0; ispecies < ns; ispecies++) {
131 // Skip it if there are no species or there is no Zeeman
132 if (not abs_species[ispecies].nelem() or
133 not abs_species[ispecies].Zeeman() or
134 not abs_lines_per_species[ispecies].nelem())
135 continue;
136 if (select_abs_species.nelem() and
137 select_abs_species not_eq abs_species[ispecies])
138 continue;
140 for (auto& band : abs_lines_per_species[ispecies]) {
141 LineShape::compute(com, sparse_com,
142 band, jacobian_quantities, rtp_nlte,
143 band.BroadeningSpeciesVMR(rtp_vmr, abs_species), abs_species[ispecies], rtp_vmr[ispecies],
144 isotopologue_ratios[band.Isotopologue()], rtp_pressure, rtp_temperature, X.H, sparse_limit,
145 polar, Options::LblSpeedup::None, false);
147 }
148 }
150 // Sum up the propagation matrix
151 Zeeman::sum(propmat_clearsky, com.F, pol);
153 // Sum up the Jacobian
154 for (Index j=0; j<nq; j++) {
155 auto& deriv = jacobian_quantities[j];
157 if (not deriv.propmattype()) continue;
159 if (deriv == Jacobian::Atm::MagneticU) {
160 Zeeman::dsum(dpropmat_clearsky_dx[j], com.F, com.dF(joker, j),
161 pol, dpol_dtheta, dpol_deta,
162 X.dH_du, X.dtheta_du, X.deta_du);
163 } else if (deriv == Jacobian::Atm::MagneticV) {
164 Zeeman::dsum(dpropmat_clearsky_dx[j], com.F, com.dF(joker, j),
165 pol, dpol_dtheta, dpol_deta,
166 X.dH_dv, X.dtheta_dv, X.deta_dv);
167 } else if (deriv == Jacobian::Atm::MagneticW) {
168 Zeeman::dsum(dpropmat_clearsky_dx[j], com.F, com.dF(joker, j),
169 pol, dpol_dtheta, dpol_deta,
170 X.dH_dw, X.dtheta_dw, X.deta_dw);
171 } else {
172 Zeeman::sum(dpropmat_clearsky_dx[j], com.dF(joker, j), pol);
173 }
174 }
176 if (nlte_do) {
177 // Sum up the source vector
178 Zeeman::sum(nlte_source, com.N, pol, false);
180 // Sum up the Jacobian
181 for (Index j=0; j<nq; j++) {
182 auto& deriv = jacobian_quantities[j];
184 if (not deriv.propmattype()) continue;
186 if (deriv == Jacobian::Atm::MagneticU) {
187 Zeeman::dsum(dnlte_source_dx[j], com.N, com.dN(joker, j),
188 pol, dpol_dtheta, dpol_deta,
189 X.dH_du, X.dtheta_du, X.deta_du, false);
190 } else if (deriv == Jacobian::Atm::MagneticV) {
191 Zeeman::dsum(dnlte_source_dx[j], com.N, com.dN(joker, j),
192 pol, dpol_dtheta, dpol_deta,
193 X.dH_dv, X.dtheta_dv, X.deta_dv, false);
194 } else if (deriv == Jacobian::Atm::MagneticW) {
195 Zeeman::dsum(dnlte_source_dx[j], com.N, com.dN(joker, j),
196 pol, dpol_dtheta, dpol_deta,
197 X.dH_dw, X.dtheta_dw, X.deta_dw, false);
198 } else {
199 Zeeman::sum(dnlte_source_dx[j], com.dN(joker, j), pol, false);
200 }
201 }
202 }
203 }
Common ARTS conversions.
This can be used to make arrays out of anything.
Definition: array.h:48
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
Index NumberOfFrequencies() const
The number of frequencies of the propagation matrix.
Index StokesDimensions() const
The stokes dimension of the propagation matrix.
Stokes vector is as Propagation matrix but only has 4 possible values.
The Vector class.
Definition: matpackI.h:910
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
Constains various line scaling functions.
constexpr bool any_negative(const MatpackType &var) noexcept
Checks for negative values.
Definition: math_funcs.h:147
The type to use for all floating point numbers.
Definition: matpack.h:33
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
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:568
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.
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.
Some molecular constants.
void check_abs_species(const ArrayOfArrayOfSpeciesTag &abs_species)
Main computational data for the line shape and strength calculations.
Definition: lineshape.h:793
ComplexVector N
Definition: lineshape.h:794
ComplexVector F
Definition: lineshape.h:794
ComplexMatrix dF
Definition: lineshape.h:795
ComplexMatrix dN
Definition: lineshape.h:795
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.