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