ARTS 2.5.11 (git: 725533f0)
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#include <vector>
7
8#include "arts_conversions.h"
9#include "debug.h"
10#include "lin_alg.h"
11#include "linemixing.h"
12#include "lineshape.h"
13#include "matpack_complex.h"
14#include "matpack_data.h"
15#include "matpack_eigen.h"
16#include "messages.h"
17#include "minimize.h"
18#include "physics_funcs.h"
19#include "quantum_numbers.h"
20#include "rational.h"
21#include "species.h"
22#include "wigner_functions.h"
23
24#if DO_FAST_WIGNER
25#define WIGNER3 fw3jja6
26#define WIGNER6 fw6jja
27#else
28#define WIGNER3 wig3jj
29#define WIGNER6 wig6jj
30#endif
31
32Numeric wig3(const Rational& a,
33 const Rational& b,
34 const Rational& c,
35 const Rational& d,
36 const Rational& e,
37 const Rational& f) noexcept {
38 return WIGNER3(
39 a.toInt(2), b.toInt(2), c.toInt(2), d.toInt(2), e.toInt(2), f.toInt(2));
40}
41
42Numeric wig6(const Rational& a,
43 const Rational& b,
44 const Rational& c,
45 const Rational& d,
46 const Rational& e,
47 const Rational& f) noexcept {
48 return WIGNER6(
49 a.toInt(2), b.toInt(2), c.toInt(2), d.toInt(2), e.toInt(2), f.toInt(2));
50}
51
52#undef WIGNER3
53#undef WIGNER6
54
56EquivalentLines::EquivalentLines(const ComplexMatrix& W,
57 const Vector& pop,
58 const Vector& dip) noexcept :
59 val(pop.nelem(), 0), str(pop.nelem(), 0) {
60 /* FIXME: (Added 2021-01-19; Richard Larsson)
61 *
62 * This function cannot easily be used with partial derivatives
63 *
64 * Doing so would allow the entire ECS approach to be analytical
65 * in its derivatives.
66 *
67 * The problem in short:
68 * There is an Eigenvalue decomposition happening (W = V diag(e) V^-1)
69 *
70 * We use V and e in a strange way. We do not know how to get the
71 * partial derivatives of these two variables
72 *
73 * The question:
74 * How do we compute dV and de? We can get dW somewhat easily.
75 */
76
77 const Index n = pop.nelem();
78
79 // Compute values
80 ComplexMatrix V(n, n);
81
82 // Main computations
83 diagonalize(V, val, W);
84
85 // Do the matrix forward multiplication
86 for (Index i=0; i<n; i++) {
87 for (Index j=0; j<n; j++) {
88 str[i] += dip[j] * V(j, i);
89 }
90 }
91
92 // Do the matrix backward multiplication
93 inv(V, V);
94 for (Index i=0; i<n; i++) {
95 Complex z(0, 0);
96 for (Index j=0; j<n; j++) {
97 z += pop[j] * dip[j] * V(i, j);
98 }
99 str[i] *= z;
100 }
101}
102
103
104void EquivalentLines::sort_by_frequency(Vector& f, const ArrayOfIndex& sorting) {
105 const Index N = val.nelem();
106
107 // First sort by frequency and assume this sorting
108 // is a shifted version of the line centers
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]);
114 }
115 }
116 }
117
118 const Vector fc = f;
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]];
125 }
126}
127
128std::ostream& operator<<(std::ostream& os, const EquivalentLines& eqv) {
129 const Index n = eqv.str.nelem();
130 for (Index i=0; i<n; i++) {
131 if (i) os << '\n';
132 os << eqv.str[i] << ' ' << eqv.val[i];
133 }
134 return os;
135}
136
137std::ostream& operator<<(std::ostream& os, const SpeciesErrorCorrectedSuddenData& srbd) {
138 os << Species::toShortName(srbd.spec);
139 os << ' ' << srbd.scaling;
140 os << ' ' << srbd.beta;
141 os << ' ' << srbd.lambda;
142 os << ' ' << srbd.collisional_distance;
143 os << ' ' << srbd.mass;
144 return os;
145}
146
147std::istream& operator>>(std::istream& is, SpeciesErrorCorrectedSuddenData& srbd) {
148 std::string spec_name;
149 is >> spec_name;
150 is >> srbd.scaling;
151 is >> srbd.beta;
152 is >> srbd.lambda;
153 is >> srbd.collisional_distance;
154 is >> double_imanip() >> srbd.mass;
155 srbd.spec = Species::fromShortName(spec_name);
156 ARTS_USER_ERROR_IF(not good_enum(srbd.spec), "Cannot recognize species: ", spec_name)
157 return is;
158}
159
160namespace Makarov2020etal {
168Numeric reduced_dipole(const Rational Ju, const Rational Jl, const Rational N) {
169 return (iseven(Jl + N) ? 1 : -1) * sqrt(6 * (2*Jl + 1) * (2*Ju + 1)) * wigner6j(1, 1, 1, Jl, Ju, N);
170};
171} // namespace Makarov2020etal
172
173
175 const AbsorptionLines& band) :
176 pop(band.NumLines()), dip(band.NumLines()) {
177 const Index n = band.NumLines();
178
179 const Numeric QT = single_partition_function(T, band.Isotopologue());
180 const Numeric QT0 = single_partition_function(band.T0, band.Isotopologue());
181 const Numeric ratiopart = QT0 / QT;
182
183 for (Index i=0; i<n; i++) {
184 const Numeric pop0 = (band.lines[i].gupp / QT0) * boltzman_factor(band.T0, band.lines[i].E0);
185 pop[i] = pop0 * ratiopart * boltzman_ratio(T, band.T0, band.lines[i].E0);
186 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)));
187
188 // Adjust the sign depending on type
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];
192 dip[i] *= std::signbit(Makarov2020etal::reduced_dipole(J.upp(), J.low(), N.upp())) ? - 1 : 1;
193 } else if (band.population == Absorption::PopulationType::ByRovibLinearDipoleLineMixing) {
194 // pass
195 }
196 }
197}
198
199
201 const Index N = pop.nelem();
202
203 // List that starts as [0, 1, ... N-2, N-1]
204 ArrayOfIndex out(N);
205 std::iota(out.begin(), out.end(), 0);
206
207 // Strength
208 Vector s(N);
209 for (Index i=0; i<N; i++) {
210 s[i] = band.lines[i].F0 * pop[i] * Math::pow2(dip[i]);
211 }
212
213 // Sorter
214 for (Index i=0; i<N; i++) {
215 for (Index j=i+1; j<N; j++) {
216 if (s[j] > s[i]) {
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]);
221 }
222 }
223 }
224
225 return out;
226}
227
228
229void PopulationAndDipole::sort(const ArrayOfIndex& presorting) noexcept {
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]];
235 }
236}
237
238
245std::pair<ArrayOfIndex, PopulationAndDipole> sorted_population_and_dipole(const Numeric T,
246 const AbsorptionLines& band) {
247 PopulationAndDipole tp(T, band);
248 return {tp.sort(band), tp};
249}
250
251
260 const ArrayOfIndex& presorting,
261 const AbsorptionLines& band) {
262 PopulationAndDipole tp(T, band);
263 tp.sort(presorting);
264 return tp;
265}
266
267Numeric SpeciesErrorCorrectedSuddenData::Q(const Rational J,
268 const Numeric T,
269 const Numeric T0,
270 const Numeric energy) const {
271 return std::exp(- beta.at(T, T0) * energy / (Constant::k * T)) * scaling.at(T, T0) / pow(J * (J+1), lambda.at(T, T0));
272}
273
275 const Numeric T0,
276 const Numeric other_mass,
277 const Numeric energy_x,
278 const Numeric energy_xm2) const {
279 using Constant::h;
280 using Constant::k;
281 using Constant::pi;
282 using Constant::m_u;
283 using Constant::h_bar;
284 using Math::pow2;
285
286 // Constants for the expression
287 constexpr Numeric fac = 8 * k / (m_u * pi);
288
289 const Numeric wnnm2 = (energy_x - energy_xm2) / h_bar;
290
291 const Numeric inv_eff_mass = 1 / mass + 1 / other_mass;
292 const Numeric v_bar_pow2 = fac*T*inv_eff_mass;
293 const Numeric tauc_pow2 = pow2(collisional_distance.at(T, T0)) / v_bar_pow2;
294
295 return 1.0 / pow2(1 + pow2(wnnm2) * tauc_pow2 / 24.0);
296}
297
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;
311
312 if constexpr (rescale_pure_rotational) {
313 return erot<false>(N, J) - erot<false>(1, 0);
314 } else {
316 using Math::pow2;
317 using Math::pow3;
318
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;
328
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);
334
335 if (J < N) {
336 if (N == 1) // erot<false>(1, 0)
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));
339 }
340 if (J > N)
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);
343 }
344}
345
357 MatrixView W,
358 const AbsorptionLines& band,
359 const ArrayOfIndex& sorting,
360 const SpeciesErrorCorrectedSuddenData& rovib_data,
361 const Numeric T) {
363
364 const auto bk = [](const Rational& r) -> Numeric { return sqrt(2 * r + 1); };
365
366 const Index n = band.NumLines();
367 if (n == 0) return;
368
369 auto& S = band.quantumidentity.val[QuantumNumberType::S];
370 const Rational Si = S.upp();
371 const Rational Sf = S.low();
372
373 Vector dipr(n);
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];
377
378 dipr[i] = reduced_dipole(J.upp(), J.low(), N.upp());
379 }
380
381 const auto maxL = temp_init_size(band.max(QuantumNumberType::J), band.max(QuantumNumberType::N));
382
383 const auto Om = [&]() {
384 std::vector<Numeric> out(maxL);
385 for (Index i = 0; i < maxL; i++)
386 out[i] = rovib_data.Omega(
387 T, band.T0, band.SpeciesMass(), erot(i), erot(i - 2));
388 return out;
389 }();
390
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));
395 return out;
396 }();
397
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];
402
403 const Rational Ji = J.upp();
404 const Rational Jf = J.low();
405 const Rational Ni = N.upp();
406 const Rational Nf = N.low();
407
408 for (Index j = 0; j < n; j++) {
409 if (i == j) continue;
410
411 auto& J_p = band.lines[sorting[j]].localquanta.val[QuantumNumberType::J];
412 auto& N_p = band.lines[sorting[j]].localquanta.val[QuantumNumberType::N];
413
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();
418
419 if (Jf_p > Jf) continue;
420
421 // Tran etal 2006 symbol with modifications:
422 // 1) [Ji] * [Ji_p] instead of [Ji_p] ^ 2 in partial accordance with Makarov etal 2013
423 Numeric sum = 0;
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) *
426 bk(Ji_p);
427 const auto [L0, L1] =
428 wigner_limits(wigner3j_limits<3>(Ni_p, Ni),
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()];
437 }
438 sum *= scl * Om[Ni.toIndex()];
439
440 // Add to W and rescale to upwards element by the populations
441 W(i, j) = sum;
442 W(j, i) = sum * std::exp((erot(Nf_p) - erot(Nf)) / kelvin2joule(T));
443 }
444 }
445 wig_temp_free();
446
447 ARTS_USER_ERROR_IF(errno == EDOM, "Cannot compute the wigner symbols")
448
449 // Sum rule correction
450 for (Index i = 0; i < n; i++) {
451 Numeric sumlw = 0.0;
452 Numeric sumup = 0.0;
453
454 for (Index j = 0; j < n; j++) {
455 if (j > i) {
456 sumlw += dipr[j] * W(j, i);
457 } else {
458 sumup += dipr[j] * W(j, i);
459 }
460 }
461
462 const Rational Ni =
463 band.lines[sorting[i]].localquanta.val[QuantumNumberType::N].low();
464 for (Index j = i + 1; j < n; j++) {
465 const Rational Nj =
466 band.lines[sorting[j]].localquanta.val[QuantumNumberType::N].low();
467 if (sumlw == 0) {
468 W(j, i) = 0.0;
469 W(i, j) = 0.0;
470 } else {
471 W(j, i) *= -sumup / sumlw;
472 W(i, j) = W(j, i) * std::exp((erot(Ni) - erot(Nj)) /
473 kelvin2joule(T)); // This gives LTE
474 }
475 }
476 }
477}
478} // namespace Makarov2020etal
479
480
489namespace LinearRovibErrorCorrectedSudden {
490
492 if (isot.spec == Species::Species::CarbonDioxide and isot.isotname == "626") {
493 return [](const Rational J) -> Numeric {return Conversion::kaycm2joule(0.39021) * Numeric(J * (J + 1));};
494 }
495
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();};
498}
499
501 const AbsorptionLines& band,
502 const ArrayOfIndex& sorting,
503 const SpeciesErrorCorrectedSuddenData& rovib_data,
504 const Numeric T) {
506
507 const Index n = band.NumLines();
508 if (not n) return;
509
510 // These are constant for a band
511 auto& l2 = band.quantumidentity.val[QuantumNumberType::l2];
512 Rational li = l2.upp();
513 Rational lf = l2.low();
514
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;
519
520 const EnergyFunction erot = erot_selection(band.Isotopologue());
521
522 Vector dipr(n);
523 for (Index i=0; i<n; i++) {
524 auto& J = band.lines[i].localquanta.val[QuantumNumberType::J];
525 dipr[i] = Absorption::reduced_rovibrational_dipole(J.upp(), J.low(), lf, li);
526 }
527
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);
533
534 for (Index j=0; j<n; j++) {
535 if(i == j) continue;
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);
540
541 // Select upper quantum number
542 if (Jf_p > Jf) continue;
543
544 Index L = std::max(std::abs((Ji-Ji_p).toIndex()), std::abs((Jf-Jf_p).toIndex()));
545 L += L % 2;
546 const Index Lf = std::min((Ji+Ji_p).toIndex(), (Jf+Jf_p).toIndex());
547
548 Numeric sum=0;
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;
556 }
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));
559 sum *= scl;
560
561 // Add to W and rescale to upwards element by the populations
562 W(j, i) = sum;
563 W(i, j) = sum * std::exp((erot(Jf_p) - erot(Jf)) / kelvin2joule(T));
564 }
565 }
566
567 // Undocumented negative absolute sign
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)
571 W(i, j) *= -1;
572
573 // Sum rule correction
574 for (Index i=0; i<n; i++) {
575 Numeric sumlw = 0.0;
576 Numeric sumup = 0.0;
577
578 for (Index j=0; j<n; j++) {
579 if (j > i) {
580 sumlw += std::abs(dipr[j]) * W(j, i); // Undocumented abs-sign
581 } else {
582 sumup += std::abs(dipr[j]) * W(j, i); // Undocumented abs-sign
583 }
584 }
585
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();
589 if (sumlw == 0) {
590 W(j, i) = 0.0;
591 W(i, j) = 0.0;
592 } else {
593 W(j, i) *= - sumup / sumlw;
594 W(i, j) = W(j, i) * std::exp((erot(Ji) - erot(Jj)) / kelvin2joule(T)); // This gives LTE
595 }
596 }
597 }
598}
599} // namespace LinearRovibErrorCorrectedSudden
600
614 const ArrayOfIndex& sorting,
615 const Numeric T,
616 const Numeric P,
617 const SpeciesErrorCorrectedSuddenData& species_ecs_data,
618 const Index species_pos) {
619 const Index N = band.NumLines();
620
621 // Allocate the matrix
622 ComplexMatrix W(N, N, 0);
623
624 // Fill diagonal keeping track of the reshuffle (sorting)
625 for (Index I=0; I<N; I++) {
626 const Index i = sorting[I];
627 auto shape = band.ShapeParameters(i, T, P, species_pos);
628 W(I, I) = Complex(shape.D0, shape.G0); // nb. This must be fixed if SDVP ???
629 }
630
631 // Set the off-diagonal part of the matrix for this broadener
632 switch (band.population) {
633 case PopulationType::ByMakarovFullRelmat:
634 Makarov2020etal::relaxation_matrix_offdiagonal(W.imag(), band, sorting, species_ecs_data, T);
635 break;
636 case PopulationType::ByRovibLinearDipoleLineMixing: {
637 LinearRovibErrorCorrectedSudden::relaxation_matrix_offdiagonal(W.imag(), band, sorting, species_ecs_data, T);
638 } break;
639 default:
640 ARTS_ASSERT(false, "Bad type, we don't support band population type: ", band.population,
641 "\nin this code. It must either be added or computations aborted earlier");
642 break;
643 }
644
645 return W;
646}
647
648
660ComplexMatrix ecs_relaxation_matrix(const Numeric T,
661 const Numeric P,
662 const Vector& vmrs,
663 const ErrorCorrectedSuddenData& ecs_data,
664 const AbsorptionLines& band,
665 const ArrayOfIndex& sorting,
666 const Numeric frenorm) {
667 const Index N = band.NumLines();
668 const Index M = vmrs.nelem();
669
670 // Create output
671 ComplexMatrix W(N, N, 0);
672
673 // Loop over all the broadeners
674 for (Index k=0; k<M; k++) {
675 if (vmrs[k] == 0) continue;
676
677 // Sum up all atmospheric components
678 matpack::eigen::mat(W).noalias() += vmrs[k] * matpack::eigen::mat(
679 single_species_ecs_relaxation_matrix(band, sorting, T, P, ecs_data[band.broadeningspecies[k]], k));
680 }
681
682 // Deal with line frequency and its re-normalization
683 for (Index i=0; i<N; i++) {
684 real_val(W(i, i)) += band.lines[sorting[i]].F0 - frenorm;
685 }
686
687 return W;
688}
689
690
691std::pair<ComplexVector, bool> ecs_absorption_impl(const Numeric T,
692 const Numeric H,
693 const Numeric P,
694 const Numeric this_vmr,
695 const Vector& vmrs,
696 const ErrorCorrectedSuddenData& ecs_data,
697 const Vector& f_grid,
698 const Zeeman::Polarization zeeman_polarization,
699 const AbsorptionLines& band) {
700 constexpr Numeric sq_ln2pi = Constant::sqrt_ln_2 / Constant::sqrt_pi;
701
702 // Weighted center of the band
703 const Numeric frenorm = band.F_mean();
704
705 // Band Doppler broadening constant
706 const Numeric GD_div_F0 = band.DopplerConstant(T);
707
708 // Sorted population
709 const auto [sorting, tp] = sorted_population_and_dipole(T, band);
710
711 // Relaxation matrix
712 const ComplexMatrix W = ecs_relaxation_matrix(T, P, vmrs, ecs_data, band, sorting, frenorm);
713
714 // Equivalent lines computations
715 const EquivalentLines eqv(W, tp.pop, tp.dip);
716
717 // Return the absorption and if this works
718 std::pair<ComplexVector, bool> retval{ComplexVector(f_grid.nelem(), 0), false};
719 auto& [absorption,works] = retval;
720
721 // Add the lines
722 for (Index i=0; i<band.NumLines(); i++) {
723 // Zeeman lines if necessary
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);
728
729 if (band.lineshapetype == LineShape::Type::LP) {
730 for (Index iv=0; iv<f_grid.nelem(); iv++) {
731 absorption[iv] -= 1i * Sz * ((eqv.str[i] / (f_grid[iv] - frenorm - dzeeman - eqv.val[i]))) / (Constant::sqrt_ln_2 * Constant::sqrt_pi);
732 }
733 } else /*if (band.LineShapeType() == LineShape::Type::VP)*/ {
734 const Numeric gamd = GD_div_F0 * (eqv.val[i].real() + frenorm + dzeeman);
735 const Numeric cte = Constant::sqrt_ln_2 / gamd;
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;
740 }
741 }
742 }
743 }
744
745 // Adjust by frequency and number density
746 const Numeric numdens = this_vmr * number_density(P, T);
747 for (Index iv=0; iv<f_grid.nelem(); iv++) {
748 const Numeric f = f_grid[iv];
749 const Numeric fact = - f * std::expm1(- (Constant::h * f) / (Constant::k * T));
750 absorption[iv] *= fact * numdens * sq_ln2pi;
751
752 // correct for too low value
753 if (auto& a = real_val(absorption[iv]); not std::isnormal(a) or a < 0) {
754 absorption[iv] = 0;
755 works = false;
756 }
757 }
758
759 return retval;
760}
761
762
763EcsReturn ecs_absorption(const Numeric T,
764 const Numeric H,
765 const Numeric P,
766 const Numeric this_vmr,
767 const Vector& vmrs,
768 const ErrorCorrectedSuddenData& ecs_data,
769 const Vector& f_grid,
770 const Zeeman::Polarization zeeman_polarization,
771 const AbsorptionLines& band,
772 const ArrayOfRetrievalQuantity& jacobian_quantities) {
773 auto [absorption, work] = ecs_absorption_impl(T, H, P, this_vmr, vmrs, ecs_data, f_grid, zeeman_polarization, band);
774
775 // Start as original, so remove new and divide with the negative to get forward derivative
776 ArrayOfComplexVector jacobian(jacobian_quantities.nelem(), absorption);
777
778 for (Index i=0; i<jacobian_quantities.nelem(); i++) {
779 auto& vec = jacobian[i];
780 auto& target = jacobian_quantities[i].Target();
781
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);
785 vec -= dabs;
786 vec /= -dT;
787 work &= dwork;
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);
791 vec -= dabs;
792 vec /= -dH;
793 work &= dwork;
794 } else if (target.isWind()) {
795 const Numeric df = target.perturbation;
796 Vector f_grid_copy = f_grid;
797 f_grid_copy += df;
798 const auto [dabs, dwork] = ecs_absorption_impl(T, H, P, this_vmr, vmrs, ecs_data, f_grid_copy, zeeman_polarization, band);
799 vec -= dabs;
800 vec /= df;
801 work &= dwork;
802 } else if (target == Jacobian::Line::VMR) {
803 if (band.DoVmrDerivative(target.qid)) {
804 Vector vmrs_copy = vmrs;
805 Numeric this_vmr_copy = this_vmr;
806 const Numeric dvmr = target.perturbation;
807
808 // Alter the VMRs for self
809 if (band.Isotopologue() == target.qid.Isotopologue()) {
810 this_vmr_copy += dvmr;
811 if (band.selfbroadening) vmrs_copy[0] += dvmr; // First value is self if band has self broadener
812 } else {
813 for (Index j=band.selfbroadening; j<band.broadeningspecies.nelem()-band.bathbroadening; j++) {
814 if (band.broadeningspecies[j] == target.qid.Species()) {
815 vmrs_copy[j] += dvmr;
816 }
817 }
818 }
819
820 // Computations
821 const auto [dabs, dwork] = ecs_absorption_impl(T, H, P, this_vmr_copy, vmrs_copy, ecs_data, f_grid, zeeman_polarization, band);
822 vec -= dabs;
823 vec /= -dvmr;
824 work &= dwork;
825 }
826 } else if (target.needQuantumIdentity()) {
827 Numeric d=1e-6;
828
829 for (Index iline=0; iline<band.NumLines(); iline++) {
830 if (Quantum::Number::StateMatch(target.qid, band.lines[iline].localquanta, band.quantumidentity) == Quantum::Number::StateMatchType::Full) {
831 AbsorptionLines band_copy = band;
832
833 const Index pos = band.BroadeningSpeciesPosition(target.qid.Species());
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;
838 break;
839 case Jacobian::Line::ShapeG0X1:
840 d *= band.lines[iline].lineshape[pos].G0().X1;
841 band_copy.lines[iline].lineshape[pos].G0().X1 += d;
842 break;
843 case Jacobian::Line::ShapeG0X2:
844 d *= band.lines[iline].lineshape[pos].G0().X2;
845 band_copy.lines[iline].lineshape[pos].G0().X2 += d;
846 break;
847 case Jacobian::Line::ShapeG0X3:
848 d *= band.lines[iline].lineshape[pos].G0().X3;
849 band_copy.lines[iline].lineshape[pos].G0().X3 += d;
850 break;
851 case Jacobian::Line::ShapeD0X0:
852 d *= band.lines[iline].lineshape[pos].D0().X0;
853 band_copy.lines[iline].lineshape[pos].D0().X0 += d;
854 break;
855 case Jacobian::Line::ShapeD0X1:
856 d *= band.lines[iline].lineshape[pos].D0().X1;
857 band_copy.lines[iline].lineshape[pos].D0().X1 += d;
858 break;
859 case Jacobian::Line::ShapeD0X2:
860 d *= band.lines[iline].lineshape[pos].D0().X2;
861 band_copy.lines[iline].lineshape[pos].D0().X2 += d;
862 break;
863 case Jacobian::Line::ShapeD0X3:
864 d *= band.lines[iline].lineshape[pos].D0().X3;
865 band_copy.lines[iline].lineshape[pos].D0().X3 += d;
866 break;
867 case Jacobian::Line::ShapeG2X0:
868 d *= band.lines[iline].lineshape[pos].G2().X0;
869 band_copy.lines[iline].lineshape[pos].G2().X0 += d;
870 break;
871 case Jacobian::Line::ShapeG2X1:
872 d *= band.lines[iline].lineshape[pos].G2().X1;
873 band_copy.lines[iline].lineshape[pos].G2().X1 += d;
874 break;
875 case Jacobian::Line::ShapeG2X2:
876 d *= band.lines[iline].lineshape[pos].G2().X2;
877 band_copy.lines[iline].lineshape[pos].G2().X2 += d;
878 break;
879 case Jacobian::Line::ShapeG2X3:
880 d *= band.lines[iline].lineshape[pos].G2().X3;
881 band_copy.lines[iline].lineshape[pos].G2().X3 += d;
882 break;
883 case Jacobian::Line::ShapeD2X0:
884 d *= band.lines[iline].lineshape[pos].D2().X0;
885 band_copy.lines[iline].lineshape[pos].D2().X0 += d;
886 break;
887 case Jacobian::Line::ShapeD2X1:
888 d *= band.lines[iline].lineshape[pos].D2().X1;
889 band_copy.lines[iline].lineshape[pos].D2().X1 += d;
890 break;
891 case Jacobian::Line::ShapeD2X2:
892 d *= band.lines[iline].lineshape[pos].D2().X2;
893 band_copy.lines[iline].lineshape[pos].D2().X2 += d;
894 break;
895 case Jacobian::Line::ShapeD2X3:
896 d *= band.lines[iline].lineshape[pos].D2().X3;
897 band_copy.lines[iline].lineshape[pos].D2().X3 += d;
898 break;
899 case Jacobian::Line::ShapeFVCX0:
900 d *= band.lines[iline].lineshape[pos].FVC().X0;
901 band_copy.lines[iline].lineshape[pos].FVC().X0 += d;
902 break;
903 case Jacobian::Line::ShapeFVCX1:
904 d *= band.lines[iline].lineshape[pos].FVC().X1;
905 band_copy.lines[iline].lineshape[pos].FVC().X1 += d;
906 break;
907 case Jacobian::Line::ShapeFVCX2:
908 d *= band.lines[iline].lineshape[pos].FVC().X2;
909 band_copy.lines[iline].lineshape[pos].FVC().X2 += d;
910 break;
911 case Jacobian::Line::ShapeFVCX3:
912 d *= band.lines[iline].lineshape[pos].FVC().X3;
913 band_copy.lines[iline].lineshape[pos].FVC().X3 += d;
914 break;
915 case Jacobian::Line::ShapeETAX0:
916 d *= band.lines[iline].lineshape[pos].ETA().X0;
917 band_copy.lines[iline].lineshape[pos].ETA().X0 += d;
918 break;
919 case Jacobian::Line::ShapeETAX1:
920 d *= band.lines[iline].lineshape[pos].ETA().X1;
921 band_copy.lines[iline].lineshape[pos].ETA().X1 += d;
922 break;
923 case Jacobian::Line::ShapeETAX2:
924 d *= band.lines[iline].lineshape[pos].ETA().X2;
925 band_copy.lines[iline].lineshape[pos].ETA().X2 += d;
926 break;
927 case Jacobian::Line::ShapeETAX3:
928 d *= band.lines[iline].lineshape[pos].ETA().X3;
929 band_copy.lines[iline].lineshape[pos].ETA().X3 += d;
930 break;
931 case Jacobian::Line::Center:
932 d *= band.lines[iline].F0;
933 band_copy.lines[iline].F0 += d;
934 break;
935 case Jacobian::Line::Strength:
936 d *= band.lines[iline].I0;
937 band_copy.lines[iline].I0 += d;
938 break;
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: {
970 /* do nothing */
971 }
972 }
973
974 // Perform calculations and estimate derivative
975 const auto [dabs, dwork] = ecs_absorption_impl(T, H, P, this_vmr, vmrs, ecs_data, f_grid, zeeman_polarization, band_copy);
976 vec -= dabs;
977 vec /= -d;
978 work &= dwork;
979 } else if (Quantum::Number::StateMatch(target.qid, band.quantumidentity) == Quantum::Number::StateMatchType::Full) {
980 ErrorCorrectedSuddenData ecs_data_copy = ecs_data;
981
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;
988 break;
989 case Jacobian::Line::ECS_SCALINGX1:
990 d *= ecs_data_copy[spec].scaling.X1;
991 ecs_data_copy[spec].scaling.X1 += d;
992 break;
993 case Jacobian::Line::ECS_SCALINGX2:
994 d *= ecs_data_copy[spec].scaling.X2;
995 ecs_data_copy[spec].scaling.X2 += d;
996 break;
997 case Jacobian::Line::ECS_SCALINGX3:
998 d *= ecs_data_copy[spec].scaling.X3;
999 ecs_data_copy[spec].scaling.X3 += d;
1000 break;
1001 case Jacobian::Line::ECS_BETAX0:
1002 d *= ecs_data_copy[spec].beta.X0;
1003 ecs_data_copy[spec].beta.X0 += d;
1004 break;
1005 case Jacobian::Line::ECS_BETAX1:
1006 d *= ecs_data_copy[spec].beta.X1;
1007 ecs_data_copy[spec].beta.X1 += d;
1008 break;
1009 case Jacobian::Line::ECS_BETAX2:
1010 d *= ecs_data_copy[spec].beta.X2;
1011 ecs_data_copy[spec].beta.X2 += d;
1012 break;
1013 case Jacobian::Line::ECS_BETAX3:
1014 d *= ecs_data_copy[spec].beta.X3;
1015 ecs_data_copy[spec].beta.X3 += d;
1016 break;
1017 case Jacobian::Line::ECS_LAMBDAX0:
1018 d *= ecs_data_copy[spec].lambda.X0;
1019 ecs_data_copy[spec].lambda.X0 += d;
1020 break;
1021 case Jacobian::Line::ECS_LAMBDAX1:
1022 d *= ecs_data_copy[spec].lambda.X1;
1023 ecs_data_copy[spec].lambda.X1 += d;
1024 break;
1025 case Jacobian::Line::ECS_LAMBDAX2:
1026 d *= ecs_data_copy[spec].lambda.X1;
1027 ecs_data_copy[spec].lambda.X1 += d;
1028 break;
1029 case Jacobian::Line::ECS_LAMBDAX3:
1030 d *= ecs_data_copy[spec].lambda.X0;
1031 ecs_data_copy[spec].lambda.X0 += d;
1032 break;
1033 case Jacobian::Line::ECS_DCX0:
1034 d *= ecs_data_copy[spec].collisional_distance.X0;
1035 ecs_data_copy[spec].collisional_distance.X0 += d;
1036 break;
1037 case Jacobian::Line::ECS_DCX1:
1038 d *= ecs_data_copy[spec].collisional_distance.X1;
1039 ecs_data_copy[spec].collisional_distance.X1 += d;
1040 break;
1041 case Jacobian::Line::ECS_DCX2:
1042 d *= ecs_data_copy[spec].collisional_distance.X2;
1043 ecs_data_copy[spec].collisional_distance.X2 += d;
1044 break;
1045 case Jacobian::Line::ECS_DCX3:
1046 d *= ecs_data_copy[spec].collisional_distance.X3;
1047 ecs_data_copy[spec].collisional_distance.X3 += d;
1048 break;
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: {
1090 /* do nothing */
1091 }
1092 }
1093
1094 // Perform calculations and estimate derivative
1095 const auto [dabs, dwork] = ecs_absorption_impl(T, H, P, this_vmr, vmrs, ecs_data_copy, f_grid, zeeman_polarization, band);
1096 vec -= dabs;
1097 vec /= -d;
1098 work &= dwork;
1099 break; // Leave early because it is a band derivative
1100 }
1101 }
1102 } else {
1103 vec *= 0; // No derivative, so don't mess around and remove everything
1104 }
1105 }
1106
1107 return {std::move(absorption), std::move(jacobian), not work};
1108}
1109
1111 AbsorptionLines& band,
1112 const Tensor4& tempdata,
1113 const Vector& temperatures,
1114 const Numeric P0,
1115 const Index ord)
1116{
1117 // Sizes
1118 const Index S = band.NumBroadeners();
1119 const Index N = band.NumLines();
1120
1121 for (Index i=0; i<S; i++) {
1122 // Allocate Vector for fitting the data
1123 Vector targ(N);
1124
1125 // Rosenkranz 1st order parameter
1126 for (Index j=0; j<N; j++) {
1127 targ = tempdata(1, j, i, joker); // Get values
1128
1129 // Assume linear pressure dependency
1130 targ /= P0;
1131
1132 // Fit to a polynomial
1133 auto [fit, c] = Minimize::curve_fit<Minimize::Polynom>(temperatures, targ, LineShape::ModelParameters::N - 1);
1134 if (not fit) return EXIT_FAILURE;
1135 band.lines[j].lineshape[i].Y() = LineShape::ModelParameters(LineShape::TemperatureModel::POLY, c);
1136 }
1137
1138 // Rosenkranz 2nd order frequency parameter
1139 if (ord > 1) {
1140 for (Index j=0; j<N; j++) {
1141 targ = tempdata(2, j, i, joker); // Get values
1142
1143 // Assume squared pressure dependency
1144 targ /= P0 * P0;
1145
1146 // Fit to a polynomial
1147 auto [fit, c] = Minimize::curve_fit<Minimize::Polynom>(temperatures, targ, LineShape::ModelParameters::N - 1);
1148 if (not fit) return EXIT_FAILURE;
1149 band.lines[j].lineshape[i].DV() = LineShape::ModelParameters(LineShape::TemperatureModel::POLY, c);
1150 }
1151 }
1152
1153 // Rosenkranz 2nd order strength parameter
1154 if (ord > 1) {
1155 for (Index j=0; j<N; j++) {
1156 targ = tempdata(0, j, i, joker); // Get values
1157
1158 // Assume squared pressure dependency
1159 targ /= P0 * P0;
1160
1161 // Fit to a polynomial
1162 auto [fit, c] = Minimize::curve_fit<Minimize::Polynom>(temperatures, targ, LineShape::ModelParameters::N - 1);
1163 if (not fit) return EXIT_FAILURE;
1164 band.lines[j].lineshape[i].G() = LineShape::ModelParameters(LineShape::TemperatureModel::POLY, c);
1165 }
1166 }
1167
1168 // 3rd order broadening parameter [[FIXME: UNTESTED]]
1169 if (ord > 2) {
1170 for (Index j=0; j<N; j++) {
1171 targ = tempdata(3, j, i, joker);
1172
1173 // Assume cubic pressure dependency
1174 targ /= P0 * P0 * P0;
1175
1176 // Fit to a polynomial
1177 auto [fit, c] = Minimize::curve_fit<Minimize::Polynom>(temperatures, targ, LineShape::ModelParameters::N - 1);
1178 if (not fit) return EXIT_FAILURE;
1179 ARTS_USER_ERROR("Not yet implemented: 3rd order line mixing")
1180 // band.lines[j].lineshape[i].DG() = LineShape::ModelParameters(LineShape::TemperatureModel::POLY, c);
1181 }
1182 }
1183 }
1184
1185 // If we reach here, we have to set the band population type
1186 // to LTE and renormalize the strength to physical frequency
1187 band.normalization = Absorption::NormalizationType::SFS;
1188 band.population = Absorption::PopulationType::LTE;
1189
1190 return EXIT_SUCCESS;
1191}
1192
1193
1201Vector RosenkranzY(const Vector& dip,
1202 const ConstMatrixView& W,
1203 const AbsorptionLines& band) {
1204 const Index N = dip.nelem();
1205
1206 // Output
1207 Vector Y(N, 0);
1208
1209 // Rosenkranz coefficients
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);
1214 }
1215 }
1216
1217 return Y;
1218}
1219
1220
1228Vector RosenkranzG(const Vector& dip,
1229 const ConstMatrixView& W,
1230 const AbsorptionLines& band) {
1231 const Index N = dip.nelem();
1232
1233 // Output
1234 Vector G(N, 0);
1235
1236 // Rosenkranz coefficients
1237 for (Index k=0; k<N; k++) {
1238 for (Index j=0; j<N; j++) {
1239 if (k == j) continue;
1240 G[k] += W(k, j) * W(j, k) / Math::pow2(band.lines[j].F0 - band.lines[k].F0);
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));
1246
1247 }
1248 }
1249 }
1250
1251 return G;
1252}
1253
1254
1262Vector RosenkranzDV(const Vector& dip,
1263 const ConstMatrixView& W,
1264 const AbsorptionLines& band) {
1265 const Index N = dip.nelem();
1266
1267 // Output
1268 Vector DV(N, 0);
1269
1270 // Rosenkranz coefficients
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);
1275 }
1276 }
1277
1278 return DV;
1279}
1280
1281
1283 const ComplexMatrix& W,
1284 const Vector& pop,
1285 const Vector& dip,
1286 const AbsorptionLines& band,
1287 const Numeric frenorm,
1288 const Numeric T,
1289 const Numeric P,
1290 const Numeric QT,
1291 const Numeric QT0,
1292 const Index broadener) {
1293 // Compute and adapt values for the Voigt adaptation
1294 EquivalentLines eig(W, pop, dip);
1295
1296 // Compute all the shifted line centers in band-order
1297 Vector f0(band.NumLines());
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);
1300 }
1301
1302 // Find how the band-orders have shifted line order at the current temp/pres
1303 ArrayOfIndex sorting(band.NumLines());
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]);
1310 }
1311 }
1312 }
1313
1314 // Sort the eigen-values and strengths in the same order as the band
1315 eig.sort_by_frequency(f0, sorting);
1316
1317 // Eigenvalue should now be (F0 + P * D0(T) + i * P * G0(T)) + (P * P * DV(T) + i * P * P * P * DG(T)),
1318 // or in form (1) + (2). We only want to keep (2), since it is new from the line mixing
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));
1321 }
1322
1323 // The strength term should now be (1 + i y(T) + g(T)) * d**2 * rho(T)
1324 // d**2 * rho(T) * F0 * (1 - exp(-hF0/kT)) should be I0(T)
1325 // We only want to keep (i y(T) + g(T)).
1326 // So we divide by d**2 * rho(T) through the use of I0(T) and remove 1 from the real component
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;
1329 eig.str[i] *= - band.lines[i].F0 * std::expm1(- (Constant::h * band.lines[i].F0) / (Constant::k * T)) / i0;
1330 }
1331 eig.str -= 1;
1332
1333 return eig;
1334}
1335
1336
1352 const Vector& temperatures,
1353 const ErrorCorrectedSuddenData& ecs_data,
1354 const Numeric P) {
1355 const Index N = band.NumLines();
1356 const Index M = band.NumBroadeners();
1357 const Index K = temperatures.nelem();
1358
1359 // Weighted center of the band
1360 const Numeric frenorm = band.F_mean();
1361
1362 // Need sorting to put weak lines last, but we need the sorting constant or the output jumps
1363 const ArrayOfIndex sorting = sorted_population_and_dipole(band.T0, band).first;
1364 const Numeric QT0 = single_partition_function(band.T0, band.Isotopologue());
1365
1366 // Output
1367 Tensor4 out(4, N, M, K);
1368
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];
1373 const Numeric QT = single_partition_function(T, band.Isotopologue());
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[sorting[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 const auto eig = eigenvalue_adaptation_of_relmat(W, pop, dip, band, frenorm, T, P, QT, QT0, m);
1385
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();
1390 }
1391 }
1392
1393 return out;
1394}
1395
1396
1412 const Vector& temperatures,
1413 const ErrorCorrectedSuddenData& ecs_data,
1414 const Numeric P) {
1415 const Index N = band.NumLines();
1416 const Index M = band.NumBroadeners();
1417 const Index K = temperatures.nelem();
1418
1419 // Weighted center of the band
1420 const Numeric frenorm = band.F_mean();
1421
1422 // Need sorting to put weak lines last, but we need the sorting constant or the output jumps
1423 const ArrayOfIndex sorting = sorted_population_and_dipole(band.T0, band).first;
1424
1425 // Output
1426 Tensor4 out(4, N, M, K);
1427
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];
1432
1433 // Relaxation matrix of T0 sorting at T
1434 ComplexMatrix W = single_species_ecs_relaxation_matrix(band, sorting, T, P, ecs_data[band.broadeningspecies[m]], m);
1435 for (Index n=0; n<N; n++) {
1436 W(n, n) += band.lines[n].F0 - frenorm;
1437 }
1438
1439 // Populations and dipoles of T0 sorting at T
1440 const auto [pop, dip] = presorted_population_and_dipole(T, sorting, band);
1441
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;
1446 }
1447 }
1448
1449 return out;
1450}
1451
1453 const Vector& temperatures,
1454 const ErrorCorrectedSuddenData& ecs_data,
1455 const Numeric P0,
1456 const Index ord,
1457 const bool robust,
1458 const bool rosenkranz_adaptation,
1459 const Verbosity& verbosity) {
1461 ARTS_USER_ERROR_IF(P0 <= 0, P0, " Pa is not possible")
1462
1463 ARTS_USER_ERROR_IF(not is_sorted(temperatures),
1464 "The temperature list [",
1465 temperatures,
1466 "] K\n"
1467 "must be fully sorted from low to high")
1468
1470 ord < 1 or ord > 3, "Order not in list [1, 2, 3], is: ", ord)
1471
1473 band,
1474 rosenkranz_adaptation not_eq 0
1475 ? rosenkranz_approximation(band, temperatures, ecs_data, P0)
1476 : ecs_eigenvalue_approximation(band, temperatures, ecs_data, P0),
1477 temperatures,
1478 P0,
1479 ord)) {
1480 ARTS_USER_ERROR_IF(not robust,
1481 "Bad eigenvalue adaptation for band: ",
1482 band.quantumidentity,
1483 '\n')
1484 out3 << "Bad eigenvalue adaptation for band: " << band.quantumidentity
1485 << '\n';
1486
1487 band.normalization = Absorption::NormalizationType::SFS;
1488 band.population = Absorption::PopulationType::LTE;
1489 for (auto& line : band.lines) {
1490 for (auto& sm : line.lineshape.Data()) {
1491 sm.Y() = LineShape::ModelParameters(LineShape::TemperatureModel::None);
1492 sm.G() = LineShape::ModelParameters(LineShape::TemperatureModel::None);
1493 sm.DV() = LineShape::ModelParameters(LineShape::TemperatureModel::None);
1494 }
1495 }
1496 }
1497}
1498
1499
1501 const Vector& temperatures,
1502 const ErrorCorrectedSuddenData& ecs_data,
1503 const Vector& pressures) {
1504 const Index N = band.NumLines();
1505 const Index M = band.NumBroadeners();
1506 const Index K = temperatures.nelem();
1507 const Index L = pressures.size();
1508
1509 Tensor5 out(4, N, M, K, L);
1510 for (Index l=0; l<L; l++) {
1511 out(joker, joker, joker, joker, l) = ecs_eigenvalue_approximation(band, temperatures, ecs_data, pressures[l]);
1512 }
1513 return out;
1514}
1515
1516std::ostream& operator<<(std::ostream& os,
1517 const ErrorCorrectedSuddenData& rbd) {
1518 for (Index i = 0; i < rbd.data.nelem(); i++) {
1519 if (i) os << '\n';
1520 os << rbd.data[i];
1521 }
1522 return os;
1523}
1524
1525std::istream& operator>>(std::istream& is, ErrorCorrectedSuddenData& rbd) {
1526 for (auto& x : rbd.data) is >> x;
1527 return is;
1528}
1529
1530Index ErrorCorrectedSuddenData::pos(Species::Species spec) const {
1531 return std::distance(data.begin(),
1532 std::find(data.cbegin(), data.cend(), spec));
1533}
1534
1536 Species::Species spec) const {
1537 if (auto ptr = std::find(data.cbegin(), data.cend(), spec);
1538 ptr not_eq data.cend())
1539 return *ptr;
1540 return *std::find(data.cbegin(), data.cend(), Species::Species::Bath);
1541}
1542
1544 Species::Species spec) {
1545 if (auto ptr = std::find(data.begin(), data.end(), spec);
1546 ptr not_eq data.end())
1547 return *ptr;
1548 return data.emplace_back(spec);
1549}
1550
1552 const QuantumIdentifier& id) {
1553 if (auto ptr = std::find(begin(), end(), id); ptr not_eq end()) return *ptr;
1554 return emplace_back(id);
1555}
1556
1558 const QuantumIdentifier& id) const {
1559 if (auto ptr = std::find(cbegin(), cend(), id); ptr not_eq cend()) {
1560 return *ptr;
1561 }
1562 ARTS_USER_ERROR("Cannot find data for QuantumIdentifier:\n",
1563 id,
1564 '\n',
1565 "Available data:\n",
1566 *this);
1567 return front(); // To get rid of potential warnings...
1568}
1569
1570std::ostream& operator<<(std::ostream& os,
1572 std::for_each(m.cbegin(), m.cend(), [&](auto& x) { os << x << '\n'; });
1573 return os;
1574}
1575} // namespace Absorption::LineMixing
Common ARTS conversions.
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
Input manipulator class for doubles to enable nan and inf parsing.
Helper macros for debugging.
#define ARTS_ASSERT(condition,...)
Definition debug.h:86
#define ARTS_USER_ERROR(...)
Definition debug.h:153
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
constexpr bool good_enum(EnumType x) noexcept
Checks if the enum number is good.
Definition enums.h:21
#define WIGNER6
Definition linemixing.cc:29
#define WIGNER3
Definition linemixing.cc:28
Numeric wig6(const Rational &a, const Rational &b, const Rational &c, const Rational &d, const Rational &e, const Rational &f) noexcept
Definition linemixing.cc:42
Numeric wig3(const Rational &a, const Rational &b, const Rational &c, const Rational &d, const Rational &e, const Rational &f) noexcept
Definition linemixing.cc:32
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.
Definition linescaling.cc:4
Numeric fac(const Index n)
fac
Definition math_funcs.cc:46
Declarations having to do with the four output streams.
#define CREATE_OUT3
Definition messages.h:189
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
Definition linemixing.h:70
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.
Definition zeemandata.h:27
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)
Definition linemixing.h:18
void sort_by_frequency(Vector &f, const ArrayOfIndex &sorting)
Rovibrational line mixing data following the ideas of Collisional Effects On Molecular Spectra by J.
Definition linemixing.h:118
ArrayOfSpeciesErrorCorrectedSuddenData data
Data of species data The program is considered ill-formed if data does not contain a Bath-species,...
Definition linemixing.h:124
Index pos(Species::Species spec) const
const SpeciesErrorCorrectedSuddenData & operator[](Species::Species spec) const
ErrorCorrectedSuddenData & operator[](const QuantumIdentifier &id)
Contains the population distribution and dipole.
Definition linemixing.h:44
ArrayOfIndex sort(const AbsorptionLines &band) noexcept
PopulationAndDipole(const Numeric T, const AbsorptionLines &band)
Numeric Omega(const Numeric T, const Numeric T0, const Numeric other_mass, const Numeric energy_x, const Numeric energy_xm2) const
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
static constexpr Index N
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.
#define d
#define w
#define a
#define c
#define b
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