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