ARTS 2.5.4 (git: bcd8c674)
linemixing.cc
Go to the documentation of this file.
1#include <iomanip>
2#include <iostream>
3#include <numeric>
4
5#include <Faddeeva/Faddeeva.hh>
6
7#include "constants.h"
8#include "debug.h"
9#include "lin_alg.h"
10#include "linemixing.h"
11#include "lineshape.h"
12#include "matpack.h"
13#include "matpack_complex.h"
14#include "messages.h"
15#include "minimize.h"
16#include "physics_funcs.h"
17#include "rational.h"
18#include "species.h"
19
22 const Vector& pop,
23 const Vector& dip) noexcept :
24 val(pop.nelem(), 0), str(pop.nelem(), 0) {
25 /* FIXME: (Added 2021-01-19; Richard Larsson)
26 *
27 * This function cannot easily be used with partial derivatives
28 *
29 * Doing so would allow the entire ECS approach to be analytical
30 * in its derivatives.
31 *
32 * The problem in short:
33 * There is an Eigenvalue decomposition happening (W = V diag(e) V^-1)
34 *
35 * We use V and e in a strange way. We do not know how to get the
36 * partial derivatives of these two variables
37 *
38 * The question:
39 * How do we compute dV and de? We can get dW somewhat easily.
40 */
41
42 const Index n = pop.nelem();
43
44 // Compute values
45 ComplexMatrix V(n, n);
46
47 // Main computations
48 diagonalize(V, val, W);
49
50 // Do the matrix forward multiplication
51 for (Index i=0; i<n; i++) {
52 for (Index j=0; j<n; j++) {
53 str[i] += dip[j] * V(j, i);
54 }
55 }
56
57 // Do the matrix backward multiplication
58 auto& invV=V.inv();
59 for (Index i=0; i<n; i++) {
60 Complex z(0, 0);
61 for (Index j=0; j<n; j++) {
62 z += pop[j] * dip[j] * invV(i, j);
63 }
64 str[i] *= z;
65 }
66}
67
68
70 const Index N = val.nelem();
71
72 // First sort by frequency and assume this sorting
73 // is a shifted version of the line centers
74 for (Index i=0; i<N; i++) {
75 for (Index j=i+1; j<N; j++) {
76 if (val.real()[j] < val.real()[i]) {
77 std::swap(val[i], val[j]);
78 std::swap(str[i], str[j]);
79 }
80 }
81 }
82
83 const Vector fc = f;
84 const ComplexVector strc = str;
85 const ComplexVector valc = val;
86 for (Index i=0; i<N; i++) {
87 f[i] = fc[sorting[i]];
88 str[i] = strc[sorting[i]];
89 val[i] = valc[sorting[i]];
90 }
91}
92
93std::ostream& operator<<(std::ostream& os, const EquivalentLines& eqv) {
94 const Index n = eqv.str.nelem();
95 for (Index i=0; i<n; i++) {
96 if (i) os << '\n';
97 os << eqv.str[i] << ' ' << eqv.val[i];
98 }
99 return os;
100}
101
102std::ostream& operator<<(std::ostream& os, const SpeciesErrorCorrectedSuddenData& srbd) {
103 os << Species::toShortName(srbd.spec);
104 os << ' ' << srbd.scaling;
105 os << ' ' << srbd.beta;
106 os << ' ' << srbd.lambda;
107 os << ' ' << srbd.collisional_distance;
108 os << ' ' << srbd.mass;
109 return os;
110}
111
112std::istream& operator>>(std::istream& is, SpeciesErrorCorrectedSuddenData& srbd) {
113 std::string spec_name;
114 is >> spec_name;
115 is >> srbd.scaling;
116 is >> srbd.beta;
117 is >> srbd.lambda;
118 is >> srbd.collisional_distance;
119 is >> double_imanip() >> srbd.mass;
120 srbd.spec = Species::fromShortName(spec_name);
121 ARTS_USER_ERROR_IF(not good_enum(srbd.spec), "Cannot recognize species: ", spec_name)
122 return is;
123}
124
125namespace Makarov2020etal {
133Numeric reduced_dipole(const Rational Ju, const Rational Jl, const Rational N) {
134 return (iseven(Jl + N) ? 1 : -1) * sqrt(6 * (2*Jl + 1) * (2*Ju + 1)) * wigner6j(1, 1, 1, Jl, Ju, N);
135};
136} // namespace Makarov2020etal
137
138
140 const AbsorptionLines& band) noexcept :
141 pop(band.NumLines()), dip(band.NumLines()) {
142 const Index n = band.NumLines();
143
144 const Numeric QT = single_partition_function(T, band.Isotopologue());
145 const Numeric QT0 = single_partition_function(band.T0, band.Isotopologue());
146 const Numeric ratiopart = QT0 / QT;
147
148 for (Index i=0; i<n; i++) {
149 const Numeric pop0 = (band.lines[i].gupp / QT0) * boltzman_factor(band.T0, band.lines[i].E0);
150 pop[i] = pop0 * ratiopart * boltzman_ratio(T, band.T0, band.lines[i].E0);
151 dip[i] = std::sqrt(- band.lines[i].I0/(pop0 * band.lines[i].F0 * std::expm1(- (Constant::h * band.lines[i].F0)) / (Constant::k * band.T0)));
152
153 // Adjust the sign depending on type
154 if (band.population == Absorption::PopulationType::ByMakarovFullRelmat) {
155 auto& J = band.lines[i].localquanta.val[QuantumNumberType::J];
156 auto& N = band.lines[i].localquanta.val[QuantumNumberType::N];
157 dip[i] *= std::signbit(Makarov2020etal::reduced_dipole(J.upp(), J.low(), N.upp())) ? - 1 : 1;
158 } else if (band.population == Absorption::PopulationType::ByRovibLinearDipoleLineMixing) {
159 // pass
160 }
161 }
162}
163
164
166 const Index N = pop.nelem();
167
168 // List that starts as [0, 1, ... N-2, N-1]
169 ArrayOfIndex out(N);
170 std::iota(out.begin(), out.end(), 0);
171
172 // Strength
173 Vector s(N);
174 for (Index i=0; i<N; i++) {
175 s[i] = band.lines[i].F0 * pop[i] * Constant::pow2(dip[i]);
176 }
177
178 // Sorter
179 for (Index i=0; i<N; i++) {
180 for (Index j=i+1; j<N; j++) {
181 if (s[j] > s[i]) {
182 std::swap(out[i], out[j]);
183 std::swap(dip[i], dip[j]);
184 std::swap(pop[i], pop[j]);
185 std::swap(s[i], s[j]);
186 }
187 }
188 }
189
190 return out;
191}
192
193
194void PopulationAndDipole::sort(const ArrayOfIndex& presorting) noexcept {
195 const Index N = presorting.size();
196 Vector dipcopy(dip), popcopy(pop);
197 for (Index i=0; i<N; i++) {
198 dip[i] = dipcopy[presorting[i]];
199 pop[i] = popcopy[presorting[i]];
200 }
201}
202
203
210std::pair<ArrayOfIndex, PopulationAndDipole> sorted_population_and_dipole(const Numeric T,
211 const AbsorptionLines& band) {
212 PopulationAndDipole tp(T, band);
213 return {tp.sort(band), tp};
214}
215
216
225 const ArrayOfIndex& presorting,
226 const AbsorptionLines& band) {
227 PopulationAndDipole tp(T, band);
228 tp.sort(presorting);
229 return tp;
230}
231
233 const Numeric T,
234 const Numeric T0,
235 const Numeric energy) const noexcept {
236 return std::exp(- beta.at(T, T0) * energy / (Constant::k * T)) * scaling.at(T, T0) / pow(J * (J+1), lambda.at(T, T0));
237}
238
240 const Numeric T0,
241 const Numeric other_mass,
242 const Numeric energy_x,
243 const Numeric energy_xm2) const noexcept
244{
245 using Constant::h;
246 using Constant::k;
247 using Constant::pi;
248 using Constant::m_u;
249 using Constant::h_bar;
250 using Constant::pow2;
251
252 // Constants for the expression
253 constexpr Numeric fac = 8 * k / (m_u * pi);
254
255 const Numeric wnnm2 = (energy_x - energy_xm2) / h_bar;
256
257 const Numeric inv_eff_mass = 1 / mass + 1 / other_mass;
258 const Numeric v_bar_pow2 = fac*T*inv_eff_mass;
259 const Numeric tauc_pow2 = pow2(collisional_distance.at(T, T0)) / v_bar_pow2;
260
261 return 1.0 / pow2(1 + pow2(wnnm2) * tauc_pow2 / 24.0);
262}
263
264namespace Makarov2020etal {
274template<bool rescale_pure_rotational=true> constexpr
275Numeric erot(const Rational N, const Rational j=-1) {
276 const Rational J = j<0 ? N : j;
277
278 if constexpr (rescale_pure_rotational) {
279 return erot<false>(N, J) - erot<false>(1, 0);
280 } else {
282 using Constant::pow2;
283 using Constant::pow3;
284
285 constexpr Numeric B0=43100.4425e0;
286 constexpr Numeric D0=.145123e0;
287 constexpr Numeric H0=3.8e-08;
288 constexpr Numeric xl0=59501.3435e0;
289 constexpr Numeric xg0=-252.58633e0;
290 constexpr Numeric xl1=0.058369e0;
291 constexpr Numeric xl2=2.899e-07;
292 constexpr Numeric xg1=-2.4344e-04;
293 constexpr Numeric xg2=-1.45e-09;
294
295 const auto XN = Numeric(N);
296 const Numeric XX = XN * (XN + 1);
297 const Numeric xlambda=xl0+xl1*XX+xl2*pow2(XX);
298 const Numeric xgama=xg0+xg1*XX+xg2*pow2(XX);
299 const Numeric C1=B0*XX-D0*pow2(XX)+H0*pow3(XX);
300
301 if (J < N) {
302 if (N == 1) // erot<false>(1, 0)
303 return mhz2joule(C1 - (xlambda+B0*(2.*XN-1.)+xgama*XN));
304 return mhz2joule(C1 - (xlambda+B0*(2.*XN-1.)+xgama*XN) + std::sqrt(pow2(B0*(2.*XN-1.))+pow2(xlambda)-2.*B0*xlambda));
305 }
306 if (J > N)
307 return mhz2joule(C1 - (xlambda-B0*(2.*XN+3.)-xgama*(XN+1.)) - std::sqrt(pow2(B0*(2.*XN+3.))+pow2(xlambda)-2.*B0*xlambda));
308 return mhz2joule(C1);
309 }
310}
311
323 const AbsorptionLines& band,
324 const ArrayOfIndex& sorting,
325 const SpeciesErrorCorrectedSuddenData& rovib_data,
326 const Numeric T) ARTS_NOEXCEPT {
328
329 const auto bk = [](const Rational& r) -> Numeric {return sqrt(2*r + 1);};
330
331 const Index n = band.NumLines();
332 if (n == 0) return;
333
334 auto& S = band.quantumidentity.val[QuantumNumberType::S];
335 const Rational Si = S.upp();
336 const Rational Sf = S.low();
337 ARTS_ASSERT(Si.isDefined() and Sf.isDefined(), "Need S for band selection")
338
339 Vector dipr(n);
340 for (Index i=0; i<n; i++) {
341 auto& J = band.lines[sorting[i]].localquanta.val[QuantumNumberType::J];
342 auto& N = band.lines[sorting[i]].localquanta.val[QuantumNumberType::N];
343
344 dipr[i] = reduced_dipole(J.upp(), J.low(), N.upp());
345 }
346
347 for (Index i = 0; i < n; i++) {
348 auto& J = band.lines[sorting[i]].localquanta.val[QuantumNumberType::J];
349 auto& N = band.lines[sorting[i]].localquanta.val[QuantumNumberType::N];
350
351 const Rational Ji = J.upp();
352 const Rational Jf = J.low();
353 const Rational Ni = N.upp();
354 const Rational Nf = N.low();
355
356 for (Index j = 0; j < n; j++) {
357 if (i == j) continue;
358
359 auto& J_p = band.lines[sorting[j]].localquanta.val[QuantumNumberType::J];
360 auto& N_p = band.lines[sorting[j]].localquanta.val[QuantumNumberType::N];
361
362 const Rational Ji_p = J_p.upp();
363 const Rational Jf_p = J_p.low();
364 const Rational Ni_p = N_p.upp();
365 const Rational Nf_p = N_p.low();
366
367 if (Jf_p > Jf) continue;
368
369 // Tran etal 2006 symbol with modifications:
370 // 1) [Ji] * [Ji_p] instead of [Ji_p] ^ 2 in partial accordance with Makarov etal 2013
371 Numeric sum=0;
372 const Numeric scl = (iseven(Ji_p + Ji + 1) ? 1 : -1) * bk(Ni) * bk(Nf) * bk(Nf_p) * bk(Ni_p) * bk(Jf) * bk(Jf_p) * bk(Ji) * bk(Ji_p);
373 const auto [L0, L1] = wigner_limits(wigner3j_limits<3>(Ni_p, Ni), {Rational(2), std::numeric_limits<Index>::max()});
374 for (Rational L=L0; L<=L1; L+=2) {
375 const Numeric a = wigner3j(Ni_p, Ni, L, 0, 0, 0);
376 const Numeric b = wigner3j(Nf_p, Nf, L, 0, 0, 0);
377 const Numeric c = wigner6j(L, Ji, Ji_p, Si, Ni_p, Ni);
378 const Numeric d = wigner6j(L, Jf, Jf_p, Sf, Nf_p, Nf);
379 const Numeric e = wigner6j(L, Ji, Ji_p, 1, Jf_p, Jf);
380 sum += a * b * c * d * e * Numeric(2*L + 1) * rovib_data.Q(L, T, band.T0, erot(L)) / rovib_data.Omega(T, band.T0, band.SpeciesMass(), erot(L), erot(L-2));
381 }
382 sum *= scl * rovib_data.Omega(T, band.T0, band.SpeciesMass(), erot(Ni), erot(Ni - 2));
383
384 // Add to W and rescale to upwards element by the populations
385 W(i, j) = sum;
386 W(j, i) = sum * std::exp((erot(Nf_p, Jf_p) - erot(Nf, Jf)) / kelvin2joule(T));
387 }
388 }
389
390 // Sum rule correction
391 for (Index i=0; i<n; i++) {
392 Numeric sumlw = 0.0;
393 Numeric sumup = 0.0;
394
395 for (Index j=0; j<n; j++) {
396 if (j > i) {
397 sumlw += dipr[j] * W(j, i);
398 } else {
399 sumup += dipr[j] * W(j, i);
400 }
401 }
402
403 const Rational Ji = band.lines[sorting[i]].localquanta.val[QuantumNumberType::J].low();
404 const Rational Ni = band.lines[sorting[i]].localquanta.val[QuantumNumberType::N].low();
405 for (Index j = i + 1; j < n; j++) {
406 const Rational Jj = band.lines[sorting[j]].localquanta.val[QuantumNumberType::J].low();
407 const Rational Nj = band.lines[sorting[j]].localquanta.val[QuantumNumberType::N].low();
408 if (sumlw == 0) {
409 W(j, i) = 0.0;
410 W(i, j) = 0.0;
411 } else {
412 W(j, i) *= - sumup / sumlw;
413 W(i, j) = W(j, i) * std::exp((erot(Ni, Ji) - erot(Nj, Jj)) / kelvin2joule(T)); // This gives LTE
414 }
415 }
416 }
417}
418} // namespace Makarov2020etal
419
420
429namespace LinearRovibErrorCorrectedSudden {
430
432{
433 if (isot.spec == Species::Species::CarbonDioxide and isot.isotname == "626") {
434 return [](const Rational J) -> Numeric {return Conversion::kaycm2joule(0.39021) * Numeric(J * (J + 1));};
435 }
436
437 ARTS_USER_ERROR(isot.FullName(), " has no rotational energies in ARTS")
438 return [](const Rational J) -> Numeric {return Numeric(J) * std::numeric_limits<Numeric>::signaling_NaN();};
439}
440
442 const AbsorptionLines& band,
443 const ArrayOfIndex& sorting,
444 const SpeciesErrorCorrectedSuddenData& rovib_data,
445 const Numeric T) ARTS_NOEXCEPT
446{
448
449 const Index n = band.NumLines();
450 if (not n) return;
451
452 // These are constant for a band
453 auto& l2 = band.quantumidentity.val[QuantumNumberType::l2];
454 Rational li = l2.upp();
455 Rational lf = l2.low();
456
457 const bool swap_order = li > lf;
458 if (swap_order) swap(li, lf);
459 const int sgn = iseven(li + lf + 1) ? -1 : 1;
460 if (abs(li - lf) > 1) return;
461
462 const EnergyFunction erot = erot_selection(band.Isotopologue());
463
464 Vector dipr(n);
465 for (Index i=0; i<n; i++) {
466 auto& J = band.lines[i].localquanta.val[QuantumNumberType::J];
467 dipr[i] = Absorption::reduced_rovibrational_dipole(J.upp(), J.low(), lf, li);
468 }
469
470 for (Index i=0; i<n; i++) {
471 auto& J = band.lines[sorting[i]].localquanta.val[QuantumNumberType::J];
472 Rational Ji = J.upp();
473 Rational Jf = J.low();
474 if (swap_order) swap(Ji, Jf);
475
476 for (Index j=0; j<n; j++) {
477 if(i == j) continue;
478 auto& J_p = band.lines[sorting[j]].localquanta.val[QuantumNumberType::J];
479 Rational Ji_p = J_p.upp();
480 Rational Jf_p = J_p.low();
481 if (swap_order) swap(Ji_p, Jf_p);
482
483 // Select upper quantum number
484 if (Jf_p > Jf) continue;
485
486 Index L = std::max(std::abs((Ji-Ji_p).toIndex()), std::abs((Jf-Jf_p).toIndex()));
487 L += L % 2;
488 const Index Lf = std::min((Ji+Ji_p).toIndex(), (Jf+Jf_p).toIndex());
489
490 Numeric sum=0;
491 for (; L<=Lf; L+=2) {
492 const Numeric a = wigner3j(Ji_p, L, Ji, li, 0, -li);
493 const Numeric b = wigner3j(Jf_p, L, Jf, lf, 0, -lf);
494 const Numeric c = wigner6j(Ji, Jf, 1, Jf_p, Ji_p, L);
495 const Numeric QL = rovib_data.Q(L, T, band.T0, erot(L));
496 const Numeric ECS = rovib_data.Omega(T, band.T0, band.SpeciesMass(), erot(L), erot(L-2));
497 sum += a * b * c * Numeric(2 * L + 1) * QL / ECS;
498 }
499 const Numeric ECS = rovib_data.Omega(T, band.T0, band.SpeciesMass(), erot(Ji), erot(Ji-2));
500 const Numeric scl = sgn * ECS * Numeric(2 * Ji_p + 1) * sqrt((2 * Jf + 1) * (2 * Jf_p + 1));
501 sum *= scl;
502
503 // Add to W and rescale to upwards element by the populations
504 W(j, i) = sum;
505 W(i, j) = sum * std::exp((erot(Jf_p) - erot(Jf)) / kelvin2joule(T));
506 }
507 }
508
509 // Undocumented negative absolute sign
510 for (Index i=0; i<n; i++)
511 for (Index j=0; j<n; j++)
512 if (j not_eq i and W(i, j) > 0)
513 W(i, j) *= -1;
514
515 // Sum rule correction
516 for (Index i=0; i<n; i++) {
517 Numeric sumlw = 0.0;
518 Numeric sumup = 0.0;
519
520 for (Index j=0; j<n; j++) {
521 if (j > i) {
522 sumlw += std::abs(dipr[j]) * W(j, i); // Undocumented abs-sign
523 } else {
524 sumup += std::abs(dipr[j]) * W(j, i); // Undocumented abs-sign
525 }
526 }
527
528 const Rational Ji = band.lines[sorting[i]].localquanta.val[QuantumNumberType::J].low();
529 for (Index j=i+1; j<n; j++) {
530 const Rational Jj = band.lines[sorting[j]].localquanta.val[QuantumNumberType::J].low();
531 if (sumlw == 0) {
532 W(j, i) = 0.0;
533 W(i, j) = 0.0;
534 } else {
535 W(j, i) *= - sumup / sumlw;
536 W(i, j) = W(j, i) * std::exp((erot(Ji) - erot(Jj)) / kelvin2joule(T)); // This gives LTE
537 }
538 }
539 }
540}
541} // namespace LinearRovibErrorCorrectedSudden
542
556 const ArrayOfIndex& sorting,
557 const Numeric T,
558 const Numeric P,
559 const SpeciesErrorCorrectedSuddenData& species_ecs_data,
560 const Index species_pos) {
561 const Index N = band.NumLines();
562
563 // Allocate the matrix
564 ComplexMatrix W(N, N, 0);
565
566 // Fill diagonal keeping track of the reshuffle (sorting)
567 for (Index I=0; I<N; I++) {
568 const Index i = sorting[I];
569 auto shape = band.ShapeParameters(i, T, P, species_pos);
570 W(I, I) = Complex(shape.D0, shape.G0); // nb. This must be fixed if SDVP ???
571 }
572
573 // Set the off-diagonal part of the matrix for this broadener
574 switch (band.population) {
575 case PopulationType::ByMakarovFullRelmat:
576 Makarov2020etal::relaxation_matrix_offdiagonal(W.imag(), band, sorting, species_ecs_data, T);
577 break;
578 case PopulationType::ByRovibLinearDipoleLineMixing: {
579 LinearRovibErrorCorrectedSudden::relaxation_matrix_offdiagonal(W.imag(), band, sorting, species_ecs_data, T);
580 } break;
581 default:
582 ARTS_ASSERT(false, "Bad type, we don't support band population type: ", band.population,
583 "\nin this code. It must either be added or computations aborted earlier");
584 break;
585 }
586
587 return W;
588}
589
590
603 const Numeric P,
604 const Vector& vmrs,
605 const ErrorCorrectedSuddenData& ecs_data,
606 const AbsorptionLines& band,
607 const ArrayOfIndex& sorting,
608 const Numeric frenorm) {
609 const Index N = band.NumLines();
610 const Index M = vmrs.nelem();
611
612 // Create output
613 ComplexMatrix W(N, N, 0);
614
615 // Loop over all the broadeners
616 for (Index k=0; k<M; k++) {
617 if (vmrs[k] == 0) continue;
618
619 // Sum up all atmospheric components
620 MapToEigen(W).noalias() += vmrs[k] * MapToEigen(
621 single_species_ecs_relaxation_matrix(band, sorting, T, P, ecs_data[band.broadeningspecies[k]], k));
622 }
623
624 // Deal with line frequency and its re-normalization
625 for (Index i=0; i<N; i++) {
626 real_val(W(i, i)) += band.lines[sorting[i]].F0 - frenorm;
627 }
628
629 return W;
630}
631
632
633std::pair<ComplexVector, bool> ecs_absorption_impl(const Numeric T,
634 const Numeric H,
635 const Numeric P,
636 const Numeric this_vmr,
637 const Vector& vmrs,
638 const ErrorCorrectedSuddenData& ecs_data,
639 const Vector& f_grid,
640 const Zeeman::Polarization zeeman_polarization,
641 const AbsorptionLines& band) {
642 constexpr Numeric sq_ln2pi = Constant::sqrt_ln_2 / Constant::sqrt_pi;
643
644 // Weighted center of the band
645 const Numeric frenorm = band.F_mean();
646
647 // Band Doppler broadening constant
648 const Numeric GD_div_F0 = band.DopplerConstant(T);
649
650 // Sorted population
651 const auto [sorting, tp] = sorted_population_and_dipole(T, band);
652
653 // Relaxation matrix
654 const ComplexMatrix W = ecs_relaxation_matrix(T, P, vmrs, ecs_data, band, sorting, frenorm);
655
656 // Equivalent lines computations
657 const EquivalentLines eqv(W, tp.pop, tp.dip);
658
659 // Return the absorption and if this works
660 std::pair<ComplexVector, bool> retval{ComplexVector(f_grid.nelem(), 0), false};
661 auto& [absorption,works] = retval;
662
663 // Add the lines
664 for (Index i=0; i<band.NumLines(); i++) {
665 // Zeeman lines if necessary
666 const Index nz = band.ZeemanCount(i, zeeman_polarization);
667 for (Index j=0; j<nz; j++) {
668 const Numeric Sz = band.ZeemanStrength(i, zeeman_polarization, j);
669 const Numeric dzeeman = H * band.ZeemanSplitting(i, zeeman_polarization, j);
670
671 if (band.lineshapetype == LineShape::Type::LP) {
672 for (Index iv=0; iv<f_grid.nelem(); iv++) {
673 absorption[iv] -= 1i * Sz * ((eqv.str[i] / (f_grid[iv] - frenorm - dzeeman - eqv.val[i]))) / (Constant::sqrt_ln_2 * Constant::sqrt_pi);
674 }
675 } else /*if (band.LineShapeType() == LineShape::Type::VP)*/ {
676 const Numeric gamd = GD_div_F0 * (eqv.val[i].real() + frenorm + dzeeman);
677 const Numeric cte = Constant::sqrt_ln_2 / gamd;
678 for (Index iv=0; iv<f_grid.nelem(); iv++) {
679 const Complex z = (eqv.val[i] + frenorm + dzeeman - f_grid[iv]) * cte;
680 const Complex w = Faddeeva::w(z);
681 absorption[iv] += Sz * eqv.str[i] * w / gamd;
682 }
683 }
684 }
685 }
686
687 // Adjust by frequency and number density
688 const Numeric numdens = this_vmr * number_density(P, T);
689 for (Index iv=0; iv<f_grid.nelem(); iv++) {
690 const Numeric f = f_grid[iv];
691 const Numeric fact = - f * std::expm1(- (Constant::h * f) / (Constant::k * T));
692 absorption[iv] *= fact * numdens * sq_ln2pi;
693
694 // correct for too low value
695 if (auto& a = real_val(absorption[iv]); not std::isnormal(a) or a < 0) {
696 absorption[iv] = 0;
697 works = false;
698 }
699 }
700
701 return retval;
702}
703
704
706 const Numeric H,
707 const Numeric P,
708 const Numeric this_vmr,
709 const Vector& vmrs,
710 const ErrorCorrectedSuddenData& ecs_data,
711 const Vector& f_grid,
712 const Zeeman::Polarization zeeman_polarization,
713 const AbsorptionLines& band,
714 const ArrayOfRetrievalQuantity& jacobian_quantities) {
715 auto [absorption, work] = ecs_absorption_impl(T, H, P, this_vmr, vmrs, ecs_data, f_grid, zeeman_polarization, band);
716
717 // Start as original, so remove new and divide with the negative to get forward derivative
718 ArrayOfComplexVector jacobian(jacobian_quantities.nelem(), absorption);
719
720 for (Index i=0; i<jacobian_quantities.nelem(); i++) {
721 auto& vec = jacobian[i];
722 auto& target = jacobian_quantities[i].Target();
723
724 if (target == Jacobian::Atm::Temperature) {
725 const Numeric dT = target.perturbation;
726 const auto [dabs, dwork] = ecs_absorption_impl(T+dT, H, P, this_vmr, vmrs, ecs_data, f_grid, zeeman_polarization, band);
727 vec -= dabs;
728 vec /= -dT;
729 work &= dwork;
730 } else if (target.isMagnetic()) {
731 const Numeric dH = target.perturbation;
732 const auto [dabs, dwork] = ecs_absorption_impl(T, H+dH, P, this_vmr, vmrs, ecs_data, f_grid, zeeman_polarization, band);
733 vec -= dabs;
734 vec /= -dH;
735 work &= dwork;
736 } else if (target.isWind()) {
737 const Numeric df = target.perturbation;
738 Vector f_grid_copy = f_grid;
739 f_grid_copy += df;
740 const auto [dabs, dwork] = ecs_absorption_impl(T, H, P, this_vmr, vmrs, ecs_data, f_grid_copy, zeeman_polarization, band);
741 vec -= dabs;
742 vec /= df;
743 work &= dwork;
744 } else if (target == Jacobian::Line::VMR) {
745 if (band.DoVmrDerivative(target.qid)) {
746 Vector vmrs_copy = vmrs;
747 Numeric this_vmr_copy = this_vmr;
748 const Numeric dvmr = target.perturbation;
749
750 // Alter the VMRs for self
751 if (band.Isotopologue() == target.qid.Isotopologue()) {
752 this_vmr_copy += dvmr;
753 if (band.selfbroadening) vmrs_copy[0] += dvmr; // First value is self if band has self broadener
754 } else {
755 for (Index j=band.selfbroadening; j<band.broadeningspecies.nelem()-band.bathbroadening; j++) {
756 if (band.broadeningspecies[j] == target.qid.Species()) {
757 vmrs_copy[j] += dvmr;
758 }
759 }
760 }
761
762 // Computations
763 const auto [dabs, dwork] = ecs_absorption_impl(T, H, P, this_vmr_copy, vmrs_copy, ecs_data, f_grid, zeeman_polarization, band);
764 vec -= dabs;
765 vec /= -dvmr;
766 work &= dwork;
767 }
768 } else if (target.needQuantumIdentity()) {
769 Numeric d=1e-6;
770
771 for (Index iline=0; iline<band.NumLines(); iline++) {
772 if (Quantum::Number::StateMatch(target.qid, band.lines[iline].localquanta, band.quantumidentity) == Quantum::Number::StateMatchType::Full) {
773 AbsorptionLines band_copy = band;
774
775 const Index pos = band.BroadeningSpeciesPosition(target.qid.Species());
776 switch (target.line) {
777 case Jacobian::Line::ShapeG0X0:
778 d *= band.lines[iline].lineshape[pos].G0().X0;
779 band_copy.lines[iline].lineshape[pos].G0().X0 += d;
780 break;
781 case Jacobian::Line::ShapeG0X1:
782 d *= band.lines[iline].lineshape[pos].G0().X1;
783 band_copy.lines[iline].lineshape[pos].G0().X1 += d;
784 break;
785 case Jacobian::Line::ShapeG0X2:
786 d *= band.lines[iline].lineshape[pos].G0().X2;
787 band_copy.lines[iline].lineshape[pos].G0().X2 += d;
788 break;
789 case Jacobian::Line::ShapeG0X3:
790 d *= band.lines[iline].lineshape[pos].G0().X3;
791 band_copy.lines[iline].lineshape[pos].G0().X3 += d;
792 break;
793 case Jacobian::Line::ShapeD0X0:
794 d *= band.lines[iline].lineshape[pos].D0().X0;
795 band_copy.lines[iline].lineshape[pos].D0().X0 += d;
796 break;
797 case Jacobian::Line::ShapeD0X1:
798 d *= band.lines[iline].lineshape[pos].D0().X1;
799 band_copy.lines[iline].lineshape[pos].D0().X1 += d;
800 break;
801 case Jacobian::Line::ShapeD0X2:
802 d *= band.lines[iline].lineshape[pos].D0().X2;
803 band_copy.lines[iline].lineshape[pos].D0().X2 += d;
804 break;
805 case Jacobian::Line::ShapeD0X3:
806 d *= band.lines[iline].lineshape[pos].D0().X3;
807 band_copy.lines[iline].lineshape[pos].D0().X3 += d;
808 break;
809 case Jacobian::Line::ShapeG2X0:
810 d *= band.lines[iline].lineshape[pos].G2().X0;
811 band_copy.lines[iline].lineshape[pos].G2().X0 += d;
812 break;
813 case Jacobian::Line::ShapeG2X1:
814 d *= band.lines[iline].lineshape[pos].G2().X1;
815 band_copy.lines[iline].lineshape[pos].G2().X1 += d;
816 break;
817 case Jacobian::Line::ShapeG2X2:
818 d *= band.lines[iline].lineshape[pos].G2().X2;
819 band_copy.lines[iline].lineshape[pos].G2().X2 += d;
820 break;
821 case Jacobian::Line::ShapeG2X3:
822 d *= band.lines[iline].lineshape[pos].G2().X3;
823 band_copy.lines[iline].lineshape[pos].G2().X3 += d;
824 break;
825 case Jacobian::Line::ShapeD2X0:
826 d *= band.lines[iline].lineshape[pos].D2().X0;
827 band_copy.lines[iline].lineshape[pos].D2().X0 += d;
828 break;
829 case Jacobian::Line::ShapeD2X1:
830 d *= band.lines[iline].lineshape[pos].D2().X1;
831 band_copy.lines[iline].lineshape[pos].D2().X1 += d;
832 break;
833 case Jacobian::Line::ShapeD2X2:
834 d *= band.lines[iline].lineshape[pos].D2().X2;
835 band_copy.lines[iline].lineshape[pos].D2().X2 += d;
836 break;
837 case Jacobian::Line::ShapeD2X3:
838 d *= band.lines[iline].lineshape[pos].D2().X3;
839 band_copy.lines[iline].lineshape[pos].D2().X3 += d;
840 break;
841 case Jacobian::Line::ShapeFVCX0:
842 d *= band.lines[iline].lineshape[pos].FVC().X0;
843 band_copy.lines[iline].lineshape[pos].FVC().X0 += d;
844 break;
845 case Jacobian::Line::ShapeFVCX1:
846 d *= band.lines[iline].lineshape[pos].FVC().X1;
847 band_copy.lines[iline].lineshape[pos].FVC().X1 += d;
848 break;
849 case Jacobian::Line::ShapeFVCX2:
850 d *= band.lines[iline].lineshape[pos].FVC().X2;
851 band_copy.lines[iline].lineshape[pos].FVC().X2 += d;
852 break;
853 case Jacobian::Line::ShapeFVCX3:
854 d *= band.lines[iline].lineshape[pos].FVC().X3;
855 band_copy.lines[iline].lineshape[pos].FVC().X3 += d;
856 break;
857 case Jacobian::Line::ShapeETAX0:
858 d *= band.lines[iline].lineshape[pos].ETA().X0;
859 band_copy.lines[iline].lineshape[pos].ETA().X0 += d;
860 break;
861 case Jacobian::Line::ShapeETAX1:
862 d *= band.lines[iline].lineshape[pos].ETA().X1;
863 band_copy.lines[iline].lineshape[pos].ETA().X1 += d;
864 break;
865 case Jacobian::Line::ShapeETAX2:
866 d *= band.lines[iline].lineshape[pos].ETA().X2;
867 band_copy.lines[iline].lineshape[pos].ETA().X2 += d;
868 break;
869 case Jacobian::Line::ShapeETAX3:
870 d *= band.lines[iline].lineshape[pos].ETA().X3;
871 band_copy.lines[iline].lineshape[pos].ETA().X3 += d;
872 break;
873 case Jacobian::Line::Center:
874 d *= band.lines[iline].F0;
875 band_copy.lines[iline].F0 += d;
876 break;
877 case Jacobian::Line::Strength:
878 d *= band.lines[iline].I0;
879 band_copy.lines[iline].I0 += d;
880 break;
881 case Jacobian::Line::ShapeYX0:
882 case Jacobian::Line::ShapeYX1:
883 case Jacobian::Line::ShapeYX2:
884 case Jacobian::Line::ShapeYX3:
885 case Jacobian::Line::ShapeGX0:
886 case Jacobian::Line::ShapeGX1:
887 case Jacobian::Line::ShapeGX2:
888 case Jacobian::Line::ShapeGX3:
889 case Jacobian::Line::ShapeDVX0:
890 case Jacobian::Line::ShapeDVX1:
891 case Jacobian::Line::ShapeDVX2:
892 case Jacobian::Line::ShapeDVX3:
893 case Jacobian::Line::NLTE:
894 case Jacobian::Line::VMR:
895 case Jacobian::Line::ECS_SCALINGX0:
896 case Jacobian::Line::ECS_SCALINGX1:
897 case Jacobian::Line::ECS_SCALINGX2:
898 case Jacobian::Line::ECS_SCALINGX3:
899 case Jacobian::Line::ECS_BETAX0:
900 case Jacobian::Line::ECS_BETAX1:
901 case Jacobian::Line::ECS_BETAX2:
902 case Jacobian::Line::ECS_BETAX3:
903 case Jacobian::Line::ECS_LAMBDAX0:
904 case Jacobian::Line::ECS_LAMBDAX1:
905 case Jacobian::Line::ECS_LAMBDAX2:
906 case Jacobian::Line::ECS_LAMBDAX3:
907 case Jacobian::Line::ECS_DCX0:
908 case Jacobian::Line::ECS_DCX1:
909 case Jacobian::Line::ECS_DCX2:
910 case Jacobian::Line::ECS_DCX3:
911 case Jacobian::Line::FINAL: {
912 /* do nothing */
913 }
914 }
915
916 // Perform calculations and estimate derivative
917 const auto [dabs, dwork] = ecs_absorption_impl(T, H, P, this_vmr, vmrs, ecs_data, f_grid, zeeman_polarization, band_copy);
918 vec -= dabs;
919 vec /= -d;
920 work &= dwork;
921 } else if (Quantum::Number::StateMatch(target.qid, band.quantumidentity) == Quantum::Number::StateMatchType::Full) {
922 ErrorCorrectedSuddenData ecs_data_copy = ecs_data;
923
924 const auto spec = target.qid.Species();
925 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)
926 switch (target.line) {
927 case Jacobian::Line::ECS_SCALINGX0:
928 d *= ecs_data_copy[spec].scaling.X0;
929 ecs_data_copy[spec].scaling.X0 += d;
930 break;
931 case Jacobian::Line::ECS_SCALINGX1:
932 d *= ecs_data_copy[spec].scaling.X1;
933 ecs_data_copy[spec].scaling.X1 += d;
934 break;
935 case Jacobian::Line::ECS_SCALINGX2:
936 d *= ecs_data_copy[spec].scaling.X2;
937 ecs_data_copy[spec].scaling.X2 += d;
938 break;
939 case Jacobian::Line::ECS_SCALINGX3:
940 d *= ecs_data_copy[spec].scaling.X3;
941 ecs_data_copy[spec].scaling.X3 += d;
942 break;
943 case Jacobian::Line::ECS_BETAX0:
944 d *= ecs_data_copy[spec].beta.X0;
945 ecs_data_copy[spec].beta.X0 += d;
946 break;
947 case Jacobian::Line::ECS_BETAX1:
948 d *= ecs_data_copy[spec].beta.X1;
949 ecs_data_copy[spec].beta.X1 += d;
950 break;
951 case Jacobian::Line::ECS_BETAX2:
952 d *= ecs_data_copy[spec].beta.X2;
953 ecs_data_copy[spec].beta.X2 += d;
954 break;
955 case Jacobian::Line::ECS_BETAX3:
956 d *= ecs_data_copy[spec].beta.X3;
957 ecs_data_copy[spec].beta.X3 += d;
958 break;
959 case Jacobian::Line::ECS_LAMBDAX0:
960 d *= ecs_data_copy[spec].lambda.X0;
961 ecs_data_copy[spec].lambda.X0 += d;
962 break;
963 case Jacobian::Line::ECS_LAMBDAX1:
964 d *= ecs_data_copy[spec].lambda.X1;
965 ecs_data_copy[spec].lambda.X1 += d;
966 break;
967 case Jacobian::Line::ECS_LAMBDAX2:
968 d *= ecs_data_copy[spec].lambda.X1;
969 ecs_data_copy[spec].lambda.X1 += d;
970 break;
971 case Jacobian::Line::ECS_LAMBDAX3:
972 d *= ecs_data_copy[spec].lambda.X0;
973 ecs_data_copy[spec].lambda.X0 += d;
974 break;
975 case Jacobian::Line::ECS_DCX0:
976 d *= ecs_data_copy[spec].collisional_distance.X0;
977 ecs_data_copy[spec].collisional_distance.X0 += d;
978 break;
979 case Jacobian::Line::ECS_DCX1:
980 d *= ecs_data_copy[spec].collisional_distance.X1;
981 ecs_data_copy[spec].collisional_distance.X1 += d;
982 break;
983 case Jacobian::Line::ECS_DCX2:
984 d *= ecs_data_copy[spec].collisional_distance.X2;
985 ecs_data_copy[spec].collisional_distance.X2 += d;
986 break;
987 case Jacobian::Line::ECS_DCX3:
988 d *= ecs_data_copy[spec].collisional_distance.X3;
989 ecs_data_copy[spec].collisional_distance.X3 += d;
990 break;
991 case Jacobian::Line::ShapeG0X0:
992 case Jacobian::Line::ShapeG0X1:
993 case Jacobian::Line::ShapeG0X2:
994 case Jacobian::Line::ShapeG0X3:
995 case Jacobian::Line::ShapeD0X0:
996 case Jacobian::Line::ShapeD0X1:
997 case Jacobian::Line::ShapeD0X2:
998 case Jacobian::Line::ShapeD0X3:
999 case Jacobian::Line::ShapeG2X0:
1000 case Jacobian::Line::ShapeG2X1:
1001 case Jacobian::Line::ShapeG2X2:
1002 case Jacobian::Line::ShapeG2X3:
1003 case Jacobian::Line::ShapeD2X0:
1004 case Jacobian::Line::ShapeD2X1:
1005 case Jacobian::Line::ShapeD2X2:
1006 case Jacobian::Line::ShapeD2X3:
1007 case Jacobian::Line::ShapeFVCX0:
1008 case Jacobian::Line::ShapeFVCX1:
1009 case Jacobian::Line::ShapeFVCX2:
1010 case Jacobian::Line::ShapeFVCX3:
1011 case Jacobian::Line::ShapeETAX0:
1012 case Jacobian::Line::ShapeETAX1:
1013 case Jacobian::Line::ShapeETAX2:
1014 case Jacobian::Line::ShapeETAX3:
1015 case Jacobian::Line::Center:
1016 case Jacobian::Line::Strength:
1017 case Jacobian::Line::ShapeYX0:
1018 case Jacobian::Line::ShapeYX1:
1019 case Jacobian::Line::ShapeYX2:
1020 case Jacobian::Line::ShapeYX3:
1021 case Jacobian::Line::ShapeGX0:
1022 case Jacobian::Line::ShapeGX1:
1023 case Jacobian::Line::ShapeGX2:
1024 case Jacobian::Line::ShapeGX3:
1025 case Jacobian::Line::ShapeDVX0:
1026 case Jacobian::Line::ShapeDVX1:
1027 case Jacobian::Line::ShapeDVX2:
1028 case Jacobian::Line::ShapeDVX3:
1029 case Jacobian::Line::NLTE:
1030 case Jacobian::Line::VMR:
1031 case Jacobian::Line::FINAL: {
1032 /* do nothing */
1033 }
1034 }
1035
1036 // Perform calculations and estimate derivative
1037 const auto [dabs, dwork] = ecs_absorption_impl(T, H, P, this_vmr, vmrs, ecs_data_copy, f_grid, zeeman_polarization, band);
1038 vec -= dabs;
1039 vec /= -d;
1040 work &= dwork;
1041 break; // Leave early because it is a band derivative
1042 }
1043 }
1044 } else {
1045 vec *= 0; // No derivative, so don't mess around and remove everything
1046 }
1047 }
1048
1049 return EcsReturn(std::move(absorption), std::move(jacobian), not work);
1050}
1051
1053 AbsorptionLines& band,
1054 const Tensor4& tempdata,
1055 const Vector& temperatures,
1056 const Numeric P0,
1057 const Index ord)
1058{
1059 // Sizes
1060 const Index S = band.NumBroadeners();
1061 const Index N = band.NumLines();
1062
1063 for (Index i=0; i<S; i++) {
1064 // Allocate Vector for fitting the data
1065 Vector targ(N);
1066
1067 // Rosenkranz 1st order parameter
1068 for (Index j=0; j<N; j++) {
1069 targ = tempdata(1, j, i, joker); // Get values
1070
1071 // Assume linear pressure dependency
1072 targ /= P0;
1073
1074 // Fit to a polynomial
1075 auto [fit, c] = Minimize::curve_fit<Minimize::Polynom>(temperatures, targ, LineShape::ModelParameters::N - 1);
1076 if (not fit) return EXIT_FAILURE;
1077 band.lines[j].lineshape[i].Y() = LineShape::ModelParameters(LineShape::TemperatureModel::POLY, c);
1078 }
1079
1080 // Rosenkranz 2nd order frequency parameter
1081 if (ord > 1) {
1082 for (Index j=0; j<N; j++) {
1083 targ = tempdata(2, j, i, joker); // Get values
1084
1085 // Assume squared pressure dependency
1086 targ /= P0 * P0;
1087
1088 // Fit to a polynomial
1089 auto [fit, c] = Minimize::curve_fit<Minimize::Polynom>(temperatures, targ, LineShape::ModelParameters::N - 1);
1090 if (not fit) return EXIT_FAILURE;
1091 band.lines[j].lineshape[i].DV() = LineShape::ModelParameters(LineShape::TemperatureModel::POLY, c);
1092 }
1093 }
1094
1095 // Rosenkranz 2nd order strength parameter
1096 if (ord > 1) {
1097 for (Index j=0; j<N; j++) {
1098 targ = tempdata(0, j, i, joker); // Get values
1099
1100 // Assume squared pressure dependency
1101 targ /= P0 * P0;
1102
1103 // Fit to a polynomial
1104 auto [fit, c] = Minimize::curve_fit<Minimize::Polynom>(temperatures, targ, LineShape::ModelParameters::N - 1);
1105 if (not fit) return EXIT_FAILURE;
1106 band.lines[j].lineshape[i].G() = LineShape::ModelParameters(LineShape::TemperatureModel::POLY, c);
1107 }
1108 }
1109
1110 // 3rd order broadening parameter [[FIXME: UNTESTED]]
1111 if (ord > 2) {
1112 for (Index j=0; j<N; j++) {
1113 targ = tempdata(3, j, i, joker);
1114
1115 // Assume cubic pressure dependency
1116 targ /= P0 * P0 * P0;
1117
1118 // Fit to a polynomial
1119 auto [fit, c] = Minimize::curve_fit<Minimize::Polynom>(temperatures, targ, LineShape::ModelParameters::N - 1);
1120 if (not fit) return EXIT_FAILURE;
1121 ARTS_USER_ERROR("Not yet implemented: 3rd order line mixing")
1122 // band.lines[j].lineshape[i].DG() = LineShape::ModelParameters(LineShape::TemperatureModel::POLY, c);
1123 }
1124 }
1125 }
1126
1127 // If we reach here, we have to set the band population type
1128 // to LTE and renormalize the strength to physical frequency
1129 band.normalization = Absorption::NormalizationType::SFS;
1130 band.population = Absorption::PopulationType::LTE;
1131
1132 return EXIT_SUCCESS;
1133}
1134
1135
1144 const ConstMatrixView& W,
1145 const AbsorptionLines& band) {
1146 const Index N = dip.nelem();
1147
1148 // Output
1149 Vector Y(N, 0);
1150
1151 // Rosenkranz coefficients
1152 for (Index k=0; k<N; k++) {
1153 for (Index j=0; j<N; j++) {
1154 if (k == j) continue;
1155 Y[k] += 2 * std::abs(dip[j] / dip[k]) * W(j, k) / (band.lines[k].F0 - band.lines[j].F0);
1156 }
1157 }
1158
1159 return Y;
1160}
1161
1162
1171 const ConstMatrixView& W,
1172 const AbsorptionLines& band) {
1173 const Index N = dip.nelem();
1174
1175 // Output
1176 Vector G(N, 0);
1177
1178 // Rosenkranz coefficients
1179 for (Index k=0; k<N; k++) {
1180 for (Index j=0; j<N; j++) {
1181 if (k == j) continue;
1182 G[k] += W(k, j) * W(j, k) / Constant::pow2(band.lines[j].F0 - band.lines[k].F0);
1183 G[k] += Constant::pow2(std::abs(dip[j] / dip[k]) * W(j, k) / (band.lines[j].F0 - band.lines[k].F0));
1184 G[k] += 2 * std::abs(dip[j] / dip[k]) * W(j, k) * W(k, k) / Constant::pow2(band.lines[j].F0 - band.lines[k].F0);
1185 for (Index l=0; l<N; l++) {
1186 if (l == k or l == j) continue;
1187 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));
1188
1189 }
1190 }
1191 }
1192
1193 return G;
1194}
1195
1196
1205 const ConstMatrixView& W,
1206 const AbsorptionLines& band) {
1207 const Index N = dip.nelem();
1208
1209 // Output
1210 Vector DV(N, 0);
1211
1212 // Rosenkranz coefficients
1213 for (Index k=0; k<N; k++) {
1214 for (Index j=0; j<N; j++) {
1215 if (k == j) continue;
1216 DV[k] += W(k, j) * W(j, k) / (band.lines[j].F0 - band.lines[k].F0);
1217 }
1218 }
1219
1220 return DV;
1221}
1222
1223
1225 const ComplexMatrix& W,
1226 const Vector& pop,
1227 const Vector& dip,
1228 const AbsorptionLines& band,
1229 const Numeric frenorm,
1230 const Numeric T,
1231 const Numeric P,
1232 const Numeric QT,
1233 const Numeric QT0,
1234 const Index broadener) {
1235 // Compute and adapt values for the Voigt adaptation
1236 EquivalentLines eig(W, pop, dip);
1237
1238 // Compute all the shifted line centers in band-order
1239 Vector f0(band.NumLines());
1240 for (Index i=0; i<pop.size(); i++) {
1241 f0[i] = band.lines[i].F0 - frenorm + P * band.lines[i].lineshape[broadener].D0().at(T, band.T0);
1242 }
1243
1244 // Find how the band-orders have shifted line order at the current temp/pres
1245 ArrayOfIndex sorting(band.NumLines());
1246 std::iota(sorting.begin(), sorting.end(), 0);
1247 for (Index i=0; i<band.NumLines(); i++) {
1248 for (Index j=i+1; j<band.NumLines(); j++) {
1249 if (f0[j] < f0[i]) {
1250 std::swap(sorting[i], sorting[j]);
1251 std::swap(f0[i], f0[j]);
1252 }
1253 }
1254 }
1255
1256 // Sort the eigen-values and strengths in the same order as the band
1257 eig.sort_by_frequency(f0, sorting);
1258
1259 // Eigenvalue should now be (F0 + P * D0(T) + i * P * G0(T)) + (P * P * DV(T) + i * P * P * P * DG(T)),
1260 // or in form (1) + (2). We only want to keep (2), since it is new from the line mixing
1261 for (Index i=0; i<pop.size(); i++) {
1262 eig.val[i] -= Complex(f0[i], P * band.lines[i].lineshape[broadener].G0().at(T, band.T0));
1263 }
1264
1265 // The strength term should now be (1 + i y(T) + g(T)) * d**2 * rho(T)
1266 // d**2 * rho(T) * F0 * (1 - exp(-hF0/kT)) should be I0(T)
1267 // We only want to keep (i y(T) + g(T)).
1268 // So we divide by d**2 * rho(T) through the use of I0(T) and remove 1 from the real component
1269 for (Index i=0; i<pop.size(); i++) {
1270 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;
1271 eig.str[i] *= - band.lines[i].F0 * std::expm1(- (Constant::h * band.lines[i].F0) / (Constant::k * T)) / i0;
1272 }
1273 eig.str -= 1;
1274
1275 return eig;
1276}
1277
1278
1294 const Vector& temperatures,
1295 const ErrorCorrectedSuddenData& ecs_data,
1296 const Numeric P) {
1297 const Index N = band.NumLines();
1298 const Index M = band.NumBroadeners();
1299 const Index K = temperatures.nelem();
1300
1301 // Weighted center of the band
1302 const Numeric frenorm = band.F_mean();
1303
1304 // Need sorting to put weak lines last, but we need the sorting constant or the output jumps
1305 const ArrayOfIndex sorting = sorted_population_and_dipole(band.T0, band).first;
1306 const Numeric QT0 = single_partition_function(band.T0, band.Isotopologue());
1307
1308 // Output
1309 Tensor4 out(4, N, M, K);
1310
1311 #pragma omp parallel for collapse(2) if (!arts_omp_in_parallel())
1312 for (Index m=0; m<M; m++) {
1313 for (Index k=0; k<K; k++) {
1314 const Numeric T = temperatures[k];
1315 const Numeric QT = single_partition_function(T, band.Isotopologue());
1316
1317 // Relaxation matrix of T0 sorting at T
1318 ComplexMatrix W = single_species_ecs_relaxation_matrix(band, sorting, T, P, ecs_data[band.broadeningspecies[m]], m);
1319 for (Index n=0; n<N; n++) {
1320 W(n, n) += band.lines[sorting[n]].F0 - frenorm;
1321 }
1322
1323 // Populations and dipoles of T0 sorting at T
1324 const auto [pop, dip] = presorted_population_and_dipole(T, sorting, band);
1325
1326 const auto eig = eigenvalue_adaptation_of_relmat(W, pop, dip, band, frenorm, T, P, QT, QT0, m);
1327
1328 out(0, joker, m, k) = eig.str.real();
1329 out(1, joker, m, k) = eig.str.imag();
1330 out(2, joker, m, k) = eig.val.real();
1331 out(3, joker, m, k) = eig.val.imag();
1332 }
1333 }
1334
1335 return out;
1336}
1337
1338
1354 const Vector& temperatures,
1355 const ErrorCorrectedSuddenData& ecs_data,
1356 const Numeric P) {
1357 const Index N = band.NumLines();
1358 const Index M = band.NumBroadeners();
1359 const Index K = temperatures.nelem();
1360
1361 // Weighted center of the band
1362 const Numeric frenorm = band.F_mean();
1363
1364 // Need sorting to put weak lines last, but we need the sorting constant or the output jumps
1365 const ArrayOfIndex sorting = sorted_population_and_dipole(band.T0, band).first;
1366
1367 // Output
1368 Tensor4 out(4, N, M, K);
1369
1370 #pragma omp parallel for collapse(2) if (!arts_omp_in_parallel())
1371 for (Index m=0; m<M; m++) {
1372 for (Index k=0; k<K; k++) {
1373 const Numeric T = temperatures[k];
1374
1375 // Relaxation matrix of T0 sorting at T
1376 ComplexMatrix W = single_species_ecs_relaxation_matrix(band, sorting, T, P, ecs_data[band.broadeningspecies[m]], m);
1377 for (Index n=0; n<N; n++) {
1378 W(n, n) += band.lines[n].F0 - frenorm;
1379 }
1380
1381 // Populations and dipoles of T0 sorting at T
1382 const auto [pop, dip] = presorted_population_and_dipole(T, sorting, band);
1383
1384 out(0, joker, m, k) = RosenkranzG(dip, W.imag(), band);
1385 out(1, joker, m, k) = RosenkranzY(dip, W.imag(), band);
1386 out(2, joker, m, k) = RosenkranzDV(dip, W.imag(), band);
1387 out(3, joker, m, k) = 0;
1388 }
1389 }
1390
1391 return out;
1392}
1393
1395 const Vector& temperatures,
1396 const ErrorCorrectedSuddenData& ecs_data,
1397 const Numeric P0,
1398 const Index ord,
1399 const bool robust,
1400 const bool rosenkranz_adaptation,
1401 const Verbosity& verbosity) {
1403 ARTS_USER_ERROR_IF(P0 <= 0, P0, " Pa is not possible")
1404
1405 ARTS_USER_ERROR_IF(not is_sorted(temperatures),
1406 "The temperature list [",
1407 temperatures,
1408 "] K\n"
1409 "must be fully sorted from low to high")
1410
1412 ord < 1 or ord > 3, "Order not in list [1, 2, 3], is: ", ord)
1413
1415 band,
1416 rosenkranz_adaptation not_eq 0
1417 ? rosenkranz_approximation(band, temperatures, ecs_data, P0)
1418 : ecs_eigenvalue_approximation(band, temperatures, ecs_data, P0),
1419 temperatures,
1420 P0,
1421 ord)) {
1422 ARTS_USER_ERROR_IF(not robust,
1423 "Bad eigenvalue adaptation for band: ",
1424 band.quantumidentity,
1425 '\n')
1426 out3 << "Bad eigenvalue adaptation for band: " << band.quantumidentity
1427 << '\n';
1428
1429 band.normalization = Absorption::NormalizationType::SFS;
1430 band.population = Absorption::PopulationType::LTE;
1431 for (auto& line : band.lines) {
1432 for (auto& sm : line.lineshape.Data()) {
1433 sm.Y() = LineShape::ModelParameters(LineShape::TemperatureModel::None);
1434 sm.G() = LineShape::ModelParameters(LineShape::TemperatureModel::None);
1435 sm.DV() = LineShape::ModelParameters(LineShape::TemperatureModel::None);
1436 }
1437 }
1438 }
1439}
1440
1441
1443 const Vector& temperatures,
1444 const ErrorCorrectedSuddenData& ecs_data,
1445 const Vector& pressures) {
1446 const Index N = band.NumLines();
1447 const Index M = band.NumBroadeners();
1448 const Index K = temperatures.nelem();
1449 const Index L = pressures.size();
1450
1451 Tensor5 out(4, N, M, K, L);
1452 for (Index l=0; l<L; l++) {
1453 out(joker, joker, joker, joker, l) = ecs_eigenvalue_approximation(band, temperatures, ecs_data, pressures[l]);
1454 }
1455 return out;
1456}
1457} // namespace Absorption::LineMixing
1458
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
MatrixView imag()
Get a view of the imaginary parts of the matrix.
The ComplexMatrix class.
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.
The ComplexVector class.
Index nelem() const noexcept
A constant view of a Matrix.
Definition: matpackI.h:1064
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:554
Index size() const noexcept
Definition: matpackI.h:557
The MatrixView class.
Definition: matpackI.h:1182
Implements rational numbers to work with other ARTS types.
Definition: rational.h:43
constexpr bool isDefined() const noexcept
Is the object defined.
Definition: rational.h:100
The Tensor4 class.
Definition: matpackIV.h:429
The Tensor5 class.
Definition: matpackV.h:516
The Vector class.
Definition: matpackI.h:922
Input manipulator class for doubles to enable nan and inf parsing.
Definition: double_imanip.h:42
Constants of physical expressions as constexpr.
Helper macros for debugging.
#define ARTS_NOEXCEPT
Definition: debug.h:80
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
constexpr bool good_enum(EnumType x) noexcept
Checks if the enum number is good.
Definition: enums.h:22
#define abs(x)
#define beta
#define min(a, b)
#define dabs(x)
#define max(a, b)
void diagonalize(MatrixView P, VectorView WR, VectorView WI, ConstMatrixView A)
Matrix Diagonalization.
Definition: lin_alg.cc:241
MatrixViewMap MapToEigen(Tensor3View &A, const Index &i)
Definition: lin_alg.cc:684
Linear algebra functions.
Numeric boltzman_factor(Numeric T, Numeric E0)
Computes exp(- E0/kT)
Definition: linescaling.cc:114
Numeric boltzman_ratio(const Numeric &T, const Numeric &T0, const Numeric &E0)
Computes exp(E0/c (T - T0) / (T * T0))
Definition: linescaling.cc:97
Numeric single_partition_function(const Numeric &T, const Species::IsotopeRecord &ir)
Computes the partition function at one temperature.
Definition: linescaling.cc:23
bool is_sorted(ConstVectorView x)
Checks if a vector is sorted in ascending order.
Definition: logic.cc:202
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:63
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
Numeric & real_val(Complex &c) noexcept
Return a reference to the real value of c.
const Joker joker
std::complex< Numeric > Complex
Declarations having to do with the four output streams.
#define CREATE_OUT3
Definition: messages.h:206
void relaxation_matrix_offdiagonal(MatrixView W, const AbsorptionLines &band, const ArrayOfIndex &sorting, const SpeciesErrorCorrectedSuddenData &rovib_data, const Numeric T) ARTS_NOEXCEPT
Definition: linemixing.cc:441
EnergyFunction erot_selection(const SpeciesIsotopeRecord &isot)
Definition: linemixing.cc:431
Numeric reduced_dipole(const Rational Ju, const Rational Jl, const Rational N)
Definition: linemixing.cc:133
constexpr Numeric erot(const Rational N, const Rational j=-1)
Definition: linemixing.cc:275
void relaxation_matrix_offdiagonal(MatrixView W, const AbsorptionLines &band, const ArrayOfIndex &sorting, const SpeciesErrorCorrectedSuddenData &rovib_data, const Numeric T) ARTS_NOEXCEPT
Definition: linemixing.cc:322
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...
Definition: linemixing.cc:1224
std::ostream & operator<<(std::ostream &os, const EquivalentLines &eqv)
Definition: linemixing.cc:93
Tensor4 ecs_eigenvalue_approximation(const AbsorptionLines &band, const Vector &temperatures, const ErrorCorrectedSuddenData &ecs_data, const Numeric P)
Definition: linemixing.cc:1293
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)
Definition: linemixing.cc:633
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)
Definition: linemixing.cc:1394
std::pair< ArrayOfIndex, PopulationAndDipole > sorted_population_and_dipole(const Numeric T, const AbsorptionLines &band)
Definition: linemixing.cc:210
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)
Definition: linemixing.cc:705
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)
Definition: linemixing.cc:555
Vector RosenkranzY(const Vector &dip, const ConstMatrixView &W, const AbsorptionLines &band)
Definition: linemixing.cc:1143
Vector RosenkranzG(const Vector &dip, const ConstMatrixView &W, const AbsorptionLines &band)
Definition: linemixing.cc:1170
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.
Definition: linemixing.cc:1052
Vector RosenkranzDV(const Vector &dip, const ConstMatrixView &W, const AbsorptionLines &band)
Definition: linemixing.cc:1204
Tensor5 ecs_eigenvalue_adaptation_test(const AbsorptionLines &band, const Vector &temperatures, const ErrorCorrectedSuddenData &ecs_data, const Vector &pressures)
Definition: linemixing.cc:1442
std::function< Numeric(const Rational)> EnergyFunction
Definition: linemixing.h:82
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)
Definition: linemixing.cc:602
std::istream & operator>>(std::istream &is, SpeciesErrorCorrectedSuddenData &srbd)
Definition: linemixing.cc:112
PopulationAndDipole presorted_population_and_dipole(const Numeric T, const ArrayOfIndex &presorting, const AbsorptionLines &band)
Definition: linemixing.cc:224
Tensor4 rosenkranz_approximation(const AbsorptionLines &band, const Vector &temperatures, const ErrorCorrectedSuddenData &ecs_data, const Numeric P)
Definition: linemixing.cc:1353
Numeric reduced_rovibrational_dipole(Rational Jf, Rational Ji, Rational lf, Rational li, Rational k=Rational(1))
Compute the reduced rovibrational dipole moment.
constexpr auto pow2(T x) noexcept -> decltype(x *x)
power of two
Definition: constants.h:66
constexpr auto pow3(T x) noexcept -> decltype(pow2(x) *x)
power of three
Definition: constants.h:72
constexpr auto kaycm2joule(T x) noexcept -> decltype(x *kaycm2freq(h))
Conversion from cm-1 to Joule.
Definition: constants.h:537
constexpr auto mhz2joule(T x) noexcept -> decltype(hz2joule(x) *1e6)
Conversion from MHz to Joule.
Definition: constants.h:549
constexpr auto kelvin2joule(T x) noexcept -> decltype(x *k)
Conversion from Kelvin to Joule.
Definition: constants.h:555
constexpr Numeric l(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Index j, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
Temperature
Definition: jacobian.h:58
Vector mass(const ConstVectorView &atmospheric_vmrs, const ArrayOfArrayOfSpeciesTag &atmospheric_species, const ArrayOfSpecies &lineshape_species, const SpeciesIsotopologueRatios &ir) ARTS_NOEXCEPT
Returns a mass vector for this model's main calculations.
Vector vmrs(const ConstVectorView &atmospheric_vmrs, const ArrayOfArrayOfSpeciesTag &atmospheric_species, const ArrayOfSpecies &lineshape_species) ARTS_NOEXCEPT
Returns a VMR vector for this model's main calculations.
Polarization
Zeeman polarization selection.
Definition: zeemandata.h:44
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:431
This file contains declerations of functions of physical character.
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Definition: physics_funcs.h:70
#define d
#define w
#define a
#define c
#define b
Contains the rational class definition.
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:736
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:728
constexpr bool iseven(const Rational r) noexcept
Returns true if even integer.
Definition: rational.h:925
#define N
Definition: rng.cc:164
#define M
Definition: rng.cc:165
Contains recomputed equivalent lines (sorting is unknown)
Definition: linemixing.h:18
void sort_by_frequency(Vector &f, const ArrayOfIndex &sorting)
Definition: linemixing.cc:69
Rovibrational line mixing data following the ideas of Collisional Effects On Molecular Spectra by J.
Definition: linemixing.h:134
Index pos(Species::Species spec) const noexcept
Definition: linemixing.h:162
Contains the population distribution and dipole.
Definition: linemixing.h:56
ArrayOfIndex sort(const AbsorptionLines &band) noexcept
Definition: linemixing.cc:165
PopulationAndDipole(const Numeric T, const AbsorptionLines &band) noexcept
Definition: linemixing.cc:139
Numeric Q(const Rational J, const Numeric T, const Numeric T0, const Numeric energy) const noexcept
Definition: linemixing.cc:232
Numeric Omega(const Numeric T, const Numeric T0, const Numeric other_mass, const Numeric energy_x, const Numeric energy_xm2) const noexcept
Definition: linemixing.cc:239
Numeric ZeemanSplitting(size_t k, Zeeman::Polarization type, Index i) const noexcept
Returns the splitting of a Zeeman split line.
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.
Array< SingleLine > lines
A list of individual lines.
Index NumLines() const noexcept
Number of lines.
Index NumBroadeners() const noexcept
Number of broadening species.
Numeric DopplerConstant(Numeric T) const noexcept
NormalizationType normalization
Line normalization type.
LineShape::Output ShapeParameters(size_t k, Numeric T, Numeric P, const Vector &vmrs) const noexcept
Line shape parameters.
Species::IsotopeRecord Isotopologue() const noexcept
Isotopologue Index.
LineShape::Type lineshapetype
Type of line shape.
ArrayOfSpecies broadeningspecies
A list of broadening species.
Index ZeemanCount(size_t k, Zeeman::Polarization type) const noexcept
Returns the number of Zeeman split lines.
bool selfbroadening
Does the line broadening have self broadening.
QuantumIdentifier quantumidentity
Catalog ID.
Numeric ZeemanStrength(size_t k, Zeeman::Polarization type, Index i) const noexcept
Returns the strength of a Zeeman split line.
bool DoVmrDerivative(const QuantumIdentifier &qid) const noexcept
Coefficients and temperature model for SingleSpeciesModel.
static constexpr Index N
StateMatchType operates so that a check less than a level should be 'better', bar None.
Struct containing all information needed about one isotope.
Definition: isotopologues.h:16
String FullName() const noexcept
Definition: isotopologues.h:49
Species spec
Species type as defined in species.h.
Definition: isotopologues.h:18
std::string_view isotname
A custom name that is unique for this Species type.
Definition: isotopologues.h:21
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)