5#include <Faddeeva/Faddeeva.hh>
25#define WIGNER3 fw3jja6
39 a.toInt(2),
b.toInt(2),
c.toInt(2),
d.toInt(2), e.toInt(2), f.toInt(2));
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();
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] * invV(i, j);
109 for (
Index i=0; i<
N; i++) {
110 for (
Index j=i+1; j<
N; j++) {
112 std::swap(
val[i],
val[j]);
113 std::swap(
str[i],
str[j]);
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]];
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};
278 const Numeric energy_xm2)
const {
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
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;
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;
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);
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];
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];
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) {
436 sum +=
a *
b *
c *
d * e *
Numeric(2 * L + 1) * Q[L.toIndex()] / Om[L.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") {
497 return [](
const Rational J) ->
Numeric {
return Numeric(J) * std::numeric_limits<Numeric>::signaling_NaN();};
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];
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];
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) {
553 const Numeric QL = rovib_data.
Q(L, T, band.
T0, erot(L));
555 sum +=
a *
b *
c *
Numeric(2 * L + 1) * QL / ECS;
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) {
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");
674 for (
Index k=0; k<
M; k++) {
675 if (vmrs[k] == 0)
continue;
683 for (
Index i=0; i<
N; i++) {
719 auto& [absorption,works] = retval;
725 for (
Index j=0; j<nz; j++) {
734 const Numeric gamd = GD_div_F0 * (eqv.
val[i].
real() + frenorm + dzeeman);
737 const Complex z = (eqv.
val[i] + frenorm + dzeeman - f_grid[iv]) * cte;
739 absorption[iv] += Sz * eqv.
str[i] *
w / gamd;
750 absorption[iv] *= fact * numdens * sq_ln2pi;
753 if (
auto&
a =
real_val(absorption[iv]); not std::isnormal(
a) or
a < 0) {
773 auto [absorption, work] =
ecs_absorption_impl(T, H, P, this_vmr, vmrs, ecs_data, f_grid, zeeman_polarization, band);
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) {
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()) {
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();
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};
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;
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);
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));
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);
1292 const Index broadener) {
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);
1307 if (f0[j] < f0[i]) {
1308 std::swap(sorting[i], sorting[j]);
1309 std::swap(f0[i], f0[j]);
1320 eig.
val[i] -=
Complex(f0[i], P * band.
lines[i].lineshape[broadener].G0().at(T, band.
T0));
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,
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,
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;
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) {
1510 for (
Index l=0; l<L; l++) {
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
MatrixView imag()
Get a view of the imaginary parts of the matrix.
ComplexMatrix & inv(const lapack_help::Inverse< Complex > &help=lapack_help::Inverse< Complex >{0})
VectorView real()
Get a view of the real part of the vector.
Index nelem() const noexcept
A constant view of a Matrix.
Index nelem() const noexcept
Returns the number of elements.
Index size() const 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.
void diagonalize(MatrixView P, VectorView WR, VectorView WI, ConstMatrixView A)
Matrix Diagonalization.
Linear algebra functions.
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.
bool is_sorted(ConstVectorView x)
Checks if a vector is sorted in ascending order.
Numeric fac(const Index n)
fac
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
NUMERIC Numeric
The type to use for all floating point numbers.
INDEX Index
The type to use for all integer numbers and indices.
void swap(ComplexVector &v1, ComplexVector &v2) noexcept
Swaps two objects.
Numeric & real_val(Complex &c) noexcept
Return a reference to the real value of c.
std::complex< Numeric > Complex
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.
auto mat(matpack::matrix auto &&x)
Map the input to a non-owning const-correct Eigen Map representing a row matrix.
This file contains declerations of functions of physical character.
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Contains the rational class definition.
Numeric pow(const Rational base, Numeric exp)
Power of.
Numeric sqrt(const Rational r)
Square root.
constexpr bool iseven(const Rational r) noexcept
Returns true if even integer.
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.
Implements rational numbers to work with other ARTS types.
constexpr Index toIndex(int n=1) const noexcept
Converts the value to index by n-scaled division.
Struct containing all information needed about one isotope.
String FullName() const noexcept
Species spec
Species type as defined in species.h.
std::string_view isotname
A custom name that is unique for this Species type.
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