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