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
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;
82} GamVT0AIR;
83
85 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWSDVT0AIR;
86 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> rHWT0AIR;
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;
98} GamVT0CO2;
99
101 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWSDVT0SELF;
102 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> rHWT0SELF;
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;
114} GamVT0H2O;
115
117 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWSDVT0H2O;
118 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> rHWT0H2O;
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;
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;
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;
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;
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;
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;
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
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;
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
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
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
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;
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;
964};
965
966
967struct Sorter {
968 Vector dop0;
969 std::vector<Rational> Ji;
970 std::vector<Rational> Ju;
971};
972
973
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
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
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
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
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
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;
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
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
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;
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
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,
2091 cmn.DTGAMAIR.BHWAIR(j, i)};
2092 const LineShape::ModelParameters G0_sdvp_h2o{LineShape::TemperatureModel::T1,
2094 cmn.DTGAMH2O.BHWH2O(j, i)};
2095 const LineShape::ModelParameters G0_sdvp_co2{LineShape::TemperatureModel::T1,
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,
2107 cmn.DTGAMAIR.BHWAIR(j, i)};
2108 const LineShape::ModelParameters G0_vp_h2o{LineShape::TemperatureModel::T1,
2110 cmn.DTGAMH2O.BHWH2O(j, i)};
2111 const LineShape::ModelParameters G0_vp_co2{LineShape::TemperatureModel::T1,
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)),
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
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
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
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
2222 not is_sorted(temperatures),
2223 "The temperature list [", temperatures, "] K\n"
2224 "must be fully sorted from low to high"
2225 )
2226
2228 hitran_lm_eigenvalue_approximation(band, temperatures, hitran, P0),
2229 temperatures, P0, ord)) {
2230 ARTS_USER_ERROR("Bad eigenvalue adaptation")
2231 }
2232}
2233
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
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
This can be used to make arrays out of anything.
Definition array.h:31
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
Main line shape model class.
Compute the line shape parameters for a single broadening species.
void set(Value v)
Sets the value if it exists or adds it otherwise.
Main Zeeman Model.
Definition zeemandata.h:279
Index nelem() const
Definition mystring.h:172
#define ARTS_ASSERT(condition,...)
Definition debug.h:86
#define ARTS_USER_ERROR(...)
Definition debug.h:153
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
Namespace and functions to deal with HITRAN linemixing.
Numeric boltzman_factor(Numeric T, Numeric E0)
Computes exp(- E0/kT)
Numeric boltzman_ratio(const Numeric &T, const Numeric &T0, const Numeric &E0)
Computes exp(E0/c (T - T0) / (T * T0))
Numeric single_partition_function(const Numeric &T, const Species::IsotopeRecord &ir)
Computes the partition function at one temperature.
Definition linescaling.cc:4
Contains the line shape namespace.
my_basic_string< char > String
The String type for ARTS.
Definition mystring.h:199
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...
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.
Numeric reduced_rovibrational_dipole(Rational Jf, Rational Ji, Rational lf, Rational li, Rational k=Rational(1))
Compute the reduced rovibrational dipole moment.
constexpr Numeric sqrt_ln_2
Square root of natural logarithm of 2.
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric k
Boltzmann constant convenience name [J/K].
constexpr Numeric sqrt_pi
Square root of pi.
constexpr Numeric inv_pi
Inverse of pi.
constexpr Numeric h
Planck constant convenience name [J s].
constexpr auto kaycm2joule(auto x) noexcept
Conversion from cm-1 to Joule.
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.
constexpr auto kaycm_per_atm2hz_per_pa(auto x) noexcept
Conversion from cm-1 per atmosphere to Hz per Pascal.
constexpr auto freq2kaycm(auto x) noexcept
Conversion from Hz to Kayser wavenumber.
constexpr auto kaycm2freq(auto x) noexcept
Conversion from Kayser wavenumber to Hz.
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.
constexpr auto pow3(auto x) noexcept
power of three
constexpr auto pow2(auto x) noexcept
power of two
void eqvlines(CommonBlock &cmn, const Index &iband, const Index &n, const Numeric &sigmoy)
void hitran_lm_eigenvalue_adaptation(AbsorptionLines &band, const Vector &temperatures, const HitranRelaxationMatrixData &hitran, const Numeric P0, const Index ord)
Complex qsdv_si(const Numeric F0, const Numeric gamd, const Numeric gam0, const Numeric gam2, const Numeric shift0, const Numeric shift2, const Numeric f)
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.
void calcw(CommonBlock &cmn, const Index &n, const Index &iband, const Numeric &temp)
Tensor5 hitran_lm_eigenvalue_adaptation_test(const AbsorptionLines &band, const Vector &temperatures, const HitranRelaxationMatrixData &hitran, const Vector &pressures)
constexpr bool typeFull(ModeOfLineMixing x)
Sorter sorter_calcw(ConvTPOut &out, const AbsorptionLines &band, const bool at_t0)
void qt_co2(const Numeric &t, const Index &iso, Numeric &gsi, Numeric &qt)
void readw(CommonBlock &cmn, const String &basedir="data_new/")
Vector compute(const Numeric p, const Numeric t, const Numeric xco2, const Numeric xh2o, const Vector &invcm_grid, const Numeric stotmax, const calctype type)
void detband(CommonBlock &cmn, const Numeric &sgminr, const Numeric &sgmaxr, const Numeric &stotmax, const String &basedir="data_new/")
Rational toRationalSum(char a, char b=' ')
Numeric atob(const Numeric &aa, const std::array< Numeric, NT > &a, const std::array< Numeric, NT > &b)
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)
void readlines(CommonBlock &cmn, const String &basedir="data_new/")
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)
constexpr bool typeLP(ModeOfLineMixing x)
Tensor4 hitran_lm_eigenvalue_approximation(const AbsorptionLines &band, const Vector &temperatures, const HitranRelaxationMatrixData &hitran, const Numeric P)
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)
Parameters parameters
Holds the command line parameters.
Definition parameters.cc:24
This file contains declerations of functions of physical character.
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Species::Tag SpeciesTag
Contains recomputed equivalent lines (sorting is unknown)
Definition linemixing.h:18
Index NumBroadeners() const ARTS_NOEXCEPT
Number of broadening species.
Numeric T0
Reference temperature for all parameters of the lines.
PopulationType population
Line population distribution.
Numeric F_mean(Numeric T=0) const noexcept
Mean frequency by weight of line strength.
Array< SingleLine > lines
A list of individual lines.
Index NumLines() const noexcept
Number of lines.
Species::IsotopeRecord Isotopologue() const noexcept
Isotopologue Index.
QuantumIdentifier quantumidentity
Catalog ID.
Coefficients and temperature model for SingleSpeciesModel.
A logical struct for global quantum numbers with species identifiers.
A logical struct for local quantum numbers.
A complete quantum number value with type information.
std::string_view isotname
A custom name that is unique for this Species type.
std::array< Index, parameters::nBmx > lf
std::array< Index, parameters::nBmx > Isot
std::array< String, parameters::nBmx > BandFile
std::array< Index, parameters::nBmx > li
std::array< Index, parameters::nBmx > nLines
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > BHWAIR
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > BHWSELF
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > BHWH2O
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > OpI
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > OpR
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > Dipo0
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > DipoT
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > E
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > HWSDVT0AIR
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > rHWT0AIR
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > rHWT0SELF
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > HWSDVT0SELF
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > HWSDVT0H2O
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > rHWT0H2O
std::array< Numeric, parameters::nLmx > HWT
std::array< Numeric, parameters::nLmx > HWSDV2T
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > HWVT0AIR
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > HWVT0SELF
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > HWVT0H2O
Eigen::Matrix< Index, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > Jf
Eigen::Matrix< Index, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > Ji
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > Sig
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > PopuT0
std::array< Numeric, parameters::nLmx > PopuT
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > shft0
std::array< Numeric, parameters::nLmx > shft
std::array< Rational, parameters::nBmx > fl2
std::array< Rational, parameters::nBmx > fv1
std::array< Rational, parameters::nBmx > iv1
std::array< Rational, parameters::nBmx > iv2
std::array< Rational, parameters::nBmx > ir
std::array< Rational, parameters::nBmx > fv2
std::array< Rational, parameters::nBmx > iv3
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > intens
std::array< Rational, parameters::nBmx > fr
std::array< Rational, parameters::nBmx > fv3
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > eina
std::array< Rational, parameters::nBmx > il2
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > W
std::array< Numeric, parameters::nLmx > YT
std::array< Complex, parameters::nLmx > ZS
struct lm_hitran_2017::CommonBlock::Zaa Zaa
struct lm_hitran_2017::CommonBlock::DTGAMAIR DTGAMAIR
struct lm_hitran_2017::CommonBlock::Wfittedp Wfittedp
struct lm_hitran_2017::CommonBlock::GamVT0CO2 GamVT0CO2
struct lm_hitran_2017::CommonBlock::Energy Energy
struct lm_hitran_2017::CommonBlock::GamSDVT0CO2 GamSDVT0CO2
struct lm_hitran_2017::CommonBlock::Wfittedr Wfittedr
struct lm_hitran_2017::CommonBlock::Bfittedq Bfittedq
struct lm_hitran_2017::CommonBlock::PopTrf PopTrf
struct lm_hitran_2017::CommonBlock::Bfittedr Bfittedr
struct lm_hitran_2017::CommonBlock::Wfittedq Wfittedq
struct lm_hitran_2017::CommonBlock::YLT YLT
struct lm_hitran_2017::CommonBlock::Bands Bands
struct lm_hitran_2017::CommonBlock::GamT GamT
struct lm_hitran_2017::CommonBlock::SHIFT0 SHIFT0
struct lm_hitran_2017::CommonBlock::DipoRigid DipoRigid
struct lm_hitran_2017::CommonBlock::DTGAMCO2 DTGAMCO2
struct lm_hitran_2017::CommonBlock::SHIFT SHIFT
struct lm_hitran_2017::CommonBlock::Wmatrix Wmatrix
struct lm_hitran_2017::CommonBlock::DiagnI DiagnI
struct lm_hitran_2017::CommonBlock::UnusedBandParams UnusedBandParams
struct lm_hitran_2017::CommonBlock::DipoTcm DipoTcm
struct lm_hitran_2017::CommonBlock::Jfln Jfln
struct lm_hitran_2017::CommonBlock::GamSDVT0AIR GamSDVT0AIR
struct lm_hitran_2017::CommonBlock::DiagnR DiagnR
struct lm_hitran_2017::CommonBlock::Bfittedp Bfittedp
struct lm_hitran_2017::CommonBlock::Zss Zss
struct lm_hitran_2017::CommonBlock::GamVT0H2O GamVT0H2O
struct lm_hitran_2017::CommonBlock::GamVT0AIR GamVT0AIR
struct lm_hitran_2017::CommonBlock::DTGAMH2O DTGAMH2O
struct lm_hitran_2017::CommonBlock::LineSg LineSg
struct lm_hitran_2017::CommonBlock::GamSDVT0H2O GamSDVT0H2O
struct lm_hitran_2017::CommonBlock::PopuT PopuT
struct lm_hitran_2017::CommonBlock::Jiln Jiln
ConvTPOut & operator=(ConvTPOut &&)=default
ConvTPOut(ConvTPOut &&)=default
ConvTPOut & operator=(const ConvTPOut &)=delete
ConvTPOut(Index n=0) noexcept
ConvTPOut(const ConvTPOut &)=delete
std::vector< Rational > Ju
std::vector< Rational > Ji
#define w
#define a
#define b