ARTS 2.5.0 (git: 9ee3ac6c)
zeeman.cc
Go to the documentation of this file.
1/* Copyright (C) 2014
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. */
30#include "zeeman.h"
31#include "constants.h"
32#include "linescaling.h"
33#include "species_info.h"
34
35#include "lineshape.h"
36
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) {
58
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();
63
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.Data()), "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")
92
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]));
103
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);
110
111 // Deal with sparse computational grid
112 const Vector f_grid_sparse(0);
113 const Numeric sparse_limit = 0;
114
115 for (auto polar : {Zeeman::Polarization::SigmaMinus,
118
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);
122
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);
128
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;
133
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 true, polar, Options::LblSpeedup::None);
140
141 }
142 }
143
144 // Sum up the propagation matrix
145 Zeeman::sum(propmat_clearsky, com.F, pol);
146
147 // Sum up the Jacobian
148 for (Index j=0; j<nq; j++) {
149 auto& deriv = jacobian_quantities[j];
150
151 if (not deriv.propmattype()) continue;
152
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 }
169
170 if (nlte_do) {
171 // Sum up the source vector
172 Zeeman::sum(nlte_source, com.N, pol, false);
173
174 // Sum up the Jacobian
175 for (Index j=0; j<nq; j++) {
176 auto& deriv = jacobian_quantities[j];
177
178 if (not deriv.propmattype()) continue;
179
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 }
198}
This can be used to make arrays out of anything.
Definition: array.h:107
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
const Tensor4 & Data() const noexcept
Energy level type.
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:876
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:135
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
const Joker joker
Index nelem(const Lines &l)
Number of lines.
constexpr auto deg2rad(T x) -> decltype(x *one_degree_in_radians)
Converts degrees to radians.
Definition: constants.h:328
SpeciesIsotopologueRatios isotopologue_ratios(Type type)
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 bool do_zeeman, const Zeeman::Polarization zeeman_polarization, const Options::LblSpeedup speedup_type) ARTS_NOEXCEPT
Compute the line shape in its entirety.
Definition: lineshape.cc:3116
X3 LineCenter None
Definition: constants.h:576
Implements Zeeman modeling.
Definition: zeemandata.cc:281
const PolarizationVector & SelectPolarization(const AllPolarizationVectors &data, Polarization type) noexcept
Selects the polarization vector depending on polarization type.
Definition: zeemandata.cc:362
constexpr Derived FromPreDerived(Numeric H, Numeric theta, Numeric eta) noexcept
Sets Derived from predefined Derived parameters.
Definition: zeemandata.h:676
Derived FromGrids(Numeric u, Numeric v, Numeric w, Numeric z, Numeric a) noexcept
Computes the derived plane from ARTS grids.
Definition: zeemandata.cc:227
AllPolarizationVectors AllPolarization_deta(Numeric theta, Numeric eta) noexcept
The derivative of AllPolarization wrt eta.
Definition: zeemandata.cc:344
AllPolarizationVectors AllPolarization_dtheta(Numeric theta, const Numeric eta) noexcept
The derivative of AllPolarization wrt theta.
Definition: zeemandata.cc:327
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.
Definition: zeemandata.cc:389
void sum(PropagationMatrix &pm, const ComplexVectorView &abs, const PolarizationVector &polvec, const bool do_phase) ARTS_NOEXCEPT
Sums the Zeeman components into a propagation matrix.
Definition: zeemandata.cc:375
AllPolarizationVectors AllPolarization(Numeric theta, Numeric eta) noexcept
Computes the polarization of each polarization type.
Definition: zeemandata.cc:311
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)
ComplexVector N
Definition: lineshape.h:687
ComplexVector F
Definition: lineshape.h:687
ComplexMatrix dF
Definition: lineshape.h:688
ComplexMatrix dN
Definition: lineshape.h:688
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.
Definition: zeeman.cc:37