5#include <Faddeeva/Faddeeva.hh>
13#include "matpack_complex.h"
14#include "matpack_data.h"
15#include "matpack_eigen.h"
25#define WIGNER3 fw3jja6
37 const Rational& f)
noexcept {
39 a.toInt(2),
b.toInt(2),
c.toInt(2),
d.toInt(2), e.toInt(2), f.toInt(2));
47 const Rational& f)
noexcept {
49 a.toInt(2),
b.toInt(2),
c.toInt(2),
d.toInt(2), e.toInt(2), f.toInt(2));
58 const Vector& dip) noexcept :
59 val(pop.nelem(), 0), str(pop.nelem(), 0) {
77 const Index n = pop.nelem();
80 ComplexMatrix V(n, n);
83 diagonalize(V, val, W);
86 for (Index i=0; i<n; i++) {
87 for (Index j=0; j<n; j++) {
88 str[i] += dip[j] * V(j, i);
94 for (Index i=0; i<n; i++) {
96 for (Index j=0; j<n; j++) {
97 z += pop[j] * dip[j] * V(i, j);
105 const Index N =
val.nelem();
109 for (Index i=0; i<N; i++) {
110 for (Index j=i+1; j<N; j++) {
111 if (
val.real()[j] <
val.real()[i]) {
112 std::swap(
val[i],
val[j]);
113 std::swap(
str[i],
str[j]);
119 const ComplexVector strc =
str;
120 const ComplexVector valc =
val;
121 for (Index i=0; i<N; i++) {
122 f[i] = fc[sorting[i]];
123 str[i] = strc[sorting[i]];
124 val[i] = valc[sorting[i]];
129 const Index n = eqv.
str.nelem();
130 for (Index i=0; i<n; i++) {
132 os << eqv.
str[i] <<
' ' << eqv.
val[i];
138 os << Species::toShortName(srbd.
spec);
140 os <<
' ' << srbd.
beta;
143 os <<
' ' << srbd.
mass;
148 std::string spec_name;
155 srbd.
spec = Species::fromShortName(spec_name);
160namespace Makarov2020etal {
169 return (iseven(Jl + N) ? 1 : -1) * sqrt(6 * (2*Jl + 1) * (2*Ju + 1)) *
wigner6j(1, 1, 1, Jl, Ju, N);
176 pop(band.NumLines()), dip(band.NumLines()) {
181 const Numeric ratiopart = QT0 / QT;
183 for (Index i=0; i<n; i++) {
189 if (band.
population == Absorption::PopulationType::ByMakarovFullRelmat) {
190 auto& J = band.
lines[i].localquanta.val[QuantumNumberType::J];
191 auto& N = band.
lines[i].localquanta.val[QuantumNumberType::N];
193 }
else if (band.
population == Absorption::PopulationType::ByRovibLinearDipoleLineMixing) {
201 const Index N = pop.
nelem();
205 std::iota(out.begin(), out.end(), 0);
209 for (Index i=0; i<N; i++) {
210 s[i] = band.lines[i].F0 * pop[i] *
Math::pow2(dip[i]);
214 for (Index i=0; i<N; i++) {
215 for (Index j=i+1; j<N; j++) {
217 std::swap(out[i], out[j]);
218 std::swap(dip[i], dip[j]);
219 std::swap(pop[i], pop[j]);
220 std::swap(s[i], s[j]);
230 const Index N = presorting.size();
231 Vector dipcopy(dip), popcopy(pop);
232 for (Index i=0; i<N; i++) {
233 dip[i] = dipcopy[presorting[i]];
234 pop[i] = popcopy[presorting[i]];
248 return {tp.
sort(band), tp};
270 const Numeric energy)
const {
276 const Numeric other_mass,
277 const Numeric energy_x,
278 const Numeric energy_xm2)
const {
287 constexpr Numeric
fac = 8 * k / (m_u * pi);
289 const Numeric wnnm2 = (energy_x - energy_xm2) / h_bar;
291 const Numeric inv_eff_mass = 1 /
mass + 1 / other_mass;
292 const Numeric v_bar_pow2 =
fac*T*inv_eff_mass;
295 return 1.0 / pow2(1 + pow2(wnnm2) * tauc_pow2 / 24.0);
298namespace Makarov2020etal {
308template<
bool rescale_pure_rotational=true>
constexpr
309Numeric
erot(
const Rational N,
const Rational j=-1) {
310 const Rational J = j<0 ? N : j;
312 if constexpr (rescale_pure_rotational) {
313 return erot<false>(N, J) - erot<false>(1, 0);
319 constexpr Numeric B0=43100.4425e0;
320 constexpr Numeric D0=.145123e0;
321 constexpr Numeric H0=3.8e-08;
322 constexpr Numeric xl0=59501.3435e0;
323 constexpr Numeric xg0=-252.58633e0;
324 constexpr Numeric xl1=0.058369e0;
325 constexpr Numeric xl2=2.899e-07;
326 constexpr Numeric xg1=-2.4344e-04;
327 constexpr Numeric xg2=-1.45e-09;
329 const auto XN = Numeric(N);
330 const Numeric XX = XN * (XN + 1);
331 const Numeric xlambda=xl0+xl1*XX+xl2*pow2(XX);
332 const Numeric xgama=xg0+xg1*XX+xg2*pow2(XX);
333 const Numeric C1=B0*XX-D0*pow2(XX)+H0*pow3(XX);
337 return mhz2joule(C1 - (xlambda+B0*(2.*XN-1.)+xgama*XN));
338 return mhz2joule(C1 - (xlambda+B0*(2.*XN-1.)+xgama*XN) + std::sqrt(pow2(B0*(2.*XN-1.))+pow2(xlambda)-2.*B0*xlambda));
341 return mhz2joule(C1 - (xlambda-B0*(2.*XN+3.)-xgama*(XN+1.)) - std::sqrt(pow2(B0*(2.*XN+3.))+pow2(xlambda)-2.*B0*xlambda));
342 return mhz2joule(C1);
364 const auto bk = [](
const Rational& r) -> Numeric {
return sqrt(2 * r + 1); };
370 const Rational Si = S.upp();
371 const Rational Sf = S.low();
374 for (Index i = 0; i < n; i++) {
375 auto& J = band.
lines[sorting[i]].localquanta.val[QuantumNumberType::J];
376 auto& N = band.
lines[sorting[i]].localquanta.val[QuantumNumberType::N];
381 const auto maxL =
temp_init_size(band.
max(QuantumNumberType::J), band.
max(QuantumNumberType::N));
383 const auto Om = [&]() {
384 std::vector<Numeric> out(maxL);
385 for (Index i = 0; i < maxL; i++)
386 out[i] = rovib_data.
Omega(
391 const auto Q = [&]() {
392 std::vector<Numeric> out(maxL);
393 for (Index i = 0; i < maxL; i++)
394 out[i] = rovib_data.
Q(i, T, band.
T0,
erot(i));
398 wig_thread_temp_init(maxL);
399 for (Index i = 0; i < n; i++) {
400 auto& J = band.
lines[sorting[i]].localquanta.val[QuantumNumberType::J];
401 auto& N = band.
lines[sorting[i]].localquanta.val[QuantumNumberType::N];
403 const Rational Ji = J.upp();
404 const Rational Jf = J.low();
405 const Rational Ni = N.upp();
406 const Rational Nf = N.low();
408 for (Index j = 0; j < n; j++) {
409 if (i == j)
continue;
411 auto& J_p = band.
lines[sorting[j]].localquanta.val[QuantumNumberType::J];
412 auto& N_p = band.
lines[sorting[j]].localquanta.val[QuantumNumberType::N];
414 const Rational Ji_p = J_p.upp();
415 const Rational Jf_p = J_p.low();
416 const Rational Ni_p = N_p.upp();
417 const Rational Nf_p = N_p.low();
419 if (Jf_p > Jf)
continue;
424 const Numeric scl = (iseven(Ji_p + Ji + 1) ? 1 : -1) * bk(Ni) * bk(Nf) *
425 bk(Nf_p) * bk(Ni_p) * bk(Jf) * bk(Jf_p) * bk(Ji) *
427 const auto [L0, L1] =
429 {Rational(2), std::numeric_limits<Index>::max()});
430 for (Rational L = L0; L <= L1; L += 2) {
431 const Numeric
a =
wig3(Ni_p, Ni, L, 0, 0, 0);
432 const Numeric
b =
wig3(Nf_p, Nf, L, 0, 0, 0);
433 const Numeric
c =
wig6(L, Ji, Ji_p, Si, Ni_p, Ni);
434 const Numeric
d =
wig6(L, Jf, Jf_p, Sf, Nf_p, Nf);
435 const Numeric e =
wig6(L, Ji, Ji_p, 1, Jf_p, Jf);
436 sum +=
a *
b *
c *
d * e * Numeric(2 * L + 1) * Q[L.toIndex()] / Om[L.toIndex()];
438 sum *= scl * Om[Ni.toIndex()];
442 W(j, i) = sum * std::exp((
erot(Nf_p) -
erot(Nf)) / kelvin2joule(T));
450 for (Index i = 0; i < n; i++) {
454 for (Index j = 0; j < n; j++) {
456 sumlw += dipr[j] * W(j, i);
458 sumup += dipr[j] * W(j, i);
463 band.
lines[sorting[i]].localquanta.val[QuantumNumberType::N].low();
464 for (Index j = i + 1; j < n; j++) {
466 band.
lines[sorting[j]].localquanta.val[QuantumNumberType::N].low();
471 W(j, i) *= -sumup / sumlw;
472 W(i, j) = W(j, i) * std::exp((
erot(Ni) -
erot(Nj)) /
489namespace LinearRovibErrorCorrectedSudden {
492 if (isot.spec == Species::Species::CarbonDioxide and isot.isotname ==
"626") {
496 ARTS_USER_ERROR(isot.FullName(),
" has no rotational energies in ARTS")
497 return [](
const Rational J) -> Numeric {
return Numeric(J) * std::numeric_limits<Numeric>::signaling_NaN();};
512 Rational li = l2.upp();
513 Rational lf = l2.low();
515 const bool swap_order = li > lf;
516 if (swap_order) swap(li, lf);
517 const int sgn = iseven(li + lf + 1) ? -1 : 1;
518 if (abs(li - lf) > 1)
return;
523 for (Index i=0; i<n; i++) {
524 auto& J = band.
lines[i].localquanta.val[QuantumNumberType::J];
528 for (Index i=0; i<n; i++) {
529 auto& J = band.
lines[sorting[i]].localquanta.val[QuantumNumberType::J];
530 Rational Ji = J.upp();
531 Rational Jf = J.low();
532 if (swap_order) swap(Ji, Jf);
534 for (Index j=0; j<n; j++) {
536 auto& J_p = band.
lines[sorting[j]].localquanta.val[QuantumNumberType::J];
537 Rational Ji_p = J_p.upp();
538 Rational Jf_p = J_p.low();
539 if (swap_order) swap(Ji_p, Jf_p);
542 if (Jf_p > Jf)
continue;
544 Index L = std::max(std::abs((Ji-Ji_p).toIndex()), std::abs((Jf-Jf_p).toIndex()));
546 const Index Lf = std::min((Ji+Ji_p).toIndex(), (Jf+Jf_p).toIndex());
549 for (; L<=Lf; L+=2) {
550 const Numeric
a =
wigner3j(Ji_p, L, Ji, li, 0, -li);
551 const Numeric
b =
wigner3j(Jf_p, L, Jf, lf, 0, -lf);
552 const Numeric
c =
wigner6j(Ji, Jf, 1, Jf_p, Ji_p, L);
553 const Numeric QL = rovib_data.
Q(L, T, band.
T0, erot(L));
554 const Numeric ECS = rovib_data.
Omega(T, band.
T0, band.
SpeciesMass(), erot(L), erot(L-2));
555 sum +=
a *
b *
c * Numeric(2 * L + 1) * QL / ECS;
557 const Numeric ECS = rovib_data.
Omega(T, band.
T0, band.
SpeciesMass(), erot(Ji), erot(Ji-2));
558 const Numeric scl = sgn * ECS * Numeric(2 * Ji_p + 1) * sqrt((2 * Jf + 1) * (2 * Jf_p + 1));
563 W(i, j) = sum * std::exp((erot(Jf_p) - erot(Jf)) / kelvin2joule(T));
568 for (Index i=0; i<n; i++)
569 for (Index j=0; j<n; j++)
570 if (j not_eq i and W(i, j) > 0)
574 for (Index i=0; i<n; i++) {
578 for (Index j=0; j<n; j++) {
580 sumlw += std::abs(dipr[j]) * W(j, i);
582 sumup += std::abs(dipr[j]) * W(j, i);
586 const Rational Ji = band.
lines[sorting[i]].localquanta.val[QuantumNumberType::J].low();
587 for (Index j=i+1; j<n; j++) {
588 const Rational Jj = band.
lines[sorting[j]].localquanta.val[QuantumNumberType::J].low();
593 W(j, i) *= - sumup / sumlw;
594 W(i, j) = W(j, i) * std::exp((erot(Ji) - erot(Jj)) / kelvin2joule(T));
618 const Index species_pos) {
622 ComplexMatrix W(N, N, 0);
625 for (Index I=0; I<N; I++) {
626 const Index i = sorting[I];
628 W(I, I) = Complex(shape.D0, shape.G0);
633 case PopulationType::ByMakarovFullRelmat:
636 case PopulationType::ByRovibLinearDipoleLineMixing: {
641 "\nin this code. It must either be added or computations aborted earlier");
666 const Numeric frenorm) {
668 const Index M = vmrs.nelem();
671 ComplexMatrix W(N, N, 0);
674 for (Index k=0; k<M; k++) {
675 if (vmrs[k] == 0)
continue;
678 matpack::eigen::mat(W).noalias() += vmrs[k] * matpack::eigen::mat(
683 for (Index i=0; i<N; i++) {
684 real_val(W(i, i)) += band.
lines[sorting[i]].F0 - frenorm;
694 const Numeric this_vmr,
697 const Vector& f_grid,
703 const Numeric frenorm = band.
F_mean();
718 std::pair<ComplexVector, bool> retval{ComplexVector(f_grid.nelem(), 0),
false};
719 auto& [absorption,works] = retval;
722 for (Index i=0; i<band.
NumLines(); i++) {
724 const Index nz = band.
ZeemanCount(i, zeeman_polarization);
725 for (Index j=0; j<nz; j++) {
726 const Numeric Sz = band.
ZeemanStrength(i, zeeman_polarization, j);
727 const Numeric dzeeman = H * band.
ZeemanSplitting(i, zeeman_polarization, j);
730 for (Index iv=0; iv<f_grid.nelem(); iv++) {
734 const Numeric gamd = GD_div_F0 * (eqv.
val[i].real() + frenorm + dzeeman);
736 for (Index iv=0; iv<f_grid.nelem(); iv++) {
737 const Complex z = (eqv.
val[i] + frenorm + dzeeman - f_grid[iv]) * cte;
738 const Complex
w = Faddeeva::w(z);
739 absorption[iv] += Sz * eqv.
str[i] *
w / gamd;
747 for (Index iv=0; iv<f_grid.nelem(); iv++) {
748 const Numeric f = f_grid[iv];
750 absorption[iv] *= fact * numdens * sq_ln2pi;
753 if (
auto&
a = real_val(absorption[iv]); not std::isnormal(
a) or
a < 0) {
766 const Numeric this_vmr,
769 const Vector& f_grid,
773 auto [absorption, work] =
ecs_absorption_impl(T, H, P, this_vmr, vmrs, ecs_data, f_grid, zeeman_polarization, band);
776 ArrayOfComplexVector jacobian(jacobian_quantities.
nelem(), absorption);
778 for (Index i=0; i<jacobian_quantities.
nelem(); i++) {
779 auto& vec = jacobian[i];
780 auto& target = jacobian_quantities[i].Target();
782 if (target == Jacobian::Atm::Temperature) {
783 const Numeric dT = target.perturbation;
784 const auto [dabs, dwork] =
ecs_absorption_impl(T+dT, H, P, this_vmr, vmrs, ecs_data, f_grid, zeeman_polarization, band);
788 }
else if (target.isMagnetic()) {
789 const Numeric dH = target.perturbation;
790 const auto [dabs, dwork] =
ecs_absorption_impl(T, H+dH, P, this_vmr, vmrs, ecs_data, f_grid, zeeman_polarization, band);
794 }
else if (target.isWind()) {
795 const Numeric df = target.perturbation;
796 Vector f_grid_copy = f_grid;
798 const auto [dabs, dwork] =
ecs_absorption_impl(T, H, P, this_vmr, vmrs, ecs_data, f_grid_copy, zeeman_polarization, band);
802 }
else if (target == Jacobian::Line::VMR) {
804 Vector vmrs_copy = vmrs;
805 Numeric this_vmr_copy = this_vmr;
806 const Numeric dvmr = target.perturbation;
810 this_vmr_copy += dvmr;
815 vmrs_copy[j] += dvmr;
821 const auto [dabs, dwork] =
ecs_absorption_impl(T, H, P, this_vmr_copy, vmrs_copy, ecs_data, f_grid, zeeman_polarization, band);
826 }
else if (target.needQuantumIdentity()) {
829 for (Index iline=0; iline<band.
NumLines(); iline++) {
834 switch (target.line) {
835 case Jacobian::Line::ShapeG0X0:
836 d *= band.
lines[iline].lineshape[pos].G0().X0;
837 band_copy.
lines[iline].lineshape[pos].G0().X0 +=
d;
839 case Jacobian::Line::ShapeG0X1:
840 d *= band.
lines[iline].lineshape[pos].G0().X1;
841 band_copy.
lines[iline].lineshape[pos].G0().X1 +=
d;
843 case Jacobian::Line::ShapeG0X2:
844 d *= band.
lines[iline].lineshape[pos].G0().X2;
845 band_copy.
lines[iline].lineshape[pos].G0().X2 +=
d;
847 case Jacobian::Line::ShapeG0X3:
848 d *= band.
lines[iline].lineshape[pos].G0().X3;
849 band_copy.
lines[iline].lineshape[pos].G0().X3 +=
d;
851 case Jacobian::Line::ShapeD0X0:
852 d *= band.
lines[iline].lineshape[pos].D0().X0;
853 band_copy.
lines[iline].lineshape[pos].D0().X0 +=
d;
855 case Jacobian::Line::ShapeD0X1:
856 d *= band.
lines[iline].lineshape[pos].D0().X1;
857 band_copy.
lines[iline].lineshape[pos].D0().X1 +=
d;
859 case Jacobian::Line::ShapeD0X2:
860 d *= band.
lines[iline].lineshape[pos].D0().X2;
861 band_copy.
lines[iline].lineshape[pos].D0().X2 +=
d;
863 case Jacobian::Line::ShapeD0X3:
864 d *= band.
lines[iline].lineshape[pos].D0().X3;
865 band_copy.
lines[iline].lineshape[pos].D0().X3 +=
d;
867 case Jacobian::Line::ShapeG2X0:
868 d *= band.
lines[iline].lineshape[pos].G2().X0;
869 band_copy.
lines[iline].lineshape[pos].G2().X0 +=
d;
871 case Jacobian::Line::ShapeG2X1:
872 d *= band.
lines[iline].lineshape[pos].G2().X1;
873 band_copy.
lines[iline].lineshape[pos].G2().X1 +=
d;
875 case Jacobian::Line::ShapeG2X2:
876 d *= band.
lines[iline].lineshape[pos].G2().X2;
877 band_copy.
lines[iline].lineshape[pos].G2().X2 +=
d;
879 case Jacobian::Line::ShapeG2X3:
880 d *= band.
lines[iline].lineshape[pos].G2().X3;
881 band_copy.
lines[iline].lineshape[pos].G2().X3 +=
d;
883 case Jacobian::Line::ShapeD2X0:
884 d *= band.
lines[iline].lineshape[pos].D2().X0;
885 band_copy.
lines[iline].lineshape[pos].D2().X0 +=
d;
887 case Jacobian::Line::ShapeD2X1:
888 d *= band.
lines[iline].lineshape[pos].D2().X1;
889 band_copy.
lines[iline].lineshape[pos].D2().X1 +=
d;
891 case Jacobian::Line::ShapeD2X2:
892 d *= band.
lines[iline].lineshape[pos].D2().X2;
893 band_copy.
lines[iline].lineshape[pos].D2().X2 +=
d;
895 case Jacobian::Line::ShapeD2X3:
896 d *= band.
lines[iline].lineshape[pos].D2().X3;
897 band_copy.
lines[iline].lineshape[pos].D2().X3 +=
d;
899 case Jacobian::Line::ShapeFVCX0:
900 d *= band.
lines[iline].lineshape[pos].FVC().X0;
901 band_copy.
lines[iline].lineshape[pos].FVC().X0 +=
d;
903 case Jacobian::Line::ShapeFVCX1:
904 d *= band.
lines[iline].lineshape[pos].FVC().X1;
905 band_copy.
lines[iline].lineshape[pos].FVC().X1 +=
d;
907 case Jacobian::Line::ShapeFVCX2:
908 d *= band.
lines[iline].lineshape[pos].FVC().X2;
909 band_copy.
lines[iline].lineshape[pos].FVC().X2 +=
d;
911 case Jacobian::Line::ShapeFVCX3:
912 d *= band.
lines[iline].lineshape[pos].FVC().X3;
913 band_copy.
lines[iline].lineshape[pos].FVC().X3 +=
d;
915 case Jacobian::Line::ShapeETAX0:
916 d *= band.
lines[iline].lineshape[pos].ETA().X0;
917 band_copy.
lines[iline].lineshape[pos].ETA().X0 +=
d;
919 case Jacobian::Line::ShapeETAX1:
920 d *= band.
lines[iline].lineshape[pos].ETA().X1;
921 band_copy.
lines[iline].lineshape[pos].ETA().X1 +=
d;
923 case Jacobian::Line::ShapeETAX2:
924 d *= band.
lines[iline].lineshape[pos].ETA().X2;
925 band_copy.
lines[iline].lineshape[pos].ETA().X2 +=
d;
927 case Jacobian::Line::ShapeETAX3:
928 d *= band.
lines[iline].lineshape[pos].ETA().X3;
929 band_copy.
lines[iline].lineshape[pos].ETA().X3 +=
d;
931 case Jacobian::Line::Center:
932 d *= band.
lines[iline].F0;
933 band_copy.
lines[iline].F0 +=
d;
935 case Jacobian::Line::Strength:
936 d *= band.
lines[iline].I0;
937 band_copy.
lines[iline].I0 +=
d;
939 case Jacobian::Line::ShapeYX0:
940 case Jacobian::Line::ShapeYX1:
941 case Jacobian::Line::ShapeYX2:
942 case Jacobian::Line::ShapeYX3:
943 case Jacobian::Line::ShapeGX0:
944 case Jacobian::Line::ShapeGX1:
945 case Jacobian::Line::ShapeGX2:
946 case Jacobian::Line::ShapeGX3:
947 case Jacobian::Line::ShapeDVX0:
948 case Jacobian::Line::ShapeDVX1:
949 case Jacobian::Line::ShapeDVX2:
950 case Jacobian::Line::ShapeDVX3:
951 case Jacobian::Line::NLTE:
952 case Jacobian::Line::VMR:
953 case Jacobian::Line::ECS_SCALINGX0:
954 case Jacobian::Line::ECS_SCALINGX1:
955 case Jacobian::Line::ECS_SCALINGX2:
956 case Jacobian::Line::ECS_SCALINGX3:
957 case Jacobian::Line::ECS_BETAX0:
958 case Jacobian::Line::ECS_BETAX1:
959 case Jacobian::Line::ECS_BETAX2:
960 case Jacobian::Line::ECS_BETAX3:
961 case Jacobian::Line::ECS_LAMBDAX0:
962 case Jacobian::Line::ECS_LAMBDAX1:
963 case Jacobian::Line::ECS_LAMBDAX2:
964 case Jacobian::Line::ECS_LAMBDAX3:
965 case Jacobian::Line::ECS_DCX0:
966 case Jacobian::Line::ECS_DCX1:
967 case Jacobian::Line::ECS_DCX2:
968 case Jacobian::Line::ECS_DCX3:
969 case Jacobian::Line::FINAL: {
975 const auto [dabs, dwork] =
ecs_absorption_impl(T, H, P, this_vmr, vmrs, ecs_data, f_grid, zeeman_polarization, band_copy);
982 const auto spec = target.qid.Species();
983 ARTS_USER_ERROR_IF(
const Index pos = ecs_data.
pos(spec); pos == ecs_data.
size(),
"No data for species ", spec,
" in ecs_data:\n", ecs_data)
984 switch (target.line) {
985 case Jacobian::Line::ECS_SCALINGX0:
986 d *= ecs_data_copy[spec].scaling.X0;
987 ecs_data_copy[spec].scaling.X0 +=
d;
989 case Jacobian::Line::ECS_SCALINGX1:
990 d *= ecs_data_copy[spec].scaling.X1;
991 ecs_data_copy[spec].scaling.X1 +=
d;
993 case Jacobian::Line::ECS_SCALINGX2:
994 d *= ecs_data_copy[spec].scaling.X2;
995 ecs_data_copy[spec].scaling.X2 +=
d;
997 case Jacobian::Line::ECS_SCALINGX3:
998 d *= ecs_data_copy[spec].scaling.X3;
999 ecs_data_copy[spec].scaling.X3 +=
d;
1001 case Jacobian::Line::ECS_BETAX0:
1002 d *= ecs_data_copy[spec].beta.X0;
1003 ecs_data_copy[spec].beta.X0 +=
d;
1005 case Jacobian::Line::ECS_BETAX1:
1006 d *= ecs_data_copy[spec].beta.X1;
1007 ecs_data_copy[spec].beta.X1 +=
d;
1009 case Jacobian::Line::ECS_BETAX2:
1010 d *= ecs_data_copy[spec].beta.X2;
1011 ecs_data_copy[spec].beta.X2 +=
d;
1013 case Jacobian::Line::ECS_BETAX3:
1014 d *= ecs_data_copy[spec].beta.X3;
1015 ecs_data_copy[spec].beta.X3 +=
d;
1017 case Jacobian::Line::ECS_LAMBDAX0:
1018 d *= ecs_data_copy[spec].lambda.X0;
1019 ecs_data_copy[spec].lambda.X0 +=
d;
1021 case Jacobian::Line::ECS_LAMBDAX1:
1022 d *= ecs_data_copy[spec].lambda.X1;
1023 ecs_data_copy[spec].lambda.X1 +=
d;
1025 case Jacobian::Line::ECS_LAMBDAX2:
1026 d *= ecs_data_copy[spec].lambda.X1;
1027 ecs_data_copy[spec].lambda.X1 +=
d;
1029 case Jacobian::Line::ECS_LAMBDAX3:
1030 d *= ecs_data_copy[spec].lambda.X0;
1031 ecs_data_copy[spec].lambda.X0 +=
d;
1033 case Jacobian::Line::ECS_DCX0:
1034 d *= ecs_data_copy[spec].collisional_distance.X0;
1035 ecs_data_copy[spec].collisional_distance.X0 +=
d;
1037 case Jacobian::Line::ECS_DCX1:
1038 d *= ecs_data_copy[spec].collisional_distance.X1;
1039 ecs_data_copy[spec].collisional_distance.X1 +=
d;
1041 case Jacobian::Line::ECS_DCX2:
1042 d *= ecs_data_copy[spec].collisional_distance.X2;
1043 ecs_data_copy[spec].collisional_distance.X2 +=
d;
1045 case Jacobian::Line::ECS_DCX3:
1046 d *= ecs_data_copy[spec].collisional_distance.X3;
1047 ecs_data_copy[spec].collisional_distance.X3 +=
d;
1049 case Jacobian::Line::ShapeG0X0:
1050 case Jacobian::Line::ShapeG0X1:
1051 case Jacobian::Line::ShapeG0X2:
1052 case Jacobian::Line::ShapeG0X3:
1053 case Jacobian::Line::ShapeD0X0:
1054 case Jacobian::Line::ShapeD0X1:
1055 case Jacobian::Line::ShapeD0X2:
1056 case Jacobian::Line::ShapeD0X3:
1057 case Jacobian::Line::ShapeG2X0:
1058 case Jacobian::Line::ShapeG2X1:
1059 case Jacobian::Line::ShapeG2X2:
1060 case Jacobian::Line::ShapeG2X3:
1061 case Jacobian::Line::ShapeD2X0:
1062 case Jacobian::Line::ShapeD2X1:
1063 case Jacobian::Line::ShapeD2X2:
1064 case Jacobian::Line::ShapeD2X3:
1065 case Jacobian::Line::ShapeFVCX0:
1066 case Jacobian::Line::ShapeFVCX1:
1067 case Jacobian::Line::ShapeFVCX2:
1068 case Jacobian::Line::ShapeFVCX3:
1069 case Jacobian::Line::ShapeETAX0:
1070 case Jacobian::Line::ShapeETAX1:
1071 case Jacobian::Line::ShapeETAX2:
1072 case Jacobian::Line::ShapeETAX3:
1073 case Jacobian::Line::Center:
1074 case Jacobian::Line::Strength:
1075 case Jacobian::Line::ShapeYX0:
1076 case Jacobian::Line::ShapeYX1:
1077 case Jacobian::Line::ShapeYX2:
1078 case Jacobian::Line::ShapeYX3:
1079 case Jacobian::Line::ShapeGX0:
1080 case Jacobian::Line::ShapeGX1:
1081 case Jacobian::Line::ShapeGX2:
1082 case Jacobian::Line::ShapeGX3:
1083 case Jacobian::Line::ShapeDVX0:
1084 case Jacobian::Line::ShapeDVX1:
1085 case Jacobian::Line::ShapeDVX2:
1086 case Jacobian::Line::ShapeDVX3:
1087 case Jacobian::Line::NLTE:
1088 case Jacobian::Line::VMR:
1089 case Jacobian::Line::FINAL: {
1095 const auto [dabs, dwork] =
ecs_absorption_impl(T, H, P, this_vmr, vmrs, ecs_data_copy, f_grid, zeeman_polarization, band);
1107 return {std::move(absorption), std::move(jacobian), not work};
1112 const Tensor4& tempdata,
1113 const Vector& temperatures,
1121 for (Index i=0; i<S; i++) {
1126 for (Index j=0; j<N; j++) {
1127 targ = tempdata(1, j, i, joker);
1134 if (not fit)
return EXIT_FAILURE;
1140 for (Index j=0; j<N; j++) {
1141 targ = tempdata(2, j, i, joker);
1148 if (not fit)
return EXIT_FAILURE;
1155 for (Index j=0; j<N; j++) {
1156 targ = tempdata(0, j, i, joker);
1163 if (not fit)
return EXIT_FAILURE;
1170 for (Index j=0; j<N; j++) {
1171 targ = tempdata(3, j, i, joker);
1174 targ /= P0 * P0 * P0;
1178 if (not fit)
return EXIT_FAILURE;
1188 band.
population = Absorption::PopulationType::LTE;
1190 return EXIT_SUCCESS;
1202 const ConstMatrixView& W,
1204 const Index N = dip.nelem();
1210 for (Index k=0; k<N; k++) {
1211 for (Index j=0; j<N; j++) {
1212 if (k == j)
continue;
1213 Y[k] += 2 * std::abs(dip[j] / dip[k]) * W(j, k) / (band.
lines[k].F0 - band.
lines[j].F0);
1229 const ConstMatrixView& W,
1231 const Index N = dip.nelem();
1237 for (Index k=0; k<N; k++) {
1238 for (Index j=0; j<N; j++) {
1239 if (k == j)
continue;
1241 G[k] +=
Math::pow2(std::abs(dip[j] / dip[k]) * W(j, k) / (band.
lines[j].F0 - band.
lines[k].F0));
1242 G[k] += 2 * std::abs(dip[j] / dip[k]) * W(j, k) * W(k, k) /
Math::pow2(band.
lines[j].F0 - band.
lines[k].F0);
1243 for (Index l=0; l<N; l++) {
1244 if (l == k or l == j)
continue;
1245 G[k] -= 2 * std::abs(dip[j] / dip[k]) * W(j, l) * W(l, k) / ((band.
lines[j].F0 - band.
lines[k].F0) * (band.
lines[l].F0 - band.
lines[k].F0));
1263 const ConstMatrixView& W,
1265 const Index N = dip.nelem();
1271 for (Index k=0; k<N; k++) {
1272 for (Index j=0; j<N; j++) {
1273 if (k == j)
continue;
1274 DV[k] += W(k, j) * W(j, k) / (band.
lines[j].F0 - band.
lines[k].F0);
1283 const ComplexMatrix& W,
1287 const Numeric frenorm,
1292 const Index broadener) {
1298 for (Index i=0; i<pop.size(); i++) {
1299 f0[i] = band.
lines[i].F0 - frenorm + P * band.
lines[i].lineshape[broadener].D0().at(T, band.
T0);
1304 std::iota(sorting.begin(), sorting.end(), 0);
1305 for (Index i=0; i<band.
NumLines(); i++) {
1306 for (Index j=i+1; j<band.
NumLines(); j++) {
1307 if (f0[j] < f0[i]) {
1308 std::swap(sorting[i], sorting[j]);
1309 std::swap(f0[i], f0[j]);
1319 for (Index i=0; i<pop.size(); i++) {
1320 eig.
val[i] -= Complex(f0[i], P * band.
lines[i].lineshape[broadener].G0().at(T, band.
T0));
1327 for (Index i=0; i<pop.size(); i++) {
1328 const Numeric i0 =
LineShape::LocalThermodynamicEquilibrium(band.
lines[i].I0, band.
T0, T, band.
lines[i].F0, band.
lines[i].E0, QT, QT0, 0, 1, 0, 0).
S;
1352 const Vector& temperatures,
1357 const Index K = temperatures.nelem();
1360 const Numeric frenorm = band.
F_mean();
1367 Tensor4 out(4, N, M, K);
1369 #pragma omp parallel for collapse(2) if (!arts_omp_in_parallel())
1370 for (Index m=0; m<M; m++) {
1371 for (Index k=0; k<K; k++) {
1372 const Numeric T = temperatures[k];
1377 for (Index n=0; n<N; n++) {
1378 W(n, n) += band.
lines[sorting[n]].F0 - frenorm;
1386 out(0, joker, m, k) = eig.str.real();
1387 out(1, joker, m, k) = eig.str.imag();
1388 out(2, joker, m, k) = eig.val.real();
1389 out(3, joker, m, k) = eig.val.imag();
1412 const Vector& temperatures,
1417 const Index K = temperatures.nelem();
1420 const Numeric frenorm = band.
F_mean();
1426 Tensor4 out(4, N, M, K);
1428 #pragma omp parallel for collapse(2) if (!arts_omp_in_parallel())
1429 for (Index m=0; m<M; m++) {
1430 for (Index k=0; k<K; k++) {
1431 const Numeric T = temperatures[k];
1435 for (Index n=0; n<N; n++) {
1436 W(n, n) += band.
lines[n].F0 - frenorm;
1442 out(0, joker, m, k) =
RosenkranzG(dip, W.imag(), band);
1443 out(1, joker, m, k) =
RosenkranzY(dip, W.imag(), band);
1444 out(2, joker, m, k) =
RosenkranzDV(dip, W.imag(), band);
1445 out(3, joker, m, k) = 0;
1453 const Vector& temperatures,
1458 const bool rosenkranz_adaptation,
1464 "The temperature list [",
1467 "must be fully sorted from low to high")
1470 ord < 1 or ord > 3,
"Order not in list [1, 2, 3], is: ", ord)
1474 rosenkranz_adaptation not_eq 0
1481 "Bad eigenvalue adaptation for band: ",
1484 out3 <<
"Bad eigenvalue adaptation for band: " << band.
quantumidentity
1488 band.
population = Absorption::PopulationType::LTE;
1489 for (
auto& line : band.
lines) {
1490 for (
auto& sm : line.lineshape.Data()) {
1501 const Vector& temperatures,
1503 const Vector& pressures) {
1506 const Index K = temperatures.nelem();
1507 const Index L = pressures.size();
1509 Tensor5 out(4, N, M, K, L);
1510 for (Index l=0; l<L; l++) {
1518 for (Index i = 0; i < rbd.
data.
nelem(); i++) {
1526 for (
auto& x : rbd.
data) is >> x;
1531 return std::distance(
data.begin(),
1532 std::find(
data.cbegin(),
data.cend(), spec));
1536 Species::Species spec)
const {
1537 if (
auto ptr = std::find(
data.cbegin(),
data.cend(), spec);
1538 ptr not_eq
data.cend())
1540 return *std::find(
data.cbegin(),
data.cend(), Species::Species::Bath);
1544 Species::Species spec) {
1545 if (
auto ptr = std::find(
data.begin(),
data.end(), spec);
1546 ptr not_eq
data.end())
1548 return data.emplace_back(spec);
1553 if (
auto ptr = std::find(begin(), end(),
id); ptr not_eq end())
return *ptr;
1554 return emplace_back(
id);
1559 if (
auto ptr = std::find(cbegin(), cend(),
id); ptr not_eq cend()) {
1565 "Available data:\n",
1572 std::for_each(m.cbegin(), m.cend(), [&](
auto& x) { os << x <<
'\n'; });
Index nelem() const ARTS_NOEXCEPT
Input manipulator class for doubles to enable nan and inf parsing.
Helper macros for debugging.
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR(...)
#define ARTS_USER_ERROR_IF(condition,...)
constexpr bool good_enum(EnumType x) noexcept
Checks if the enum number is good.
Numeric wig6(const Rational &a, const Rational &b, const Rational &c, const Rational &d, const Rational &e, const Rational &f) noexcept
Numeric wig3(const Rational &a, const Rational &b, const Rational &c, const Rational &d, const Rational &e, const Rational &f) noexcept
Numeric boltzman_factor(Numeric T, Numeric E0)
Computes exp(- E0/kT)
Numeric boltzman_ratio(const Numeric &T, const Numeric &T0, const Numeric &E0)
Computes exp(E0/c (T - T0) / (T * T0))
Numeric single_partition_function(const Numeric &T, const Species::IsotopeRecord &ir)
Computes the partition function at one temperature.
Numeric fac(const Index n)
fac
Declarations having to do with the four output streams.
EnergyFunction erot_selection(const SpeciesIsotopeRecord &isot)
void relaxation_matrix_offdiagonal(MatrixView W, const AbsorptionLines &band, const ArrayOfIndex &sorting, const SpeciesErrorCorrectedSuddenData &rovib_data, const Numeric T)
Numeric reduced_dipole(const Rational Ju, const Rational Jl, const Rational N)
constexpr Numeric erot(const Rational N, const Rational j=-1)
void relaxation_matrix_offdiagonal(MatrixView W, const AbsorptionLines &band, const ArrayOfIndex &sorting, const SpeciesErrorCorrectedSuddenData &rovib_data, const Numeric T)
EquivalentLines eigenvalue_adaptation_of_relmat(const ComplexMatrix &W, const Vector &pop, const Vector &dip, const AbsorptionLines &band, const Numeric frenorm, const Numeric T, const Numeric P, const Numeric QT, const Numeric QT0, const Index broadener)
Adapts the relaxation matrix eigenvalues to a form where they represent additions towards the three R...
std::ostream & operator<<(std::ostream &os, const EquivalentLines &eqv)
Tensor4 ecs_eigenvalue_approximation(const AbsorptionLines &band, const Vector &temperatures, const ErrorCorrectedSuddenData &ecs_data, const Numeric P)
std::pair< ComplexVector, bool > ecs_absorption_impl(const Numeric T, const Numeric H, const Numeric P, const Numeric this_vmr, const Vector &vmrs, const ErrorCorrectedSuddenData &ecs_data, const Vector &f_grid, const Zeeman::Polarization zeeman_polarization, const AbsorptionLines &band)
void ecs_eigenvalue_adaptation(AbsorptionLines &band, const Vector &temperatures, const ErrorCorrectedSuddenData &ecs_data, const Numeric P0, const Index ord, const bool robust, const bool rosenkranz_adaptation, const Verbosity &verbosity)
std::pair< ArrayOfIndex, PopulationAndDipole > sorted_population_and_dipole(const Numeric T, const AbsorptionLines &band)
EcsReturn ecs_absorption(const Numeric T, const Numeric H, const Numeric P, const Numeric this_vmr, const Vector &vmrs, const ErrorCorrectedSuddenData &ecs_data, const Vector &f_grid, const Zeeman::Polarization zeeman_polarization, const AbsorptionLines &band, const ArrayOfRetrievalQuantity &jacobian_quantities)
ComplexMatrix single_species_ecs_relaxation_matrix(const AbsorptionLines &band, const ArrayOfIndex &sorting, const Numeric T, const Numeric P, const SpeciesErrorCorrectedSuddenData &species_ecs_data, const Index species_pos)
Vector RosenkranzY(const Vector &dip, const ConstMatrixView &W, const AbsorptionLines &band)
Vector RosenkranzG(const Vector &dip, const ConstMatrixView &W, const AbsorptionLines &band)
Index band_eigenvalue_adaptation(AbsorptionLines &band, const Tensor4 &tempdata, const Vector &temperatures, const Numeric P0, const Index ord)
Adapts the band to the temperature data.
Vector RosenkranzDV(const Vector &dip, const ConstMatrixView &W, const AbsorptionLines &band)
Tensor5 ecs_eigenvalue_adaptation_test(const AbsorptionLines &band, const Vector &temperatures, const ErrorCorrectedSuddenData &ecs_data, const Vector &pressures)
std::function< Numeric(const Rational)> EnergyFunction
ComplexMatrix ecs_relaxation_matrix(const Numeric T, const Numeric P, const Vector &vmrs, const ErrorCorrectedSuddenData &ecs_data, const AbsorptionLines &band, const ArrayOfIndex &sorting, const Numeric frenorm)
std::istream & operator>>(std::istream &is, SpeciesErrorCorrectedSuddenData &srbd)
PopulationAndDipole presorted_population_and_dipole(const Numeric T, const ArrayOfIndex &presorting, const AbsorptionLines &band)
Tensor4 rosenkranz_approximation(const AbsorptionLines &band, const Vector &temperatures, const ErrorCorrectedSuddenData &ecs_data, const Numeric P)
Numeric reduced_rovibrational_dipole(Rational Jf, Rational Ji, Rational lf, Rational li, Rational k=Rational(1))
Compute the reduced rovibrational dipole moment.
constexpr Numeric h_bar
Reduced planck constant convenience name [J s].
constexpr Numeric sqrt_ln_2
Square root of natural logarithm of 2.
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric k
Boltzmann constant convenience name [J/K].
constexpr Numeric sqrt_pi
Square root of pi.
constexpr Numeric h
Planck constant convenience name [J s].
constexpr Numeric m_u
Unified atomic mass unit convenience name [kg].
constexpr auto kelvin2joule(auto x) noexcept
Conversion from Kelvin to Joule.
constexpr auto kaycm2joule(auto x) noexcept
Conversion from cm-1 to Joule.
constexpr auto mhz2joule(auto x) noexcept
Conversion from MHz to Joule.
constexpr auto pow3(auto x) noexcept
power of three
constexpr auto pow2(auto x) noexcept
power of two
Polarization
Zeeman polarization selection.
This file contains declerations of functions of physical character.
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Contains recomputed equivalent lines (sorting is unknown)
void sort_by_frequency(Vector &f, const ArrayOfIndex &sorting)
EquivalentLines(Index n=0) noexcept
Rovibrational line mixing data following the ideas of Collisional Effects On Molecular Spectra by J.
ArrayOfSpeciesErrorCorrectedSuddenData data
Data of species data The program is considered ill-formed if data does not contain a Bath-species,...
Index size() const noexcept
Index pos(Species::Species spec) const
const SpeciesErrorCorrectedSuddenData & operator[](Species::Species spec) const
ErrorCorrectedSuddenData & operator[](const QuantumIdentifier &id)
Contains the population distribution and dipole.
ArrayOfIndex sort(const AbsorptionLines &band) noexcept
PopulationAndDipole(const Numeric T, const AbsorptionLines &band)
LineShapeModelParameters collisional_distance
LineShapeModelParameters lambda
LineShapeModelParameters beta
Numeric Omega(const Numeric T, const Numeric T0, const Numeric other_mass, const Numeric energy_x, const Numeric energy_xm2) const
LineShapeModelParameters scaling
Numeric Q(const Rational J, const Numeric T, const Numeric T0, const Numeric energy) const
Index NumBroadeners() const ARTS_NOEXCEPT
Number of broadening species.
Index BroadeningSpeciesPosition(Species::Species spec) const noexcept
Position of species if available or -1 else.
Numeric T0
Reference temperature for all parameters of the lines.
PopulationType population
Line population distribution.
Numeric F_mean(Numeric T=0) const noexcept
Mean frequency by weight of line strength.
bool bathbroadening
Does the line broadening have bath broadening.
Numeric ZeemanSplitting(size_t k, Zeeman::Polarization type, Index i) const ARTS_NOEXCEPT
Returns the splitting of a Zeeman split line.
Array< SingleLine > lines
A list of individual lines.
Index NumLines() const noexcept
Number of lines.
Numeric DopplerConstant(Numeric T) const noexcept
Numeric ZeemanStrength(size_t k, Zeeman::Polarization type, Index i) const ARTS_NOEXCEPT
Returns the strength of a Zeeman split line.
NormalizationType normalization
Line normalization type.
Species::IsotopeRecord Isotopologue() const noexcept
Isotopologue Index.
LineShape::Type lineshapetype
Type of line shape.
LineShape::Output ShapeParameters(size_t k, Numeric T, Numeric P, const Vector &vmrs) const ARTS_NOEXCEPT
Line shape parameters.
Index ZeemanCount(size_t k, Zeeman::Polarization type) const ARTS_NOEXCEPT
Returns the number of Zeeman split lines.
ArrayOfSpecies broadeningspecies
A list of broadening species.
bool selfbroadening
Does the line broadening have self broadening.
Rational max(QuantumNumberType) const
Numeric SpeciesMass() const noexcept
Mass of the molecule.
QuantumIdentifier quantumidentity
Catalog ID.
bool DoVmrDerivative(const QuantumIdentifier &qid) const noexcept
Coefficients and temperature model for SingleSpeciesModel.
Numeric at(Numeric T, Numeric T0) const noexcept
A logical struct for global quantum numbers with species identifiers.
StateMatchType operates so that a check less than a level should be 'better', bar None.
Struct containing all information needed about one isotope.
Numeric wigner6j(const Rational j1, const Rational j2, const Rational j3, const Rational l1, const Rational l2, const Rational l3)
Wigner 6J symbol.
Numeric wigner3j(const Rational j1, const Rational j2, const Rational j3, const Rational m1, const Rational m2, const Rational m3)
Wigner 3J symbol.
std::pair< Rational, Rational > wigner_limits(std::pair< Rational, Rational > a, std::pair< Rational, Rational > b)
Wigner symbol interactions.
constexpr int temp_init_size(Integer... vals) noexcept