ARTS 2.5.4 (git: 4c0d3b4d)
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"
31#include "constants.h"
32#include "linescaling.h"
33#include "species_info.h"
35#include "lineshape.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 ArrayOfRetrievalQuantity& jacobian_quantities,
44 const ArrayOfArrayOfAbsorptionLines& abs_lines_per_species,
46 const Vector& f_grid,
47 const Vector& rtp_vmr,
48 const EnergyLevelMap& rtp_nlte,
49 const Vector& rtp_mag,
50 const Vector& rtp_los,
51 const Numeric& rtp_pressure,
52 const Numeric& rtp_temperature,
53 const Index& nlte_do,
54 const Index& manual_tag,
55 const Numeric& H0,
56 const Numeric& theta0,
57 const Numeric& eta0) {
59 // Size of problem
60 const Index nf = f_grid.nelem();
61 const Index nq = jacobian_quantities.nelem();
62 const Index ns = abs_species.nelem();
64 // Possible things that can go wrong in this code (excluding line parameters)
65 check_abs_species(abs_species);
66 ARTS_USER_ERROR_IF((rtp_mag.nelem() not_eq 3) and (not manual_tag),
67 "Only for 3D *rtp_mag* or a manual magnetic field")
68 ARTS_USER_ERROR_IF(rtp_vmr.nelem() not_eq abs_species.nelem(),
69 "*rtp_vmr* must match *abs_species*")
70 ARTS_USER_ERROR_IF(propmat_clearsky.NumberOfFrequencies() not_eq nf,
71 "*f_grid* must match *propmat_clearsky*")
72 ARTS_USER_ERROR_IF(propmat_clearsky.StokesDimensions() not_eq 4,
73 "*propmat_clearsky* must have *stokes_dim* 4")
74 ARTS_USER_ERROR_IF(nlte_source.NumberOfFrequencies() not_eq nf,
75 "*f_grid* must match *nlte_source*")
76 ARTS_USER_ERROR_IF(nlte_source.StokesDimensions() not_eq 4,
77 "*nlte_source* must have *stokes_dim* 4")
78 ARTS_USER_ERROR_IF(not nq and (nq not_eq dpropmat_clearsky_dx.nelem()),
79 "*dpropmat_clearsky_dx* must match derived form of *jacobian_quantities*")
80 ARTS_USER_ERROR_IF(not nq and bad_propmat(dpropmat_clearsky_dx, f_grid, 4),
81 "*dpropmat_clearsky_dx* must have Stokes dim 4 and frequency dim same as *f_grid*")
82 ARTS_USER_ERROR_IF(nlte_do and (nq not_eq dnlte_source_dx.nelem()),
83 "*dnlte_source_dx* must match derived form of *jacobian_quantities* when non-LTE is on")
84 ARTS_USER_ERROR_IF(nlte_do and bad_propmat(dnlte_source_dx, f_grid, 4),
85 "*dnlte_source_dx* must have Stokes dim 4 and frequency dim same as *f_grid* when non-LTE is on")
86 ARTS_USER_ERROR_IF(any_negative(f_grid), "Negative frequency (at least one value).")
87 ARTS_USER_ERROR_IF(any_negative(rtp_vmr), "Negative VMR (at least one value).")
88 ARTS_USER_ERROR_IF(any_negative(rtp_nlte.value), "Negative NLTE (at least one value).")
89 ARTS_USER_ERROR_IF(rtp_temperature <= 0, "Non-positive temperature")
90 ARTS_USER_ERROR_IF(rtp_pressure <= 0, "Non-positive pressure")
91 ARTS_USER_ERROR_IF(manual_tag and H0 < 0, "Negative manual magnetic field strength")
93 // Magnetic field internals and derivatives...
94 const auto X =
95 manual_tag
98 : Zeeman::FromGrids(rtp_mag[0],
99 rtp_mag[1],
100 rtp_mag[2],
101 Conversion::deg2rad(rtp_los[0]),
102 Conversion::deg2rad(rtp_los[1]));
104 // Polarization
105 const auto polarization_scale_data = Zeeman::AllPolarization(X.theta, X.eta);
106 const auto polarization_scale_dtheta_data =
107 Zeeman::AllPolarization_dtheta(X.theta, X.eta);
108 const auto polarization_scale_deta_data =
109 Zeeman::AllPolarization_deta(X.theta, X.eta);
111 // Deal with sparse computational grid
112 const Vector f_grid_sparse(0);
113 const Numeric sparse_limit = 0;
115 for (auto polar : {Zeeman::Polarization::SigmaMinus,
119 // Calculations data
120 LineShape::ComputeData com(f_grid, jacobian_quantities, nlte_do);
121 LineShape::ComputeData sparse_com(f_grid_sparse, jacobian_quantities, nlte_do);
123 auto& pol = Zeeman::SelectPolarization(polarization_scale_data, polar);
124 auto& dpol_dtheta =
125 Zeeman::SelectPolarization(polarization_scale_dtheta_data, polar);
126 auto& dpol_deta =
127 Zeeman::SelectPolarization(polarization_scale_deta_data, polar);
129 for (Index ispecies = 0; ispecies < ns; ispecies++) {
130 // Skip it if there are no species or there is no Zeeman
131 if (not abs_species[ispecies].nelem() or not abs_species[ispecies].Zeeman() or not abs_lines_per_species[ispecies].nelem())
132 continue;
134 for (auto& band : abs_lines_per_species[ispecies]) {
135 LineShape::compute(com, sparse_com,
136 band, jacobian_quantities, rtp_nlte,
137 band.BroadeningSpeciesVMR(rtp_vmr, abs_species), abs_species[ispecies], rtp_vmr[ispecies],
138 isotopologue_ratios[band.Isotopologue()], rtp_pressure, rtp_temperature, X.H, sparse_limit,
139 polar, Options::LblSpeedup::None, false);
141 }
142 }
144 // Sum up the propagation matrix
145 Zeeman::sum(propmat_clearsky, com.F, pol);
147 // Sum up the Jacobian
148 for (Index j=0; j<nq; j++) {
149 auto& deriv = jacobian_quantities[j];
151 if (not deriv.propmattype()) continue;
153 if (deriv == Jacobian::Atm::MagneticU) {
154 Zeeman::dsum(dpropmat_clearsky_dx[j], com.F, com.dF(joker, j),
155 pol, dpol_dtheta, dpol_deta,
156 X.dH_du, X.dtheta_du, X.deta_du);
157 } else if (deriv == Jacobian::Atm::MagneticV) {
158 Zeeman::dsum(dpropmat_clearsky_dx[j], com.F, com.dF(joker, j),
159 pol, dpol_dtheta, dpol_deta,
160 X.dH_dv, X.dtheta_dv, X.deta_dv);
161 } else if (deriv == Jacobian::Atm::MagneticW) {
162 Zeeman::dsum(dpropmat_clearsky_dx[j], com.F, com.dF(joker, j),
163 pol, dpol_dtheta, dpol_deta,
164 X.dH_dw, X.dtheta_dw, X.deta_dw);
165 } else {
166 Zeeman::sum(dpropmat_clearsky_dx[j], com.dF(joker, j), pol);
167 }
168 }
170 if (nlte_do) {
171 // Sum up the source vector
172 Zeeman::sum(nlte_source, com.N, pol, false);
174 // Sum up the Jacobian
175 for (Index j=0; j<nq; j++) {
176 auto& deriv = jacobian_quantities[j];
178 if (not deriv.propmattype()) continue;
180 if (deriv == Jacobian::Atm::MagneticU) {
181 Zeeman::dsum(dnlte_source_dx[j], com.N, com.dN(joker, j),
182 pol, dpol_dtheta, dpol_deta,
183 X.dH_du, X.dtheta_du, X.deta_du, false);
184 } else if (deriv == Jacobian::Atm::MagneticV) {
185 Zeeman::dsum(dnlte_source_dx[j], com.N, com.dN(joker, j),
186 pol, dpol_dtheta, dpol_deta,
187 X.dH_dv, X.dtheta_dv, X.deta_dv, false);
188 } else if (deriv == Jacobian::Atm::MagneticW) {
189 Zeeman::dsum(dnlte_source_dx[j], com.N, com.dN(joker, j),
190 pol, dpol_dtheta, dpol_deta,
191 X.dH_dw, X.dtheta_dw, X.deta_dw, false);
192 } else {
193 Zeeman::sum(dnlte_source_dx[j], com.dN(joker, j), pol, false);
194 }
195 }
196 }
197 }
This can be used to make arrays out of anything.
Definition: array.h:108
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:197
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:541
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:908
Constants of physical expressions as constexpr.
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
#define ns
Constains various line scaling functions.
constexpr bool any_negative(const MatpackType &var) noexcept
Checks for negative values.
Definition: math_funcs.h:138
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
Index nelem(const Lines &l)
Number of lines.
constexpr auto deg2rad(T x) noexcept -> decltype(x *one_degree_in_radians)
Converts degrees to radians.
Definition: constants.h:331
SpeciesIsotopologueRatios isotopologue_ratios()
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:614
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:781
ComplexVector N
Definition: lineshape.h:782
ComplexVector F
Definition: lineshape.h:782
ComplexMatrix dF
Definition: lineshape.h:783
ComplexMatrix dN
Definition: lineshape.h:783
void zeeman_on_the_fly(PropagationMatrix &propmat_clearsky, StokesVector &nlte_source, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, ArrayOfStokesVector &dnlte_source_dx, const ArrayOfArrayOfSpeciesTag &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.