15#include "matpack_data.h"
16#include "matpack_eigen.h"
22 const Numeric& GL)
noexcept {
24 qns.has(QuantumNumberType::Omega,
26 QuantumNumberType::Lambda,
27 QuantumNumberType::S)) {
28 auto& Omega = qns[QuantumNumberType::Omega];
29 auto& J = qns[QuantumNumberType::J];
30 auto& Lambda = qns[QuantumNumberType::Lambda];
31 auto& S = qns[QuantumNumberType::S];
33 Omega.upp(), J.upp(), Lambda.upp(), S.upp(), GS, GL),
35 Omega.low(), J.low(), Lambda.low(), S.low(), GS, GL)};
39 qns.has(QuantumNumberType::N,
41 QuantumNumberType::Lambda,
42 QuantumNumberType::S)) {
43 auto& N = qns[QuantumNumberType::N];
44 auto& J = qns[QuantumNumberType::J];
45 auto& Lambda = qns[QuantumNumberType::Lambda];
46 auto& S = qns[QuantumNumberType::S];
59 return SimpleG(qid.val, GS, GL);
77 using std::atan2, std::cos, std::sin;
79 if (J.isUndefined() or N.isUndefined())
return NAN;
82 auto nom = (lB + lD * (J * J + J + 1) + lH * pow2(J * J + J + 1)) *
83 (2 * sqrt(J * J + J) / (2 * J + 1));
86 B * J * (J - 1) - D * pow2(J * (J - 1)) + H * pow3(J * (J - 1)) +
87 (gB + gD * J * (J - 1) + gH * pow2(J * (J - 1))) * (J - 1) +
88 (lB + lD * J * (J - 1) + lH * pow2(J * (J - 1))) *
89 (2. / 3. - 2 * J / (2 * J + 1)) -
90 (B * (J + 2) * (J + 1) - D * pow2((J + 2) * (J + 1)) +
91 H * pow3((J + 2) * (J + 1)) -
92 (gB + gD * (J + 2) * (J + 1) + gH * pow2((J + 2) * (J + 1))) * (J + 2) +
93 (lB + lD * (J + 2) * (J + 1) + lH * pow2((J + 2) * (J + 1))) *
94 (2. / 3. - 2 * (J + 1) / (2 * J + 1)));
96 auto phi = atan2(2 * nom, denom) / 2;
98 if (J == N)
return (GS + GR) / (J * (J + 1)) - GR;
100 return (GS + GR) * (pow2(cos(phi)) / J - pow2(sin(phi)) / (J + 1)) +
101 2 * GLE * cos(2 * phi) / (2 * J + 1) - GR;
102 return (GS + GR) * (pow2(sin(phi)) / J - pow2(cos(phi)) / (J + 1)) -
103 2 * GLE * cos(2 * phi) / (2 * J + 1) - GR;
111 if (k.isUndefined() or j.isUndefined() or j == 0)
return 0;
112 return gperp + (gperp + gpara) * (pow2(k) / (j * (j + 1)));
117 if (qid.Isotopologue() ==
"O2-66") {
118 if (qid.val.has(QuantumNumberType::J,
119 QuantumNumberType::N,
120 QuantumNumberType::v1)) {
121 if (qid.val[QuantumNumberType::v1].low() == 0 and
122 qid.val[QuantumNumberType::v1].upp() == 0) {
123 constexpr Numeric GS = 2.002084;
124 constexpr Numeric GLE = 2.77e-3;
125 constexpr Numeric GR = -1.16e-4;
126 constexpr Numeric B = 43100.44276e6;
127 constexpr Numeric D = 145.1271e3;
128 constexpr Numeric H = 49e-3;
129 constexpr Numeric lB = 59501.3438e6;
130 constexpr Numeric lD = 58.3680e3;
131 constexpr Numeric lH = 290.8e-3;
132 constexpr Numeric gB = -252.58634e6;
133 constexpr Numeric gD = -243.42;
134 constexpr Numeric gH = -1.46e-3;
136 auto JU = qid.val[QuantumNumberType::J].upp();
137 auto NU = qid.val[QuantumNumberType::N].upp();
139 JU, NU, GS, GR, GLE, B, D, H, gB, gD, gH, lB, lD, lH);
140 auto JL = qid.val[QuantumNumberType::J].low();
141 auto NL = qid.val[QuantumNumberType::N].low();
143 JL, NL, GS, GR, GLE, B, D, H, gB, gD, gH, lB, lD, lH);
147 }
else if (qid.Isotopologue() ==
"O2-68") {
148 if (qid.val.has(QuantumNumberType::J,
149 QuantumNumberType::N,
150 QuantumNumberType::v1)) {
151 if (qid.val[QuantumNumberType::v1].low() == 0 and
152 qid.val[QuantumNumberType::v1].upp() == 0) {
153 constexpr Numeric GS = 2.002025;
154 constexpr Numeric GLE = 2.813e-3;
155 constexpr Numeric GR = -1.26e-4;
156 constexpr Numeric B = 40707.38657e6;
157 constexpr Numeric D = 129.4142e3;
158 constexpr Numeric H = 0;
159 constexpr Numeric lB = 59499.0375e6;
160 constexpr Numeric lD = 54.9777e3;
161 constexpr Numeric lH = 272.1e-3;
162 constexpr Numeric gB = -238.51530e6;
163 constexpr Numeric gD = -217.77;
164 constexpr Numeric gH = -1.305e-3;
166 auto JU = qid.val[QuantumNumberType::J].upp();
167 auto NU = qid.val[QuantumNumberType::N].upp();
169 JU, NU, GS, GR, GLE, B, D, H, gB, gD, gH, lB, lD, lH);
170 auto JL = qid.val[QuantumNumberType::J].low();
171 auto NL = qid.val[QuantumNumberType::N].low();
173 JL, NL, GS, GR, GLE, B, D, H, gB, gD, gH, lB, lD, lH);
177 }
else if (qid.Isotopologue() ==
"CO-26") {
178 constexpr Numeric gperp =
182 return {gperp, gperp};
183 }
else if (qid.Isotopologue() ==
"OCS-622") {
184 constexpr Numeric gperp =
187 constexpr Numeric gpara =
190 if (qid.val.has(QuantumNumberType::J, QuantumNumberType::Ka)) {
191 auto JU = qid.val[QuantumNumberType::J].upp();
192 auto KU = qid.val[QuantumNumberType::Ka].upp();
193 auto JL = qid.val[QuantumNumberType::J].low();
194 auto KL = qid.val[QuantumNumberType::Ka].low();
199 }
else if (qid.Isotopologue() ==
"OCS-624") {
200 constexpr Numeric gperp =
203 constexpr Numeric gpara =
207 if (qid.val.has(QuantumNumberType::J, QuantumNumberType::Ka)) {
208 auto JU = qid.val[QuantumNumberType::J].upp();
209 auto KU = qid.val[QuantumNumberType::Ka].upp();
210 auto JL = qid.val[QuantumNumberType::J].low();
211 auto KL = qid.val[QuantumNumberType::Ka].low();
216 }
else if (qid.Isotopologue() ==
"CO2-626") {
217 constexpr Numeric gperp =
220 constexpr Numeric gpara =
224 if (qid.val.has(QuantumNumberType::J, QuantumNumberType::Ka)) {
225 auto JU = qid.val[QuantumNumberType::J].upp();
226 auto KU = qid.val[QuantumNumberType::Ka].upp();
227 auto JL = qid.val[QuantumNumberType::J].low();
228 auto KL = qid.val[QuantumNumberType::Ka].low();
246 return Eigen::Vector3d(
v,
u,
w).normalized();
252 return Eigen::Vector3d(cos(
a) * sin(z), sin(
a) * sin(z), cos(z)).normalized();
258 return Eigen::Vector3d(cos(
a) * cos(z), sin(
a) * cos(z), -sin(z)).normalized();
262 Numeric
u, Numeric
v, Numeric
w, Numeric z, Numeric
a)
noexcept {
266 output.
H = std::hypot(
u,
v,
w);
276 output.
dH_dv = nH[0];
277 output.
dH_du = nH[1];
278 output.
dH_dw = nH[2];
281 const Numeric cos_theta = n.dot(nH);
282 const Numeric sin_theta = std::sqrt(1 -
Math::pow2(cos_theta));
283 output.
theta = std::acos(cos_theta);
284 if (sin_theta not_eq 0) {
285 const Eigen::Vector3d dtheta =
286 (nH * cos_theta - n) / (output.
H * sin_theta);
297 const Eigen::Vector3d inplane = nH - nH.dot(n) * n;
298 const Numeric y = ev.cross(inplane).dot(n);
299 const Numeric x = ev.dot(inplane);
300 output.
eta = std::atan2(y, x);
301 if (x not_eq 0 or y not_eq 0) {
302 const Eigen::Vector3d deta =
325 auto ml =
Ml(Ju, Jl, type, n);
326 auto mu =
Mu(Ju, Jl, type, n);
327 auto dm = Rational(
dM(type));
329 pow2(
wigner3j(Jl, Rational(1), Ju, ml, -dm, -mu));
353 const Numeric ST = std::sin(theta), CT = std::cos(theta), ST2 = ST * ST,
354 CT2 = CT * CT, ST2C2E = ST2 * std::cos(2 * eta),
355 ST2S2E = ST2 * std::sin(2 * eta);
359 1 + CT2, ST2C2E, ST2S2E, 2 * CT, 4 * CT, 2 * ST2S2E, -2 * ST2C2E);
363 1 + CT2, ST2C2E, ST2S2E, -2 * CT, -4 * CT, 2 * ST2S2E, -2 * ST2C2E);
368 const Numeric eta)
noexcept {
369 const Numeric ST = std::sin(theta), CT = std::cos(theta),
370 C2E = std::cos(2 * eta), S2E = std::sin(2 * eta), dST = CT,
371 dST2 = 2 * ST * dST, dCT = -ST, dST2C2E = dST2 * C2E,
372 dST2S2E = dST2 * S2E, dCT2 = 2 * CT * dCT;
376 dCT2, dST2C2E, dST2S2E, 2 * dCT, 4 * dCT, 2 * dST2S2E, -2 * dST2C2E);
378 dST2, -dST2C2E, -dST2S2E, 0, 0, -2 * dST2S2E, 2 * dST2C2E);
380 dCT2, dST2C2E, dST2S2E, -2 * dCT, -4 * dCT, 2 * dST2S2E, -2 * dST2C2E);
385 Numeric eta)
noexcept {
386 const Numeric ST = std::sin(theta), ST2 = ST * ST, C2E = std::cos(2 * eta),
387 S2E = std::sin(2 * eta), dST2C2E = -2 * ST2 * S2E,
388 dST2S2E = 2 * ST2 * C2E;
394 0, -dST2C2E, -dST2S2E, 0, 0, -2 * dST2S2E, 2 * dST2C2E);
400#pragma GCC diagnostic push
401#pragma GCC diagnostic ignored "-Wreturn-type"
415#pragma GCC diagnostic pop
417void sum(PropagationMatrix& pm,
418 const ComplexVectorView& abs,
423 ARTS_ASSERT(pm.NumberOfFrequencies() == abs.nelem())
424 ARTS_ASSERT(do_phase ? pm.NumberOfNeededVectors() == 7
425 : pm.NumberOfNeededVectors() == 4)
427 const ExhaustiveConstVectorView pol_real(polvec.att);
428 const ExhaustiveConstVectorView pol_imag(polvec.dis);
430 MatrixView out = pm.Data()(0, 0, joker, joker);
431 matpack::eigen::mat(out).leftCols<4>().noalias() += matpack::eigen::row_vec(abs.real()) * matpack::eigen::col_vec(pol_real);
433 matpack::eigen::mat(out).rightCols<3>().noalias() +=
434 matpack::eigen::row_vec(abs.imag()) * matpack::eigen::col_vec(pol_imag);
437void dsum(PropagationMatrix& pm,
438 const ComplexVectorView& abs,
439 const ComplexVectorView& dabs,
449 ARTS_ASSERT(pm.NumberOfFrequencies() == abs.nelem())
450 ARTS_ASSERT(do_phase ? pm.NumberOfNeededVectors() == 7
451 : pm.NumberOfNeededVectors() == 4)
453 const ExhaustiveConstVectorView pol_real(polvec.att);
454 const ExhaustiveConstVectorView pol_imag(polvec.dis);
455 const ExhaustiveConstVectorView dpolvec_dtheta_real(dpolvec_dtheta.att);
456 const ExhaustiveConstVectorView dpolvec_dtheta_imag(dpolvec_dtheta.dis);
457 const ExhaustiveConstVectorView dpolvec_deta_real(dpolvec_deta.att);
458 const ExhaustiveConstVectorView dpolvec_deta_imag(dpolvec_deta.dis);
461 (dt * matpack::eigen::col_vec(dpolvec_dtheta_real) + de * matpack::eigen::col_vec(dpolvec_deta_real));
463 (dt * matpack::eigen::col_vec(dpolvec_dtheta_imag) + de * matpack::eigen::col_vec(dpolvec_deta_imag));
465 MatrixView out = pm.Data()(0, 0, joker, joker);
466 matpack::eigen::mat(out).leftCols<4>().noalias() +=
467 dH * matpack::eigen::row_vec(dabs.real()) * matpack::eigen::col_vec(pol_real) + matpack::eigen::row_vec(abs.real()) * da_r;
469 matpack::eigen::mat(out).rightCols<3>().noalias() +=
470 dH * matpack::eigen::row_vec(dabs.imag()) * matpack::eigen::col_vec(pol_imag) + matpack::eigen::row_vec(abs.imag()) * da_i;
A list of many quantum numbers. Should always remain sorted.
Numeric Strength(Rational Ju, Rational Jl, Polarization type, Index n) const ARTS_NOEXCEPT
Gives the strength of one subline of a given polarization.
bool empty() const noexcept
Returns true if the Model represents no Zeeman effect.
constexpr Model(SplittingData gs={NAN, NAN}) noexcept
Default copy/init of Model from its only private variable.
Binary output file stream class.
Binary output file stream class.
Input manipulator class for doubles to enable nan and inf parsing.
Helper macros for debugging.
#define ARTS_ASSERT(condition,...)
constexpr Numeric mass_ratio_electrons_per_proton
Mass ratio of electrons to protons [-] From: https://physics.nist.gov/cgi-bin/cuu/Value?...
constexpr auto pow3(auto x) noexcept
power of three
constexpr auto pow2(auto x) noexcept
power of two
bool vamdcCheck(const ValueList &l, VAMDC type) ARTS_NOEXCEPT
Implements Zeeman modeling.
std::ostream & operator<<(std::ostream &os, const Model &m)
constexpr Index dM(Polarization type) noexcept
Gives the change of M given a polarization type.
constexpr Numeric SimpleGCaseB(Rational N, Rational J, Rational Lambda, Rational S, Numeric GS, Numeric GL) noexcept
Computes the Zeeman splitting coefficient.
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.
constexpr Numeric SimpleGCaseA(Rational Omega, Rational J, Rational Lambda, Rational Sigma, Numeric GS, Numeric GL) noexcept
Computes the Zeeman splitting coefficient.
constexpr Numeric PolarizationFactor(Polarization type) noexcept
The renormalization factor of a polarization type.
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.
Polarization
Zeeman polarization selection.
Model GetAdvancedModel(const QuantumIdentifier &qid) ARTS_NOEXCEPT
Returns an advanced Zeeman model.
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.
std::istream & operator>>(std::istream &is, Model &m)
Model GetSimpleModel(const QuantumIdentifier &qid) ARTS_NOEXCEPT
Returns a simple Zeeman model.
constexpr Rational Ml(Rational Ju, Rational Jl, Polarization type, Index n) noexcept
Gives the lower state M value at an index.
AllPolarizationVectors AllPolarization(Numeric theta, Numeric eta) noexcept
Computes the polarization of each polarization type.
constexpr Rational Mu(Rational Ju, Rational Jl, Polarization type, Index n) noexcept
Gives the upper state M value at an index.
Numeric get_lande_spin_constant(const Species::Species species) noexcept
Get the Lande spin constant.
Numeric get_lande_lambda_constant() noexcept
Get the Lande Lambda constant.
Some molecular constants.
A logical struct for global quantum numbers with species identifiers.
PolarizationVector for each Polarization.
Contains derived values useful for Zeeman calculations.
Polarization vector for Zeeman Propagation Matrix.
Main storage for Zeeman splitting coefficients.
Numeric wigner3j(const Rational j1, const Rational j2, const Rational j3, const Rational m1, const Rational m2, const Rational m3)
Wigner 3J symbol.
Wigner symbol interactions.
Eigen::Vector3d los_xyz_by_za_local(Numeric z, Numeric a)
Numeric case_b_g_coefficient_o2(Rational J, Rational N, Numeric GS, Numeric GR, Numeric GLE, Numeric B, Numeric D, Numeric H, Numeric gB, Numeric gD, Numeric gH, Numeric lB, Numeric lD, Numeric lH)
Zeeman::SplittingData SimpleG(const Quantum::Number::ValueList &qns, const Numeric &GS, const Numeric &GL) noexcept
constexpr Numeric closed_shell_trilinear(Rational k, Rational j, Numeric gperp, Numeric gpara)
Eigen::Vector3d ev_xyz_by_za_local(Numeric z, Numeric a)
Eigen::Vector3d los_xyz_by_uvw_local(Numeric u, Numeric v, Numeric w)
Headers and class definition of Zeeman modeling.