10#include <Faddeeva/Faddeeva.hh>
19#include "matpack_data.h"
20#include "matpack_math.h"
23#pragma GCC diagnostic push
24#pragma GCC diagnostic ignored "-Wshadow"
25#pragma GCC diagnostic ignored "-Wconversion"
26#pragma GCC diagnostic ignored "-Wfloat-conversion"
28#pragma GCC diagnostic ignored "-Wdeprecated-copy-with-dtor"
29#pragma GCC diagnostic ignored "-Wdeprecated-copy-with-user-provided-dtor"
31#pragma GCC diagnostic ignored "-Wclass-memaccess"
36#pragma GCC diagnostic pop
40 static constexpr Index nBmx=7'000;
41 static constexpr Index nLmx=700;
42 static constexpr Index Nlifmax=10;
43 static constexpr Index Jmax=131;
45 static constexpr Numeric Ct=1.4387686e0;
46 static constexpr Numeric T0=296;
47 static constexpr Numeric CtGamD=1.1325e-08;
48 static constexpr Numeric aMolAtm=7.33889e+21;
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);
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;
65 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
Sig;
70 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
Dipo0;
75 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
E;
80 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
HWVT0AIR;
85 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
HWSDVT0AIR;
86 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
rHWT0AIR;
91 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
BHWAIR;
96 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
HWVT0SELF;
101 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
HWSDVT0SELF;
102 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
rHWT0SELF;
107 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
BHWSELF;
112 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
HWVT0H2O;
117 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
HWSDVT0H2O;
118 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
rHWT0H2O;
123 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
BHWH2O;
128 std::array<Numeric, parameters::nLmx>
HWT;
129 std::array<Numeric, parameters::nLmx>
HWSDV2T;
133 std::array<Numeric, parameters::nLmx>
shft;
137 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
shft0;
142 std::array<Numeric, parameters::nLmx>
PopuT;
146 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
PopuT0;
151 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
DipoT;
156 Eigen::Matrix<Index, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
Ji;
161 Eigen::Matrix<Index, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
Jf;
166 std::array<Complex, parameters::nLmx>
ZS;
171 Zaa() noexcept :
ZA(parameters::nLmx) {}
175 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
W;
240 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
OpR;
245 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
OpI;
250 std::array<Numeric, parameters::nLmx>
YT;
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;
274 if (
b ==
' ' and
a ==
' ')
277 return Rational(
a-
'0');
279 return Rational(
b-
'0');
280 return Rational(10*(
a-
'0') +
b-
'0');
288 for (Index iband=0; iband<cmn.
Bands.
nBand; iband++) {
289 std::ifstream fortranfile;
291 fortranfile.open(fname.c_str());
293 if (fortranfile.is_open()) {
295 getline(fortranfile, line);
298 while (fortranfile.good()) {
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;
307 "%c%c" "%1" PRId64
"%12lf" "%10lf"
308 "%10lf" "%5lf" "%5lf" "%4lf"
309 "%5lf" "%5lf" "%4lf" "%10lf"
312 "%c%c" "%c%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"
335 &iv11, &iv12, &iv21, &iv22, &il21, &il22, &iv31, &iv32, &ir1, &fv32, &fr1,
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,
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,
348 getline(fortranfile, line);
374 "Bad read, bands do not have the same global quantum numbers...");
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');
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);
393 cmn.
Jfln.
Jf(nliner, iband) = cmn.
Jiln.
Ji(nliner, iband) + 1;
410 const std::array<Numeric, NT>&
a,
411 const std::array<Numeric, NT>&
b)
413 for (
size_t i=1; i<NT; i++) {
415 if (i < 2 or i == NT) {
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;
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);
434 return A0*
b[j-2] + A1*
b[j-1] + A2*
b[j];
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;
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);
461 return A0*
b[j-2] + A1*
b[j-1] + A2*
b[j] + A3*
b[j+1];
465 return std::numeric_limits<Numeric>::quiet_NaN();
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.,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
776 const auto& q = qoft[iso-1];
778 if (t < 70 or t > 3000)
781 qt =
atob(t, tdat, q);
789 for (Index i=0; i<n; i++) {
791 for (Index j=0; j<n; j++) {
797 for (Index i=0; i<n; i++) {
804 for (Index i=0; i<n; i++) {
808 for (Index i=0; i<n; i++) {
809 for (Index j=i+1; j<n; j++) {
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));
823 std::swap(s[i], s[j]);
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]);
833 for (Index i=0; i<n; i++) {
836 jji=cmn.
Jiln.
Ji(i, iband);
837 jjf=cmn.
Jfln.
Jf(i, iband);
839 jji=cmn.
Jfln.
Jf(i, iband);
840 jjf=cmn.
Jiln.
Ji(i, iband);
843 for (Index j=0; j<n; j++) {
846 jjip=cmn.
Jiln.
Ji(j, iband);
847 jjfp=cmn.
Jfln.
Jf(j, iband);
849 jjip=cmn.
Jfln.
Jf(j, iband);
850 jjfp=cmn.
Jiln.
Ji(j, iband);
853 if (jjip > jji)
continue;
858 if (jji > jjf and jjip > jjfp) {
861 }
else if (jji > jjf and jjip == jjfp) {
864 }
else if (jji > jjf and jjip < jjfp) {
867 }
else if (jji < jjf and jjip > jjfp) {
870 }
else if (jji < jjf and jjip == jjfp) {
873 }
else if (jji < jjf and jjip < jjfp) {
876 }
else if (jji == jjf and jjip > jjfp) {
879 }
else if (jji == jjf and jjip == jjfp) {
882 }
else if (jji == jjf and jjip < jjfp) {
887 const Numeric ycal = std::exp(w0 - b0*dlgt0t);
894 for (Index i=0; i<n; i++) {
895 for (Index j=0; j<n; j++) {
903 for (Index i=0; i<n; i++) {
908 for (Index i=0; i<n; i++) {
912 for (Index j=0; j<n; j++) {
922 for (Index j=i+1; j<n; j++) {
934 for (Index i=0; i<n; i++) {
936 for (Index j=0; j<n; j++) {
937 if (i == j)
continue;
942 if (std::abs(deltasig) < 1e-4 ) deltasig = 1e-4 ;
947 cmn.
YLT.
YT[i] = sum0;
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) {}
969 std::vector<Rational>
Ji;
970 std::vector<Rational>
Ju;
978 Vector& g0 = out.
hwt;
979 Vector& g2 = out.
hwt2;
980 Vector& d0 = out.
shft;
982 Vector& pop = out.
pop;
983 Vector& dip = out.
dip;
986 const Index n=Y.nelem();
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();
1001 for (Index i=0; i<n; i++) {
1003 s[i] = band.
lines[i].I0;
1009 for (Index i=0; i<n; i++) {
1010 for (Index j=i+1; j<n; j++) {
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]);
1029 return {dip0, Ji, Jf};
1037 const bool at_t0=
false) {
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());
1046 const Index n=Y.nelem();
1054 const auto [dip0, Ji, Jf] =
sorter_calcw(out, band, at_t0);
1056 if (li > 8 or abs(li - lf) > 1) {
1057 for (Index i=0; i<n; i++) {
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);
1066 for (Index i=0; i<n; i++) {
1076 for (Index j=0; j<n; j++) {
1077 Rational jjip, jjfp;
1086 if (jjip > jji)
continue;
1093 (abs(Ji[i] - Ji[j]) % 2) not_eq 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());
1128 const Numeric ycal = std::exp(w0 - b0*dlgt0t);
1130 W(i, j) = ycal * pop[i] / pop[j];
1135 for (Index i=0; i<n; i++) {
1136 for (Index j=0; j<n; j++) {
1138 W(i, j) = - std::abs(W(i, j));
1144 for (Index i=0; i<n; i++) {
1149 for (Index i=0; i<n; i++) {
1150 Numeric sumlw = 0.0;
1151 Numeric sumup = 0.0;
1153 for (Index j=0; j<n; j++) {
1160 (abs(Ji[i] - Ji[j]) % 2) not_eq 0;
1165 sumlw += std::abs(dip0[j]) * W(j, i);
1167 sumup += std::abs(dip0[j]) * W(j, i);
1171 for (Index j=i+1; j<n; j++) {
1176 W(j, i) *= - sumup / sumlw;
1177 W(i, j) = W(j, i) * pop[i] / pop[j];
1183 for (Index i=0; i<n; i++) {
1185 for (Index j=0; j<n; j++) {
1186 if (i == j)
continue;
1193 (abs(Ji[i] - Ji[j]) % 2) not_eq 0;
1197 Numeric deltasig = f0[i] - f0[j];
1201 sum0 += 2 * std::abs(dip[j]) / std::abs(dip[i]) * W(j, i) / deltasig;
1212 const Numeric& sigmoy)
1214 ComplexMatrix zop(n, n), zvec(n, n);
1215 ComplexMatrixView inv_zvec = zop;
1216 ComplexVectorView zval = cmn.
Zaa.
ZA[Range(0, n)];
1217 for (Index i=0; i<n; i++) {
1218 for (Index j=0; j<n; j++) {
1224 diagonalize(zvec, zval, zop);
1225 inv(inv_zvec, zvec);
1231 for (Index i=0; i<n; i++) {
1234 for (Index j=0; j<n; j++) {
1245 const Numeric& fmean)
1248 const Index n = pop.nelem();
1251 ComplexMatrix zvec(n, n);
1255 diagonalize(zvec, out.
val, W);
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);
1266 for (Index i=0; i<n; i++) {
1268 for (Index j=0; j<n; j++) {
1269 z += pop[j] * dip[j] * zvec(i, j);
1284 const Index& nlinec,
1285 const Numeric& xh2o,
1286 const Numeric& xco2,
1287 const Numeric& temp,
1288 const Numeric& ptot,
1290 Numeric& gamd_gam0_mx,
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;
1300 const Numeric sqrtm = std::sqrt(temp / parameters::aMass[isotc-1]);
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));
1323 sigmoy += cmn.
LineSg.
Sig(iline,iband) * wgt;
1327 calcw(cmn, nlinec, iband, temp);
1330 for (Index iline=0; iline<nlinec; iline++) {
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]);
1339 for (Index iline=0; iline<nlinec; iline++) {
1340 for (Index ilinep=0; ilinep<nlinec; ilinep++) {
1342 if (iline == ilinep) {
1350 eqvlines(cmn, iband, nlinec, sigmoy);
1364 const Numeric ratiopart = QT0 / QT;
1369 for (Index i=0; i<n; i++) {
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);
1377 out.
hwt2[i] = band.
lines[i].lineshape.G2(T, band.
T0, 1, vmrs);
1382 calcw(out, hitran, band, T);
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;
1404 const Numeric& gamd,
1405 const Numeric& gam0,
1406 const Numeric& gam2,
1407 const Numeric& shift0,
1408 const Numeric& shift2,
1414 Complex x , y, csqrty, z1, z2, w1, w2, aterm;
1415 Numeric xz1, yz1, xz2, yz2;
1420 const Complex iz = Complex(0, 1);
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;
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;
1432 if (std::abs(y) <= 1.e-15*std::abs(x))
goto label140;
1433 z1 = std::sqrt(x + y) - csqrty;
1434 z2 = z1 + 2 * csqrty;
1441 w1 = Faddeeva::w(Complex(xz1 ,yz1));
1442 w2 = Faddeeva::w(Complex(xz2 ,yz2));
1443 aterm = rpi * cte * (w1 - w2);
1447 z1 = (iz*(sg0-sg)+c0t)*cte;
1451 w1 = Faddeeva::w(Complex(xz1 ,yz1));
1452 aterm = rpi*cte* w1;
1457 z1 = (iz*(sg0-sg)+c0t)*cte;
1458 z2 = std::sqrt(x+y)+csqrty;
1464 w1 = Faddeeva::w(Complex(xz1 ,yz1));
1465 w2 = Faddeeva::w(Complex(xz2 ,yz2));
1466 aterm = rpi * cte * (w1-w2);
1471 if (std::abs(std::sqrt(x)) <= 4.e3) {
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);
1477 aterm = (1 / c2t) * (1 / x - 1.5 / pow2(x));
1482 const Complex ls_qsdv = (1 / pi) * aterm;
1484 ls_qsdv_r = ls_qsdv.real();
1485 ls_qsdv_i = ls_qsdv.imag();
1494 const Numeric shift0,
1495 const Numeric shift2,
1499 Complex x , y, csqrty, z1, z2, w1, w2, aterm;
1500 Numeric xz1, yz1, xz2, yz2;
1505 const Complex iz = Complex(0, 1);
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;
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;
1526 w1 = Faddeeva::w(Complex(xz1 ,yz1));
1527 w2 = Faddeeva::w(Complex(xz2 ,yz2));
1528 aterm = rpi * cte * (w1 - w2);
1532 z1 = (iz*(F0 - f) + c0t)*cte;
1536 w1 = Faddeeva::w(Complex(xz1 ,yz1));
1537 aterm = rpi*cte* w1;
1542 z1 = (iz*(F0 - f) + c0t)*cte;
1543 z2 = std::sqrt(x+y)+csqrty;
1549 w1 = Faddeeva::w(Complex(xz1 ,yz1));
1550 w2 = Faddeeva::w(Complex(xz2 ,yz2));
1551 aterm = rpi * cte * (w1-w2);
1556 if (std::abs(std::sqrt(x)) <= 4.e3) {
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);
1562 aterm = (1 / c2t) * (1 / x - 1.5 / pow2(x));
1567 return (1 / pi) * aterm;
1572 const Numeric& temp,
1573 const Numeric& ptot,
1574 const Numeric& xco2,
1575 const Numeric& xh2o,
1576 const Vector& invcm_grid,
1586 constexpr Numeric rdmult = 30;
1589 const Index nsig = invcm_grid.nelem();
1598 const Numeric dens = xco2 * ptot * parameters::aMolAtm / temp;
1600 constexpr Numeric u_sqln2pi = 1 / sq_ln2pi;
1602 for (Index iband=0; iband<cmn.
Bands.
nBand; iband++) {
1606 xh2o, xco2, temp, ptot, sigmoy, gamdmx, mixfull, mixsdv);
1607 const Numeric sqrtm = std::sqrt(temp / parameters::aMass[cmn.
Bands.
Isot[iband]-1]);
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));
1612 for (Index isig=0; isig<nsig; isig++) {
1613 const Numeric sigc = invcm_grid[isig];
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;
1619 const Numeric dsigc = sigc - cmn.
LineSg.
Sig(iline, iband) - cmn.
SHIFT.
shft[iline] * ptot;
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))
1631 const Numeric order2i = dsigc / (pow2(cmn.
GamT.
HWT[iline]) + pow2(dsigc)) +
1635 absv[isig] += popudipo * u_sqln2pi * u_pi * order2r;
1637 absy[isig] += popudipo * u_sqln2pi * u_pi * (order2r + cmn.
YLT.
YT[iline]*order2i);
1639 const Numeric shft0 = cmn.
SHIFT.
shft[iline] * ptot;
1644 absv[isig] += popudipo * wr * u_sqln2pi;
1646 absy[isig] += popudipo * u_sqln2pi * (wr - cmn.
YLT.
YT[iline] * wi);
1649 if (std::abs(cmn.
LineSg.
Sig(iline, iband)-sigc) > (rdmult*gamd)) {
1650 absv[isig] += popudipo * u_sqln2pi * u_pi * cmn.
GamT.
HWT[iline] / (pow2(cmn.
GamT.
HWT[iline]) + pow2(dsigc));
1652 absy[isig] += popudipo * u_sqln2pi * u_pi * (cmn.
GamT.
HWT[iline] + cmn.
YLT.
YT[iline] * dsigc) / (pow2(cmn.
GamT.
HWT[iline]) + pow2(dsigc));
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));
1658 absv[isig] += popudipo *
w.real() / gamd;
1660 absy[isig] += popudipo * (
w.real() - cmn.
YLT.
YT[iline] *
w.imag()) / gamd;
1672 absw[isig] += u_sqln2pi * u_pi * (cmn.
Zss.
ZS[iline] / (sigc - cmn.
Zaa.
ZA[iline])).imag();
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;
1694 const Vector &f_grid)
1700 const Index nf = f_grid.nelem();
1703 Vector absorption(nf, 0);
1709 constexpr Numeric u_sqln2pi = 1 / sq_ln2pi;
1711 for (Index iband=0; iband<bands.
nelem(); iband++) {
1712 if (not bands[iband].DoLineMixing(P))
continue;
1714 const Numeric rat_isot = isotopologue_ratio[bands[iband].Isotopologue()];
1716 auto tp =
convtp(vmrs, hitran, bands[iband], T, P);
1717 const Numeric GD_div_F0 = bands[iband].DopplerConstant(T);
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;
1724 for (Index iv=0; iv<nf; iv++) {
1725 const Numeric f = f_grid[iv];
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]);
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();
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;
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;
1748 a += u_sqln2pi * u_pi * (tp.eqv.str[iline] / (f - tp.eqv.val[iline])).imag();
1750 ARTS_USER_ERROR (
"Cannot understand the combination of calculations requested...");
1755 if (std::isnormal(
a) and
a > 0)
1756 absorption[iv] += rat_isot *
a;
1760 for (Index iv=0; iv<nf; iv++) {
1761 const Numeric f = f_grid[iv];
1763 absorption[iv] *= fact * dens * sq_ln2pi;
1770 const Numeric& sgminr,
1771 const Numeric& sgmaxr,
1772 const Numeric& stotmax,
1773 const String& basedir=
"data_new/")
1776 std::ifstream fortranfile;
1777 fortranfile.open(basedir +
"/BandInfo.dat");
1780 "Cannot read the file. Please make sure you point at BandInfo.dat basedir and have the right to read.");
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(),
1790 "%c%c" "%1" PRId64
"%c%c"
1791 "%c%c" "%1" PRId64
"%c%c"
1797 "%4" PRId64
"%4" PRId64
"%4" PRId64,
1799 &c11, &c12, &lfr, &c21, &c22,
1800 &c31, &c32, &lir, &c41, &c42,
1805 &x, &x, &x, &x, &x, &x, &x, &x,
1806 &jmxp, &jmxq, &jmxr);
1808 getline(fortranfile, line);
1812 if ((sgminr < sgmax) and (sgmaxr > sgmin)) {
1819 std::array<char, 15> name;
1820 snprintf(name.data(),
1822 "S%" PRId64
"%c%c%" PRId64
"%c%c%c%c%" PRId64
"%c%c%c%c",
1842 fortranfile.close();
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;
1852 const String cr = std::to_string(lli) + std::to_string(llf);
1854 std::ifstream fortranfile;
1856 fortranfile.open(fname.c_str());
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(),
1867 "%4" PRId64
"%4" PRId64
"%4" PRId64
"%4" PRId64,
1868 sw0r, sb0r, &dmaxdt, &wtmax, &jic, &jfc, &jipc, &jfpc);
1870 std::replace(ssw0r.begin(), ssw0r.end(),
'D',
'E');
1871 const Numeric w0r = std::stod(ssw0r);
1873 std::replace(ssb0r.begin(), ssb0r.end(),
'D',
'E');
1874 const Numeric b0r = std::stod(ssb0r);
1876 getline(fortranfile, line);
1882 "Out of bounds in reading...");
1884 if (jic > jfc and jipc > jfpc) {
1887 }
else if(jic > jfc and jipc == jfpc) {
1890 }
else if(jic > jfc and jipc < jfpc) {
1893 }
else if(jic < jfc and jipc > jfpc) {
1896 }
else if(jic < jfc and jipc == jfpc) {
1899 }
else if(jic < jfc and jipc < jfpc) {
1902 }
else if(jic == jfc and jipc > jfpc) {
1905 }
else if(jic == jfc and jipc == jfpc) {
1908 }
else if(jic == jfc and jipc < jfpc) {
1917Vector
compute(
const Numeric p,
const Numeric t,
const Numeric xco2,
const Numeric xh2o,
const Vector& invcm_grid,
const Numeric stotmax,
const calctype type)
1919 const Index n = invcm_grid.nelem();
1920 Vector absorption(n);
1926 detband(cmn, invcm_grid[0], invcm_grid[n-1], stotmax);
1934 compabs(cmn, t,p, xco2, xh2o, invcm_grid,
false,
true, absv, absy, absw);
1936 compabs(cmn, t,p, xco2, xh2o, invcm_grid,
true,
false, absv, absy, absw);
1938 compabs(cmn, t,p, xco2, xh2o, invcm_grid,
false,
false, absv, absy, absw);
1944 for (Index i=0; i<n; i++)
1945 absorption[i] = absv[i];
1950 for (Index i=0; i<n; i++)
1951 absorption[i] = absy[i];
1956 for (Index i=0; i<n; i++)
1957 absorption[i] = absw[i];
1970 const Vector& f_grid)
1972 return f_grid.nelem() ?
compabs(T, P, hitran, bands, isotopologue_ratio, vmrs, f_grid) : Vector(0);
1979 const Numeric linemixinglimit,
1983 const ModeOfLineMixing mode)
1985 String newbase = basedir;
1986 if (newbase.
nelem() == 0)
1991 readw(cmn, newbase);
1994 Numeric linemixinglimit_internal=0;
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.");
2005 const auto lstype =
typeLP(mode) ?
2006 LineShape::Type::LP :
2008 LineShape::Type::VP :
2009 LineShape::Type::SDVP) ;
2010 const auto poptype =
typeFull(mode) ?
2011 Absorption::PopulationType::ByHITRANFullRelmat :
2012 Absorption::PopulationType::ByHITRANRosenkranzRelmat ;
2069 Absorption::CutoffType::None,
2070 Absorption::MirroringType::None,
2072 Absorption::NormalizationType::None,
2076 linemixinglimit_internal,
2078 {Species::fromShortName(
"CO2"), Species::fromShortName(
"H2O"), Species::fromShortName(
"AIR")}};
2080 const Numeric rat_isot = isotopologue_ratio[bands[i].Isotopologue()];
2086 Rational(cmn.
Jiln.
Ji(j, i)),
2087 Rational(cmn.
Jfln.
Jf(j, i))));
2121 const auto lsmodel = typeVP(mode) ?
2125 Numeric qt0_co2, gsi0;
2132 gsi0*Numeric(cmn.
Jfln.
Jf(j, i) * 2 + 1),
2133 gsi0*Numeric(cmn.
Jiln.
Ji(j, i) * 2 + 1),
2143 const Vector& temperatures,
2148 const Index K = temperatures.nelem();
2152 const Numeric fmean = band.
F_mean();
2155 Tensor4 out(4, N, M, K, 0);
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++) {
2162 const Numeric T = temperatures[k];
2165 const Numeric ratiopart = QT0 / QT;
2167 for (Index i=0; i<N; i++) {
2170 calc.
f0[i] = band.
lines[i].F0;
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);
2175 calc.
hwt2[i] = P * band.
lines[i].lineshape[m].G2().at(T, band.
T0);
2179 calcw(calc, hitran, band, T,
true);
2182 if (Absorption::PopulationType::ByHITRANRosenkranzRelmat == band.
population) {
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]);
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;
2200 calc.
W, calc.
pop, calc.
dip, band, fmean, T, P, QT, QT0, m);
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();
2214 const Vector& temperatures,
2222 not is_sorted(temperatures),
2223 "The temperature list [", temperatures,
"] K\n"
2224 "must be fully sorted from low to high"
2229 temperatures, P0, ord)) {
2235 const Vector& temperatures,
2237 const Vector& pressures) {
2240 const Index K = temperatures.nelem();
2241 const Index L = pressures.size();
2243 Tensor5 out(4, N, M, K, L);
2244 for (Index l=0; l<L; l++) {
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)
This can be used to make arrays out of anything.
Index nelem() const ARTS_NOEXCEPT
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.
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR(...)
#define ARTS_USER_ERROR_IF(condition,...)
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.
Contains the line shape namespace.
my_basic_string< char > String
The String type for ARTS.
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.
This file contains declerations of functions of physical character.
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Contains recomputed equivalent lines (sorting is unknown)
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