ARTS 2.5.11 (git: 725533f0)
linemixing_hitran.cc
Go to the documentation of this file.
1
9#include "linemixing_hitran.h"
10#include <Faddeeva/Faddeeva.hh>
11#include <cinttypes>
12#include <fstream>
13#include <functional>
14#include <numeric>
15
16#include "lin_alg.h"
17#include "linemixing.h"
18#include "lineshapemodel.h"
19#include "matpack_data.h"
20#include "matpack_math.h"
21#include "physics_funcs.h"
22
23#pragma GCC diagnostic push
24#pragma GCC diagnostic ignored "-Wshadow"
25#pragma GCC diagnostic ignored "-Wconversion"
26#pragma GCC diagnostic ignored "-Wfloat-conversion"
27#if defined(__clang__)
28#pragma GCC diagnostic ignored "-Wdeprecated-copy-with-dtor"
29#pragma GCC diagnostic ignored "-Wdeprecated-copy-with-user-provided-dtor"
30#else
31#pragma GCC diagnostic ignored "-Wclass-memaccess"
32#endif
33
34#include <Eigen/Dense>
35
36#pragma GCC diagnostic pop
37
38namespace lm_hitran_2017 {
39namespace parameters {
40 static constexpr Index nBmx=7'000; // Max Number of Bands
41 static constexpr Index nLmx=700; // Max Number of Lines per Band
42 static constexpr Index Nlifmax=10; // Max number of l values
43 static constexpr Index Jmax=131; // Max number of j values
44
45 static constexpr Numeric Ct=1.4387686e0; // Constant
46 static constexpr Numeric T0=296; // Constant
47 static constexpr Numeric CtGamD=1.1325e-08; // Constant
48 static constexpr Numeric aMolAtm=7.33889e+21; // Constant
49
50 static constexpr auto aMass = stdarrayify(44.e-3,45.e-3,46.e-3,45.e-3,47.e-3,46.e-3,48.e-3,47.e-3,46.e-3,49.e-3); // Constant
51
52} // namespace parameters
53
54struct CommonBlock {
55struct Bands {
56 Index nBand;
57 std::array<Index, parameters::nBmx> Isot;
58 std::array<Index, parameters::nBmx> nLines;
59 std::array<Index, parameters::nBmx> li;
60 std::array<Index, parameters::nBmx> lf;
61 std::array<String, parameters::nBmx> BandFile;
62} Bands;
63
64struct LineSg {
65 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> Sig;
66 LineSg() : Sig(parameters::nLmx, parameters::nBmx) {};
67} LineSg;
68
69struct DipoRigid {
70 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> Dipo0;
71 DipoRigid() : Dipo0(parameters::nLmx, parameters::nBmx) {};
72} DipoRigid;
73
74struct Energy {
75 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> E;
76 Energy() : E(parameters::nLmx, parameters::nBmx) {};
77} Energy;
78
79struct GamVT0AIR {
80 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWVT0AIR;
81 GamVT0AIR() : HWVT0AIR(parameters::nLmx, parameters::nBmx) {};
82} GamVT0AIR;
83
84struct GamSDVT0AIR {
85 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWSDVT0AIR;
86 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> rHWT0AIR;
87 GamSDVT0AIR() : HWSDVT0AIR(parameters::nLmx, parameters::nBmx), rHWT0AIR(parameters::nLmx, parameters::nBmx) {};
88} GamSDVT0AIR;
89
90struct DTGAMAIR {
91 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> BHWAIR;
92 DTGAMAIR() : BHWAIR(parameters::nLmx, parameters::nBmx) {};
93} DTGAMAIR;
94
95struct GamVT0CO2 {
96 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWVT0SELF;
97 GamVT0CO2() : HWVT0SELF(parameters::nLmx, parameters::nBmx) {};
98} GamVT0CO2;
99
100struct GamSDVT0CO2 {
101 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWSDVT0SELF;
102 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> rHWT0SELF;
103 GamSDVT0CO2() : HWSDVT0SELF(parameters::nLmx, parameters::nBmx), rHWT0SELF(parameters::nLmx, parameters::nBmx) {};
104} GamSDVT0CO2;
105
106struct DTGAMCO2 {
107 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> BHWSELF;
108 DTGAMCO2() : BHWSELF(parameters::nLmx, parameters::nBmx) {};
109} DTGAMCO2;
110
111struct GamVT0H2O {
112 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWVT0H2O;
113 GamVT0H2O() : HWVT0H2O(parameters::nLmx, parameters::nBmx) {};
114} GamVT0H2O;
115
116struct GamSDVT0H2O {
117 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWSDVT0H2O;
118 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> rHWT0H2O;
119 GamSDVT0H2O() : HWSDVT0H2O(parameters::nLmx, parameters::nBmx), rHWT0H2O(parameters::nLmx, parameters::nBmx) {};
120} GamSDVT0H2O;
121
122struct DTGAMH2O {
123 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> BHWH2O;
124 DTGAMH2O() : BHWH2O(parameters::nLmx, parameters::nBmx) {};
125} DTGAMH2O;
126
127struct GamT {
128 std::array<Numeric, parameters::nLmx> HWT;
129 std::array<Numeric, parameters::nLmx> HWSDV2T;
130} GamT;
131
132struct SHIFT {
133 std::array<Numeric, parameters::nLmx> shft;
134} SHIFT;
135
136struct SHIFT0 {
137 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> shft0;
138 SHIFT0() : shft0(parameters::nLmx, parameters::nBmx) {};
139} SHIFT0;
140
141struct PopuT {
142 std::array<Numeric, parameters::nLmx> PopuT;
143} PopuT;
144
145struct PopTrf {
146 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> PopuT0;
147 PopTrf() : PopuT0(parameters::nLmx, parameters::nBmx) {};
148} PopTrf;
149
150struct DipoTcm {
151 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> DipoT;
152 DipoTcm() : DipoT(parameters::nLmx, parameters::nBmx) {};
153} DipoTcm;
154
155struct Jiln {
156 Eigen::Matrix<Index, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> Ji;
157 Jiln() : Ji(parameters::nLmx, parameters::nBmx) {};
158} Jiln;
159
160struct Jfln {
161 Eigen::Matrix<Index, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> Jf;
162 Jfln() : Jf(parameters::nLmx, parameters::nBmx) {};
163} Jfln;
164
165struct Zss {
166 std::array<Complex, parameters::nLmx> ZS;
167} Zss;
168
169struct Zaa {
170 ComplexVector ZA;
171 Zaa() noexcept : ZA(parameters::nLmx) {}
172} Zaa;
173
174struct Wmatrix {
175 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> W;
176 Wmatrix() : W(parameters::nLmx, parameters::nLmx) {};
177} Wmatrix;
178
179struct Wfittedp {
180 Tensor4 W0pp;
181 Tensor4 W0pq;
182 Tensor4 W0pr;
183 Wfittedp() :
184 W0pp(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
185 W0pq(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
186 W0pr(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax) {}
187} Wfittedp;
188
189struct Wfittedq {
190 Tensor4 W0qp;
191 Tensor4 W0qq;
192 Tensor4 W0qr;
193 Wfittedq() :
194 W0qp(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
195 W0qq(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
196 W0qr(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax) {}
197} Wfittedq;
198
199struct Wfittedr {
200 Tensor4 W0rp;
201 Tensor4 W0rq;
202 Tensor4 W0rr;
203 Wfittedr() :
204 W0rp(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
205 W0rq(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
206 W0rr(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax) {}
207} Wfittedr;
208
209struct Bfittedp {
210 Tensor4 B0pp;
211 Tensor4 B0pq;
212 Tensor4 B0pr;
213 Bfittedp() :
214 B0pp(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
215 B0pq(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
216 B0pr(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax) {}
217} Bfittedp;
218
219struct Bfittedq {
220 Tensor4 B0qp;
221 Tensor4 B0qq;
222 Tensor4 B0qr;
223 Bfittedq() :
224 B0qp(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
225 B0qq(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
226 B0qr(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax) {}
227} Bfittedq;
228
229struct Bfittedr {
230 Tensor4 B0rp;
231 Tensor4 B0rq;
232 Tensor4 B0rr;
233 Bfittedr() :
234 B0rp(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
235 B0rq(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
236 B0rr(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax) {}
237} Bfittedr;
238
239struct DiagnR {
240 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> OpR;
241 DiagnR() : OpR(parameters::nLmx, parameters::nLmx) {};
242} DiagnR;
243
244struct DiagnI {
245 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> OpI;
246 DiagnI() : OpI(parameters::nLmx, parameters::nLmx) {};
247} DiagnI;
248
249struct YLT {
250 std::array<Numeric, parameters::nLmx> YT;
251} YLT;
252
253struct UnusedBandParams {
254 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> intens;
255 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> eina;
256 std::array<Rational, parameters::nBmx> iv1;
257 std::array<Rational, parameters::nBmx> iv2;
258 std::array<Rational, parameters::nBmx> il2;
259 std::array<Rational, parameters::nBmx> iv3;
260 std::array<Rational, parameters::nBmx> ir;
261 std::array<Rational, parameters::nBmx> fv1;
262 std::array<Rational, parameters::nBmx> fv2;
263 std::array<Rational, parameters::nBmx> fl2;
264 std::array<Rational, parameters::nBmx> fv3;
265 std::array<Rational, parameters::nBmx> fr;
266 UnusedBandParams() :
267 intens(parameters::nLmx, parameters::nBmx),
268 eina(parameters::nLmx, parameters::nBmx) {};
269} UnusedBandParams;
270}; // CommonBlock
271
272Rational toRationalSum(char a, char b=' ')
273{
274 if (b == ' ' and a == ' ')
275 return Rational();
276 if (b == ' ')
277 return Rational(a-'0');
278 if (a == ' ')
279 return Rational(b-'0');
280 return Rational(10*(a-'0') + b-'0');
281}
282
283void readlines(CommonBlock& cmn, const String& basedir="data_new/")
284{
285 ARTS_USER_ERROR_IF (cmn.Bands.nBand > parameters::nBmx,
286 "Too many bands");
287
288 for (Index iband=0; iband<cmn.Bands.nBand; iband++) {
289 std::ifstream fortranfile;
290 const String fname = basedir + String("/") + cmn.Bands.BandFile[iband] + String(".dat");
291 fortranfile.open(fname.c_str());
292
293 if (fortranfile.is_open()) {
294 String line;
295 getline(fortranfile, line);
296
297 Index nliner=0;
298 while (fortranfile.good()) {
299 ARTS_USER_ERROR_IF (nliner >= parameters::nLmx,
300 "Too many lines");
301
302 char tpline, x;
303 char sDipoRigid[21], sPopTrf[21];
304 char iv11, iv12, iv21, iv22, il21, il22, iv31, iv32, ir1, fv11, fv12, fv21, fv22, fl21, fl22, fv31, fv32, fr1;
305
306 sscanf(line.c_str(),
307 "%c%c" "%1" PRId64 "%12lf" "%10lf"
308 "%10lf" "%5lf" "%5lf" "%4lf"
309 "%5lf" "%5lf" "%4lf" "%10lf"
310 "%4lf" "%4lf" "%8lf"
311 "%c%c%c%c%c%c"
312 "%c%c" "%c%c" "%c%c" "%c%c" "%c"
313 "%c%c%c%c%c%c"
314 "%c%c" "%c%c" "%c%c" "%c%c" "%c"
315 "%c%c%c%c%c%c%c%c%c%c%c%c%c%c%c%c%c%c%c%c"
316 "%c" "%3" PRId64 "%c%c%c%c%c%c%c%c%c%c%c%c%c%c%c%c%c%c%c%c%c"
317 "%5lf" "%5lf" "%4lf" "%5lf"
318 "%20s" "%20s",
319 &x,&x,
320 &cmn.Bands.Isot[iband],
321 &cmn.LineSg.Sig(nliner, iband),
322 &cmn.UnusedBandParams.intens(nliner, iband),
323 &cmn.UnusedBandParams.eina(nliner, iband),
324 &cmn.GamVT0AIR.HWVT0AIR(nliner, iband),
325 &cmn.GamSDVT0AIR.HWSDVT0AIR(nliner, iband),
326 &cmn.GamSDVT0AIR.rHWT0AIR(nliner, iband),
327 &cmn.GamVT0CO2.HWVT0SELF(nliner, iband),
328 &cmn.GamSDVT0CO2.HWSDVT0SELF(nliner, iband),
329 &cmn.GamSDVT0CO2.rHWT0SELF(nliner, iband),
330 &cmn.Energy.E(nliner, iband),
331 &cmn.DTGAMAIR.BHWAIR(nliner, iband),
332 &cmn.DTGAMCO2.BHWSELF(nliner, iband),
333 &cmn.SHIFT0.shft0(nliner, iband),
334 &x,&x,&x,&x,&x,&x,
335 &iv11, &iv12, &iv21, &iv22, &il21, &il22, &iv31, &iv32, &ir1, &fv32, &fr1,
336 &x,&x,&x,&x,&x,&x,
337 &fv11, &fv12, &fv21, &fv22, &fl21, &fl22, &fv31,
338 &x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,
339 &tpline,
340 &cmn.Jiln.Ji(nliner, iband),
341 &x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,&x,
342 &cmn.GamVT0H2O.HWVT0H2O(nliner, iband),
343 &cmn.GamSDVT0H2O.HWSDVT0H2O(nliner, iband),
344 &cmn.GamSDVT0H2O.rHWT0H2O(nliner, iband),
345 &cmn.DTGAMH2O.BHWH2O(nliner, iband),
346 sDipoRigid,
347 sPopTrf);
348 getline(fortranfile, line);
349
350 if (nliner==0) {
351 cmn.UnusedBandParams.iv1[iband]=Rational(toRationalSum(iv11, iv12));
352 cmn.UnusedBandParams.iv2[iband]=Rational(toRationalSum(iv21, iv22));
353 cmn.UnusedBandParams.il2[iband]=Rational(toRationalSum(il21, il22));
354 cmn.UnusedBandParams.iv3[iband]=Rational(toRationalSum(iv31, iv32));
355 cmn.UnusedBandParams.ir[iband]=Rational(toRationalSum(ir1));
356 cmn.UnusedBandParams.fv1[iband]=Rational(toRationalSum(fv11, fv12));
357 cmn.UnusedBandParams.fv2[iband]=Rational(toRationalSum(fv21, fv22));
358 cmn.UnusedBandParams.fl2[iband]=Rational(toRationalSum(fl21, fl22));
359 cmn.UnusedBandParams.fv3[iband]=Rational(toRationalSum(fv31, fv32));
360 cmn.UnusedBandParams.fr[iband]=Rational(toRationalSum(fr1));
361 }
362
363 ARTS_USER_ERROR_IF(not (
364 cmn.UnusedBandParams.iv1[iband]==Rational(toRationalSum(iv11, iv12)) and
365 cmn.UnusedBandParams.iv2[iband]==Rational(toRationalSum(iv21, iv22)) and
366 cmn.UnusedBandParams.il2[iband]==Rational(toRationalSum(il21, il22)) and
367 cmn.UnusedBandParams.iv3[iband]==Rational(toRationalSum(iv31, iv32)) and
368 cmn.UnusedBandParams.ir[iband]==Rational(toRationalSum(ir1)) and
369 cmn.UnusedBandParams.fv1[iband]==Rational(toRationalSum(fv11, fv12)) and
370 cmn.UnusedBandParams.fv2[iband]==Rational(toRationalSum(fv21, fv22)) and
371 cmn.UnusedBandParams.fl2[iband]==Rational(toRationalSum(fl21, fl22)) and
372 cmn.UnusedBandParams.fv3[iband]==Rational(toRationalSum(fv31, fv32)) and
373 cmn.UnusedBandParams.fr[iband]==Rational(toRationalSum(fr1))),
374 "Bad read, bands do not have the same global quantum numbers...");
375
376 // Fix...
377 String ssDipoRigid = sDipoRigid;
378 String ssPopTrf = sPopTrf;
379 std::replace(ssDipoRigid.begin(), ssDipoRigid.end(), 'D', 'E');
380 std::replace(ssPopTrf.begin(), ssPopTrf.end(), 'D', 'E');
381 cmn.DipoRigid.Dipo0(nliner, iband) = std::stod(ssDipoRigid);
382 cmn.PopTrf.PopuT0(nliner, iband) = std::stod(ssPopTrf);
383
384 // Dipole at temperature
385 cmn.DipoTcm.DipoT(nliner, iband) = std::sqrt(cmn.UnusedBandParams.intens(nliner, iband)/(cmn.PopTrf.PopuT0(nliner,iband) * cmn.LineSg.Sig(nliner,iband) * (1-std::exp(-1.4388*cmn.LineSg.Sig(nliner,iband)/296))));
386
387 // Fix Js
388 if (tpline == 'P')
389 cmn.Jfln.Jf(nliner, iband) = cmn.Jiln.Ji(nliner, iband) - 1;
390 else if (tpline == 'Q')
391 cmn.Jfln.Jf(nliner, iband) = cmn.Jiln.Ji(nliner, iband);
392 else
393 cmn.Jfln.Jf(nliner, iband) = cmn.Jiln.Ji(nliner, iband) + 1;
394
395 // Fix isotologue
396 if (cmn.Bands.Isot[iband] == 0)
397 cmn.Bands.Isot[iband] = 10;
398
399 nliner++;
400 }
401 cmn.Bands.nLines[iband] = nliner;
402 }
403
404 fortranfile.close();
405 }
406}
407
408template<size_t NT>
409Numeric atob(const Numeric& aa,
410 const std::array<Numeric, NT>& a,
411 const std::array<Numeric, NT>& b)
412{
413 for (size_t i=1; i<NT; i++) {
414 if (a[i] >= aa) {
415 if (i < 2 or i == NT) {
416 size_t j=i;
417
418 if (i < 2)
419 j = 2;
420 if (i == NT-1)
421 j = NT-1;
422
423 Numeric A0D1=a[j-2]-a[j-1]; if (A0D1 == 0) A0D1=0.0001;
424 Numeric A0D2=a[j-2]-a[j]; if (A0D2 == 0) A0D2=0.0001;
425 Numeric A1D1=a[j-1]-a[j-2]; if (A1D1 == 0) A1D1=0.0001;
426 Numeric A1D2=a[j-1]-a[j]; if (A1D2 == 0) A1D2=0.0001;
427 Numeric A2D1=a[j]-a[j-2]; if (A2D1 == 0) A2D1=0.0001;
428 Numeric A2D2=a[j]-a[j-1]; if (A2D2 == 0) A2D2=0.0001;
429
430 const Numeric A0=(aa-a[j-1])*(aa-a[j])/(A0D1*A0D2);
431 const Numeric A1=(aa-a[j-2])*(aa-a[j])/(A1D1*A1D2);
432 const Numeric A2=(aa-a[j-2])*(aa-a[j-1])/(A2D1*A2D2);
433
434 return A0*b[j-2] + A1*b[j-1] + A2*b[j];
435 }
436 size_t j = i;
437
438 Numeric A0D1=a[j-2]-a[j-1]; if (A0D1 == 0) A0D1=0.0001;
439 Numeric A0D2=a[j-2]-a[j]; if (A0D2 == 0) A0D2=0.0001;
440 Numeric A0D3=a[j-2]-a[j+1]; if (A0D3 == 0) A0D3=0.0001;
441 Numeric A1D1=a[j-1]-a[j-2]; if (A1D1 == 0) A1D1=0.0001;
442 Numeric A1D2=a[j-1]-a[j]; if (A1D2 == 0) A1D2=0.0001;
443 Numeric A1D3=a[j-1]-a[j+1]; if (A1D3 == 0) A1D3=0.0001;
444 Numeric A2D1=a[j]-a[j-2]; if (A2D1 == 0) A2D1=0.0001;
445 Numeric A2D2=a[j]-a[j-1]; if (A2D2 == 0) A2D2=0.0001;
446 Numeric A2D3=a[j]-a[j+1]; if (A2D3 == 0) A2D3=0.0001;
447 Numeric A3D1=a[j+1]-a[j-2]; if (A3D1 == 0) A3D1=0.0001;
448 Numeric A3D2=a[j+1]-a[j-1]; if (A3D2 == 0) A3D2=0.0001;
449 Numeric A3D3=a[j+1]-a[j]; if (A3D3 == 0) A3D3=0.0001;
450
451
452 Numeric A0=(aa-a[j-1])*(aa-a[j])*(aa-a[j+1]);
453 A0=A0/(A0D1*A0D2*A0D3);
454 Numeric A1=(aa-a[j-2])*(aa-a[j])*(aa-a[j+1]);
455 A1=A1/(A1D1*A1D2*A1D3);
456 Numeric A2=(aa-a[j-2])*(aa-a[j-1])*(aa-a[j+1]);
457 A2=A2/(A2D1*A2D2*A2D3);
458 Numeric A3=(aa-a[j-2])*(aa-a[j-1])*(aa-a[j]);
459 A3=A3/(A3D1*A3D2*A3D3);
460
461 return A0*b[j-2] + A1*b[j-1] + A2*b[j] + A3*b[j+1];
462 }
463 }
464
465 return std::numeric_limits<Numeric>::quiet_NaN();
466}
467
468void qt_co2(const Numeric& t,
469 const Index& iso,
470 Numeric& gsi,
471 Numeric& qt)
472{
473 constexpr auto xgj = stdarrayify(1.,2.,1.,6.,2.,12.,1.,6.,1.,2.,12.);
474 constexpr auto tdat = stdarrayify(60., 85., 110., 135., 160., 185., 210., 235.,
475 260., 285., 310., 335., 360., 385., 410., 435., 460., 485.,
476 510., 535., 560., 585., 610., 635., 660., 685., 710., 735.,
477 760., 785., 810., 835., 860., 885., 910., 935., 960., 985.,
478 1010.,1035.,1060.,1085.,1110.,1135.,1160.,1185.,1210.,1235.,
479 1260.,1285.,1310.,1335.,1360.,1385.,1410.,1435.,1460.,1485.,
480 1510.,1535.,1560.,1585.,1610.,1635.,1660.,1685.,1710.,1735.,
481 1760.,1785.,1810.,1835.,1860.,1885.,1910.,1935.,1960.,1985.,
482 2010.,2035.,2060.,2085.,2110.,2135.,2160.,2185.,2210.,2235.,
483 2260.,2285.,2310.,2335.,2360.,2385.,2410.,2435.,2460.,2485.,
484 2510.,2535.,2560.,2585.,2610.,2635.,2660.,2685.,2710.,2735.,
485 2760.,2785.,2810.,2835.,2860.,2885.,2910.,2935.,2960.,2985.,
486 3010.);
487 constexpr auto qoft = stdarrayify(
488// c... -- 626
489 stdarrayify( 0.53642E+02, 0.75947E+02, 0.98292E+02,
490 + 0.12078E+03, 0.14364E+03, 0.16714E+03, 0.19160E+03, 0.21731E+03,
491 + 0.24454E+03, 0.27355E+03, 0.30456E+03, 0.33778E+03, 0.37343E+03,
492 + 0.41170E+03, 0.45280E+03, 0.49692E+03, 0.54427E+03, 0.59505E+03,
493 + 0.64948E+03, 0.70779E+03, 0.77019E+03, 0.83693E+03, 0.90825E+03,
494 + 0.98440E+03, 0.10656E+04, 0.11522E+04, 0.12445E+04, 0.13427E+04,
495 + 0.14471E+04, 0.15580E+04, 0.16759E+04, 0.18009E+04, 0.19334E+04,
496 + 0.20739E+04, 0.22225E+04, 0.23798E+04, 0.25462E+04, 0.27219E+04,
497 + 0.29074E+04, 0.31032E+04, 0.33097E+04, 0.35272E+04, 0.37564E+04,
498 + 0.39976E+04, 0.42514E+04, 0.45181E+04, 0.47985E+04, 0.50929E+04,
499 + 0.54019E+04, 0.57260E+04, 0.60659E+04, 0.64221E+04, 0.67952E+04,
500 + 0.71859E+04, 0.75946E+04, 0.80222E+04, 0.84691E+04, 0.89362E+04,
501 + 0.94241E+04, 0.99335E+04, 0.10465E+05, 0.11020E+05, 0.11598E+05,
502 + 0.12201E+05, 0.12828E+05, 0.13482E+05, 0.14163E+05, 0.14872E+05,
503 + 0.15609E+05, 0.16376E+05, 0.17173E+05, 0.18001E+05, 0.18861E+05,
504 + 0.19754E+05, 0.20682E+05, 0.21644E+05, 0.22643E+05, 0.23678E+05,
505 + 0.24752E+05, 0.25865E+05, 0.27018E+05, 0.28212E+05, 0.29449E+05,
506 + 0.30730E+05, 0.32055E+05, 0.33426E+05, 0.34845E+05, 0.36312E+05,
507 + 0.37828E+05, 0.39395E+05, 0.41015E+05, 0.42688E+05, 0.44416E+05,
508 + 0.46199E+05, 0.48041E+05, 0.49942E+05, 0.51902E+05, 0.53925E+05,
509 + 0.56011E+05, 0.58162E+05, 0.60379E+05, 0.62664E+05, 0.65019E+05,
510 + 0.67444E+05, 0.69942E+05, 0.72515E+05, 0.75163E+05, 0.77890E+05,
511 + 0.80695E+05, 0.83582E+05, 0.86551E+05, 0.89605E+05, 0.92746E+05,
512 + 0.95975E+05, 0.99294E+05, 0.10271E+06, 0.10621E+06, 0.10981E+06,
513 + 0.11351E+06),
514// c... -- 636
515 stdarrayify( 0.10728E+03, 0.15189E+03, 0.19659E+03,
516 + 0.24164E+03, 0.28753E+03, 0.33486E+03, 0.38429E+03, 0.43643E+03,
517 + 0.49184E+03, 0.55104E+03, 0.61449E+03, 0.68263E+03, 0.75589E+03,
518 + 0.83468E+03, 0.91943E+03, 0.10106E+04, 0.11085E+04, 0.12137E+04,
519 + 0.13266E+04, 0.14477E+04, 0.15774E+04, 0.17163E+04, 0.18649E+04,
520 + 0.20237E+04, 0.21933E+04, 0.23743E+04, 0.25673E+04, 0.27729E+04,
521 + 0.29917E+04, 0.32245E+04, 0.34718E+04, 0.37345E+04, 0.40132E+04,
522 + 0.43087E+04, 0.46218E+04, 0.49533E+04, 0.53041E+04, 0.56749E+04,
523 + 0.60668E+04, 0.64805E+04, 0.69171E+04, 0.73774E+04, 0.78626E+04,
524 + 0.83736E+04, 0.89114E+04, 0.94772E+04, 0.10072E+05, 0.10697E+05,
525 + 0.11353E+05, 0.12042E+05, 0.12765E+05, 0.13523E+05, 0.14317E+05,
526 + 0.15148E+05, 0.16019E+05, 0.16930E+05, 0.17883E+05, 0.18879E+05,
527 + 0.19920E+05, 0.21008E+05, 0.22143E+05, 0.23328E+05, 0.24563E+05,
528 + 0.25852E+05, 0.27195E+05, 0.28594E+05, 0.30051E+05, 0.31568E+05,
529 + 0.33146E+05, 0.34788E+05, 0.36496E+05, 0.38271E+05, 0.40115E+05,
530 + 0.42031E+05, 0.44021E+05, 0.46086E+05, 0.48230E+05, 0.50453E+05,
531 + 0.52759E+05, 0.55150E+05, 0.57628E+05, 0.60195E+05, 0.62854E+05,
532 + 0.65608E+05, 0.68459E+05, 0.71409E+05, 0.74461E+05, 0.77618E+05,
533 + 0.80883E+05, 0.84258E+05, 0.87746E+05, 0.91350E+05, 0.95073E+05,
534 + 0.98918E+05, 0.10289E+06, 0.10698E+06, 0.11121E+06, 0.11558E+06,
535 + 0.12008E+06, 0.12472E+06, 0.12950E+06, 0.13443E+06, 0.13952E+06,
536 + 0.14475E+06, 0.15015E+06, 0.15571E+06, 0.16143E+06, 0.16732E+06,
537 + 0.17338E+06, 0.17962E+06, 0.18604E+06, 0.19264E+06, 0.19943E+06,
538 + 0.20642E+06, 0.21360E+06, 0.22098E+06, 0.22856E+06, 0.23636E+06,
539 + 0.24436E+06),
540// c... -- 628
541 stdarrayify( 0.11368E+03, 0.16096E+03, 0.20833E+03,
542 + 0.25603E+03, 0.30452E+03, 0.35442E+03, 0.40640E+03, 0.46110E+03,
543 + 0.51910E+03, 0.58093E+03, 0.64709E+03, 0.71804E+03, 0.79422E+03,
544 + 0.87607E+03, 0.96402E+03, 0.10585E+04, 0.11600E+04, 0.12689E+04,
545 + 0.13857E+04, 0.15108E+04, 0.16449E+04, 0.17883E+04, 0.19416E+04,
546 + 0.21054E+04, 0.22803E+04, 0.24668E+04, 0.26655E+04, 0.28770E+04,
547 + 0.31021E+04, 0.33414E+04, 0.35956E+04, 0.38654E+04, 0.41516E+04,
548 + 0.44549E+04, 0.47761E+04, 0.51160E+04, 0.54755E+04, 0.58555E+04,
549 + 0.62568E+04, 0.66804E+04, 0.71273E+04, 0.75982E+04, 0.80944E+04,
550 + 0.86169E+04, 0.91666E+04, 0.97446E+04, 0.10352E+05, 0.10990E+05,
551 + 0.11660E+05, 0.12363E+05, 0.13101E+05, 0.13874E+05, 0.14683E+05,
552 + 0.15531E+05, 0.16418E+05, 0.17347E+05, 0.18317E+05, 0.19332E+05,
553 + 0.20392E+05, 0.21499E+05, 0.22654E+05, 0.23859E+05, 0.25116E+05,
554 + 0.26426E+05, 0.27792E+05, 0.29214E+05, 0.30695E+05, 0.32236E+05,
555 + 0.33840E+05, 0.35508E+05, 0.37242E+05, 0.39045E+05, 0.40917E+05,
556 + 0.42862E+05, 0.44881E+05, 0.46977E+05, 0.49152E+05, 0.51407E+05,
557 + 0.53746E+05, 0.56171E+05, 0.58683E+05, 0.61286E+05, 0.63981E+05,
558 + 0.66772E+05, 0.69661E+05, 0.72650E+05, 0.75742E+05, 0.78940E+05,
559 + 0.82246E+05, 0.85664E+05, 0.89196E+05, 0.92845E+05, 0.96613E+05,
560 + 0.10050E+06, 0.10452E+06, 0.10867E+06, 0.11295E+06, 0.11736E+06,
561 + 0.12191E+06, 0.12661E+06, 0.13145E+06, 0.13643E+06, 0.14157E+06,
562 + 0.14687E+06, 0.15232E+06, 0.15794E+06, 0.16372E+06, 0.16968E+06,
563 + 0.17580E+06, 0.18211E+06, 0.18859E+06, 0.19526E+06, 0.20213E+06,
564 + 0.20918E+06, 0.21643E+06, 0.22388E+06, 0.23154E+06, 0.23941E+06,
565 + 0.24750E+06),
566// c... -- 627
567 stdarrayify( 0.66338E+03, 0.93923E+03, 0.12156E+04,
568 + 0.14938E+04, 0.17766E+04, 0.20676E+04, 0.23705E+04, 0.26891E+04,
569 + 0.30267E+04, 0.33866E+04, 0.37714E+04, 0.41839E+04, 0.46267E+04,
570 + 0.51023E+04, 0.56132E+04, 0.61618E+04, 0.67508E+04, 0.73827E+04,
571 + 0.80603E+04, 0.87863E+04, 0.95636E+04, 0.10395E+05, 0.11284E+05,
572 + 0.12233E+05, 0.13246E+05, 0.14326E+05, 0.15477E+05, 0.16702E+05,
573 + 0.18005E+05, 0.19390E+05, 0.20861E+05, 0.22422E+05, 0.24077E+05,
574 + 0.25832E+05, 0.27689E+05, 0.29655E+05, 0.31734E+05, 0.33931E+05,
575 + 0.36250E+05, 0.38698E+05, 0.41280E+05, 0.44002E+05, 0.46869E+05,
576 + 0.49886E+05, 0.53062E+05, 0.56400E+05, 0.59909E+05, 0.63594E+05,
577 + 0.67462E+05, 0.71521E+05, 0.75777E+05, 0.80238E+05, 0.84911E+05,
578 + 0.89804E+05, 0.94925E+05, 0.10028E+06, 0.10588E+06, 0.11173E+06,
579 + 0.11785E+06, 0.12423E+06, 0.13090E+06, 0.13785E+06, 0.14510E+06,
580 + 0.15265E+06, 0.16053E+06, 0.16873E+06, 0.17727E+06, 0.18615E+06,
581 + 0.19540E+06, 0.20501E+06, 0.21501E+06, 0.22540E+06, 0.23619E+06,
582 + 0.24740E+06, 0.25904E+06, 0.27112E+06, 0.28365E+06, 0.29664E+06,
583 + 0.31012E+06, 0.32409E+06, 0.33856E+06, 0.35356E+06, 0.36908E+06,
584 + 0.38516E+06, 0.40180E+06, 0.41902E+06, 0.43683E+06, 0.45525E+06,
585 + 0.47429E+06, 0.49397E+06, 0.51431E+06, 0.53532E+06, 0.55702E+06,
586 + 0.57943E+06, 0.60256E+06, 0.62644E+06, 0.65107E+06, 0.67648E+06,
587 + 0.70269E+06, 0.72972E+06, 0.75758E+06, 0.78629E+06, 0.81588E+06,
588 + 0.84636E+06, 0.87775E+06, 0.91008E+06, 0.94337E+06, 0.97763E+06,
589 + 0.10129E+07, 0.10492E+07, 0.10865E+07, 0.11249E+07, 0.11644E+07,
590 + 0.12050E+07, 0.12467E+07, 0.12896E+07, 0.13337E+07, 0.13789E+07,
591 + 0.14255E+07),
592// c... -- 638
593 stdarrayify( 0.22737E+03, 0.32194E+03, 0.41671E+03,
594 + 0.51226E+03, 0.60963E+03, 0.71017E+03, 0.81528E+03, 0.92628E+03,
595 + 0.10444E+04, 0.11707E+04, 0.13061E+04, 0.14518E+04, 0.16085E+04,
596 + 0.17772E+04, 0.19588E+04, 0.21542E+04, 0.23644E+04, 0.25903E+04,
597 + 0.28330E+04, 0.30934E+04, 0.33726E+04, 0.36717E+04, 0.39918E+04,
598 + 0.43342E+04, 0.47001E+04, 0.50907E+04, 0.55074E+04, 0.59515E+04,
599 + 0.64244E+04, 0.69276E+04, 0.74626E+04, 0.80310E+04, 0.86344E+04,
600 + 0.92744E+04, 0.99528E+04, 0.10671E+05, 0.11432E+05, 0.12236E+05,
601 + 0.13086E+05, 0.13984E+05, 0.14932E+05, 0.15932E+05, 0.16985E+05,
602 + 0.18096E+05, 0.19265E+05, 0.20495E+05, 0.21788E+05, 0.23148E+05,
603 + 0.24576E+05, 0.26075E+05, 0.27648E+05, 0.29298E+05, 0.31027E+05,
604 + 0.32839E+05, 0.34736E+05, 0.36721E+05, 0.38798E+05, 0.40970E+05,
605 + 0.43240E+05, 0.45611E+05, 0.48087E+05, 0.50671E+05, 0.53368E+05,
606 + 0.56180E+05, 0.59111E+05, 0.62165E+05, 0.65347E+05, 0.68659E+05,
607 + 0.72107E+05, 0.75694E+05, 0.79425E+05, 0.83303E+05, 0.87334E+05,
608 + 0.91522E+05, 0.95872E+05, 0.10039E+06, 0.10507E+06, 0.10994E+06,
609 + 0.11498E+06, 0.12021E+06, 0.12563E+06, 0.13125E+06, 0.13707E+06,
610 + 0.14309E+06, 0.14933E+06, 0.15579E+06, 0.16247E+06, 0.16938E+06,
611 + 0.17653E+06, 0.18392E+06, 0.19156E+06, 0.19946E+06, 0.20761E+06,
612 + 0.21604E+06, 0.22473E+06, 0.23371E+06, 0.24298E+06, 0.25254E+06,
613 + 0.26240E+06, 0.27258E+06, 0.28307E+06, 0.29388E+06, 0.30502E+06,
614 + 0.31651E+06, 0.32834E+06, 0.34052E+06, 0.35307E+06, 0.36599E+06,
615 + 0.37929E+06, 0.39298E+06, 0.40706E+06, 0.42155E+06, 0.43645E+06,
616 + 0.45178E+06, 0.46753E+06, 0.48373E+06, 0.50038E+06, 0.51748E+06,
617 + 0.53506E+06),
618// c... -- 637
619 stdarrayify( 0.13267E+04, 0.18785E+04, 0.24314E+04,
620 + 0.29888E+04, 0.35566E+04, 0.41426E+04, 0.47550E+04, 0.54013E+04,
621 + 0.60886E+04, 0.68232E+04, 0.76109E+04, 0.84574E+04, 0.93678E+04,
622 + 0.10348E+05, 0.11402E+05, 0.12536E+05, 0.13755E+05, 0.15065E+05,
623 + 0.16471E+05, 0.17980E+05, 0.19598E+05, 0.21330E+05, 0.23184E+05,
624 + 0.25166E+05, 0.27283E+05, 0.29543E+05, 0.31953E+05, 0.34521E+05,
625 + 0.37256E+05, 0.40164E+05, 0.43256E+05, 0.46541E+05, 0.50026E+05,
626 + 0.53723E+05, 0.57641E+05, 0.61790E+05, 0.66180E+05, 0.70823E+05,
627 + 0.75729E+05, 0.80910E+05, 0.86378E+05, 0.92145E+05, 0.98224E+05,
628 + 0.10463E+06, 0.11137E+06, 0.11846E+06, 0.12592E+06, 0.13375E+06,
629 + 0.14198E+06, 0.15062E+06, 0.15969E+06, 0.16920E+06, 0.17916E+06,
630 + 0.18959E+06, 0.20052E+06, 0.21196E+06, 0.22392E+06, 0.23642E+06,
631 + 0.24949E+06, 0.26314E+06, 0.27740E+06, 0.29227E+06, 0.30779E+06,
632 + 0.32398E+06, 0.34085E+06, 0.35842E+06, 0.37673E+06, 0.39579E+06,
633 + 0.41563E+06, 0.43626E+06, 0.45772E+06, 0.48003E+06, 0.50322E+06,
634 + 0.52730E+06, 0.55232E+06, 0.57829E+06, 0.60524E+06, 0.63320E+06,
635 + 0.66219E+06, 0.69226E+06, 0.72342E+06, 0.75571E+06, 0.78916E+06,
636 + 0.82380E+06, 0.85966E+06, 0.89678E+06, 0.93518E+06, 0.97490E+06,
637 + 0.10160E+07, 0.10585E+07, 0.11023E+07, 0.11477E+07, 0.11946E+07,
638 + 0.12430E+07, 0.12929E+07, 0.13445E+07, 0.13977E+07, 0.14526E+07,
639 + 0.15093E+07, 0.15677E+07, 0.16280E+07, 0.16901E+07, 0.17541E+07,
640 + 0.18200E+07, 0.18880E+07, 0.19579E+07, 0.20300E+07, 0.21042E+07,
641 + 0.21805E+07, 0.22591E+07, 0.23400E+07, 0.24232E+07, 0.25087E+07,
642 + 0.25967E+07, 0.26871E+07, 0.27801E+07, 0.28757E+07, 0.29739E+07,
643 + 0.30747E+07),
644// c... -- 828
645 stdarrayify( 0.60334E+02, 0.85430E+02, 0.11058E+03,
646 + 0.13590E+03, 0.16167E+03, 0.18821E+03, 0.21588E+03, 0.24502E+03,
647 + 0.27595E+03, 0.30896E+03, 0.34431E+03, 0.38225E+03, 0.42301E+03,
648 + 0.46684E+03, 0.51397E+03, 0.56464E+03, 0.61907E+03, 0.67753E+03,
649 + 0.74027E+03, 0.80753E+03, 0.87961E+03, 0.95676E+03, 0.10393E+04,
650 + 0.11275E+04, 0.12217E+04, 0.13222E+04, 0.14293E+04, 0.15434E+04,
651 + 0.16648E+04, 0.17940E+04, 0.19312E+04, 0.20769E+04, 0.22315E+04,
652 + 0.23954E+04, 0.25691E+04, 0.27529E+04, 0.29474E+04, 0.31530E+04,
653 + 0.33702E+04, 0.35995E+04, 0.38414E+04, 0.40965E+04, 0.43654E+04,
654 + 0.46484E+04, 0.49464E+04, 0.52598E+04, 0.55892E+04, 0.59353E+04,
655 + 0.62988E+04, 0.66803E+04, 0.70804E+04, 0.74998E+04, 0.79394E+04,
656 + 0.83998E+04, 0.88817E+04, 0.93859E+04, 0.99132E+04, 0.10464E+05,
657 + 0.11040E+05, 0.11642E+05, 0.12270E+05, 0.12925E+05, 0.13609E+05,
658 + 0.14321E+05, 0.15064E+05, 0.15838E+05, 0.16643E+05, 0.17482E+05,
659 + 0.18355E+05, 0.19263E+05, 0.20207E+05, 0.21188E+05, 0.22208E+05,
660 + 0.23267E+05, 0.24366E+05, 0.25508E+05, 0.26692E+05, 0.27921E+05,
661 + 0.29195E+05, 0.30516E+05, 0.31886E+05, 0.33304E+05, 0.34773E+05,
662 + 0.36294E+05, 0.37869E+05, 0.39499E+05, 0.41185E+05, 0.42929E+05,
663 + 0.44732E+05, 0.46596E+05, 0.48522E+05, 0.50513E+05, 0.52569E+05,
664 + 0.54692E+05, 0.56884E+05, 0.59146E+05, 0.61481E+05, 0.63890E+05,
665 + 0.66375E+05, 0.68937E+05, 0.71578E+05, 0.74301E+05, 0.77107E+05,
666 + 0.79998E+05, 0.82976E+05, 0.86043E+05, 0.89201E+05, 0.92452E+05,
667 + 0.95799E+05, 0.99242E+05, 0.10278E+06, 0.10643E+06, 0.11018E+06,
668 + 0.11403E+06, 0.11799E+06, 0.12206E+06, 0.12625E+06, 0.13055E+06,
669 + 0.13497E+06),
670// c... -- 827
671 stdarrayify( 0.70354E+03, 0.99615E+03, 0.12893E+04,
672 + 0.15846E+04, 0.18848E+04, 0.21940E+04, 0.25162E+04, 0.28554E+04,
673 + 0.32152E+04, 0.35991E+04, 0.40099E+04, 0.44507E+04, 0.49242E+04,
674 + 0.54332E+04, 0.59802E+04, 0.65681E+04, 0.71996E+04, 0.78776E+04,
675 + 0.86050E+04, 0.93847E+04, 0.10220E+05, 0.11114E+05, 0.12070E+05,
676 + 0.13091E+05, 0.14182E+05, 0.15345E+05, 0.16585E+05, 0.17906E+05,
677 + 0.19311E+05, 0.20805E+05, 0.22393E+05, 0.24078E+05, 0.25865E+05,
678 + 0.27760E+05, 0.29768E+05, 0.31893E+05, 0.34140E+05, 0.36516E+05,
679 + 0.39025E+05, 0.41674E+05, 0.44469E+05, 0.47416E+05, 0.50520E+05,
680 + 0.53789E+05, 0.57229E+05, 0.60847E+05, 0.64650E+05, 0.68645E+05,
681 + 0.72840E+05, 0.77242E+05, 0.81859E+05, 0.86699E+05, 0.91770E+05,
682 + 0.97081E+05, 0.10264E+06, 0.10846E+06, 0.11454E+06, 0.12090E+06,
683 + 0.12754E+06, 0.13447E+06, 0.14171E+06, 0.14927E+06, 0.15715E+06,
684 + 0.16536E+06, 0.17392E+06, 0.18284E+06, 0.19213E+06, 0.20179E+06,
685 + 0.21185E+06, 0.22231E+06, 0.23319E+06, 0.24450E+06, 0.25625E+06,
686 + 0.26845E+06, 0.28112E+06, 0.29427E+06, 0.30791E+06, 0.32206E+06,
687 + 0.33674E+06, 0.35196E+06, 0.36772E+06, 0.38406E+06, 0.40098E+06,
688 + 0.41850E+06, 0.43663E+06, 0.45539E+06, 0.47480E+06, 0.49488E+06,
689 + 0.51564E+06, 0.53710E+06, 0.55928E+06, 0.58219E+06, 0.60586E+06,
690 + 0.63029E+06, 0.65553E+06, 0.68157E+06, 0.70844E+06, 0.73616E+06,
691 + 0.76476E+06, 0.79424E+06, 0.82464E+06, 0.85597E+06, 0.88826E+06,
692 + 0.92153E+06, 0.95580E+06, 0.99108E+06, 0.10274E+07, 0.10648E+07,
693 + 0.11033E+07, 0.11429E+07, 0.11837E+07, 0.12256E+07, 0.12687E+07,
694 + 0.13131E+07, 0.13586E+07, 0.14055E+07, 0.14536E+07, 0.15031E+07,
695 + 0.15539E+07),
696// c... -- 727
697 stdarrayify( 0.20518E+04, 0.29051E+04, 0.37601E+04,
698 + 0.46209E+04, 0.54961E+04, 0.63969E+04, 0.73353E+04, 0.83227E+04,
699 + 0.93698E+04, 0.10486E+05, 0.11681E+05, 0.12962E+05, 0.14337E+05,
700 + 0.15815E+05, 0.17403E+05, 0.19110E+05, 0.20942E+05, 0.22909E+05,
701 + 0.25018E+05, 0.27278E+05, 0.29699E+05, 0.32290E+05, 0.35060E+05,
702 + 0.38019E+05, 0.41177E+05, 0.44545E+05, 0.48135E+05, 0.51957E+05,
703 + 0.56023E+05, 0.60346E+05, 0.64938E+05, 0.69812E+05, 0.74981E+05,
704 + 0.80461E+05, 0.86264E+05, 0.92406E+05, 0.98902E+05, 0.10577E+06,
705 + 0.11302E+06, 0.12067E+06, 0.12875E+06, 0.13726E+06, 0.14622E+06,
706 + 0.15566E+06, 0.16559E+06, 0.17604E+06, 0.18702E+06, 0.19855E+06,
707 + 0.21066E+06, 0.22336E+06, 0.23669E+06, 0.25065E+06, 0.26528E+06,
708 + 0.28061E+06, 0.29664E+06, 0.31342E+06, 0.33096E+06, 0.34930E+06,
709 + 0.36845E+06, 0.38845E+06, 0.40933E+06, 0.43111E+06, 0.45383E+06,
710 + 0.47751E+06, 0.50219E+06, 0.52790E+06, 0.55466E+06, 0.58252E+06,
711 + 0.61151E+06, 0.64166E+06, 0.67300E+06, 0.70558E+06, 0.73943E+06,
712 + 0.77458E+06, 0.81108E+06, 0.84896E+06, 0.88827E+06, 0.92904E+06,
713 + 0.97131E+06, 0.10151E+07, 0.10605E+07, 0.11076E+07, 0.11563E+07,
714 + 0.12068E+07, 0.12590E+07, 0.13130E+07, 0.13689E+07, 0.14267E+07,
715 + 0.14865E+07, 0.15483E+07, 0.16121E+07, 0.16781E+07, 0.17462E+07,
716 + 0.18165E+07, 0.18892E+07, 0.19641E+07, 0.20415E+07, 0.21213E+07,
717 + 0.22036E+07, 0.22884E+07, 0.23759E+07, 0.24661E+07, 0.25590E+07,
718 + 0.26547E+07, 0.27533E+07, 0.28549E+07, 0.29594E+07, 0.30670E+07,
719 + 0.31778E+07, 0.32918E+07, 0.34090E+07, 0.35296E+07, 0.36536E+07,
720 + 0.37812E+07, 0.39123E+07, 0.40470E+07, 0.41855E+07, 0.43278E+07,
721 + 0.44739E+07),
722// c... -- 838
723 stdarrayify( 0.12066E+03, 0.17085E+03, 0.22116E+03,
724 + 0.27190E+03, 0.32364E+03, 0.37711E+03, 0.43305E+03, 0.49219E+03,
725 + 0.55516E+03, 0.62256E+03, 0.69492E+03, 0.77276E+03, 0.85657E+03,
726 + 0.94685E+03, 0.10441E+04, 0.11488E+04, 0.12614E+04, 0.13826E+04,
727 + 0.15127E+04, 0.16525E+04, 0.18024E+04, 0.19630E+04, 0.21351E+04,
728 + 0.23191E+04, 0.25158E+04, 0.27260E+04, 0.29502E+04, 0.31892E+04,
729 + 0.34438E+04, 0.37148E+04, 0.40031E+04, 0.43094E+04, 0.46346E+04,
730 + 0.49797E+04, 0.53455E+04, 0.57331E+04, 0.61434E+04, 0.65775E+04,
731 + 0.70364E+04, 0.75212E+04, 0.80330E+04, 0.85730E+04, 0.91424E+04,
732 + 0.97423E+04, 0.10374E+05, 0.11039E+05, 0.11738E+05, 0.12474E+05,
733 + 0.13246E+05, 0.14057E+05, 0.14908E+05, 0.15801E+05, 0.16737E+05,
734 + 0.17717E+05, 0.18744E+05, 0.19819E+05, 0.20944E+05, 0.22120E+05,
735 + 0.23349E+05, 0.24634E+05, 0.25975E+05, 0.27376E+05, 0.28837E+05,
736 + 0.30361E+05, 0.31950E+05, 0.33605E+05, 0.35330E+05, 0.37126E+05,
737 + 0.38996E+05, 0.40942E+05, 0.42965E+05, 0.45069E+05, 0.47256E+05,
738 + 0.49528E+05, 0.51888E+05, 0.54338E+05, 0.56882E+05, 0.59521E+05,
739 + 0.62259E+05, 0.65097E+05, 0.68040E+05, 0.71090E+05, 0.74249E+05,
740 + 0.77522E+05, 0.80910E+05, 0.84417E+05, 0.88046E+05, 0.91801E+05,
741 + 0.95684E+05, 0.99699E+05, 0.10385E+06, 0.10814E+06, 0.11257E+06,
742 + 0.11715E+06, 0.12187E+06, 0.12675E+06, 0.13179E+06, 0.13699E+06,
743 + 0.14235E+06, 0.14788E+06, 0.15358E+06, 0.15946E+06, 0.16552E+06,
744 + 0.17176E+06, 0.17819E+06, 0.18482E+06, 0.19164E+06, 0.19867E+06,
745 + 0.20590E+06, 0.21335E+06, 0.22101E+06, 0.22889E+06, 0.23699E+06,
746 + 0.24533E+06, 0.25390E+06, 0.26271E+06, 0.27177E+06, 0.28108E+06,
747 + 0.29064E+06),
748// c... -- 837
749 stdarrayify( 0.14071E+04, 0.19923E+04, 0.25789E+04,
750 + 0.31704E+04, 0.37733E+04, 0.43962E+04, 0.50477E+04, 0.57360E+04,
751 + 0.64687E+04, 0.72525E+04, 0.80938E+04, 0.89984E+04, 0.99723E+04,
752 + 0.11021E+05, 0.12150E+05, 0.13366E+05, 0.14673E+05, 0.16079E+05,
753 + 0.17589E+05, 0.19211E+05, 0.20949E+05, 0.22812E+05, 0.24807E+05,
754 + 0.26940E+05, 0.29221E+05, 0.31656E+05, 0.34254E+05, 0.37023E+05,
755 + 0.39972E+05, 0.43111E+05, 0.46449E+05, 0.49996E+05, 0.53762E+05,
756 + 0.57756E+05, 0.61991E+05, 0.66477E+05, 0.71226E+05, 0.76249E+05,
757 + 0.81558E+05, 0.87167E+05, 0.93088E+05, 0.99334E+05, 0.10592E+06,
758 + 0.11286E+06, 0.12016E+06, 0.12785E+06, 0.13594E+06, 0.14444E+06,
759 + 0.15337E+06, 0.16274E+06, 0.17258E+06, 0.18290E+06, 0.19371E+06,
760 + 0.20504E+06, 0.21691E+06, 0.22933E+06, 0.24233E+06, 0.25592E+06,
761 + 0.27012E+06, 0.28496E+06, 0.30046E+06, 0.31663E+06, 0.33351E+06,
762 + 0.35111E+06, 0.36946E+06, 0.38858E+06, 0.40850E+06, 0.42924E+06,
763 + 0.45083E+06, 0.47329E+06, 0.49666E+06, 0.52095E+06, 0.54620E+06,
764 + 0.57243E+06, 0.59967E+06, 0.62796E+06, 0.65732E+06, 0.68778E+06,
765 + 0.71938E+06, 0.75214E+06, 0.78611E+06, 0.82131E+06, 0.85777E+06,
766 + 0.89553E+06, 0.93463E+06, 0.97511E+06, 0.10170E+07, 0.10603E+07,
767 + 0.11051E+07, 0.11514E+07, 0.11993E+07, 0.12488E+07, 0.12999E+07,
768 + 0.13527E+07, 0.14073E+07, 0.14636E+07, 0.15217E+07, 0.15816E+07,
769 + 0.16435E+07, 0.17072E+07, 0.17730E+07, 0.18408E+07, 0.19107E+07,
770 + 0.19827E+07, 0.20569E+07, 0.21334E+07, 0.22121E+07, 0.22931E+07,
771 + 0.23765E+07, 0.24624E+07, 0.25507E+07, 0.26416E+07, 0.27351E+07,
772 + 0.28312E+07, 0.29301E+07, 0.30317E+07, 0.31361E+07, 0.32434E+07,
773 + 0.33537E+07));
774
775 gsi = xgj[iso-1];
776 const auto& q = qoft[iso-1];
777
778 if (t < 70 or t > 3000)
779 qt = -1.0;
780 else
781 qt = atob(t, tdat, q);
782}
783
784void calcw(CommonBlock& cmn,
785 const Index& n,
786 const Index& iband,
787 const Numeric& temp)
788{
789 for (Index i=0; i<n; i++) {
790 cmn.YLT.YT[i] = 0;
791 for (Index j=0; j<n; j++) {
792 cmn.Wmatrix.W(i, j) = 0;
793 }
794 }
795
796 if (cmn.Bands.li[iband] > 8 or std::abs(cmn.Bands.li[iband] - cmn.Bands.lf[iband]) > 1) {
797 for (Index i=0; i<n; i++) {
798 cmn.Wmatrix.W(i, i) = cmn.GamT.HWT[i];
799 }
800 return;
801 }
802
803 Vector s(n);
804 for (Index i=0; i<n; i++) {
805 s[i] = cmn.LineSg.Sig(i, iband) * cmn.PopuT.PopuT[i] * Math::pow2(cmn.DipoTcm.DipoT(i, iband));
806 }
807
808 for (Index i=0; i<n; i++) {
809 for (Index j=i+1; j<n; j++) {
810 if (s[j] > s[i]) {
811 // Switch switches 2 elements in a vector :: std::swap(A(i), A(j));
812 // SwitchTwo switches two elements in a Matrix :: std::swap(A(i, k), A(j, k));
813 // SwitchInt is as SwitchTwo but for Integers :: std::swap(A(i, k), A(j, k));
814 std::swap(cmn.DipoRigid.Dipo0(i, iband), cmn.DipoRigid.Dipo0(j, iband));
815 std::swap(cmn.DipoTcm.DipoT(i, iband), cmn.DipoTcm.DipoT(j, iband));
816 std::swap(cmn.LineSg.Sig(i, iband), cmn.LineSg.Sig(j, iband));
817 std::swap(cmn.Jiln.Ji(i, iband), cmn.Jiln.Ji(j, iband));
818 std::swap(cmn.Jfln.Jf(i, iband), cmn.Jfln.Jf(j, iband));
819 std::swap(cmn.GamT.HWSDV2T[i], cmn.GamT.HWSDV2T[j]);
820 std::swap(cmn.PopuT.PopuT[i], cmn.PopuT.PopuT[j]);
821 std::swap(cmn.SHIFT.shft[i], cmn.SHIFT.shft[j]);
822 std::swap(cmn.GamT.HWT[i], cmn.GamT.HWT[j]);
823 std::swap(s[i], s[j]);
824 }
825 }
826 }
827
828 const Numeric dlgt0t = std::log(parameters::T0 / temp);
829 const Index lli = std::min(cmn.Bands.li[iband], cmn.Bands.lf[iband]);
830 const Index llf = std::max(cmn.Bands.li[iband], cmn.Bands.lf[iband]);
831
832 // Set off-diagonal elements
833 for (Index i=0; i<n; i++) {
834 Index jji, jjf;
835 if (cmn.Bands.li[iband] <= cmn.Bands.lf[iband]) {
836 jji=cmn.Jiln.Ji(i, iband);
837 jjf=cmn.Jfln.Jf(i, iband);
838 } else {
839 jji=cmn.Jfln.Jf(i, iband);
840 jjf=cmn.Jiln.Ji(i, iband);
841 }
842
843 for (Index j=0; j<n; j++) {
844 Index jjip, jjfp;
845 if (cmn.Bands.li[iband] <= cmn.Bands.lf[iband]) {
846 jjip=cmn.Jiln.Ji(j, iband);
847 jjfp=cmn.Jfln.Jf(j, iband);
848 } else {
849 jjip=cmn.Jfln.Jf(j, iband);
850 jjfp=cmn.Jiln.Ji(j, iband);
851 }
852
853 if (jjip > jji) continue;
854
855 if (cmn.Bands.Isot[iband] > 2 and cmn.Bands.Isot[iband] not_eq 7 and cmn.Bands.Isot[iband] not_eq 10 and (std::abs(cmn.Jiln.Ji(i, iband)-cmn.Jiln.Ji(j, iband)) % 2) not_eq 0) continue;
856
857 Numeric w0=0, b0=0;
858 if (jji > jjf and jjip > jjfp) {
859 w0 = cmn.Wfittedp.W0pp(lli, llf, jji, jjip);
860 b0 = cmn.Bfittedp.B0pp(lli, llf, jji, jjip);
861 } else if (jji > jjf and jjip == jjfp) {
862 w0 = cmn.Wfittedp.W0pq(lli, llf, jji, jjip);
863 b0 = cmn.Bfittedp.B0pq(lli, llf, jji, jjip);
864 } else if (jji > jjf and jjip < jjfp) {
865 w0 = cmn.Wfittedp.W0pr(lli, llf, jji, jjip);
866 b0 = cmn.Bfittedp.B0pr(lli, llf, jji, jjip);
867 } else if (jji < jjf and jjip > jjfp) {
868 w0 = cmn.Wfittedr.W0rp(lli, llf, jji, jjip);
869 b0 = cmn.Bfittedr.B0rp(lli, llf, jji, jjip);
870 } else if (jji < jjf and jjip == jjfp) {
871 w0 = cmn.Wfittedr.W0rq(lli, llf, jji, jjip);
872 b0 = cmn.Bfittedr.B0rq(lli, llf, jji, jjip);
873 } else if (jji < jjf and jjip < jjfp) {
874 w0 = cmn.Wfittedr.W0rr(lli, llf, jji, jjip);
875 b0 = cmn.Bfittedr.B0rr(lli, llf, jji, jjip);
876 } else if (jji == jjf and jjip > jjfp) {
877 w0 = cmn.Wfittedq.W0qp(lli, llf, jji, jjip);
878 b0 = cmn.Bfittedq.B0qp(lli, llf, jji, jjip);
879 } else if (jji == jjf and jjip == jjfp) {
880 w0 = cmn.Wfittedq.W0qq(lli, llf, jji, jjip);
881 b0 = cmn.Bfittedq.B0qq(lli, llf, jji, jjip);
882 } else if (jji == jjf and jjip < jjfp) {
883 w0 = cmn.Wfittedq.W0qr(lli, llf, jji, jjip);
884 b0 = cmn.Bfittedq.B0qr(lli, llf, jji, jjip);
885 }
886
887 const Numeric ycal = std::exp(w0 - b0*dlgt0t);
888 cmn.Wmatrix.W(j, i) = ycal;
889 cmn.Wmatrix.W(i, j) = ycal * cmn.PopuT.PopuT[i] / cmn.PopuT.PopuT[j];
890 }
891 }
892
893 // Weird undocumented minus sign
894 for (Index i=0; i<n; i++) {
895 for (Index j=0; j<n; j++) {
896 if (j not_eq i) {
897 cmn.Wmatrix.W(i, j) = - std::abs(cmn.Wmatrix.W(i, j));
898 }
899 }
900 }
901
902 // Set diagonal to measured broadening
903 for (Index i=0; i<n; i++) {
904 cmn.Wmatrix.W(i, i) = cmn.GamT.HWT[i];
905 }
906
907 // Sum rule correction
908 for (Index i=0; i<n; i++) {
909 Numeric sumlw = 0.0;
910 Numeric sumup = 0.0;
911
912 for (Index j=0; j<n; j++) {
913 if (cmn.Bands.Isot[iband] > 2 and cmn.Bands.Isot[iband] not_eq 7 and cmn.Bands.Isot[iband] not_eq 10 and (std::abs(cmn.Jiln.Ji(i, iband)-cmn.Jiln.Ji(j, iband)) % 2) not_eq 0) continue;
914
915 if (j > i) {
916 sumlw += std::abs(cmn.DipoRigid.Dipo0(j, iband)) * cmn.Wmatrix.W(j, i);
917 } else {
918 sumup += std::abs(cmn.DipoRigid.Dipo0(j, iband)) * cmn.Wmatrix.W(j, i);
919 }
920 }
921
922 for (Index j=i+1; j<n; j++) {
923 if (sumlw == 0) {
924 cmn.Wmatrix.W(j, i) = 0.0;
925 cmn.Wmatrix.W(i, j) = 0.0;
926 } else {
927 cmn.Wmatrix.W(j, i) = cmn.Wmatrix.W(j, i) * (- sumup / sumlw);
928 cmn.Wmatrix.W(i, j) = cmn.Wmatrix.W(j, i) * cmn.PopuT.PopuT[i] / cmn.PopuT.PopuT[j];
929 }
930 }
931 }
932
933 // Rosenkranz coefficients
934 for (Index i=0; i<n; i++) {
935 Numeric sum0 = 0;
936 for (Index j=0; j<n; j++) {
937 if (i == j) continue;
938
939 if (cmn.Bands.Isot[iband] > 2 and cmn.Bands.Isot[iband] not_eq 7 and cmn.Bands.Isot[iband] not_eq 10 and (std::abs(cmn.Jiln.Ji(i, iband)-cmn.Jiln.Ji(j, iband)) % 2) not_eq 0) continue;
940
941 Numeric deltasig = cmn.LineSg.Sig(i, iband) - cmn.LineSg.Sig(j, iband);
942 if (std::abs(deltasig) < 1e-4 /*cm-1*/ ) deltasig = 1e-4 /*cm-1*/;
943
944 sum0 += 2 * std::abs(cmn.DipoTcm.DipoT(j, iband)) / std::abs(std::abs(cmn.DipoTcm.DipoT(i, iband))) * cmn.Wmatrix.W(j, i) / deltasig;
945 }
946
947 cmn.YLT.YT[i] = sum0;
948 }
949}
950
951// Equivalent lines --- class is not allowed to copy
952using EqvLinesOut = Absorption::LineMixing::EquivalentLines;
953
954// Converted Temperature-Pressure data --- class is not allowed to copy
955struct ConvTPOut {
956 Vector Y, hwt, hwt2, shft, f0, pop, dip;
957 ComplexMatrix W;
958 EqvLinesOut eqv;
959 explicit ConvTPOut(Index n=0) noexcept : Y(n, 0), hwt(n), hwt2(n), shft(n), f0(n), pop(n), dip(n), W(n, n, 0), eqv(n) {}
960 ConvTPOut(const ConvTPOut&) = delete;
961 ConvTPOut(ConvTPOut&&) = default;
962 ConvTPOut& operator=(const ConvTPOut&) = delete;
963 ConvTPOut& operator=(ConvTPOut&&) = default;
964};
965
966
967struct Sorter {
968 Vector dop0;
969 std::vector<Rational> Ji;
970 std::vector<Rational> Ju;
971};
972
973
974Sorter sorter_calcw(ConvTPOut& out,
975 const AbsorptionLines& band,
976 const bool at_t0) {
977 Vector& Y = out.Y;
978 Vector& g0 = out.hwt;
979 Vector& g2 = out.hwt2;
980 Vector& d0 = out.shft;
981 Vector& f0 = out.f0;
982 Vector& pop = out.pop;
983 Vector& dip = out.dip;
984
985 // Size of problem
986 const Index n=Y.nelem();
987
988 // Copies for local variables
989 const Rational li = band.quantumidentity.val[QuantumNumberType::l2].upp();
990 const Rational lf = band.quantumidentity.val[QuantumNumberType::l2].low();
991
992 Vector dip0(n);
993 std::vector<Rational> Ji(n), Jf(n);
994 for (Index i=0; i<n; i++) {
995 Ji[i] = band.lines[i].localquanta.val[QuantumNumberType::J].upp();
996 Jf[i] = band.lines[i].localquanta.val[QuantumNumberType::J].low();
997 dip0[i] = Absorption::reduced_rovibrational_dipole(Jf[i], Ji[i], lf, li, 1);
998 }
999
1000 Vector s(n);
1001 for (Index i=0; i<n; i++) {
1002 if (at_t0) {
1003 s[i] = band.lines[i].I0;
1004 } else {
1005 s[i] = f0[i] * pop[i] * Math::pow2(dip[i]);
1006 }
1007 }
1008
1009 for (Index i=0; i<n; i++) {
1010 for (Index j=i+1; j<n; j++) {
1011 if (s[j] > s[i]) {
1012 // Switch switches 2 elements in a vector :: std::swap(A(i), A(j));
1013 // SwitchTwo switches two elements in a Matrix :: std::swap(A(i, k), A(j, k));
1014 // SwitchInt is as SwitchTwo but for Integers :: std::swap(A(i, k), A(j, k));
1015 std::swap(Ji[i], Ji[j]);
1016 std::swap(Jf[i], Jf[j]);
1017 std::swap(g0[i], g0[j]);
1018 std::swap(d0[i], d0[j]);
1019 std::swap(g2[i], g2[j]);
1020 std::swap(f0[i], f0[j]);
1021 std::swap(dip[i], dip[j]);
1022 std::swap(pop[i], pop[j]);
1023 std::swap(dip0[i], dip0[j]);
1024 std::swap(s[i], s[j]);
1025 }
1026 }
1027 }
1028
1029 return {dip0, Ji, Jf};
1030}
1031
1032
1033void calcw(ConvTPOut& out,
1034 const HitranRelaxationMatrixData& hitran,
1035 const AbsorptionLines& band,
1036 const Numeric T,
1037 const bool at_t0=false) {
1038 Vector& Y = out.Y;
1039 Vector& g0 = out.hwt;
1040 Vector& f0 = out.f0;
1041 Vector& pop = out.pop;
1042 Vector& dip = out.dip;
1043 MatrixView W = transpose(out.W.imag()); // Transpose to fit Fortran-code
1044
1045 // Size of problem
1046 const Index n=Y.nelem();
1047
1048 const Rational li = band.quantumidentity.val[QuantumNumberType::l2].upp();
1049 const Rational lf = band.quantumidentity.val[QuantumNumberType::l2].low();
1050
1051 // Sort before the if-statement. This technicalyl goes awa from how
1052 // HITRAN code does it. It is however required for other line mixing
1053 // setup that we can select the cut off reliably
1054 const auto [dip0, Ji, Jf] = sorter_calcw(out, band, at_t0);
1055
1056 if (li > 8 or abs(li - lf) > 1) {
1057 for (Index i=0; i<n; i++) {
1058 W(i, i) = g0[i]; // nb... units Hz/Pa
1059 }
1060 } else {
1061 const Numeric dlgt0t = std::log(band.T0 / T);
1062 const Rational lli = std::min(li, lf);
1063 const Rational llf = std::max(li, lf);
1064
1065 // Set off-diagonal elements
1066 for (Index i=0; i<n; i++) {
1067 Rational jji, jjf;
1068 if (li <= lf) {
1069 jji = Ji[i];
1070 jjf = Jf[i];
1071 } else {
1072 jji = Jf[i];
1073 jjf = Ji[i];
1074 }
1075
1076 for (Index j=0; j<n; j++) {
1077 Rational jjip, jjfp;
1078 if (li <= lf) {
1079 jjip = Ji[j];
1080 jjfp = Jf[j];
1081 } else {
1082 jjip = Jf[j];
1083 jjfp = Ji[j];
1084 }
1085
1086 if (jjip > jji) continue;
1087
1088 // nb. FIX THIS FOR ARTS SPECIFIC Isotopologue
1089 const bool skip = not (band.Isotopologue().isotname == "626" or
1090 band.Isotopologue().isotname == "636" or
1091 band.Isotopologue().isotname == "828" or
1092 band.Isotopologue().isotname == "838") and
1093 (abs(Ji[i] - Ji[j]) % 2) not_eq 0;
1094
1095 if(skip) continue;
1096
1097 Numeric w0=0, b0=0;
1098 if (jji > jjf and jjip > jjfp) {
1099 w0 = hitran.W0pp(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1100 b0 = hitran.B0pp(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1101 } else if (jji > jjf and jjip == jjfp) {
1102 w0 = hitran.W0pq(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1103 b0 = hitran.B0pq(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1104 } else if (jji > jjf and jjip < jjfp) {
1105 w0 = hitran.W0pr(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1106 b0 = hitran.B0pr(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1107 } else if (jji < jjf and jjip > jjfp) {
1108 w0 = hitran.W0rp(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1109 b0 = hitran.B0rp(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1110 } else if (jji < jjf and jjip == jjfp) {
1111 w0 = hitran.W0rq(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1112 b0 = hitran.B0rq(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1113 } else if (jji < jjf and jjip < jjfp) {
1114 w0 = hitran.W0rr(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1115 b0 = hitran.B0rr(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1116 } else if (jji == jjf and jjip > jjfp) {
1117 w0 = hitran.W0qp(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1118 b0 = hitran.B0qp(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1119 } else if (jji == jjf and jjip == jjfp) {
1120 w0 = hitran.W0qq(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1121 b0 = hitran.B0qq(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1122 } else if (jji == jjf and jjip < jjfp) {
1123 w0 = hitran.W0qr(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1124 b0 = hitran.B0qr(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1125 }
1126
1127 // nb... should order be different???
1128 const Numeric ycal = std::exp(w0 - b0*dlgt0t);
1129 W(j, i) = ycal;
1130 W(i, j) = ycal * pop[i] / pop[j];
1131 }
1132 }
1133
1134 // Weird undocumented minus sign
1135 for (Index i=0; i<n; i++) {
1136 for (Index j=0; j<n; j++) {
1137 if (j not_eq i) {
1138 W(i, j) = - std::abs(W(i, j));
1139 }
1140 }
1141 }
1142
1143 // Set diagonal to measured broadening
1144 for (Index i=0; i<n; i++) {
1145 W(i, i) = g0[i]; // nb... units Hz/Pa
1146 }
1147
1148 // Sum rule correction
1149 for (Index i=0; i<n; i++) {
1150 Numeric sumlw = 0.0;
1151 Numeric sumup = 0.0;
1152
1153 for (Index j=0; j<n; j++) {
1154
1155 // nb. FIX THIS FOR ARTS SPECIFIC Isotopologue
1156 const bool skip = not (band.Isotopologue().isotname == "626" or
1157 band.Isotopologue().isotname == "636" or
1158 band.Isotopologue().isotname == "828" or
1159 band.Isotopologue().isotname == "838") and
1160 (abs(Ji[i] - Ji[j]) % 2) not_eq 0;
1161
1162 if (skip) continue;
1163
1164 if (j > i) {
1165 sumlw += std::abs(dip0[j]) * W(j, i);
1166 } else {
1167 sumup += std::abs(dip0[j]) * W(j, i);
1168 }
1169 }
1170
1171 for (Index j=i+1; j<n; j++) {
1172 if (sumlw == 0) {
1173 W(j, i) = 0.0;
1174 W(i, j) = 0.0;
1175 } else {
1176 W(j, i) *= - sumup / sumlw;
1177 W(i, j) = W(j, i) * pop[i] / pop[j];
1178 }
1179 }
1180 }
1181
1182 // Rosenkranz coefficients
1183 for (Index i=0; i<n; i++) {
1184 Numeric sum0 = 0;
1185 for (Index j=0; j<n; j++) {
1186 if (i == j) continue;
1187
1188 // nb. FIX THIS FOR ARTS SPECIFIC Isotopologue
1189 const bool skip = not (band.Isotopologue().isotname == "626" or
1190 band.Isotopologue().isotname == "636" or
1191 band.Isotopologue().isotname == "828" or
1192 band.Isotopologue().isotname == "838") and
1193 (abs(Ji[i] - Ji[j]) % 2) not_eq 0;
1194
1195 if (skip) continue;
1196
1197 Numeric deltasig = f0[i] - f0[j];
1198 if (std::abs(deltasig) < Conversion::kaycm2freq(1e-4 /*cm-1*/))
1199 deltasig = Conversion::kaycm2freq(1e-4) /*Hz*/;
1200
1201 sum0 += 2 * std::abs(dip[j]) / std::abs(dip[i]) * W(j, i) / deltasig;
1202 }
1203
1204 Y[i] = sum0;
1205 }
1206 }
1207}
1208
1209void eqvlines(CommonBlock& cmn,
1210 const Index& iband,
1211 const Index& n,
1212 const Numeric& sigmoy)
1213{
1214 ComplexMatrix zop(n, n), zvec(n, n);
1215 ComplexMatrixView inv_zvec = zop; // Rename to not confuse later operations
1216 ComplexVectorView zval = cmn.Zaa.ZA[Range(0, n)]; // Rename and rescale to right size
1217 for (Index i=0; i<n; i++) {
1218 for (Index j=0; j<n; j++) {
1219 zop(j, i) = Complex(cmn.DiagnR.OpR(i, j), cmn.DiagnI.OpI(i, j)); // nb. reverse from Fortran algorithm due to row-col issues
1220 }
1221 }
1222
1223 // Main computations
1224 diagonalize(zvec, zval, zop);
1225 inv(inv_zvec, zvec);
1226
1227 // Add average sigma
1228 zval += sigmoy;
1229
1230 // Do the matrix multiplication
1231 for (Index i=0; i<n; i++) {
1232 cmn.Zss.ZS[i] = 0;
1233 Complex z(0, 0);
1234 for (Index j=0; j<n; j++) {
1235 cmn.Zss.ZS[i] += cmn.DipoTcm.DipoT(j, iband) * zvec(j, i);
1236 z += cmn.PopuT.PopuT[j] * cmn.DipoTcm.DipoT(j, iband) * inv_zvec(i, j);
1237 }
1238 cmn.Zss.ZS[i] *= z;
1239 }
1240}
1241
1242EqvLinesOut eqvlines(const ComplexMatrix& W,
1243 const Vector& pop,
1244 const Vector& dip,
1245 const Numeric& fmean)
1246{
1247 // Size of problem
1248 const Index n = pop.nelem();
1249
1250 // Compute values
1251 ComplexMatrix zvec(n, n);
1252 EqvLinesOut out(n);
1253
1254 // Main computations
1255 diagonalize(zvec, out.val, W);
1256
1257 // Do the matrix forward multiplication
1258 for (Index i=0; i<n; i++) {
1259 for (Index j=0; j<n; j++) {
1260 out.str[i] += dip[j] * zvec(j, i);
1261 }
1262 }
1263
1264 // Do the matrix backward multiplication
1265 inv(zvec, zvec);
1266 for (Index i=0; i<n; i++) {
1267 Complex z(0, 0);
1268 for (Index j=0; j<n; j++) {
1269 z += pop[j] * dip[j] * zvec(i, j);
1270 }
1271 out.str[i] *= z;
1272 }
1273
1274 // Add the weighted frequency
1275 out.val += fmean;
1276
1277 return out;
1278}
1279
1280
1281void convtp(CommonBlock& cmn,
1282 const Index& iband,
1283 const Index& isotc,
1284 const Index& nlinec,
1285 const Numeric& xh2o,
1286 const Numeric& xco2,
1287 const Numeric& temp,
1288 const Numeric& ptot,
1289 Numeric& sigmoy,
1290 Numeric& gamd_gam0_mx,
1291 const bool mixfull,
1292 const bool mixsdv)
1293{
1294 const Numeric ratiot = parameters::T0 / temp;
1295 Numeric qt0_co2, qtt_co2, gsi0, gsit;
1296 qt_co2(parameters::T0, isotc, gsi0, qt0_co2);
1297 qt_co2(temp, isotc, gsit, qtt_co2);
1298 const Numeric ratiopart = qt0_co2 / qtt_co2;
1299
1300 const Numeric sqrtm = std::sqrt(temp / parameters::aMass[isotc-1]);
1301
1302 Numeric sumwgt=0;
1303 sigmoy=0;
1304 for (Index iline=0; iline<nlinec; iline++) {
1305 cmn.PopuT.PopuT[iline]=cmn.PopTrf.PopuT0(iline, iband)*ratiopart*std::exp(-parameters::Ct*cmn.Energy.E(iline, iband)*(1/temp-1/parameters::T0));
1306 if (mixsdv) {
1307 cmn.GamT.HWT[iline] = ((1-xh2o-xco2) * (cmn.GamSDVT0AIR.HWSDVT0AIR(iline, iband) * std::pow(ratiot, cmn.DTGAMAIR.BHWAIR(iline, iband))))
1308 + (xh2o * (cmn.GamSDVT0H2O.HWSDVT0H2O(iline, iband) * std::pow(ratiot, cmn.DTGAMH2O.BHWH2O(iline, iband))))
1309 + (xco2 * (cmn.GamSDVT0CO2.HWSDVT0SELF(iline, iband) * std::pow(ratiot, cmn.DTGAMCO2.BHWSELF(iline, iband))));
1310
1311 cmn.GamT.HWSDV2T[iline] =
1312 ((1-xh2o-xco2) * (cmn.GamSDVT0AIR.HWSDVT0AIR(iline, iband) * cmn.GamSDVT0AIR.rHWT0AIR(iline, iband) * std::pow(ratiot, cmn.DTGAMAIR.BHWAIR(iline, iband))))
1313 + (xh2o * (cmn.GamSDVT0H2O.HWSDVT0H2O(iline, iband)*cmn.GamSDVT0H2O.rHWT0H2O(iline, iband) * std::pow(ratiot, cmn.DTGAMH2O.BHWH2O(iline, iband))))
1314 + (xco2 * (cmn.GamSDVT0CO2.HWSDVT0SELF(iline, iband) * cmn.GamSDVT0CO2.rHWT0SELF(iline, iband) * std::pow(ratiot, cmn.DTGAMCO2.BHWSELF(iline, iband))));
1315 } else {
1316 cmn.GamT.HWT[iline] = ((1-xh2o-xco2) * (cmn.GamVT0AIR.HWVT0AIR(iline,iband) * std::pow(ratiot, cmn.DTGAMAIR.BHWAIR(iline, iband))))
1317 + (xh2o * (cmn.GamVT0H2O.HWVT0H2O(iline, iband) * std::pow(ratiot, cmn.DTGAMH2O.BHWH2O(iline, iband))))
1318 + (xco2 * (cmn.GamVT0CO2.HWVT0SELF(iline, iband) * std::pow(ratiot, cmn.DTGAMCO2.BHWSELF(iline, iband))));
1319 }
1320 cmn.SHIFT.shft[iline]=cmn.SHIFT0.shft0(iline, iband);
1321 const Numeric wgt = cmn.PopuT.PopuT[iline] * Math::pow2(cmn.DipoTcm.DipoT(iline, iband));
1322 sumwgt += wgt;
1323 sigmoy += cmn.LineSg.Sig(iline,iband) * wgt;
1324 }
1325 sigmoy /= sumwgt;
1326
1327 calcw(cmn, nlinec, iband, temp);
1328
1329 // Adjust for pressure
1330 for (Index iline=0; iline<nlinec; iline++) {
1331 cmn.GamT.HWT[iline] *= ptot;
1332 cmn.GamT.HWSDV2T[iline] *= ptot;
1333 cmn.YLT.YT[iline] *= ptot;
1334 const Numeric gamd = parameters::CtGamD * cmn.LineSg.Sig(iline, iband) * sqrtm;
1335 gamd_gam0_mx = std::max(gamd_gam0_mx, gamd/cmn.GamT.HWT[iline]);
1336 }
1337
1338 if (mixfull) {
1339 for (Index iline=0; iline<nlinec; iline++) {
1340 for (Index ilinep=0; ilinep<nlinec; ilinep++) {
1341 cmn.DiagnI.OpI(iline, ilinep) = ptot * cmn.Wmatrix.W(iline, ilinep);
1342 if (iline == ilinep) {
1343 cmn.DiagnR.OpR(iline, iline) = cmn.LineSg.Sig(iline, iband) - sigmoy + cmn.SHIFT.shft[iline] * ptot;
1344 } else {
1345 cmn.DiagnR.OpR(iline, ilinep) = 0;
1346 }
1347 }
1348 }
1349
1350 eqvlines(cmn, iband, nlinec, sigmoy);
1351 }
1352}
1353
1354ConvTPOut convtp(const Vector& vmrs,
1355 const HitranRelaxationMatrixData& hitran,
1356 const AbsorptionLines& band,
1357 const Numeric T,
1358 const Numeric P)
1359{
1360 const Index n = band.NumLines();
1361
1362 const Numeric QT = single_partition_function(T, band.Isotopologue());
1363 const Numeric QT0 = single_partition_function(band.T0, band.Isotopologue());
1364 const Numeric ratiopart = QT0 / QT;
1365
1366 ConvTPOut out(n);
1367
1368 Vector wgt(n);
1369 for (Index i=0; i<n; i++) {
1370 const Numeric pop0 = (band.lines[i].gupp / QT0) * boltzman_factor(band.T0, band.lines[i].E0);
1371
1372 out.f0[i] = band.lines[i].F0;
1373 out.pop[i] = pop0 * ratiopart * boltzman_ratio(T, band.T0, band.lines[i].E0);
1374 out.hwt[i] = band.lines[i].lineshape.G0(T, band.T0, 1, vmrs);
1375 out.shft[i] = band.lines[i].lineshape.D0(T, band.T0, 1, vmrs);
1376 out.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))));
1377 out.hwt2[i] = band.lines[i].lineshape.G2(T, band.T0, 1, vmrs);
1378 wgt[i] = out.pop[i] * Math::pow2(out.dip[i]);
1379 }
1380
1381 // Calculate the relaxation matrix
1382 calcw(out, hitran, band, T);
1383
1384 // Adjust for pressure
1385 out.hwt *= P;
1386 out.hwt2 *= P;
1387 out.shft *= P;
1388 out.Y *= P;
1389
1390 if (band.population == Absorption::PopulationType::ByHITRANFullRelmat) {
1391 const Numeric fmean = band.F_mean(wgt);
1392 out.W.diagonal().real() = out.f0;
1393 out.W.diagonal().real() += out.shft;
1394 out.W.diagonal().real() -= fmean;
1395 out.W.imag() *= P;
1396 out.eqv = eqvlines(out.W, out.pop, out.dip, fmean);
1397 }
1398
1399 return out;
1400}
1401
1402
1403void qsdv(const Numeric& sg0,
1404 const Numeric& gamd,
1405 const Numeric& gam0,
1406 const Numeric& gam2,
1407 const Numeric& shift0,
1408 const Numeric& shift2,
1409 const Numeric& sg,
1410 Numeric& ls_qsdv_r,
1411 Numeric& ls_qsdv_i)
1412{
1413 using Math::pow2;
1414 Complex x , y, csqrty, z1, z2, w1, w2, aterm;
1415 Numeric xz1, yz1, xz2, yz2;
1416
1417 const Numeric cte = Constant::sqrt_ln_2 / gamd;
1418 const Numeric pi = Constant::pi;
1419 const Numeric rpi = Constant::sqrt_pi;
1420 const Complex iz = Complex(0, 1);
1421
1422 const Complex c0 = Complex(gam0, shift0);
1423 const Complex c2 = Complex(gam2, shift2);
1424 const Complex c0t = (c0 - 1.5 * c2);
1425 const Complex c2t = c2;
1426
1427 if (std::abs(c2t) == 0) goto label110;
1428 x = (iz * (sg0 - sg) + c0t) / c2t;
1429 y = 1 / pow2(2 * cte * c2t);
1430 csqrty = (gam2 - iz * shift2) / (2 * cte * (pow2(gam2) + pow2(shift2)));
1431 if (std::abs(x) <= 3.e-8*std::abs(y)) goto label120; // IN cgs units
1432 if (std::abs(y) <= 1.e-15*std::abs(x)) goto label140; // IN cgs units
1433 z1 = std::sqrt(x + y) - csqrty;
1434 z2 = z1 + 2 * csqrty;
1435
1436 xz1 = - z1.imag();
1437 yz1 = z1.real();
1438 xz2 = - z2.imag();
1439 yz2 = z2.real();
1440
1441 w1 = Faddeeva::w(Complex(xz1 ,yz1));
1442 w2 = Faddeeva::w(Complex(xz2 ,yz2));
1443 aterm = rpi * cte * (w1 - w2);
1444 goto label10;
1445
1446 label110:
1447 z1 = (iz*(sg0-sg)+c0t)*cte;
1448 xz1 = - z1.imag();
1449 yz1 = z1.real();
1450
1451 w1 = Faddeeva::w(Complex(xz1 ,yz1));
1452 aterm = rpi*cte* w1;
1453 goto label10;
1454
1455 label120:
1456
1457 z1 = (iz*(sg0-sg)+c0t)*cte;
1458 z2 = std::sqrt(x+y)+csqrty;
1459 xz1 = - z1.imag();
1460 yz1 = z1.real();
1461 xz2 = - z2.imag();
1462 yz2 = z2.real();
1463
1464 w1 = Faddeeva::w(Complex(xz1 ,yz1));
1465 w2 = Faddeeva::w(Complex(xz2 ,yz2));
1466 aterm = rpi * cte * (w1-w2);
1467 goto label10;
1468
1469 label140:
1470
1471 if (std::abs(std::sqrt(x)) <= 4.e3) { // IN cgs
1472 const Numeric xxb = - std::sqrt(x).imag();
1473 const Numeric yxb = std::sqrt(x).real();
1474 const Complex wb = Faddeeva::w(Complex(xxb, yxb));
1475 aterm = (2 * rpi / c2t) * (1 / rpi - std::sqrt(x) * wb);
1476 } else {
1477 aterm = (1 / c2t) * (1 / x - 1.5 / pow2(x));
1478 }
1479
1480 label10:
1481
1482 const Complex ls_qsdv = (1 / pi) * aterm;
1483
1484 ls_qsdv_r = ls_qsdv.real();
1485 ls_qsdv_i = ls_qsdv.imag();
1486}
1487
1488
1489// nb.qsdv times sg0 is qsdv_si times F0
1490Complex qsdv_si(const Numeric F0,
1491 const Numeric gamd,
1492 const Numeric gam0,
1493 const Numeric gam2,
1494 const Numeric shift0,
1495 const Numeric shift2,
1496 const Numeric f)
1497{
1498 using Math::pow2;
1499 Complex x , y, csqrty, z1, z2, w1, w2, aterm;
1500 Numeric xz1, yz1, xz2, yz2;
1501
1502 const Numeric cte = Constant::sqrt_ln_2 / gamd;
1503 const Numeric pi = Constant::pi;
1504 const Numeric rpi = Constant::sqrt_pi;
1505 const Complex iz = Complex(0, 1);
1506
1507 const Complex c0 = Complex(gam0, shift0);
1508 const Complex c2 = Complex(gam2, shift2);
1509 const Complex c0t = (c0 - 1.5 * c2);
1510 const Complex c2t = c2;
1511
1512 if (std::abs(c2t) == 0) goto label110;
1513 x = (iz * (F0 - f) + c0t) / c2t;
1514 y = 1 / pow2(2 * cte * c2t);
1515 csqrty = (gam2 - iz * shift2) / (2 * cte * (pow2(gam2) + pow2(shift2)));
1516 if (std::abs(x) <= 3.e-8*std::abs(y)) goto label120;
1517 if (std::abs(y) <= 1.e-15*std::abs(x)) goto label140;
1518 z1 = std::sqrt(x + y) - csqrty;
1519 z2 = z1 + 2 * csqrty;
1520
1521 xz1 = - z1.imag();
1522 yz1 = z1.real();
1523 xz2 = - z2.imag();
1524 yz2 = z2.real();
1525
1526 w1 = Faddeeva::w(Complex(xz1 ,yz1));
1527 w2 = Faddeeva::w(Complex(xz2 ,yz2));
1528 aterm = rpi * cte * (w1 - w2);
1529 goto label10;
1530
1531 label110:
1532 z1 = (iz*(F0 - f) + c0t)*cte;
1533 xz1 = - z1.imag();
1534 yz1 = z1.real();
1535
1536 w1 = Faddeeva::w(Complex(xz1 ,yz1));
1537 aterm = rpi*cte* w1;
1538 goto label10;
1539
1540 label120:
1541
1542 z1 = (iz*(F0 - f) + c0t)*cte;
1543 z2 = std::sqrt(x+y)+csqrty;
1544 xz1 = - z1.imag();
1545 yz1 = z1.real();
1546 xz2 = - z2.imag();
1547 yz2 = z2.real();
1548
1549 w1 = Faddeeva::w(Complex(xz1 ,yz1));
1550 w2 = Faddeeva::w(Complex(xz2 ,yz2));
1551 aterm = rpi * cte * (w1-w2);
1552 goto label10;
1553
1554 label140:
1555
1556 if (std::abs(std::sqrt(x)) <= 4.e3) { // IN cgs
1557 const Numeric xxb = - std::sqrt(x).imag();
1558 const Numeric yxb = std::sqrt(x).real();
1559 const Complex wb = Faddeeva::w(Complex(xxb, yxb));
1560 aterm = (2 * rpi / c2t) * (1 / rpi - std::sqrt(x) * wb);
1561 } else {
1562 aterm = (1 / c2t) * (1 / x - 1.5 / pow2(x));
1563 }
1564
1565 label10:
1566
1567 return (1 / pi) * aterm;
1568}
1569
1570void compabs(
1571 CommonBlock& cmn,
1572 const Numeric& temp,
1573 const Numeric& ptot,
1574 const Numeric& xco2,
1575 const Numeric& xh2o,
1576 const Vector& invcm_grid,
1577 const bool mixsdv,
1578 const bool mixfull,
1579 VectorView absv,
1580 VectorView absy,
1581 VectorView absw)
1582{
1583 using Math::pow2;
1584 using Math::pow3;
1585
1586 constexpr Numeric rdmult = 30;
1587
1588 // Number of points to compute
1589 const Index nsig = invcm_grid.nelem();
1590
1591 // Set to zero
1592 absv = 0;
1593 absy = 0;
1594 absw = 0;
1595
1596 constexpr Numeric sq_ln2 = Constant::sqrt_ln_2;
1597 constexpr Numeric sq_ln2pi = sq_ln2 / Constant::sqrt_pi;
1598 const Numeric dens = xco2 * ptot * parameters::aMolAtm / temp;
1599 constexpr Numeric u_pi = Constant::inv_pi;
1600 constexpr Numeric u_sqln2pi = 1 / sq_ln2pi;
1601
1602 for (Index iband=0; iband<cmn.Bands.nBand; iband++) {
1603 Numeric sigmoy=0;
1604 Numeric gamdmx=0;
1605 convtp(cmn, iband, cmn.Bands.Isot[iband], cmn.Bands.nLines[iband],
1606 xh2o, xco2, temp, ptot, sigmoy, gamdmx, mixfull, mixsdv);
1607 const Numeric sqrtm = std::sqrt(temp / parameters::aMass[cmn.Bands.Isot[iband]-1]);
1608
1609 const Numeric ds0=(70.67e0+104.1e0*0.21e0*std::pow(gamdmx, 6.4))/(1.0+0.21*std::pow(gamdmx, 5.4));
1610 const Numeric ds2=(34.97e0+105.e0*9.e0*std::pow(gamdmx, 3.1))/(1.0+9.0*std::pow(gamdmx, 2.1));
1611
1612 for (Index isig=0; isig<nsig; isig++) {
1613 const Numeric sigc = invcm_grid[isig];
1614
1615 for (Index iline=0; iline<cmn.Bands.nLines[iband]; iline++) {
1616 const Numeric gamd=parameters::CtGamD*cmn.LineSg.Sig(iline,iband)*sqrtm;
1617 const Numeric cte = sq_ln2 / gamd;
1618 const Numeric popudipo = cmn.PopuT.PopuT[iline] * pow2(cmn.DipoTcm.DipoT(iline, iband));
1619 const Numeric dsigc = sigc - cmn.LineSg.Sig(iline, iband) - cmn.SHIFT.shft[iline] * ptot;
1620
1621 if (mixsdv) {
1622 if (std::abs(dsigc) / cmn.GamT.HWT[iline] > ds0) {
1623 absv[isig] += popudipo * u_sqln2pi * u_pi * cmn.GamT.HWT[iline] / (pow2(cmn.GamT.HWT[iline]) + pow2(dsigc));
1624 absy[isig] += popudipo * u_sqln2pi * u_pi * (cmn.GamT.HWT[iline] + cmn.YLT.YT[iline]*dsigc) / (pow2(cmn.GamT.HWT[iline]) + pow2(dsigc));
1625 } else if (std::abs(dsigc) / cmn.GamT.HWT[iline] > ds2 and
1626 std::abs(dsigc) / cmn.GamT.HWT[iline] < ds0) {
1627 const Numeric order2r = cmn.GamT.HWT[iline] / (pow2(cmn.GamT.HWT[iline]) + pow2(dsigc))
1628 + 3 * pow2(cmn.GamT.HWSDV2T[iline]) * (pow3(cmn.GamT.HWT[iline]) - 3 * cmn.GamT.HWT[iline]*pow2(dsigc)) /
1629 (2*Constant::pi*pow3(pow2(cmn.GamT.HWT[iline]) + pow2(dsigc)));
1630
1631 const Numeric order2i = dsigc / (pow2(cmn.GamT.HWT[iline]) + pow2(dsigc)) +
1632 3 * pow2(cmn.GamT.HWSDV2T[iline]) *(3*pow2(cmn.GamT.HWT[iline])*dsigc-pow3(dsigc)) /
1633 (2*Constant::pi*pow3(pow2(cmn.GamT.HWT[iline]) + pow2(dsigc)));
1634
1635 absv[isig] += popudipo * u_sqln2pi * u_pi * order2r;
1636
1637 absy[isig] += popudipo * u_sqln2pi * u_pi * (order2r + cmn.YLT.YT[iline]*order2i);
1638 } else {
1639 const Numeric shft0 = cmn.SHIFT.shft[iline] * ptot;
1640
1641 Numeric wr=0, wi=0;
1642 qsdv(cmn.LineSg.Sig(iline, iband), gamd, cmn.GamT.HWT[iline], cmn.GamT.HWSDV2T[iline], shft0, 0, sigc, wr, wi);
1643
1644 absv[isig] += popudipo * wr * u_sqln2pi;
1645
1646 absy[isig] += popudipo * u_sqln2pi * (wr - cmn.YLT.YT[iline] * wi);
1647 }
1648 } else {
1649 if (std::abs(cmn.LineSg.Sig(iline, iband)-sigc) > (rdmult*gamd)) { // NOTE: Removed in updated version of the code...
1650 absv[isig] += popudipo * u_sqln2pi * u_pi * cmn.GamT.HWT[iline] / (pow2(cmn.GamT.HWT[iline]) + pow2(dsigc));
1651
1652 absy[isig] += popudipo * u_sqln2pi * u_pi * (cmn.GamT.HWT[iline] + cmn.YLT.YT[iline] * dsigc) / (pow2(cmn.GamT.HWT[iline]) + pow2(dsigc));
1653 } else {
1654 const Numeric yy = cmn.GamT.HWT[iline] * cte;
1655 const Numeric xx= (cmn.LineSg.Sig(iline, iband)+cmn.SHIFT.shft[iline] * ptot - sigc) * cte;
1656 const Complex w = Faddeeva::w(Complex(xx, yy));
1657
1658 absv[isig] += popudipo * w.real() / gamd;
1659
1660 absy[isig] += popudipo * (w.real() - cmn.YLT.YT[iline] * w.imag()) / gamd;
1661 }
1662 }
1663
1664 if (mixfull) {
1665 /*const Numeric gamd_int=parameters::CtGamD*(cmn.FicLPR.AlphR[iline] + sigmoy)*sqrtm;
1666 const Numeric cte_int = sq_ln2 / gamd_int;
1667 if (std::abs(cmn.FicLPI.AlphI[iline] + sigmoy - sigc) <= (rdmult*gamd_int)) {
1668 const Complex z = Complex(cmn.FicLPR.AlphR[iline] + sigmoy - sigc, cmn.FicLPI.AlphI[iline]) * cte_int;
1669 const Complex w = Faddeeva::w(z);
1670 absw[isig] += (cmn.FicLSR.SSR[iline] * w.real() - cmn.FicLSI.SSI[iline] * w.imag()) / cte_int;
1671 }*/
1672 absw[isig] += u_sqln2pi * u_pi * (cmn.Zss.ZS[iline] / (sigc - cmn.Zaa.ZA[iline])).imag();
1673 }
1674 }
1675 }
1676 }
1677
1678 for (Index isig=0; isig<nsig; isig++) {
1679 const Numeric sigc = invcm_grid[isig];
1680 const Numeric fact = sigc * (1-std::exp(-parameters::Ct * sigc / temp));
1681 absv[isig] *= fact * dens * sq_ln2pi;
1682 absy[isig] *= fact * dens * sq_ln2pi;
1683 absw[isig] *= fact * dens * sq_ln2pi;
1684 }
1685}
1686
1687Vector compabs(
1688 const Numeric T,
1689 const Numeric P,
1690 const HitranRelaxationMatrixData& hitran,
1691 const ArrayOfAbsorptionLines& bands,
1692 const SpeciesIsotopologueRatios& isotopologue_ratio,
1693 const Vector &vmrs,
1694 const Vector &f_grid)
1695{
1696 using Math::pow2;
1697 using Math::pow3;
1698
1699 // Number of points to compute
1700 const Index nf = f_grid.nelem();
1701
1702 // Set to zero
1703 Vector absorption(nf, 0);
1704
1705 constexpr Numeric sq_ln2 = Constant::sqrt_ln_2;
1706 constexpr Numeric sq_ln2pi = sq_ln2 / Constant::sqrt_pi;
1707 const Numeric dens = vmrs[0] * number_density(P, T);
1708 constexpr Numeric u_pi = Constant::inv_pi;
1709 constexpr Numeric u_sqln2pi = 1 / sq_ln2pi;
1710
1711 for (Index iband=0; iband<bands.nelem(); iband++) {
1712 if (not bands[iband].DoLineMixing(P)) continue;
1713
1714 const Numeric rat_isot = isotopologue_ratio[bands[iband].Isotopologue()];
1715
1716 auto tp = convtp(vmrs, hitran, bands[iband], T, P);
1717 const Numeric GD_div_F0 = bands[iband].DopplerConstant(T);
1718
1719 const bool sdvp = bands[iband].lineshapetype == LineShape::Type::SDVP;
1720 const bool vp = bands[iband].lineshapetype == LineShape::Type::VP;
1721 const bool rosenkranz = bands[iband].population == Absorption::PopulationType::ByHITRANRosenkranzRelmat;
1722 const bool full = bands[iband].population == Absorption::PopulationType::ByHITRANFullRelmat;
1723
1724 for (Index iv=0; iv<nf; iv++) {
1725 const Numeric f = f_grid[iv];
1726
1727 Numeric a=0;
1728 for (Index iline=0; iline<bands[iband].NumLines(); iline++) {
1729 const Numeric gamd=GD_div_F0 * tp.f0[iline];
1730 const Numeric gamd_mod=GD_div_F0 * tp.eqv.val[iline].real();
1731 const Numeric cte = sq_ln2 / gamd;
1732 const Numeric cte_mod = sq_ln2 / gamd_mod;
1733 const Numeric popudipo = tp.pop[iline] * pow2(tp.dip[iline]);
1734
1735 if (rosenkranz and sdvp) {
1736 const Complex w = qsdv_si(tp.f0[iline], gamd, tp.hwt[iline], tp.hwt2[iline], tp.shft[iline], 0, f);
1737 a += popudipo * u_sqln2pi * (Complex(1, tp.Y[iline]) * w).real(); // NB. Changing sign on Y gave positive absorption but is not agreeing with measurements according to HITRAN data
1738 } else if (rosenkranz and vp) {
1739 const Numeric yy = tp.hwt[iline] * cte;
1740 const Numeric xx = (tp.f0[iline]+tp.shft[iline] - f) * cte;
1741 const Complex w = Faddeeva::w(Complex(xx, yy));
1742 a += popudipo * (Complex(1, tp.Y[iline]) * w).real() / gamd; // NB. Changing sign on Y gave positive absorption but is not agreeing with measurements according to HITRAN data
1743 } else if (full and vp) {
1744 const Complex z = (tp.eqv.val[iline]-f) * cte_mod;
1745 const Complex w = Faddeeva::w(z);
1746 a += (tp.eqv.str[iline] * w).real() / gamd_mod;
1747 } else if (full) {
1748 a += u_sqln2pi * u_pi * (tp.eqv.str[iline] / (f - tp.eqv.val[iline])).imag();
1749 } else {
1750 ARTS_USER_ERROR ("Cannot understand the combination of calculations requested...");
1751 }
1752 }
1753
1754 // Guard to not allow negative absorption inside a band
1755 if (std::isnormal(a) and a > 0)
1756 absorption[iv] += rat_isot * a;
1757 }
1758 }
1759
1760 for (Index iv=0; iv<nf; iv++) {
1761 const Numeric f = f_grid[iv];
1762 const Numeric fact = - f * std::expm1(- (Constant::h * f) / (Constant::k * T));
1763 absorption[iv] *= fact * dens * sq_ln2pi;
1764 }
1765
1766 return absorption;
1767}
1768
1769void detband(CommonBlock& cmn,
1770 const Numeric& sgminr,
1771 const Numeric& sgmaxr,
1772 const Numeric& stotmax,
1773 const String& basedir="data_new/")
1774{
1775 cmn.Bands.nBand=0;
1776 std::ifstream fortranfile;
1777 fortranfile.open(basedir + "/BandInfo.dat");
1778
1779 ARTS_USER_ERROR_IF (not fortranfile.is_open(),
1780 "Cannot read the file. Please make sure you point at BandInfo.dat basedir and have the right to read.");
1781
1782 String line;
1783 getline(fortranfile, line);
1784 while (fortranfile.good()) {
1785 Numeric stot, sgmin, sgmax;
1786 Index isotr, lfr, lir, jmxp, jmxq, jmxr;
1787 char c11, c12, c21, c22, c31, c32, c41, c42, c51, c52, x;
1788 sscanf(line.c_str(),
1789 "%1" PRId64
1790 "%c%c" "%1" PRId64 "%c%c"
1791 "%c%c" "%1" PRId64 "%c%c"
1792 "%c%c"
1793 "%12lf"
1794 "%c" "%12lf"
1795 "%c" "%12lf"
1796 "%c%c%c%c%c%c%c%c"
1797 "%4" PRId64 "%4" PRId64 "%4" PRId64,
1798 &isotr,
1799 &c11, &c12, &lfr, &c21, &c22,
1800 &c31, &c32, &lir, &c41, &c42,
1801 &c51, &c52,
1802 &stot,
1803 &x, &sgmin,
1804 &x, &sgmax,
1805 &x, &x, &x, &x, &x, &x, &x, &x,
1806 &jmxp, &jmxq, &jmxr);
1807
1808 getline(fortranfile, line);
1809 if (stot < stotmax)
1810 continue;
1811
1812 if ((sgminr < sgmax) and (sgmaxr > sgmin)) {
1813 cmn.Bands.Isot[cmn.Bands.nBand] = isotr;
1814 if (isotr == 0)
1815 cmn.Bands.Isot[cmn.Bands.nBand] = 10;
1816 cmn.Bands.li[cmn.Bands.nBand] = lir;
1817 cmn.Bands.lf[cmn.Bands.nBand] = lfr;
1818
1819 std::array<char, 15> name;
1820 snprintf(name.data(),
1821 name.size(),
1822 "S%" PRId64 "%c%c%" PRId64 "%c%c%c%c%" PRId64 "%c%c%c%c",
1823 isotr,
1824 c11,
1825 c12,
1826 lfr,
1827 c21,
1828 c22,
1829 c31,
1830 c32,
1831 lir,
1832 c41,
1833 c42,
1834 c51,
1835 c52);
1836 cmn.Bands.BandFile[cmn.Bands.nBand] = name.data();
1837 cmn.Bands.nBand++;
1838 ARTS_USER_ERROR_IF (cmn.Bands.nBand > parameters::nBmx,
1839 "Too many bands");
1840 }
1841 }
1842 fortranfile.close();
1843}
1844
1845void readw(CommonBlock& cmn, const String& basedir="data_new/")
1846{
1847 for (Index l=0; l<=8; l++) {
1848 for (Index ideltal=0; ideltal<=1; ideltal++) {
1849 const Index lli = l;
1850 const Index llf = l + ideltal;
1851
1852 const String cr = std::to_string(lli) + std::to_string(llf);
1853
1854 std::ifstream fortranfile;
1855 const String fname = basedir + String("/WTfit") + cr + String(".dat");
1856 fortranfile.open(fname.c_str());
1857
1858 String line;
1859 getline(fortranfile, line);
1860 while (fortranfile.good()) {
1861 char sw0r[21], sb0r[21];
1862 Numeric dmaxdt, wtmax;
1863 Index jic, jfc, jipc, jfpc;
1864 sscanf(line.c_str(),
1865 "%20s" "%20s"
1866 "%14lf" "%14lf"
1867 "%4" PRId64 "%4" PRId64 "%4" PRId64 "%4" PRId64,
1868 sw0r, sb0r, &dmaxdt, &wtmax, &jic, &jfc, &jipc, &jfpc);
1869 String ssw0r = sw0r;
1870 std::replace(ssw0r.begin(), ssw0r.end(), 'D', 'E');
1871 const Numeric w0r = std::stod(ssw0r);
1872 String ssb0r = sb0r;
1873 std::replace(ssb0r.begin(), ssb0r.end(), 'D', 'E');
1874 const Numeric b0r = std::stod(ssb0r);
1875
1876 getline(fortranfile, line);
1877
1878 ARTS_USER_ERROR_IF (lli > cmn.Wfittedp.W0pp.nbooks() or
1879 llf > cmn.Wfittedp.W0pp.npages() or
1880 jic > cmn.Wfittedp.W0pp.nrows() or
1881 jipc > cmn.Wfittedp.W0pp.ncols(),
1882 "Out of bounds in reading...");
1883
1884 if (jic > jfc and jipc > jfpc) {
1885 cmn.Wfittedp.W0pp(lli,llf,jic,jipc) = w0r;
1886 cmn.Bfittedp.B0pp(lli,llf,jic,jipc) = b0r;
1887 } else if(jic > jfc and jipc == jfpc) {
1888 cmn.Wfittedp.W0pq(lli,llf,jic,jipc) = w0r;
1889 cmn.Bfittedp.B0pq(lli,llf,jic,jipc) = b0r;
1890 } else if(jic > jfc and jipc < jfpc) {
1891 cmn.Wfittedp.W0pr(lli,llf,jic,jipc) = w0r;
1892 cmn.Bfittedp.B0pr(lli,llf,jic,jipc) = b0r;
1893 } else if(jic < jfc and jipc > jfpc) {
1894 cmn.Wfittedr.W0rp(lli,llf,jic,jipc) = w0r;
1895 cmn.Bfittedr.B0rp(lli,llf,jic,jipc) = b0r;
1896 } else if(jic < jfc and jipc == jfpc) {
1897 cmn.Wfittedr.W0rq(lli,llf,jic,jipc) = w0r;
1898 cmn.Bfittedr.B0rq(lli,llf,jic,jipc) = b0r;
1899 } else if(jic < jfc and jipc < jfpc) {
1900 cmn.Wfittedr.W0rr(lli,llf,jic,jipc) = w0r;
1901 cmn.Bfittedr.B0rr(lli,llf,jic,jipc) = b0r;
1902 } else if(jic == jfc and jipc > jfpc) {
1903 cmn.Wfittedq.W0qp(lli,llf,jic,jipc) = w0r;
1904 cmn.Bfittedq.B0qp(lli,llf,jic,jipc) = b0r;
1905 } else if(jic == jfc and jipc == jfpc) {
1906 cmn.Wfittedq.W0qq(lli,llf,jic,jipc) = w0r;
1907 cmn.Bfittedq.B0qq(lli,llf,jic,jipc) = b0r;
1908 } else if(jic == jfc and jipc < jfpc) {
1909 cmn.Wfittedq.W0qr(lli,llf,jic,jipc) = w0r;
1910 cmn.Bfittedq.B0qr(lli,llf,jic,jipc) = b0r;
1911 }
1912 }
1913 }
1914 }
1915}
1916
1917Vector compute(const Numeric p, const Numeric t, const Numeric xco2, const Numeric xh2o, const Vector& invcm_grid, const Numeric stotmax, const calctype type)
1918{
1919 const Index n = invcm_grid.nelem();
1920 Vector absorption(n);
1921 if (not n)
1922 return absorption;
1923
1924 // Common setup and main IO
1925 CommonBlock cmn;
1926 detband(cmn, invcm_grid[0], invcm_grid[n-1], stotmax);
1927 readw(cmn);
1928 readlines(cmn);
1929
1930 Vector absv(n);
1931 Vector absy(n);
1932 Vector absw(n);
1933 if (type == calctype::FullVP or type == calctype::FullRosenkranz or type == calctype::FullW)
1934 compabs(cmn, t,p, xco2, xh2o, invcm_grid, false, true, absv, absy, absw);
1935 else if (type == calctype::SDVP or type == calctype::SDRosenkranz or type == calctype::SDW)
1936 compabs(cmn, t,p, xco2, xh2o, invcm_grid, true, false, absv, absy, absw);
1937 else if (type == calctype::NoneVP or type == calctype::NoneRosenkranz or type == calctype::NoneW)
1938 compabs(cmn, t,p, xco2, xh2o, invcm_grid, false, false, absv, absy, absw);
1939
1940 switch(type) {
1941 case calctype::FullVP:
1942 case calctype::SDVP:
1943 case calctype::NoneVP:
1944 for (Index i=0; i<n; i++)
1945 absorption[i] = absv[i];
1946 break;
1947 case calctype::FullRosenkranz:
1948 case calctype::SDRosenkranz:
1949 case calctype::NoneRosenkranz:
1950 for (Index i=0; i<n; i++)
1951 absorption[i] = absy[i];
1952 break;
1953 case calctype::FullW:
1954 case calctype::SDW:
1955 case calctype::NoneW:
1956 for (Index i=0; i<n; i++)
1957 absorption[i] = absw[i];
1958 break;
1959 }
1960
1961 return absorption;
1962}
1963
1964Vector compute(const HitranRelaxationMatrixData& hitran,
1965 const ArrayOfAbsorptionLines& bands,
1966 const SpeciesIsotopologueRatios& isotopologue_ratio,
1967 const Numeric P,
1968 const Numeric T,
1969 const Vector& vmrs,
1970 const Vector& f_grid)
1971{
1972 return f_grid.nelem() ? compabs(T, P, hitran, bands, isotopologue_ratio, vmrs, f_grid) : Vector(0);
1973}
1974
1975void read(HitranRelaxationMatrixData& hitran,
1976 ArrayOfAbsorptionLines& bands,
1977 const SpeciesIsotopologueRatios& isotopologue_ratio,
1978 const String& basedir,
1979 const Numeric linemixinglimit,
1980 const Numeric fmin,
1981 const Numeric fmax,
1982 const Numeric stot,
1983 const ModeOfLineMixing mode)
1984{
1985 String newbase = basedir;
1986 if (newbase.nelem() == 0)
1987 newbase = ".";
1988
1989 CommonBlock cmn;
1990 detband(cmn, Conversion::freq2kaycm(fmin), Conversion::freq2kaycm(fmax), Conversion::hz_per_msquared2kaycm_per_cmsquared(stot), newbase);
1991 readw(cmn, newbase);
1992 readlines(cmn, newbase);
1993
1994 Numeric linemixinglimit_internal=0;
1995 switch(mode) {
1996 case ModeOfLineMixing::VP:
1997 case ModeOfLineMixing::SDVP: break;
1998 case ModeOfLineMixing::FullW:
1999 case ModeOfLineMixing::VP_W:
2000 case ModeOfLineMixing::VP_Y:
2001 case ModeOfLineMixing::SDVP_Y: linemixinglimit_internal=linemixinglimit; break;
2002 case ModeOfLineMixing::FINAL: ARTS_ASSERT(true, "Bad mode input. Must update function.");
2003 }
2004
2005 const auto lstype = typeLP(mode) ?
2006 LineShape::Type::LP :
2007 (typeVP(mode) ?
2008 LineShape::Type::VP :
2009 LineShape::Type::SDVP) ;
2010 const auto poptype = typeFull(mode) ?
2011 Absorption::PopulationType::ByHITRANFullRelmat :
2012 Absorption::PopulationType::ByHITRANRosenkranzRelmat ;
2013
2014 auto specs=stdarrayify(
2015 SpeciesTag("CO2-626"),
2016 SpeciesTag("CO2-636"),
2017 SpeciesTag("CO2-628"),
2018 SpeciesTag("CO2-627"),
2019 SpeciesTag("CO2-638"),
2020 SpeciesTag("CO2-637"),
2021 SpeciesTag("CO2-828"),
2022 SpeciesTag("CO2-827"),
2023 SpeciesTag("CO2-727"),
2024 SpeciesTag("CO2-838"),
2025 SpeciesTag("CO2-837"));
2026
2027 // Move data from Fortran-style common block to ARTS
2028 hitran.B0pp = std::move(cmn.Bfittedp.B0pp);
2029 hitran.B0pq = std::move(cmn.Bfittedp.B0pq);
2030 hitran.B0pr = std::move(cmn.Bfittedp.B0pr);
2031 hitran.W0pp = std::move(cmn.Wfittedp.W0pp);
2032 hitran.W0pq = std::move(cmn.Wfittedp.W0pq);
2033 hitran.W0pr = std::move(cmn.Wfittedp.W0pr);
2034 hitran.B0qp = std::move(cmn.Bfittedq.B0qp);
2035 hitran.B0qq = std::move(cmn.Bfittedq.B0qq);
2036 hitran.B0qr = std::move(cmn.Bfittedq.B0qr);
2037 hitran.W0qp = std::move(cmn.Wfittedq.W0qp);
2038 hitran.W0qq = std::move(cmn.Wfittedq.W0qq);
2039 hitran.W0qr = std::move(cmn.Wfittedq.W0qr);
2040 hitran.B0rp = std::move(cmn.Bfittedr.B0rp);
2041 hitran.B0rq = std::move(cmn.Bfittedr.B0rq);
2042 hitran.B0rr = std::move(cmn.Bfittedr.B0rr);
2043 hitran.W0rp = std::move(cmn.Wfittedr.W0rp);
2044 hitran.W0rq = std::move(cmn.Wfittedr.W0rq);
2045 hitran.W0rr = std::move(cmn.Wfittedr.W0rr);
2046
2047 // Reshape to band count size
2048 bands.resize(cmn.Bands.nBand);
2049 for (Index i{0}; i<cmn.Bands.nBand; i++) {
2050 QuantumIdentifier qid{specs[cmn.Bands.Isot[i]-1].spec_ind};
2051
2052 qid.val.set(Quantum::Number::Value(
2053 QuantumNumberType::l2, cmn.Bands.li[i], cmn.Bands.lf[i]));
2054 qid.val.set(Quantum::Number::Value(QuantumNumberType::v1,
2055 cmn.UnusedBandParams.iv1[i],
2056 cmn.UnusedBandParams.fv1[i]));
2057 qid.val.set(Quantum::Number::Value(QuantumNumberType::v2,
2058 cmn.UnusedBandParams.iv2[i],
2059 cmn.UnusedBandParams.fv2[i]));
2060 qid.val.set(Quantum::Number::Value(QuantumNumberType::v3,
2061 cmn.UnusedBandParams.iv3[i],
2062 cmn.UnusedBandParams.fv3[i]));
2063 qid.val.set(Quantum::Number::Value(QuantumNumberType::r,
2064 cmn.UnusedBandParams.ir[i],
2065 cmn.UnusedBandParams.fr[i]));
2066
2067 bands[i] = {true,
2068 true,
2069 Absorption::CutoffType::None,
2070 Absorption::MirroringType::None,
2071 poptype,
2072 Absorption::NormalizationType::None,
2073 lstype,
2074 296,
2075 -1,
2076 linemixinglimit_internal,
2077 qid,
2078 {Species::fromShortName("CO2"), Species::fromShortName("H2O"), Species::fromShortName("AIR")}};
2079
2080 const Numeric rat_isot = isotopologue_ratio[bands[i].Isotopologue()];
2081
2082 bands[i].lines.resize(cmn.Bands.nLines[i]);
2083 for (Index j{0}; j<cmn.Bands.nLines[i]; j++) {
2084 Quantum::Number::LocalState local_state(
2085 Quantum::Number::Value(QuantumNumberType::J,
2086 Rational(cmn.Jiln.Ji(j, i)),
2087 Rational(cmn.Jfln.Jf(j, i))));
2088
2089 const LineShape::ModelParameters G0_sdvp_air{LineShape::TemperatureModel::T1,
2090 Conversion::kaycm_per_atm2hz_per_pa(cmn.GamSDVT0AIR.HWSDVT0AIR(j, i)),
2091 cmn.DTGAMAIR.BHWAIR(j, i)};
2092 const LineShape::ModelParameters G0_sdvp_h2o{LineShape::TemperatureModel::T1,
2093 Conversion::kaycm_per_atm2hz_per_pa(cmn.GamSDVT0H2O.HWSDVT0H2O(j, i)),
2094 cmn.DTGAMH2O.BHWH2O(j, i)};
2095 const LineShape::ModelParameters G0_sdvp_co2{LineShape::TemperatureModel::T1,
2096 Conversion::kaycm_per_atm2hz_per_pa(cmn.GamSDVT0CO2.HWSDVT0SELF(j, i)),
2097 cmn.DTGAMCO2.BHWSELF(j, i)};
2098
2099 LineShape::ModelParameters G2_sdvp_air{G0_sdvp_air}; G2_sdvp_air.X0 *= cmn.GamSDVT0AIR.rHWT0AIR(j, i);
2100 LineShape::ModelParameters G2_sdvp_h2o{G0_sdvp_h2o}; G2_sdvp_h2o.X0 *= cmn.GamSDVT0H2O.rHWT0H2O(j, i);
2101 LineShape::ModelParameters G2_sdvp_co2{G0_sdvp_co2}; G2_sdvp_co2.X0 *= cmn.GamSDVT0CO2.rHWT0SELF(j, i);
2102
2103 const LineShape::ModelParameters D0{LineShape::TemperatureModel::T0, Conversion::kaycm_per_atm2hz_per_pa(cmn.SHIFT0.shft0(j, i))};
2104
2105 const LineShape::ModelParameters G0_vp_air{LineShape::TemperatureModel::T1,
2106 Conversion::kaycm_per_atm2hz_per_pa(cmn.GamVT0AIR.HWVT0AIR(j, i)),
2107 cmn.DTGAMAIR.BHWAIR(j, i)};
2108 const LineShape::ModelParameters G0_vp_h2o{LineShape::TemperatureModel::T1,
2109 Conversion::kaycm_per_atm2hz_per_pa(cmn.GamVT0H2O.HWVT0H2O(j, i)),
2110 cmn.DTGAMH2O.BHWH2O(j, i)};
2111 const LineShape::ModelParameters G0_vp_co2{LineShape::TemperatureModel::T1,
2112 Conversion::kaycm_per_atm2hz_per_pa(cmn.GamVT0CO2.HWVT0SELF(j, i)),
2113 cmn.DTGAMCO2.BHWSELF(j, i)};
2114
2115 const LineShape::SingleSpeciesModel sdvp_air{G0_sdvp_air, D0, G2_sdvp_air};
2116 const LineShape::SingleSpeciesModel sdvp_h2o{G0_sdvp_h2o, D0, G2_sdvp_h2o};
2117 const LineShape::SingleSpeciesModel sdvp_co2{G0_sdvp_co2, D0, G2_sdvp_co2};
2118 const LineShape::SingleSpeciesModel vp_air{G0_vp_air, D0};
2119 const LineShape::SingleSpeciesModel vp_h2o{G0_vp_h2o, D0};
2120 const LineShape::SingleSpeciesModel vp_co2{G0_vp_co2, D0};
2121 const auto lsmodel = typeVP(mode) ?
2122 LineShape::Model{{vp_co2,vp_h2o, vp_air}}:
2123 LineShape::Model{{sdvp_co2, sdvp_h2o, sdvp_air}};
2124
2125 Numeric qt0_co2, gsi0;
2126 qt_co2(parameters::T0, cmn.Bands.Isot[i], gsi0, qt0_co2);
2127
2128 // Should probably use the isotopologue ratio
2129 bands[i].lines[j] = {Conversion::kaycm2freq(cmn.LineSg.Sig(j, i)),
2130 Conversion::kaycm_per_cmsquared2hz_per_msquared(cmn.UnusedBandParams.intens(j, i)) / rat_isot,
2131 Conversion::kaycm2joule(cmn.Energy.E(j, i)),
2132 gsi0*Numeric(cmn.Jfln.Jf(j, i) * 2 + 1),
2133 gsi0*Numeric(cmn.Jiln.Ji(j, i) * 2 + 1),
2134 cmn.UnusedBandParams.eina(j, i),
2135 ZeemanModel{},
2136 lsmodel,
2137 local_state};
2138 }
2139 }
2140}
2141
2142Tensor4 hitran_lm_eigenvalue_approximation(const AbsorptionLines& band,
2143 const Vector& temperatures,
2144 const HitranRelaxationMatrixData& hitran,
2145 const Numeric P) {
2146 const Index N = band.NumLines();
2147 const Index M = band.NumBroadeners();
2148 const Index K = temperatures.nelem();
2149
2150 // Need sorting to put weak lines last, but we need the sorting constant or the output jumps
2151 const Numeric QT0 = single_partition_function(band.T0, band.Isotopologue());
2152 const Numeric fmean = band.F_mean();
2153
2154 // Output
2155 Tensor4 out(4, N, M, K, 0);
2156
2157 #pragma omp parallel for collapse(2) if (!arts_omp_in_parallel())
2158 for (Index m=0; m<M; m++) {
2159 for (Index k=0; k<K; k++) {
2160 ConvTPOut calc(N);
2161
2162 const Numeric T = temperatures[k];
2163
2164 const Numeric QT = single_partition_function(T, band.Isotopologue());
2165 const Numeric ratiopart = QT0 / QT;
2166
2167 for (Index i=0; i<N; i++) {
2168 const Numeric pop0 = (band.lines[i].gupp / QT0) * boltzman_factor(band.T0, band.lines[i].E0);
2169
2170 calc.f0[i] = band.lines[i].F0;
2171 calc.pop[i] = pop0 * ratiopart * boltzman_ratio(T, band.T0, band.lines[i].E0);
2172 calc.hwt[i] = P * band.lines[i].lineshape[m].G0().at(T, band.T0);
2173 calc.shft[i] = P * band.lines[i].lineshape[m].D0().at(T, band.T0);
2174 calc.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))));
2175 calc.hwt2[i] = P * band.lines[i].lineshape[m].G2().at(T, band.T0);
2176 }
2177
2178 // Calculate the relaxation matrix
2179 calcw(calc, hitran, band, T, true);
2180
2181 // Select types depending on population tag
2182 if (Absorption::PopulationType::ByHITRANRosenkranzRelmat == band.population) {
2183 // Sort the values
2184 for (Index i=0; i<calc.Y.nelem(); i++) {
2185 for (Index j=i+1; j<calc.Y.nelem(); j++) {
2186 if (calc.f0[j] < calc.f0[i]) {
2187 std::swap(calc.f0[i], calc.f0[j]);
2188 std::swap(calc.Y[i], calc.Y[j]);
2189 }
2190 }
2191 }
2192
2193 out(1, joker, m, k) = calc.Y;
2194 } else if (Absorption::PopulationType::ByHITRANFullRelmat == band.population) {
2195 calc.W.diagonal().real() = calc.f0;
2196 calc.W.diagonal().real() += calc.shft;
2197 calc.W.diagonal().real() -= fmean;
2198
2199 const auto eig = Absorption::LineMixing::eigenvalue_adaptation_of_relmat(
2200 calc.W, calc.pop, calc.dip, band, fmean, T, P, QT, QT0, m);
2201
2202 out(0, joker, m, k) = eig.str.real();
2203 out(1, joker, m, k) = eig.str.imag();
2204 out(2, joker, m, k) = eig.val.real();
2205 out(3, joker, m, k) = eig.val.imag();
2206 }
2207 }
2208 }
2209
2210 return out;
2211}
2212
2213void hitran_lm_eigenvalue_adaptation(AbsorptionLines& band,
2214 const Vector& temperatures,
2215 const HitranRelaxationMatrixData& hitran,
2216 const Numeric P0,
2217 const Index ord)
2218{
2219 ARTS_USER_ERROR_IF (P0 <= 0, P0, " Pa is not possible")
2220
2221 ARTS_USER_ERROR_IF(
2222 not is_sorted(temperatures),
2223 "The temperature list [", temperatures, "] K\n"
2224 "must be fully sorted from low to high"
2225 )
2226
2227 if (Absorption::LineMixing::band_eigenvalue_adaptation(band,
2228 hitran_lm_eigenvalue_approximation(band, temperatures, hitran, P0),
2229 temperatures, P0, ord)) {
2230 ARTS_USER_ERROR("Bad eigenvalue adaptation")
2231 }
2232}
2233
2234Tensor5 hitran_lm_eigenvalue_adaptation_test(const AbsorptionLines& band,
2235 const Vector& temperatures,
2236 const HitranRelaxationMatrixData& hitran,
2237 const Vector& pressures) {
2238 const Index N = band.NumLines();
2239 const Index M = band.NumBroadeners();
2240 const Index K = temperatures.nelem();
2241 const Index L = pressures.size();
2242
2243 Tensor5 out(4, N, M, K, L);
2244 for (Index l=0; l<L; l++) {
2245 out(joker, joker, joker, joker, l) = hitran_lm_eigenvalue_approximation(band, temperatures, hitran, pressures[l]);
2246 }
2247 return out;
2248}
2249} // namespace lm_hitran_2017
2250
stdarrayify
constexpr std::array< T, 1+sizeof...(Ts)> stdarrayify(const T &first, const Ts &... the_rest)
Make a std::array of a list of variables (must be 1-long at least)
Definition array.h:255
Array
This can be used to make arrays out of anything.
Definition array.h:31
Array::nelem
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
LineShape::Model
Main line shape model class.
Definition lineshapemodel.h:626
LineShape::SingleSpeciesModel
Compute the line shape parameters for a single broadening species.
Definition lineshapemodel.h:433
Quantum::Number::ValueList::set
void set(Value v)
Sets the value if it exists or adds it otherwise.
Definition quantum_numbers.cc:348
Zeeman::Model
Main Zeeman Model.
Definition zeemandata.h:279
my_basic_string< char >
my_basic_string::nelem
Index nelem() const
Definition mystring.h:172
ARTS_ASSERT
#define ARTS_ASSERT(condition,...)
Definition debug.h:86
ARTS_USER_ERROR
#define ARTS_USER_ERROR(...)
Definition debug.h:153
ARTS_USER_ERROR_IF
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
linemixing_hitran.h
Namespace and functions to deal with HITRAN linemixing.
boltzman_factor
Numeric boltzman_factor(Numeric T, Numeric E0)
Computes exp(- E0/kT)
Definition linescaling.cc:90
boltzman_ratio
Numeric boltzman_ratio(const Numeric &T, const Numeric &T0, const Numeric &E0)
Computes exp(E0/c (T - T0) / (T * T0))
Definition linescaling.cc:75
single_partition_function
Numeric single_partition_function(const Numeric &T, const Species::IsotopeRecord &ir)
Computes the partition function at one temperature.
Definition linescaling.cc:4
lineshapemodel.h
Contains the line shape namespace.
String
my_basic_string< char > String
The String type for ARTS.
Definition mystring.h:199
Absorption::LineMixing::eigenvalue_adaptation_of_relmat
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:1282
Absorption::LineMixing::band_eigenvalue_adaptation
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:1110
Absorption::reduced_rovibrational_dipole
Numeric reduced_rovibrational_dipole(Rational Jf, Rational Ji, Rational lf, Rational li, Rational k=Rational(1))
Compute the reduced rovibrational dipole moment.
Definition absorptionlines.cc:1902
Constant::sqrt_ln_2
constexpr Numeric sqrt_ln_2
Square root of natural logarithm of 2.
Definition arts_constants.h:99
Constant::pi
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
Definition arts_constants.h:57
Constant::k
constexpr Numeric k
Boltzmann constant convenience name [J/K].
Definition arts_constants.h:152
Constant::sqrt_pi
constexpr Numeric sqrt_pi
Square root of pi.
Definition arts_constants.h:69
Constant::inv_pi
constexpr Numeric inv_pi
Inverse of pi.
Definition arts_constants.h:60
Constant::h
constexpr Numeric h
Planck constant convenience name [J s].
Definition arts_constants.h:128
Conversion::kaycm2joule
constexpr auto kaycm2joule(auto x) noexcept
Conversion from cm-1 to Joule.
Definition arts_conversions.h:145
Conversion::kaycm_per_cmsquared2hz_per_msquared
constexpr auto kaycm_per_cmsquared2hz_per_msquared(auto x) noexcept
Conversion from cm-1 per molecule per cm^2 to Hz per molecule per m^2.
Definition arts_conversions.h:125
Conversion::kaycm_per_atm2hz_per_pa
constexpr auto kaycm_per_atm2hz_per_pa(auto x) noexcept
Conversion from cm-1 per atmosphere to Hz per Pascal.
Definition arts_conversions.h:135
Conversion::freq2kaycm
constexpr auto freq2kaycm(auto x) noexcept
Conversion from Hz to Kayser wavenumber.
Definition arts_conversions.h:53
Conversion::kaycm2freq
constexpr auto kaycm2freq(auto x) noexcept
Conversion from Kayser wavenumber to Hz.
Definition arts_conversions.h:50
Conversion::hz_per_msquared2kaycm_per_cmsquared
constexpr auto hz_per_msquared2kaycm_per_cmsquared(auto x) noexcept
Conversion from Hz per molecule per m^2 to cm-1 per molecule per cm^2.
Definition arts_conversions.h:130
Math::pow3
constexpr auto pow3(auto x) noexcept
power of three
Definition arts_constexpr_math.h:24
Math::pow2
constexpr auto pow2(auto x) noexcept
power of two
Definition arts_constexpr_math.h:21
lm_hitran_2017
Definition linemixing_hitran.cc:38
lm_hitran_2017::eqvlines
void eqvlines(CommonBlock &cmn, const Index &iband, const Index &n, const Numeric &sigmoy)
Definition linemixing_hitran.cc:1209
lm_hitran_2017::hitran_lm_eigenvalue_adaptation
void hitran_lm_eigenvalue_adaptation(AbsorptionLines &band, const Vector &temperatures, const HitranRelaxationMatrixData &hitran, const Numeric P0, const Index ord)
Definition linemixing_hitran.cc:2213
lm_hitran_2017::qsdv_si
Complex qsdv_si(const Numeric F0, const Numeric gamd, const Numeric gam0, const Numeric gam2, const Numeric shift0, const Numeric shift2, const Numeric f)
Definition linemixing_hitran.cc:1490
lm_hitran_2017::read
void read(HitranRelaxationMatrixData &hitran, ArrayOfAbsorptionLines &bands, const SpeciesIsotopologueRatios &isotopologue_ratio, const String &basedir, const Numeric linemixinglimit, const Numeric fmin, const Numeric fmax, const Numeric stot, const ModeOfLineMixing mode)
Read from HITRAN online line mixing file.
Definition linemixing_hitran.cc:1975
lm_hitran_2017::calcw
void calcw(CommonBlock &cmn, const Index &n, const Index &iband, const Numeric &temp)
Definition linemixing_hitran.cc:784
lm_hitran_2017::hitran_lm_eigenvalue_adaptation_test
Tensor5 hitran_lm_eigenvalue_adaptation_test(const AbsorptionLines &band, const Vector &temperatures, const HitranRelaxationMatrixData &hitran, const Vector &pressures)
Definition linemixing_hitran.cc:2234
lm_hitran_2017::typeFull
constexpr bool typeFull(ModeOfLineMixing x)
Definition linemixing_hitran.h:108
lm_hitran_2017::sorter_calcw
Sorter sorter_calcw(ConvTPOut &out, const AbsorptionLines &band, const bool at_t0)
Definition linemixing_hitran.cc:974
lm_hitran_2017::qt_co2
void qt_co2(const Numeric &t, const Index &iso, Numeric &gsi, Numeric &qt)
Definition linemixing_hitran.cc:468
lm_hitran_2017::readw
void readw(CommonBlock &cmn, const String &basedir="data_new/")
Definition linemixing_hitran.cc:1845
lm_hitran_2017::compute
Vector compute(const Numeric p, const Numeric t, const Numeric xco2, const Numeric xh2o, const Vector &invcm_grid, const Numeric stotmax, const calctype type)
Definition linemixing_hitran.cc:1917
lm_hitran_2017::detband
void detband(CommonBlock &cmn, const Numeric &sgminr, const Numeric &sgmaxr, const Numeric &stotmax, const String &basedir="data_new/")
Definition linemixing_hitran.cc:1769
lm_hitran_2017::toRationalSum
Rational toRationalSum(char a, char b=' ')
Definition linemixing_hitran.cc:272
lm_hitran_2017::atob
Numeric atob(const Numeric &aa, const std::array< Numeric, NT > &a, const std::array< Numeric, NT > &b)
Definition linemixing_hitran.cc:409
lm_hitran_2017::convtp
void convtp(CommonBlock &cmn, const Index &iband, const Index &isotc, const Index &nlinec, const Numeric &xh2o, const Numeric &xco2, const Numeric &temp, const Numeric &ptot, Numeric &sigmoy, Numeric &gamd_gam0_mx, const bool mixfull, const bool mixsdv)
Definition linemixing_hitran.cc:1281
lm_hitran_2017::readlines
void readlines(CommonBlock &cmn, const String &basedir="data_new/")
Definition linemixing_hitran.cc:283
lm_hitran_2017::calctype
calctype
Definition linemixing_hitran.h:48
lm_hitran_2017::calctype::FullW
@ FullW
lm_hitran_2017::calctype::NoneW
@ NoneW
lm_hitran_2017::calctype::SDVP
@ SDVP
lm_hitran_2017::calctype::NoneVP
@ NoneVP
lm_hitran_2017::calctype::NoneRosenkranz
@ NoneRosenkranz
lm_hitran_2017::calctype::SDRosenkranz
@ SDRosenkranz
lm_hitran_2017::calctype::FullRosenkranz
@ FullRosenkranz
lm_hitran_2017::calctype::SDW
@ SDW
lm_hitran_2017::calctype::FullVP
@ FullVP
lm_hitran_2017::compabs
void compabs(CommonBlock &cmn, const Numeric &temp, const Numeric &ptot, const Numeric &xco2, const Numeric &xh2o, const Vector &invcm_grid, const bool mixsdv, const bool mixfull, VectorView absv, VectorView absy, VectorView absw)
Definition linemixing_hitran.cc:1570
lm_hitran_2017::typeLP
constexpr bool typeLP(ModeOfLineMixing x)
Definition linemixing_hitran.h:103
lm_hitran_2017::hitran_lm_eigenvalue_approximation
Tensor4 hitran_lm_eigenvalue_approximation(const AbsorptionLines &band, const Vector &temperatures, const HitranRelaxationMatrixData &hitran, const Numeric P)
Definition linemixing_hitran.cc:2142
lm_hitran_2017::qsdv
void qsdv(const Numeric &sg0, const Numeric &gamd, const Numeric &gam0, const Numeric &gam2, const Numeric &shift0, const Numeric &shift2, const Numeric &sg, Numeric &ls_qsdv_r, Numeric &ls_qsdv_i)
Definition linemixing_hitran.cc:1403
parameters
Parameters parameters
Holds the command line parameters.
Definition parameters.cc:24
physics_funcs.h
This file contains declerations of functions of physical character.
number_density
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Definition physics_funcs.h:51
SpeciesTag
Species::Tag SpeciesTag
Definition species_tags.h:94
Absorption::LineMixing::EquivalentLines
Contains recomputed equivalent lines (sorting is unknown)
Definition linemixing.h:18
Absorption::Lines::NumBroadeners
Index NumBroadeners() const ARTS_NOEXCEPT
Number of broadening species.
Definition absorptionlines.cc:2296
Absorption::Lines::T0
Numeric T0
Reference temperature for all parameters of the lines.
Definition absorptionlines.h:295
Absorption::Lines::population
PopulationType population
Line population distribution.
Definition absorptionlines.h:286
Absorption::Lines::F_mean
Numeric F_mean(Numeric T=0) const noexcept
Mean frequency by weight of line strength.
Definition absorptionlines.cc:2392
Absorption::Lines::lines
Array< SingleLine > lines
A list of individual lines.
Definition absorptionlines.h:310
Absorption::Lines::NumLines
Index NumLines() const noexcept
Number of lines.
Definition absorptionlines.cc:2294
Absorption::Lines::Isotopologue
Species::IsotopeRecord Isotopologue() const noexcept
Isotopologue Index.
Definition absorptionlines.cc:2292
Absorption::Lines::quantumidentity
QuantumIdentifier quantumidentity
Catalog ID.
Definition absorptionlines.h:304
LineShape::ModelParameters
Coefficients and temperature model for SingleSpeciesModel.
Definition lineshapemodel.h:153
Quantum::Number::GlobalState
A logical struct for global quantum numbers with species identifiers.
Definition quantum_numbers.h:1045
Quantum::Number::LocalState
A logical struct for local quantum numbers.
Definition quantum_numbers.h:1010
Quantum::Number::Value
A complete quantum number value with type information.
Definition quantum_numbers.h:716
Species::IsotopeRecord::isotname
std::string_view isotname
A custom name that is unique for this Species type.
Definition isotopologues.h:21
lm_hitran_2017::CommonBlock::Bands
Definition linemixing_hitran.cc:55
lm_hitran_2017::CommonBlock::Bands::lf
std::array< Index, parameters::nBmx > lf
Definition linemixing_hitran.cc:60
lm_hitran_2017::CommonBlock::Bands::Isot
std::array< Index, parameters::nBmx > Isot
Definition linemixing_hitran.cc:57
lm_hitran_2017::CommonBlock::Bands::BandFile
std::array< String, parameters::nBmx > BandFile
Definition linemixing_hitran.cc:61
lm_hitran_2017::CommonBlock::Bands::nBand
Index nBand
Definition linemixing_hitran.cc:56
lm_hitran_2017::CommonBlock::Bands::li
std::array< Index, parameters::nBmx > li
Definition linemixing_hitran.cc:59
lm_hitran_2017::CommonBlock::Bands::nLines
std::array< Index, parameters::nBmx > nLines
Definition linemixing_hitran.cc:58
lm_hitran_2017::CommonBlock::Bfittedp
Definition linemixing_hitran.cc:209
lm_hitran_2017::CommonBlock::Bfittedp::B0pq
Tensor4 B0pq
Definition linemixing_hitran.cc:211
lm_hitran_2017::CommonBlock::Bfittedp::B0pr
Tensor4 B0pr
Definition linemixing_hitran.cc:212
lm_hitran_2017::CommonBlock::Bfittedp::Bfittedp
Bfittedp()
Definition linemixing_hitran.cc:213
lm_hitran_2017::CommonBlock::Bfittedp::B0pp
Tensor4 B0pp
Definition linemixing_hitran.cc:210
lm_hitran_2017::CommonBlock::Bfittedq
Definition linemixing_hitran.cc:219
lm_hitran_2017::CommonBlock::Bfittedq::B0qq
Tensor4 B0qq
Definition linemixing_hitran.cc:221
lm_hitran_2017::CommonBlock::Bfittedq::B0qr
Tensor4 B0qr
Definition linemixing_hitran.cc:222
lm_hitran_2017::CommonBlock::Bfittedq::B0qp
Tensor4 B0qp
Definition linemixing_hitran.cc:220
lm_hitran_2017::CommonBlock::Bfittedq::Bfittedq
Bfittedq()
Definition linemixing_hitran.cc:223
lm_hitran_2017::CommonBlock::Bfittedr
Definition linemixing_hitran.cc:229
lm_hitran_2017::CommonBlock::Bfittedr::B0rr
Tensor4 B0rr
Definition linemixing_hitran.cc:232
lm_hitran_2017::CommonBlock::Bfittedr::B0rp
Tensor4 B0rp
Definition linemixing_hitran.cc:230
lm_hitran_2017::CommonBlock::Bfittedr::B0rq
Tensor4 B0rq
Definition linemixing_hitran.cc:231
lm_hitran_2017::CommonBlock::Bfittedr::Bfittedr
Bfittedr()
Definition linemixing_hitran.cc:233
lm_hitran_2017::CommonBlock::DTGAMAIR
Definition linemixing_hitran.cc:90
lm_hitran_2017::CommonBlock::DTGAMAIR::BHWAIR
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > BHWAIR
Definition linemixing_hitran.cc:91
lm_hitran_2017::CommonBlock::DTGAMAIR::DTGAMAIR
DTGAMAIR()
Definition linemixing_hitran.cc:92
lm_hitran_2017::CommonBlock::DTGAMCO2
Definition linemixing_hitran.cc:106
lm_hitran_2017::CommonBlock::DTGAMCO2::DTGAMCO2
DTGAMCO2()
Definition linemixing_hitran.cc:108
lm_hitran_2017::CommonBlock::DTGAMCO2::BHWSELF
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > BHWSELF
Definition linemixing_hitran.cc:107
lm_hitran_2017::CommonBlock::DTGAMH2O
Definition linemixing_hitran.cc:122
lm_hitran_2017::CommonBlock::DTGAMH2O::BHWH2O
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > BHWH2O
Definition linemixing_hitran.cc:123
lm_hitran_2017::CommonBlock::DTGAMH2O::DTGAMH2O
DTGAMH2O()
Definition linemixing_hitran.cc:124
lm_hitran_2017::CommonBlock::DiagnI
Definition linemixing_hitran.cc:244
lm_hitran_2017::CommonBlock::DiagnI::DiagnI
DiagnI()
Definition linemixing_hitran.cc:246
lm_hitran_2017::CommonBlock::DiagnI::OpI
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > OpI
Definition linemixing_hitran.cc:245
lm_hitran_2017::CommonBlock::DiagnR
Definition linemixing_hitran.cc:239
lm_hitran_2017::CommonBlock::DiagnR::OpR
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > OpR
Definition linemixing_hitran.cc:240
lm_hitran_2017::CommonBlock::DiagnR::DiagnR
DiagnR()
Definition linemixing_hitran.cc:241
lm_hitran_2017::CommonBlock::DipoRigid
Definition linemixing_hitran.cc:69
lm_hitran_2017::CommonBlock::DipoRigid::Dipo0
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > Dipo0
Definition linemixing_hitran.cc:70
lm_hitran_2017::CommonBlock::DipoRigid::DipoRigid
DipoRigid()
Definition linemixing_hitran.cc:71
lm_hitran_2017::CommonBlock::DipoTcm
Definition linemixing_hitran.cc:150
lm_hitran_2017::CommonBlock::DipoTcm::DipoT
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > DipoT
Definition linemixing_hitran.cc:151
lm_hitran_2017::CommonBlock::DipoTcm::DipoTcm
DipoTcm()
Definition linemixing_hitran.cc:152
lm_hitran_2017::CommonBlock::Energy
Definition linemixing_hitran.cc:74
lm_hitran_2017::CommonBlock::Energy::Energy
Energy()
Definition linemixing_hitran.cc:76
lm_hitran_2017::CommonBlock::Energy::E
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > E
Definition linemixing_hitran.cc:75
lm_hitran_2017::CommonBlock::GamSDVT0AIR
Definition linemixing_hitran.cc:84
lm_hitran_2017::CommonBlock::GamSDVT0AIR::HWSDVT0AIR
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > HWSDVT0AIR
Definition linemixing_hitran.cc:85
lm_hitran_2017::CommonBlock::GamSDVT0AIR::GamSDVT0AIR
GamSDVT0AIR()
Definition linemixing_hitran.cc:87
lm_hitran_2017::CommonBlock::GamSDVT0AIR::rHWT0AIR
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > rHWT0AIR
Definition linemixing_hitran.cc:86
lm_hitran_2017::CommonBlock::GamSDVT0CO2
Definition linemixing_hitran.cc:100
lm_hitran_2017::CommonBlock::GamSDVT0CO2::GamSDVT0CO2
GamSDVT0CO2()
Definition linemixing_hitran.cc:103
lm_hitran_2017::CommonBlock::GamSDVT0CO2::rHWT0SELF
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > rHWT0SELF
Definition linemixing_hitran.cc:102
lm_hitran_2017::CommonBlock::GamSDVT0CO2::HWSDVT0SELF
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > HWSDVT0SELF
Definition linemixing_hitran.cc:101
lm_hitran_2017::CommonBlock::GamSDVT0H2O
Definition linemixing_hitran.cc:116
lm_hitran_2017::CommonBlock::GamSDVT0H2O::HWSDVT0H2O
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > HWSDVT0H2O
Definition linemixing_hitran.cc:117
lm_hitran_2017::CommonBlock::GamSDVT0H2O::GamSDVT0H2O
GamSDVT0H2O()
Definition linemixing_hitran.cc:119
lm_hitran_2017::CommonBlock::GamSDVT0H2O::rHWT0H2O
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > rHWT0H2O
Definition linemixing_hitran.cc:118
lm_hitran_2017::CommonBlock::GamT
Definition linemixing_hitran.cc:127
lm_hitran_2017::CommonBlock::GamT::HWT
std::array< Numeric, parameters::nLmx > HWT
Definition linemixing_hitran.cc:128
lm_hitran_2017::CommonBlock::GamT::HWSDV2T
std::array< Numeric, parameters::nLmx > HWSDV2T
Definition linemixing_hitran.cc:129
lm_hitran_2017::CommonBlock::GamVT0AIR
Definition linemixing_hitran.cc:79
lm_hitran_2017::CommonBlock::GamVT0AIR::HWVT0AIR
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > HWVT0AIR
Definition linemixing_hitran.cc:80
lm_hitran_2017::CommonBlock::GamVT0AIR::GamVT0AIR
GamVT0AIR()
Definition linemixing_hitran.cc:81
lm_hitran_2017::CommonBlock::GamVT0CO2
Definition linemixing_hitran.cc:95
lm_hitran_2017::CommonBlock::GamVT0CO2::HWVT0SELF
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > HWVT0SELF
Definition linemixing_hitran.cc:96
lm_hitran_2017::CommonBlock::GamVT0CO2::GamVT0CO2
GamVT0CO2()
Definition linemixing_hitran.cc:97
lm_hitran_2017::CommonBlock::GamVT0H2O
Definition linemixing_hitran.cc:111
lm_hitran_2017::CommonBlock::GamVT0H2O::HWVT0H2O
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > HWVT0H2O
Definition linemixing_hitran.cc:112
lm_hitran_2017::CommonBlock::GamVT0H2O::GamVT0H2O
GamVT0H2O()
Definition linemixing_hitran.cc:113
lm_hitran_2017::CommonBlock::Jfln
Definition linemixing_hitran.cc:160
lm_hitran_2017::CommonBlock::Jfln::Jf
Eigen::Matrix< Index, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > Jf
Definition linemixing_hitran.cc:161
lm_hitran_2017::CommonBlock::Jfln::Jfln
Jfln()
Definition linemixing_hitran.cc:162
lm_hitran_2017::CommonBlock::Jiln
Definition linemixing_hitran.cc:155
lm_hitran_2017::CommonBlock::Jiln::Jiln
Jiln()
Definition linemixing_hitran.cc:157
lm_hitran_2017::CommonBlock::Jiln::Ji
Eigen::Matrix< Index, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > Ji
Definition linemixing_hitran.cc:156
lm_hitran_2017::CommonBlock::LineSg
Definition linemixing_hitran.cc:64
lm_hitran_2017::CommonBlock::LineSg::LineSg
LineSg()
Definition linemixing_hitran.cc:66
lm_hitran_2017::CommonBlock::LineSg::Sig
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > Sig
Definition linemixing_hitran.cc:65
lm_hitran_2017::CommonBlock::PopTrf
Definition linemixing_hitran.cc:145
lm_hitran_2017::CommonBlock::PopTrf::PopuT0
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > PopuT0
Definition linemixing_hitran.cc:146
lm_hitran_2017::CommonBlock::PopTrf::PopTrf
PopTrf()
Definition linemixing_hitran.cc:147
lm_hitran_2017::CommonBlock::PopuT
Definition linemixing_hitran.cc:141
lm_hitran_2017::CommonBlock::PopuT::PopuT
std::array< Numeric, parameters::nLmx > PopuT
Definition linemixing_hitran.cc:142
lm_hitran_2017::CommonBlock::SHIFT0
Definition linemixing_hitran.cc:136
lm_hitran_2017::CommonBlock::SHIFT0::SHIFT0
SHIFT0()
Definition linemixing_hitran.cc:138
lm_hitran_2017::CommonBlock::SHIFT0::shft0
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > shft0
Definition linemixing_hitran.cc:137
lm_hitran_2017::CommonBlock::SHIFT
Definition linemixing_hitran.cc:132
lm_hitran_2017::CommonBlock::SHIFT::shft
std::array< Numeric, parameters::nLmx > shft
Definition linemixing_hitran.cc:133
lm_hitran_2017::CommonBlock::UnusedBandParams
Definition linemixing_hitran.cc:253
lm_hitran_2017::CommonBlock::UnusedBandParams::fl2
std::array< Rational, parameters::nBmx > fl2
Definition linemixing_hitran.cc:263
lm_hitran_2017::CommonBlock::UnusedBandParams::fv1
std::array< Rational, parameters::nBmx > fv1
Definition linemixing_hitran.cc:261
lm_hitran_2017::CommonBlock::UnusedBandParams::iv1
std::array< Rational, parameters::nBmx > iv1
Definition linemixing_hitran.cc:256
lm_hitran_2017::CommonBlock::UnusedBandParams::iv2
std::array< Rational, parameters::nBmx > iv2
Definition linemixing_hitran.cc:257
lm_hitran_2017::CommonBlock::UnusedBandParams::ir
std::array< Rational, parameters::nBmx > ir
Definition linemixing_hitran.cc:260
lm_hitran_2017::CommonBlock::UnusedBandParams::fv2
std::array< Rational, parameters::nBmx > fv2
Definition linemixing_hitran.cc:262
lm_hitran_2017::CommonBlock::UnusedBandParams::UnusedBandParams
UnusedBandParams()
Definition linemixing_hitran.cc:266
lm_hitran_2017::CommonBlock::UnusedBandParams::iv3
std::array< Rational, parameters::nBmx > iv3
Definition linemixing_hitran.cc:259
lm_hitran_2017::CommonBlock::UnusedBandParams::intens
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > intens
Definition linemixing_hitran.cc:254
lm_hitran_2017::CommonBlock::UnusedBandParams::fr
std::array< Rational, parameters::nBmx > fr
Definition linemixing_hitran.cc:265
lm_hitran_2017::CommonBlock::UnusedBandParams::fv3
std::array< Rational, parameters::nBmx > fv3
Definition linemixing_hitran.cc:264
lm_hitran_2017::CommonBlock::UnusedBandParams::eina
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > eina
Definition linemixing_hitran.cc:255
lm_hitran_2017::CommonBlock::UnusedBandParams::il2
std::array< Rational, parameters::nBmx > il2
Definition linemixing_hitran.cc:258
lm_hitran_2017::CommonBlock::Wfittedp
Definition linemixing_hitran.cc:179
lm_hitran_2017::CommonBlock::Wfittedp::Wfittedp
Wfittedp()
Definition linemixing_hitran.cc:183
lm_hitran_2017::CommonBlock::Wfittedp::W0pr
Tensor4 W0pr
Definition linemixing_hitran.cc:182
lm_hitran_2017::CommonBlock::Wfittedp::W0pp
Tensor4 W0pp
Definition linemixing_hitran.cc:180
lm_hitran_2017::CommonBlock::Wfittedp::W0pq
Tensor4 W0pq
Definition linemixing_hitran.cc:181
lm_hitran_2017::CommonBlock::Wfittedq
Definition linemixing_hitran.cc:189
lm_hitran_2017::CommonBlock::Wfittedq::W0qr
Tensor4 W0qr
Definition linemixing_hitran.cc:192
lm_hitran_2017::CommonBlock::Wfittedq::W0qp
Tensor4 W0qp
Definition linemixing_hitran.cc:190
lm_hitran_2017::CommonBlock::Wfittedq::Wfittedq
Wfittedq()
Definition linemixing_hitran.cc:193
lm_hitran_2017::CommonBlock::Wfittedq::W0qq
Tensor4 W0qq
Definition linemixing_hitran.cc:191
lm_hitran_2017::CommonBlock::Wfittedr
Definition linemixing_hitran.cc:199
lm_hitran_2017::CommonBlock::Wfittedr::W0rp
Tensor4 W0rp
Definition linemixing_hitran.cc:200
lm_hitran_2017::CommonBlock::Wfittedr::Wfittedr
Wfittedr()
Definition linemixing_hitran.cc:203
lm_hitran_2017::CommonBlock::Wfittedr::W0rr
Tensor4 W0rr
Definition linemixing_hitran.cc:202
lm_hitran_2017::CommonBlock::Wfittedr::W0rq
Tensor4 W0rq
Definition linemixing_hitran.cc:201
lm_hitran_2017::CommonBlock::Wmatrix
Definition linemixing_hitran.cc:174
lm_hitran_2017::CommonBlock::Wmatrix::Wmatrix
Wmatrix()
Definition linemixing_hitran.cc:176
lm_hitran_2017::CommonBlock::Wmatrix::W
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > W
Definition linemixing_hitran.cc:175
lm_hitran_2017::CommonBlock::YLT
Definition linemixing_hitran.cc:249
lm_hitran_2017::CommonBlock::YLT::YT
std::array< Numeric, parameters::nLmx > YT
Definition linemixing_hitran.cc:250
lm_hitran_2017::CommonBlock::Zaa
Definition linemixing_hitran.cc:169
lm_hitran_2017::CommonBlock::Zaa::Zaa
Zaa() noexcept
Definition linemixing_hitran.cc:171
lm_hitran_2017::CommonBlock::Zaa::ZA
ComplexVector ZA
Definition linemixing_hitran.cc:170
lm_hitran_2017::CommonBlock::Zss
Definition linemixing_hitran.cc:165
lm_hitran_2017::CommonBlock::Zss::ZS
std::array< Complex, parameters::nLmx > ZS
Definition linemixing_hitran.cc:166
lm_hitran_2017::CommonBlock
Definition linemixing_hitran.cc:54
lm_hitran_2017::CommonBlock::Zaa
struct lm_hitran_2017::CommonBlock::Zaa Zaa
lm_hitran_2017::CommonBlock::DTGAMAIR
struct lm_hitran_2017::CommonBlock::DTGAMAIR DTGAMAIR
lm_hitran_2017::CommonBlock::Wfittedp
struct lm_hitran_2017::CommonBlock::Wfittedp Wfittedp
lm_hitran_2017::CommonBlock::GamVT0CO2
struct lm_hitran_2017::CommonBlock::GamVT0CO2 GamVT0CO2
lm_hitran_2017::CommonBlock::Energy
struct lm_hitran_2017::CommonBlock::Energy Energy
lm_hitran_2017::CommonBlock::GamSDVT0CO2
struct lm_hitran_2017::CommonBlock::GamSDVT0CO2 GamSDVT0CO2
lm_hitran_2017::CommonBlock::Wfittedr
struct lm_hitran_2017::CommonBlock::Wfittedr Wfittedr
lm_hitran_2017::CommonBlock::Bfittedq
struct lm_hitran_2017::CommonBlock::Bfittedq Bfittedq
lm_hitran_2017::CommonBlock::PopTrf
struct lm_hitran_2017::CommonBlock::PopTrf PopTrf
lm_hitran_2017::CommonBlock::Bfittedr
struct lm_hitran_2017::CommonBlock::Bfittedr Bfittedr
lm_hitran_2017::CommonBlock::Wfittedq
struct lm_hitran_2017::CommonBlock::Wfittedq Wfittedq
lm_hitran_2017::CommonBlock::YLT
struct lm_hitran_2017::CommonBlock::YLT YLT
lm_hitran_2017::CommonBlock::Bands
struct lm_hitran_2017::CommonBlock::Bands Bands
lm_hitran_2017::CommonBlock::GamT
struct lm_hitran_2017::CommonBlock::GamT GamT
lm_hitran_2017::CommonBlock::SHIFT0
struct lm_hitran_2017::CommonBlock::SHIFT0 SHIFT0
lm_hitran_2017::CommonBlock::DipoRigid
struct lm_hitran_2017::CommonBlock::DipoRigid DipoRigid
lm_hitran_2017::CommonBlock::DTGAMCO2
struct lm_hitran_2017::CommonBlock::DTGAMCO2 DTGAMCO2
lm_hitran_2017::CommonBlock::SHIFT
struct lm_hitran_2017::CommonBlock::SHIFT SHIFT
lm_hitran_2017::CommonBlock::Wmatrix
struct lm_hitran_2017::CommonBlock::Wmatrix Wmatrix
lm_hitran_2017::CommonBlock::DiagnI
struct lm_hitran_2017::CommonBlock::DiagnI DiagnI
lm_hitran_2017::CommonBlock::UnusedBandParams
struct lm_hitran_2017::CommonBlock::UnusedBandParams UnusedBandParams
lm_hitran_2017::CommonBlock::DipoTcm
struct lm_hitran_2017::CommonBlock::DipoTcm DipoTcm
lm_hitran_2017::CommonBlock::Jfln
struct lm_hitran_2017::CommonBlock::Jfln Jfln
lm_hitran_2017::CommonBlock::GamSDVT0AIR
struct lm_hitran_2017::CommonBlock::GamSDVT0AIR GamSDVT0AIR
lm_hitran_2017::CommonBlock::DiagnR
struct lm_hitran_2017::CommonBlock::DiagnR DiagnR
lm_hitran_2017::CommonBlock::Bfittedp
struct lm_hitran_2017::CommonBlock::Bfittedp Bfittedp
lm_hitran_2017::CommonBlock::Zss
struct lm_hitran_2017::CommonBlock::Zss Zss
lm_hitran_2017::CommonBlock::GamVT0H2O
struct lm_hitran_2017::CommonBlock::GamVT0H2O GamVT0H2O
lm_hitran_2017::CommonBlock::GamVT0AIR
struct lm_hitran_2017::CommonBlock::GamVT0AIR GamVT0AIR
lm_hitran_2017::CommonBlock::DTGAMH2O
struct lm_hitran_2017::CommonBlock::DTGAMH2O DTGAMH2O
lm_hitran_2017::CommonBlock::LineSg
struct lm_hitran_2017::CommonBlock::LineSg LineSg
lm_hitran_2017::CommonBlock::GamSDVT0H2O
struct lm_hitran_2017::CommonBlock::GamSDVT0H2O GamSDVT0H2O
lm_hitran_2017::CommonBlock::PopuT
struct lm_hitran_2017::CommonBlock::PopuT PopuT
lm_hitran_2017::CommonBlock::Jiln
struct lm_hitran_2017::CommonBlock::Jiln Jiln
lm_hitran_2017::ConvTPOut
Definition linemixing_hitran.cc:955
lm_hitran_2017::ConvTPOut::operator=
ConvTPOut & operator=(ConvTPOut &&)=default
lm_hitran_2017::ConvTPOut::ConvTPOut
ConvTPOut(ConvTPOut &&)=default
lm_hitran_2017::ConvTPOut::operator=
ConvTPOut & operator=(const ConvTPOut &)=delete
lm_hitran_2017::ConvTPOut::shft
Vector shft
Definition linemixing_hitran.cc:956
lm_hitran_2017::ConvTPOut::pop
Vector pop
Definition linemixing_hitran.cc:956
lm_hitran_2017::ConvTPOut::ConvTPOut
ConvTPOut(Index n=0) noexcept
Definition linemixing_hitran.cc:959
lm_hitran_2017::ConvTPOut::dip
Vector dip
Definition linemixing_hitran.cc:956
lm_hitran_2017::ConvTPOut::hwt
Vector hwt
Definition linemixing_hitran.cc:956
lm_hitran_2017::ConvTPOut::eqv
EqvLinesOut eqv
Definition linemixing_hitran.cc:958
lm_hitran_2017::ConvTPOut::Y
Vector Y
Definition linemixing_hitran.cc:956
lm_hitran_2017::ConvTPOut::f0
Vector f0
Definition linemixing_hitran.cc:956
lm_hitran_2017::ConvTPOut::hwt2
Vector hwt2
Definition linemixing_hitran.cc:956
lm_hitran_2017::ConvTPOut::W
ComplexMatrix W
Definition linemixing_hitran.cc:957
lm_hitran_2017::ConvTPOut::ConvTPOut
ConvTPOut(const ConvTPOut &)=delete
lm_hitran_2017::Sorter
Definition linemixing_hitran.cc:967
lm_hitran_2017::Sorter::dop0
Vector dop0
Definition linemixing_hitran.cc:968
lm_hitran_2017::Sorter::Ju
std::vector< Rational > Ju
Definition linemixing_hitran.cc:970
lm_hitran_2017::Sorter::Ji
std::vector< Rational > Ji
Definition linemixing_hitran.cc:969
w
#define w
a
#define a
b
#define b