ARTS  2.3.1285(git:92a29ea9-dirty)
linemixing_hitran.cc
Go to the documentation of this file.
1 /* Copyright (C) 2020
2  * Richard Larsson <larsson@mps.mpg.de>
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 <fstream>
28 #include <Faddeeva/Faddeeva.hh>
29 #include "linemixing_hitran.h"
30 
31 #include "lin_alg.h"
32 #include "linefunctions.h"
33 #include "physics_funcs.h"
34 
35 namespace lm_hitran_2017 {
36 namespace parameters {
37  static constexpr Index nBmx=7'000; // Max Number of Bands
38  static constexpr Index nLmx=700; // Max Number of Lines per Band
39  static constexpr Index Nlifmax=10; // Max number of l values
40  static constexpr Index Jmax=131; // Max number of j values
41 
42  static constexpr Numeric Ct=1.4387686e0; // Constant
43  static constexpr Numeric T0=296; // Constant
44  static constexpr Numeric CtGamD=1.1325e-08; // Constant
45  static constexpr Numeric aMolAtm=7.33889e+21; // Constant
46 
47  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
48 
49 };
50 
51 struct CommonBlock {
52 struct Bands {
53  Index nBand;
54  std::array<Index, parameters::nBmx> Isot;
55  std::array<Index, parameters::nBmx> nLines;
56  std::array<Index, parameters::nBmx> li;
57  std::array<Index, parameters::nBmx> lf;
58  std::array<String, parameters::nBmx> BandFile;
59 } Bands;
60 
61 struct LineSg {
62  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> Sig;
63  LineSg() : Sig(parameters::nLmx, parameters::nBmx) {};
64 } LineSg;
65 
66 struct DipoRigid {
67  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> Dipo0;
68  DipoRigid() : Dipo0(parameters::nLmx, parameters::nBmx) {};
69 } DipoRigid;
70 
71 struct Energy {
72  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> E;
73  Energy() : E(parameters::nLmx, parameters::nBmx) {};
74 } Energy;
75 
76 struct GamVT0AIR {
77  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWVT0AIR;
78  GamVT0AIR() : HWVT0AIR(parameters::nLmx, parameters::nBmx) {};
79 } GamVT0AIR;
80 
81 struct GamSDVT0AIR {
82  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWSDVT0AIR;
83  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> rHWT0AIR;
84  GamSDVT0AIR() : HWSDVT0AIR(parameters::nLmx, parameters::nBmx), rHWT0AIR(parameters::nLmx, parameters::nBmx) {};
85 } GamSDVT0AIR;
86 
87 struct DTGAMAIR {
88  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> BHWAIR;
89  DTGAMAIR() : BHWAIR(parameters::nLmx, parameters::nBmx) {};
90 } DTGAMAIR;
91 
92 struct GamVT0CO2 {
93  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWVT0SELF;
94  GamVT0CO2() : HWVT0SELF(parameters::nLmx, parameters::nBmx) {};
95 } GamVT0CO2;
96 
97 struct GamSDVT0CO2 {
98  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWSDVT0SELF;
99  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> rHWT0SELF;
100  GamSDVT0CO2() : HWSDVT0SELF(parameters::nLmx, parameters::nBmx), rHWT0SELF(parameters::nLmx, parameters::nBmx) {};
101 } GamSDVT0CO2;
102 
103 struct DTGAMCO2 {
104  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> BHWSELF;
105  DTGAMCO2() : BHWSELF(parameters::nLmx, parameters::nBmx) {};
106 } DTGAMCO2;
107 
108 struct GamVT0H2O {
109  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWVT0H2O;
110  GamVT0H2O() : HWVT0H2O(parameters::nLmx, parameters::nBmx) {};
111 } GamVT0H2O;
112 
113 struct GamSDVT0H2O {
114  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> HWSDVT0H2O;
115  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> rHWT0H2O;
116  GamSDVT0H2O() : HWSDVT0H2O(parameters::nLmx, parameters::nBmx), rHWT0H2O(parameters::nLmx, parameters::nBmx) {};
117 } GamSDVT0H2O;
118 
119 struct DTGAMH2O {
120  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> BHWH2O;
121  DTGAMH2O() : BHWH2O(parameters::nLmx, parameters::nBmx) {};
122 } DTGAMH2O;
123 
124 struct GamT {
125  std::array<Numeric, parameters::nLmx> HWT;
126  std::array<Numeric, parameters::nLmx> HWSDV2T;
127 } GamT;
128 
129 struct SHIFT {
130  std::array<Numeric, parameters::nLmx> shft;
131 } SHIFT;
132 
133 struct SHIFT0 {
134  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> shft0;
135  SHIFT0() : shft0(parameters::nLmx, parameters::nBmx) {};
136 } SHIFT0;
137 
138 struct PopuT {
139  std::array<Numeric, parameters::nLmx> PopuT;
140 } PopuT;
141 
142 struct PopTrf {
143  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> PopuT0;
144  PopTrf() : PopuT0(parameters::nLmx, parameters::nBmx) {};
145 } PopTrf;
146 
147 struct DipoTcm {
148  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> DipoT;
149  DipoTcm() : DipoT(parameters::nLmx, parameters::nBmx) {};
150 } DipoTcm;
151 
152 struct Jiln {
153  Eigen::Matrix<Index, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> Ji;
154  Jiln() : Ji(parameters::nLmx, parameters::nBmx) {};
155 } Jiln;
156 
157 struct Jfln {
158  Eigen::Matrix<Index, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> Jf;
159  Jfln() : Jf(parameters::nLmx, parameters::nBmx) {};
160 } Jfln;
161 
162 struct Zss {
163  std::array<Complex, parameters::nLmx> ZS;
164 } Zss;
165 
166 struct Zaa {
167  ComplexVector ZA;
168  Zaa() noexcept : ZA(parameters::nLmx) {}
169 } Zaa;
170 
171 struct Wmatrix {
172  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> W;
173  Wmatrix() : W(parameters::nLmx, parameters::nLmx) {};
174 } Wmatrix;
175 
176 struct Wfittedp {
177  Tensor4 W0pp;
178  Tensor4 W0pq;
179  Tensor4 W0pr;
180  Wfittedp() :
181  W0pp(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
182  W0pq(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
183  W0pr(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax) {}
184 } Wfittedp;
185 
186 struct Wfittedq {
187  Tensor4 W0qp;
188  Tensor4 W0qq;
189  Tensor4 W0qr;
190  Wfittedq() :
191  W0qp(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
192  W0qq(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
193  W0qr(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax) {}
194 } Wfittedq;
195 
196 struct Wfittedr {
197  Tensor4 W0rp;
198  Tensor4 W0rq;
199  Tensor4 W0rr;
200  Wfittedr() :
201  W0rp(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
202  W0rq(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
203  W0rr(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax) {}
204 } Wfittedr;
205 
206 struct Bfittedp {
207  Tensor4 B0pp;
208  Tensor4 B0pq;
209  Tensor4 B0pr;
210  Bfittedp() :
211  B0pp(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
212  B0pq(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
213  B0pr(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax) {}
214 } Bfittedp;
215 
216 struct Bfittedq {
217  Tensor4 B0qp;
218  Tensor4 B0qq;
219  Tensor4 B0qr;
220  Bfittedq() :
221  B0qp(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
222  B0qq(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
223  B0qr(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax) {}
224 } Bfittedq;
225 
226 struct Bfittedr {
227  Tensor4 B0rp;
228  Tensor4 B0rq;
229  Tensor4 B0rr;
230  Bfittedr() :
231  B0rp(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
232  B0rq(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax),
233  B0rr(parameters::Nlifmax, parameters::Nlifmax, parameters::Jmax, parameters::Jmax) {}
234 } Bfittedr;
235 
236 struct DiagnR {
237  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> OpR;
238  DiagnR() : OpR(parameters::nLmx, parameters::nLmx) {};
239 } DiagnR;
240 
241 struct DiagnI {
242  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> OpI;
243  DiagnI() : OpI(parameters::nLmx, parameters::nLmx) {};
244 } DiagnI;
245 
246 struct YLT {
247  std::array<Numeric, parameters::nLmx> YT;
248 } YLT;
249 
250 struct UnusedBandParams {
251  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> intens;
252  Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> eina;
253  std::array<Rational, parameters::nBmx> iv1;
254  std::array<Rational, parameters::nBmx> iv2;
255  std::array<Rational, parameters::nBmx> il2;
256  std::array<Rational, parameters::nBmx> iv3;
257  std::array<Rational, parameters::nBmx> ir;
258  std::array<Rational, parameters::nBmx> fv1;
259  std::array<Rational, parameters::nBmx> fv2;
260  std::array<Rational, parameters::nBmx> fl2;
261  std::array<Rational, parameters::nBmx> fv3;
262  std::array<Rational, parameters::nBmx> fr;
263  UnusedBandParams() :
264  intens(parameters::nLmx, parameters::nBmx),
265  eina(parameters::nLmx, parameters::nBmx) {};
266 } UnusedBandParams;
267 }; // CommonBlock
268 
269 Rational toRationalSum(char a, char b=' ')
270 {
271  if (b == ' ' and a == ' ')
272  return Rational();
273  else if (b == ' ')
274  return Rational(a-'0');
275  else if (a == ' ')
276  return Rational(b-'0');
277  else
278  return Rational(10*(a-'0') + b-'0');
279 }
280 
281 void readlines(CommonBlock& cmn, const String& basedir="data_new/")
282 {
283  if (cmn.Bands.nBand > parameters::nBmx)
284  throw std::runtime_error("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  if (nliner >= parameters::nLmx)
298  throw std::runtime_error("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" "%1ld" "%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" "%3ld" "%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  } else if(not (
360  cmn.UnusedBandParams.iv1[iband]==Rational(toRationalSum(iv11, iv12)) and
361  cmn.UnusedBandParams.iv2[iband]==Rational(toRationalSum(iv21, iv22)) and
362  cmn.UnusedBandParams.il2[iband]==Rational(toRationalSum(il21, il22)) and
363  cmn.UnusedBandParams.iv3[iband]==Rational(toRationalSum(iv31, iv32)) and
364  cmn.UnusedBandParams.ir[iband]==Rational(toRationalSum(ir1)) and
365  cmn.UnusedBandParams.fv1[iband]==Rational(toRationalSum(fv11, fv12)) and
366  cmn.UnusedBandParams.fv2[iband]==Rational(toRationalSum(fv21, fv22)) and
367  cmn.UnusedBandParams.fl2[iband]==Rational(toRationalSum(fl21, fl22)) and
368  cmn.UnusedBandParams.fv3[iband]==Rational(toRationalSum(fv31, fv32)) and
369  cmn.UnusedBandParams.fr[iband]==Rational(toRationalSum(fr1)))) {
370  throw std::runtime_error("Bad read, bands do not have the same global quantum numbers...");
371  }
372 
373  // Fix...
374  String ssDipoRigid = sDipoRigid;
375  String ssPopTrf = sPopTrf;
376  std::replace(ssDipoRigid.begin(), ssDipoRigid.end(), 'D', 'E');
377  std::replace(ssPopTrf.begin(), ssPopTrf.end(), 'D', 'E');
378  cmn.DipoRigid.Dipo0(nliner, iband) = std::stod(ssDipoRigid);
379  cmn.PopTrf.PopuT0(nliner, iband) = std::stod(ssPopTrf);
380 
381  // Dipole at temperature
382  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))));
383 
384  // Fix Js
385  if (tpline == 'P')
386  cmn.Jfln.Jf(nliner, iband) = cmn.Jiln.Ji(nliner, iband) - 1;
387  else if (tpline == 'Q')
388  cmn.Jfln.Jf(nliner, iband) = cmn.Jiln.Ji(nliner, iband);
389  else
390  cmn.Jfln.Jf(nliner, iband) = cmn.Jiln.Ji(nliner, iband) + 1;
391 
392  // Fix isotologue
393  if (cmn.Bands.Isot[iband] == 0)
394  cmn.Bands.Isot[iband] = 10;
395 
396  nliner++;
397  }
398  cmn.Bands.nLines[iband] = nliner;
399  }
400 
401  fortranfile.close();
402  }
403 }
404 
405 template<size_t NT>
406 Numeric atob(const Numeric& aa,
407  const std::array<Numeric, NT>& a,
408  const std::array<Numeric, NT>& b)
409 {
410  for (size_t i=1; i<NT; i++) {
411  if (a[i] >= aa) {
412  if (i < 2 or i == NT) {
413  size_t j=i;
414 
415  if (i < 2)
416  j = 2;
417  if (i == NT-1)
418  j = NT-1;
419 
420  Numeric A0D1=a[j-2]-a[j-1]; if (A0D1 == 0) A0D1=0.0001;
421  Numeric A0D2=a[j-2]-a[j]; if (A0D2 == 0) A0D2=0.0001;
422  Numeric A1D1=a[j-1]-a[j-2]; if (A1D1 == 0) A1D1=0.0001;
423  Numeric A1D2=a[j-1]-a[j]; if (A1D2 == 0) A1D2=0.0001;
424  Numeric A2D1=a[j]-a[j-2]; if (A2D1 == 0) A2D1=0.0001;
425  Numeric A2D2=a[j]-a[j-1]; if (A2D2 == 0) A2D2=0.0001;
426 
427  const Numeric A0=(aa-a[j-1])*(aa-a[j])/(A0D1*A0D2);
428  const Numeric A1=(aa-a[j-2])*(aa-a[j])/(A1D1*A1D2);
429  const Numeric A2=(aa-a[j-2])*(aa-a[j-1])/(A2D1*A2D2);
430 
431  return A0*b[j-2] + A1*b[j-1] + A2*b[j];
432  } else {
433  size_t j = i;
434 
435  Numeric A0D1=a[j-2]-a[j-1]; if (A0D1 == 0) A0D1=0.0001;
436  Numeric A0D2=a[j-2]-a[j]; if (A0D2 == 0) A0D2=0.0001;
437  Numeric A0D3=a[j-2]-a[j+1]; if (A0D3 == 0) A0D3=0.0001;
438  Numeric A1D1=a[j-1]-a[j-2]; if (A1D1 == 0) A1D1=0.0001;
439  Numeric A1D2=a[j-1]-a[j]; if (A1D2 == 0) A1D2=0.0001;
440  Numeric A1D3=a[j-1]-a[j+1]; if (A1D3 == 0) A1D3=0.0001;
441  Numeric A2D1=a[j]-a[j-2]; if (A2D1 == 0) A2D1=0.0001;
442  Numeric A2D2=a[j]-a[j-1]; if (A2D2 == 0) A2D2=0.0001;
443  Numeric A2D3=a[j]-a[j+1]; if (A2D3 == 0) A2D3=0.0001;
444  Numeric A3D1=a[j+1]-a[j-2]; if (A3D1 == 0) A3D1=0.0001;
445  Numeric A3D2=a[j+1]-a[j-1]; if (A3D2 == 0) A3D2=0.0001;
446  Numeric A3D3=a[j+1]-a[j]; if (A3D3 == 0) A3D3=0.0001;
447 
448 
449  Numeric A0=(aa-a[j-1])*(aa-a[j])*(aa-a[j+1]);
450  A0=A0/(A0D1*A0D2*A0D3);
451  Numeric A1=(aa-a[j-2])*(aa-a[j])*(aa-a[j+1]);
452  A1=A1/(A1D1*A1D2*A1D3);
453  Numeric A2=(aa-a[j-2])*(aa-a[j-1])*(aa-a[j+1]);
454  A2=A2/(A2D1*A2D2*A2D3);
455  Numeric A3=(aa-a[j-2])*(aa-a[j-1])*(aa-a[j]);
456  A3=A3/(A3D1*A3D2*A3D3);
457 
458  return A0*b[j-2] + A1*b[j-1] + A2*b[j] + A3*b[j+1];
459  }
460  }
461  }
462 
463  return std::numeric_limits<Numeric>::quiet_NaN();
464 }
465 
466 void 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... -- 728
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 
782 void calcw(CommonBlock& cmn,
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 
950 // Class not allowed to copy
951 struct EqvLinesOut {
952  ComplexVector zval, zstr;
953  explicit EqvLinesOut(Index n=0) noexcept : zval(n, 0), zstr(n, 0) {}
954  EqvLinesOut(const EqvLinesOut&) = delete;
955  EqvLinesOut(EqvLinesOut&&) = default;
956  EqvLinesOut& operator=(const EqvLinesOut&) = delete;
957  EqvLinesOut& operator=(EqvLinesOut&&) = default;
958 };
959 
960 
961 // Class not allowed to copy
962 struct ConvTPOut {
963  Vector Y, hwt, hwt2, shft, f0, pop, dip;
964  ComplexMatrix W;
965  EqvLinesOut eqv;
966  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) {}
967  ConvTPOut(const ConvTPOut&) = delete;
968  ConvTPOut(ConvTPOut&&) = default;
969  ConvTPOut& operator=(const ConvTPOut&) = delete;
970  ConvTPOut& operator=(ConvTPOut&&) = default;
971 };
972 
973 
974 void calcw(ConvTPOut& out,
975  const HitranRelaxationMatrixData& hitran,
976  const AbsorptionLines& band,
977  const Numeric T)
978 {
979  Vector& Y = out.Y;
980  Vector& g0 = out.hwt;
981  Vector& g2 = out.hwt2;
982  Vector& d0 = out.shft;
983  Vector& f0 = out.f0;
984  Vector& pop = out.pop;
985  Vector& dip = out.dip;
986  MatrixView W = transpose(out.W.imag()); // Transpose to fit Fortran-code
987 
988  // Size of problem
989  const Index n=Y.nelem();
990 
991  // Copies for local variables
992  const Rational li = band.UpperQuantumNumber(0, QuantumNumberType::l1);
993  const Rational lf = band.LowerQuantumNumber(0, QuantumNumberType::l1);
994  Vector dip0(n);
995  std::vector<Rational> Ji(n), Jf(n);
996  for (Index i=0; i<n; i++) {
997  Ji[i] = band.UpperQuantumNumber(i, QuantumNumberType::J);
998  Jf[i] = band.LowerQuantumNumber(i, QuantumNumberType::J);
999  dip0[i] = Absorption::reduced_rovibrational_dipole(Jf[i], Ji[i], lf, li, 1_rat);
1000  }
1001 
1002  Vector s(n);
1003  if (li > 8 or abs(li - lf) > 1) {
1004  for (Index i=0; i<n; i++) {
1005  W(i, i) = g0[i]; // nb... units Hz/Pa
1006  }
1007  } else {
1008  for (Index i=0; i<n; i++) {
1009  s[i] = f0[i] * pop[i] * Constant::pow2(dip[i]);
1010  }
1011 
1012  for (Index i=0; i<n; i++) {
1013  for (Index j=i+1; j<n; j++) {
1014  if (s[j] > s[i]) {
1015  // Switch switches 2 elements in a vector :: std::swap(A(i), A(j));
1016  // SwitchTwo switches two elements in a Matrix :: std::swap(A(i, k), A(j, k));
1017  // SwitchInt is as SwitchTwo but for Integers :: std::swap(A(i, k), A(j, k));
1018  std::swap(Ji[i], Ji[j]);
1019  std::swap(Jf[i], Jf[j]);
1020  std::swap(g0[i], g0[j]);
1021  std::swap(d0[i], d0[j]);
1022  std::swap(g2[i], g2[j]);
1023  std::swap(f0[i], f0[j]);
1024  std::swap(dip[i], dip[j]);
1025  std::swap(pop[i], pop[j]);
1026  std::swap(dip0[i], dip0[j]);
1027  std::swap(s[i], s[j]);
1028  }
1029  }
1030  }
1031 
1032  const Numeric dlgt0t = std::log(band.T0() / T);
1033  const Rational lli = std::min(li, lf);
1034  const Rational llf = std::max(li, lf);
1035 
1036  // Set off-diagonal elements
1037  for (Index i=0; i<n; i++) {
1038  Rational jji, jjf;
1039  if (li <= lf) {
1040  jji = Ji[i];
1041  jjf = Jf[i];
1042  } else {
1043  jji = Jf[i];
1044  jjf = Ji[i];
1045  }
1046 
1047  for (Index j=0; j<n; j++) {
1048  Rational jjip, jjfp;
1049  if (li <= lf) {
1050  jjip = Ji[j];
1051  jjfp = Jf[j];
1052  } else {
1053  jjip = Jf[j];
1054  jjfp = Ji[j];
1055  }
1056 
1057  if (jjip > jji) continue;
1058 
1059  // nb. FIX THIS FOR ARTS SPECIFIC Isotopologue
1060  const bool skip = band.Isotopologue() > 2 and
1061  band.Isotopologue() not_eq 7 and
1062  band.Isotopologue() not_eq 10 and
1063  (abs(Ji[i] - Ji[j]) % 2) not_eq 0;
1064 
1065  if(skip) continue;
1066 
1067  Numeric w0=0, b0=0;
1068  if (jji > jjf and jjip > jjfp) {
1069  w0 = hitran.W0pp(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1070  b0 = hitran.B0pp(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1071  } else if (jji > jjf and jjip == jjfp) {
1072  w0 = hitran.W0pq(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1073  b0 = hitran.B0pq(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1074  } else if (jji > jjf and jjip < jjfp) {
1075  w0 = hitran.W0pr(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1076  b0 = hitran.B0pr(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1077  } else if (jji < jjf and jjip > jjfp) {
1078  w0 = hitran.W0rp(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1079  b0 = hitran.B0rp(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1080  } else if (jji < jjf and jjip == jjfp) {
1081  w0 = hitran.W0rq(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1082  b0 = hitran.B0rq(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1083  } else if (jji < jjf and jjip < jjfp) {
1084  w0 = hitran.W0rr(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1085  b0 = hitran.B0rr(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1086  } else if (jji == jjf and jjip > jjfp) {
1087  w0 = hitran.W0qp(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1088  b0 = hitran.B0qp(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1089  } else if (jji == jjf and jjip == jjfp) {
1090  w0 = hitran.W0qq(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1091  b0 = hitran.B0qq(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1092  } else if (jji == jjf and jjip < jjfp) {
1093  w0 = hitran.W0qr(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1094  b0 = hitran.B0qr(lli.toIndex(), llf.toIndex(), jji.toIndex(), jjip.toIndex());
1095  }
1096 
1097  // nb... should order be different???
1098  const Numeric ycal = std::exp(w0 - b0*dlgt0t);
1099  W(j, i) = ycal;
1100  W(i, j) = ycal * pop[i] / pop[j];
1101  }
1102  }
1103 
1104  // Weird undocumented minus sign
1105  for (Index i=0; i<n; i++) {
1106  for (Index j=0; j<n; j++) {
1107  if (j not_eq i) {
1108  W(i, j) = - std::abs(W(i, j));
1109  }
1110  }
1111  }
1112 
1113  // Set diagonal to measured broadening
1114  for (Index i=0; i<n; i++) {
1115  W(i, i) = g0[i]; // nb... units Hz/Pa
1116  }
1117 
1118  // Sum rule correction
1119  for (Index i=0; i<n; i++) {
1120  Numeric sumlw = 0.0;
1121  Numeric sumup = 0.0;
1122 
1123  for (Index j=0; j<n; j++) {
1124 
1125  // nb. FIX THIS FOR ARTS SPECIFIC Isotopologue
1126  const bool skip = band.Isotopologue() > 2 and
1127  band.Isotopologue() not_eq 7 and
1128  band.Isotopologue() not_eq 10 and
1129  (abs(Ji[i] - Ji[j]) % 2) not_eq 0;
1130 
1131  if (skip) continue;
1132 
1133  if (j > i) {
1134  sumlw += std::abs(dip0[j]) * W(j, i);
1135  } else {
1136  sumup += std::abs(dip0[j]) * W(j, i);
1137  }
1138  }
1139 
1140  for (Index j=i+1; j<n; j++) {
1141  if (sumlw == 0) {
1142  W(j, i) = 0.0;
1143  W(i, j) = 0.0;
1144  } else {
1145  W(j, i) *= - sumup / sumlw;
1146  W(i, j) = W(j, i) * pop[i] / pop[j];
1147  }
1148  }
1149  }
1150 
1151  // Rosenkranz coefficients
1152  for (Index i=0; i<n; i++) {
1153  Numeric sum0 = 0;
1154  for (Index j=0; j<n; j++) {
1155  if (i == j) continue;
1156 
1157  // nb. FIX THIS FOR ARTS SPECIFIC Isotopologue
1158  const bool skip = band.Isotopologue() > 2 and
1159  band.Isotopologue() not_eq 7 and
1160  band.Isotopologue() not_eq 10 and
1161  (abs(Ji[i] - Ji[j]) % 2) not_eq 0;
1162 
1163  if (skip) continue;
1164 
1165  Numeric deltasig = f0[i] - f0[j];
1166  if (std::abs(deltasig) < Conversion::kaycm2freq(1e-4 /*cm-1*/))
1167  deltasig = Conversion::kaycm2freq(1e-4) /*Hz*/;
1168 
1169  sum0 += 2 * std::abs(dip[j]) / std::abs(dip[i]) * W(j, i) / deltasig;
1170  }
1171 
1172  Y[i] = sum0;
1173  }
1174 
1175  }
1176 }
1177 
1178 void eqvlines(CommonBlock& cmn,
1179  const Index& iband,
1180  const Index& n,
1181  const Numeric& sigmoy)
1182 {
1183  ComplexMatrix zop(n, n), zvec(n, n);
1184  ComplexMatrixView inv_zvec = zop; // Rename to not confuse later operations
1185  ComplexVectorView zval = cmn.Zaa.ZA[Range(0, n)]; // Rename and rescale to right size
1186  for (Index i=0; i<n; i++) {
1187  for (Index j=0; j<n; j++) {
1188  zop(j, i) = Complex(cmn.DiagnR.OpR(i, j), cmn.DiagnI.OpI(i, j)); // nb. reverse from Fortran algorithm due to row-col issues
1189  }
1190  }
1191 
1192  // Main computations
1193  diagonalize(zvec, zval, zop);
1194  inv(inv_zvec, zvec);
1195 
1196  // Add average sigma
1197  zval += sigmoy;
1198 
1199  // Do the matrix multiplication
1200  for (Index i=0; i<n; i++) {
1201  cmn.Zss.ZS[i] = 0;
1202  Complex z(0, 0);
1203  for (Index j=0; j<n; j++) {
1204  cmn.Zss.ZS[i] += cmn.DipoTcm.DipoT(j, iband) * zvec(j, i);
1205  z += cmn.PopuT.PopuT[j] * cmn.DipoTcm.DipoT(j, iband) * inv_zvec(i, j);
1206  }
1207  cmn.Zss.ZS[i] *= z;
1208  }
1209 }
1210 
1211 EqvLinesOut eqvlines(const ConstComplexMatrixView W,
1212  const ConstVectorView pop,
1213  const ConstVectorView dip,
1214  const Numeric& fmean)
1215 {
1216  // Size of problem
1217  const Index n = pop.nelem();
1218 
1219  // Compute values
1220  ComplexMatrix zvec(n, n);
1221  EqvLinesOut out(n);
1222 
1223  // Main computations
1224  diagonalize(zvec, out.zval, W);
1225 
1226  // Do the matrix forward multiplication
1227  for (Index i=0; i<n; i++) {
1228  for (Index j=0; j<n; j++) {
1229  out.zstr[i] += dip[j] * zvec(j, i);
1230  }
1231  }
1232 
1233  // Do the matrix backward multiplication
1234  auto& inv_zvec=zvec.inv();
1235  for (Index i=0; i<n; i++) {
1236  Complex z(0, 0);
1237  for (Index j=0; j<n; j++) {
1238  z += pop[j] * dip[j] * inv_zvec(i, j);
1239  }
1240  out.zstr[i] *= z;
1241  }
1242 
1243  // Add the weighted frequency
1244  out.zval += fmean;
1245 
1246  return out;
1247 }
1248 
1249 
1250 void convtp(CommonBlock& cmn,
1251  const Index& iband,
1252  const Index& isotc,
1253  const Index& nlinec,
1254  const Numeric& xh2o,
1255  const Numeric& xco2,
1256  const Numeric& temp,
1257  const Numeric& ptot,
1258  Numeric& sigmoy,
1259  Numeric& gamd_gam0_mx,
1260  const bool mixfull,
1261  const bool mixsdv)
1262 {
1263  const Numeric ratiot = parameters::T0 / temp;
1264  Numeric qt0_co2, qtt_co2, gsi0, gsit;
1265  qt_co2(parameters::T0, isotc, gsi0, qt0_co2);
1266  qt_co2(temp, isotc, gsit, qtt_co2);
1267  const Numeric ratiopart = qt0_co2 / qtt_co2;
1268 
1269  const Numeric sqrtm = std::sqrt(temp / parameters::aMass[isotc-1]);
1270 
1271  Numeric sumwgt=0;
1272  sigmoy=0;
1273  for (Index iline=0; iline<nlinec; iline++) {
1274  cmn.PopuT.PopuT[iline]=cmn.PopTrf.PopuT0(iline, iband)*ratiopart*std::exp(-parameters::Ct*cmn.Energy.E(iline, iband)*(1/temp-1/parameters::T0));
1275  if (mixsdv) {
1276  cmn.GamT.HWT[iline] = ((1-xh2o-xco2) * (cmn.GamSDVT0AIR.HWSDVT0AIR(iline, iband) * std::pow(ratiot, cmn.DTGAMAIR.BHWAIR(iline, iband))))
1277  + (xh2o * (cmn.GamSDVT0H2O.HWSDVT0H2O(iline, iband) * std::pow(ratiot, cmn.DTGAMH2O.BHWH2O(iline, iband))))
1278  + (xco2 * (cmn.GamSDVT0CO2.HWSDVT0SELF(iline, iband) * std::pow(ratiot, cmn.DTGAMCO2.BHWSELF(iline, iband))));
1279 
1280  cmn.GamT.HWSDV2T[iline] =
1281  ((1-xh2o-xco2) * (cmn.GamSDVT0AIR.HWSDVT0AIR(iline, iband) * cmn.GamSDVT0AIR.rHWT0AIR(iline, iband) * std::pow(ratiot, cmn.DTGAMAIR.BHWAIR(iline, iband))))
1282  + (xh2o * (cmn.GamSDVT0H2O.HWSDVT0H2O(iline, iband)*cmn.GamSDVT0H2O.rHWT0H2O(iline, iband) * std::pow(ratiot, cmn.DTGAMH2O.BHWH2O(iline, iband))))
1283  + (xco2 * (cmn.GamSDVT0CO2.HWSDVT0SELF(iline, iband) * cmn.GamSDVT0CO2.rHWT0SELF(iline, iband) * std::pow(ratiot, cmn.DTGAMCO2.BHWSELF(iline, iband))));
1284  } else {
1285  cmn.GamT.HWT[iline] = ((1-xh2o-xco2) * (cmn.GamVT0AIR.HWVT0AIR(iline,iband) * std::pow(ratiot, cmn.DTGAMAIR.BHWAIR(iline, iband))))
1286  + (xh2o * (cmn.GamVT0H2O.HWVT0H2O(iline, iband) * std::pow(ratiot, cmn.DTGAMH2O.BHWH2O(iline, iband))))
1287  + (xco2 * (cmn.GamVT0CO2.HWVT0SELF(iline, iband) * std::pow(ratiot, cmn.DTGAMCO2.BHWSELF(iline, iband))));
1288  }
1289  cmn.SHIFT.shft[iline]=cmn.SHIFT0.shft0(iline, iband);
1290  const Numeric wgt = cmn.PopuT.PopuT[iline] * Constant::pow2(cmn.DipoTcm.DipoT(iline, iband));
1291  sumwgt += wgt;
1292  sigmoy += cmn.LineSg.Sig(iline,iband) * wgt;
1293  }
1294  sigmoy /= sumwgt;
1295 
1296  calcw(cmn, nlinec, iband, temp);
1297 
1298  // Adjust for pressure
1299  for (Index iline=0; iline<nlinec; iline++) {
1300  cmn.GamT.HWT[iline] *= ptot;
1301  cmn.GamT.HWSDV2T[iline] *= ptot;
1302  cmn.YLT.YT[iline] *= ptot;
1303  const Numeric gamd = parameters::CtGamD * cmn.LineSg.Sig(iline, iband) * sqrtm;
1304  gamd_gam0_mx = std::max(gamd_gam0_mx, gamd/cmn.GamT.HWT[iline]);
1305  }
1306 
1307  if (mixfull) {
1308  for (Index iline=0; iline<nlinec; iline++) {
1309  for (Index ilinep=0; ilinep<nlinec; ilinep++) {
1310  cmn.DiagnI.OpI(iline, ilinep) = ptot * cmn.Wmatrix.W(iline, ilinep);
1311  if (iline == ilinep) {
1312  cmn.DiagnR.OpR(iline, iline) = cmn.LineSg.Sig(iline, iband) - sigmoy + cmn.SHIFT.shft[iline] * ptot;
1313  } else {
1314  cmn.DiagnR.OpR(iline, ilinep) = 0;
1315  }
1316  }
1317  }
1318 
1319  eqvlines(cmn, iband, nlinec, sigmoy);
1320  }
1321 }
1322 
1323 ConvTPOut convtp(const ConstVectorView vmrs,
1324  const HitranRelaxationMatrixData& hitran,
1325  const AbsorptionLines& band,
1326  const Numeric T,
1327  const Numeric P,
1328  const SpeciesAuxData::AuxType& partition_type,
1329  const ArrayOfGriddedField1& partition_data)
1330 {
1331  const Index n = band.NumLines();
1332 
1333  const Numeric QT = single_partition_function(T, partition_type, partition_data);
1334  const Numeric QT0 = single_partition_function(band.T0(), partition_type, partition_data);
1335  const Numeric ratiopart = QT0 / QT;
1336 
1337  ConvTPOut out(n);
1338 
1339  Vector wgt(n);
1340  for (Index i=0; i<n; i++) {
1341  const Numeric pop0 = (band.g_upp(i) / QT0) * boltzman_factor(band.T0(), band.E0(i));
1342 
1343  out.f0[i] = band.F0(i);
1344  out.pop[i] = pop0 * ratiopart * boltzman_ratio(T, band.T0(), band.E0(i));
1345  out.hwt[i] = band.Line(i).LineShape().G0(T, band.T0(), 1, vmrs);
1346  out.shft[i] = band.Line(i).LineShape().D0(T, band.T0(), 1, vmrs);
1347  out.dip[i] = std::sqrt(band.I0(i)/(pop0 * band.F0(i) * (1-stimulated_emission(band.T0(), band.F0(i)))));
1348  out.hwt2[i] = band.Line(i).LineShape().G2(T, band.T0(), 1, vmrs);
1349  wgt[i] = out.pop[i] * Constant::pow2(out.dip[i]);
1350  }
1351 
1352  // Calculate the relaxation matrix
1353  calcw(out, hitran, band, T);
1354 
1355  // Adjust for pressure
1356  out.hwt *= P;
1357  out.hwt2 *= P;
1358  out.shft *= P;
1359  out.Y *= P;
1360 
1361  if (band.Population() == Absorption::PopulationType::ByHITRANFullRelmat) {
1362  const Numeric fmean = band.F_mean(wgt);
1363  out.W.diagonal().real() = out.f0;
1364  out.W.diagonal().real() += out.shft;
1365  out.W.diagonal().real() -= fmean;
1366  out.W.imag() *= P;
1367  out.eqv = eqvlines(out.W, out.pop, out.dip, fmean);
1368  }
1369 
1370  return out;
1371 }
1372 
1373 
1374 void qsdv(const Numeric& sg0,
1375  const Numeric& gamd,
1376  const Numeric& gam0,
1377  const Numeric& gam2,
1378  const Numeric& shift0,
1379  const Numeric& shift2,
1380  const Numeric& sg,
1381  Numeric& ls_qsdv_r,
1382  Numeric& ls_qsdv_i)
1383 {
1384  using Constant::pow2;
1385  Complex x , y, csqrty, z1, z2, w1, w2, aterm;
1386  Numeric xz1, yz1, xz2, yz2;
1387 
1388  const Numeric cte = Constant::sqrt_ln_2 / gamd;
1389  const Numeric pi = Constant::pi;
1390  const Numeric rpi = Constant::sqrt_pi;
1391  const Complex iz = Complex(0, 1);
1392 
1393  const Complex c0 = Complex(gam0, shift0);
1394  const Complex c2 = Complex(gam2, shift2);
1395  const Complex c0t = (c0 - 1.5 * c2);
1396  const Complex c2t = c2;
1397 
1398  if (std::abs(c2t) == 0) goto label110;
1399  x = (iz * (sg0 - sg) + c0t) / c2t;
1400  y = 1 / pow2(2 * cte * c2t);
1401  csqrty = (gam2 - iz * shift2) / (2 * cte * (pow2(gam2) + pow2(shift2)));
1402  if (std::abs(x) <= 3.e-8*std::abs(y)) goto label120; // IN cgs units
1403  if (std::abs(y) <= 1.e-15*std::abs(x)) goto label140; // IN cgs units
1404  z1 = std::sqrt(x + y) - csqrty;
1405  z2 = z1 + 2 * csqrty;
1406 
1407  xz1 = - z1.imag();
1408  yz1 = z1.real();
1409  xz2 = - z2.imag();
1410  yz2 = z2.real();
1411 
1412  w1 = Faddeeva::w(Complex(xz1 ,yz1));
1413  w2 = Faddeeva::w(Complex(xz2 ,yz2));
1414  aterm = rpi * cte * (w1 - w2);
1415  goto label10;
1416 
1417  label110:
1418  z1 = (iz*(sg0-sg)+c0t)*cte;
1419  xz1 = - z1.imag();
1420  yz1 = z1.real();
1421 
1422  w1 = Faddeeva::w(Complex(xz1 ,yz1));
1423  aterm = rpi*cte* w1;
1424  goto label10;
1425 
1426  label120:
1427 
1428  z1 = (iz*(sg0-sg)+c0t)*cte;
1429  z2 = std::sqrt(x+y)+csqrty;
1430  xz1 = - z1.imag();
1431  yz1 = z1.real();
1432  xz2 = - z2.imag();
1433  yz2 = z2.real();
1434 
1435  w1 = Faddeeva::w(Complex(xz1 ,yz1));
1436  w2 = Faddeeva::w(Complex(xz2 ,yz2));
1437  aterm = rpi * cte * (w1-w2);
1438  goto label10;
1439 
1440  label140:
1441 
1442  if (std::abs(std::sqrt(x)) <= 4.e3) { // IN cgs
1443  const Numeric xxb = - std::sqrt(x).imag();
1444  const Numeric yxb = std::sqrt(x).real();
1445  const Complex wb = Faddeeva::w(Complex(xxb, yxb));
1446  aterm = (2 * rpi / c2t) * (1 / rpi - std::sqrt(x) * wb);
1447  } else {
1448  aterm = (1 / c2t) * (1 / x - 1.5 / pow2(x));
1449  }
1450 
1451  label10:
1452 
1453  const Complex ls_qsdv = (1 / pi) * aterm;
1454 
1455  ls_qsdv_r = ls_qsdv.real();
1456  ls_qsdv_i = ls_qsdv.imag();
1457 }
1458 
1459 
1460 // nb.qsdv times sg0 is qsdv_si times F0
1461 Complex qsdv_si(const Numeric F0,
1462  const Numeric gamd,
1463  const Numeric gam0,
1464  const Numeric gam2,
1465  const Numeric shift0,
1466  const Numeric shift2,
1467  const Numeric f)
1468 {
1469  using Constant::pow2;
1470  Complex x , y, csqrty, z1, z2, w1, w2, aterm;
1471  Numeric xz1, yz1, xz2, yz2;
1472 
1473  const Numeric cte = Constant::sqrt_ln_2 / gamd;
1474  const Numeric pi = Constant::pi;
1475  const Numeric rpi = Constant::sqrt_pi;
1476  const Complex iz = Complex(0, 1);
1477 
1478  const Complex c0 = Complex(gam0, shift0);
1479  const Complex c2 = Complex(gam2, shift2);
1480  const Complex c0t = (c0 - 1.5 * c2);
1481  const Complex c2t = c2;
1482 
1483  if (std::abs(c2t) == 0) goto label110;
1484  x = (iz * (F0 - f) + c0t) / c2t;
1485  y = 1 / pow2(2 * cte * c2t);
1486  csqrty = (gam2 - iz * shift2) / (2 * cte * (pow2(gam2) + pow2(shift2)));
1487  if (std::abs(x) <= 3.e-8*std::abs(y)) goto label120;
1488  if (std::abs(y) <= 1.e-15*std::abs(x)) goto label140;
1489  z1 = std::sqrt(x + y) - csqrty;
1490  z2 = z1 + 2 * csqrty;
1491 
1492  xz1 = - z1.imag();
1493  yz1 = z1.real();
1494  xz2 = - z2.imag();
1495  yz2 = z2.real();
1496 
1497  w1 = Faddeeva::w(Complex(xz1 ,yz1));
1498  w2 = Faddeeva::w(Complex(xz2 ,yz2));
1499  aterm = rpi * cte * (w1 - w2);
1500  goto label10;
1501 
1502  label110:
1503  z1 = (iz*(F0 - f) + c0t)*cte;
1504  xz1 = - z1.imag();
1505  yz1 = z1.real();
1506 
1507  w1 = Faddeeva::w(Complex(xz1 ,yz1));
1508  aterm = rpi*cte* w1;
1509  goto label10;
1510 
1511  label120:
1512 
1513  z1 = (iz*(F0 - f) + c0t)*cte;
1514  z2 = std::sqrt(x+y)+csqrty;
1515  xz1 = - z1.imag();
1516  yz1 = z1.real();
1517  xz2 = - z2.imag();
1518  yz2 = z2.real();
1519 
1520  w1 = Faddeeva::w(Complex(xz1 ,yz1));
1521  w2 = Faddeeva::w(Complex(xz2 ,yz2));
1522  aterm = rpi * cte * (w1-w2);
1523  goto label10;
1524 
1525  label140:
1526 
1527  if (std::abs(std::sqrt(x)) <= 4.e3) { // IN cgs
1528  const Numeric xxb = - std::sqrt(x).imag();
1529  const Numeric yxb = std::sqrt(x).real();
1530  const Complex wb = Faddeeva::w(Complex(xxb, yxb));
1531  aterm = (2 * rpi / c2t) * (1 / rpi - std::sqrt(x) * wb);
1532  } else {
1533  aterm = (1 / c2t) * (1 / x - 1.5 / pow2(x));
1534  }
1535 
1536  label10:
1537 
1538  return (1 / pi) * aterm;
1539 }
1540 
1541 void compabs(
1542  CommonBlock& cmn,
1543  const Numeric& temp,
1544  const Numeric& ptot,
1545  const Numeric& xco2,
1546  const Numeric& xh2o,
1547  const ConstVectorView invcm_grid,
1548  const bool mixsdv,
1549  const bool mixfull,
1550  VectorView absv,
1551  VectorView absy,
1552  VectorView absw)
1553 {
1554  using Constant::pow2;
1555  using Constant::pow3;
1556 
1557  constexpr Numeric rdmult = 30;
1558 
1559  // Number of points to compute
1560  const Index nsig = invcm_grid.nelem();
1561 
1562  // Set to zero
1563  absv = 0;
1564  absy = 0;
1565  absw = 0;
1566 
1567  constexpr Numeric sq_ln2 = Constant::sqrt_ln_2;
1568  constexpr Numeric sq_ln2pi = sq_ln2 / Constant::sqrt_pi;
1569  const Numeric dens = xco2 * ptot * parameters::aMolAtm / temp;
1570  constexpr Numeric u_pi = Constant::inv_pi;
1571  constexpr Numeric u_sqln2pi = 1 / sq_ln2pi;
1572 
1573  for (Index iband=0; iband<cmn.Bands.nBand; iband++) {
1574  Numeric sigmoy=0;
1575  Numeric gamdmx=0;
1576  convtp(cmn, iband, cmn.Bands.Isot[iband], cmn.Bands.nLines[iband],
1577  xh2o, xco2, temp, ptot, sigmoy, gamdmx, mixfull, mixsdv);
1578  const Numeric sqrtm = std::sqrt(temp / parameters::aMass[cmn.Bands.Isot[iband]-1]);
1579 
1580  const Numeric ds0=(70.67e0+104.1e0*0.21e0*std::pow(gamdmx, 6.4))/(1.0+0.21*std::pow(gamdmx, 5.4));
1581  const Numeric ds2=(34.97e0+105.e0*9.e0*std::pow(gamdmx, 3.1))/(1.0+9.0*std::pow(gamdmx, 2.1));
1582 
1583  for (Index isig=0; isig<nsig; isig++) {
1584  const Numeric sigc = invcm_grid[isig];
1585 
1586  for (Index iline=0; iline<cmn.Bands.nLines[iband]; iline++) {
1587  const Numeric gamd=parameters::CtGamD*cmn.LineSg.Sig(iline,iband)*sqrtm;
1588  const Numeric cte = sq_ln2 / gamd;
1589  const Numeric popudipo = cmn.PopuT.PopuT[iline] * pow2(cmn.DipoTcm.DipoT(iline, iband));
1590  const Numeric dsigc = sigc - cmn.LineSg.Sig(iline, iband) - cmn.SHIFT.shft[iline] * ptot;
1591 
1592  if (mixsdv) {
1593  if (std::abs(dsigc) / cmn.GamT.HWT[iline] > ds0) {
1594  absv[isig] += popudipo * u_sqln2pi * u_pi * cmn.GamT.HWT[iline] / (pow2(cmn.GamT.HWT[iline]) + pow2(dsigc));
1595  absy[isig] += popudipo * u_sqln2pi * u_pi * (cmn.GamT.HWT[iline] + cmn.YLT.YT[iline]*dsigc) / (pow2(cmn.GamT.HWT[iline]) + pow2(dsigc));
1596  } else if (std::abs(dsigc) / cmn.GamT.HWT[iline] > ds2 and
1597  std::abs(dsigc) / cmn.GamT.HWT[iline] < ds0) {
1598  const Numeric order2r = cmn.GamT.HWT[iline] / (pow2(cmn.GamT.HWT[iline]) + pow2(dsigc))
1599  + 3 * pow2(cmn.GamT.HWSDV2T[iline]) * (pow3(cmn.GamT.HWT[iline]) - 3 * cmn.GamT.HWT[iline]*pow2(dsigc)) /
1600  (2*Constant::pi*pow3(pow2(cmn.GamT.HWT[iline]) + pow2(dsigc)));
1601 
1602  const Numeric order2i = dsigc / (pow2(cmn.GamT.HWT[iline]) + pow2(dsigc)) +
1603  3 * pow2(cmn.GamT.HWSDV2T[iline]) *(3*pow2(cmn.GamT.HWT[iline])*dsigc-pow3(dsigc)) /
1604  (2*Constant::pi*pow3(pow2(cmn.GamT.HWT[iline]) + pow2(dsigc)));
1605 
1606  absv[isig] += popudipo * u_sqln2pi * u_pi * order2r;
1607 
1608  absy[isig] += popudipo * u_sqln2pi * u_pi * (order2r + cmn.YLT.YT[iline]*order2i);
1609  } else {
1610  const Numeric shft0 = cmn.SHIFT.shft[iline] * ptot;
1611 
1612  Numeric wr=0, wi=0;
1613  qsdv(cmn.LineSg.Sig(iline, iband), gamd, cmn.GamT.HWT[iline], cmn.GamT.HWSDV2T[iline], shft0, 0, sigc, wr, wi);
1614 
1615  absv[isig] += popudipo * wr * u_sqln2pi;
1616 
1617  absy[isig] += popudipo * u_sqln2pi * (wr - cmn.YLT.YT[iline] * wi);
1618  }
1619  } else {
1620  if (std::abs(cmn.LineSg.Sig(iline, iband)-sigc) > (rdmult*gamd)) { // NOTE: Removed in updated version of the code...
1621  absv[isig] += popudipo * u_sqln2pi * u_pi * cmn.GamT.HWT[iline] / (pow2(cmn.GamT.HWT[iline]) + pow2(dsigc));
1622 
1623  absy[isig] += popudipo * u_sqln2pi * u_pi * (cmn.GamT.HWT[iline] + cmn.YLT.YT[iline] * dsigc) / (pow2(cmn.GamT.HWT[iline]) + pow2(dsigc));
1624  } else {
1625  const Numeric yy = cmn.GamT.HWT[iline] * cte;
1626  const Numeric xx= (cmn.LineSg.Sig(iline, iband)+cmn.SHIFT.shft[iline] * ptot - sigc) * cte;
1627  const Complex w = Faddeeva::w(Complex(xx, yy));
1628 
1629  absv[isig] += popudipo * w.real() / gamd;
1630 
1631  absy[isig] += popudipo * (w.real() - cmn.YLT.YT[iline] * w.imag()) / gamd;
1632  }
1633  }
1634 
1635  if (mixfull) {
1636  /*const Numeric gamd_int=parameters::CtGamD*(cmn.FicLPR.AlphR[iline] + sigmoy)*sqrtm;
1637  const Numeric cte_int = sq_ln2 / gamd_int;
1638  if (std::abs(cmn.FicLPI.AlphI[iline] + sigmoy - sigc) <= (rdmult*gamd_int)) {
1639  const Complex z = Complex(cmn.FicLPR.AlphR[iline] + sigmoy - sigc, cmn.FicLPI.AlphI[iline]) * cte_int;
1640  const Complex w = Faddeeva::w(z);
1641  absw[isig] += (cmn.FicLSR.SSR[iline] * w.real() - cmn.FicLSI.SSI[iline] * w.imag()) / cte_int;
1642  }*/
1643  absw[isig] += u_sqln2pi * u_pi * (cmn.Zss.ZS[iline] / (sigc - cmn.Zaa.ZA[iline])).imag();
1644  }
1645  }
1646  }
1647  }
1648 
1649  for (Index isig=0; isig<nsig; isig++) {
1650  const Numeric sigc = invcm_grid[isig];
1651  const Numeric fact = sigc * (1-std::exp(-parameters::Ct * sigc / temp));
1652  absv[isig] *= fact * dens * sq_ln2pi;
1653  absy[isig] *= fact * dens * sq_ln2pi;
1654  absw[isig] *= fact * dens * sq_ln2pi;
1655  }
1656 }
1657 
1658 Vector compabs(
1659  const Numeric T,
1660  const Numeric P,
1661  const HitranRelaxationMatrixData& hitran,
1662  const ArrayOfAbsorptionLines& bands,
1663  const ConstVectorView vmrs,
1664  const ConstVectorView f_grid,
1665  const SpeciesAuxData& partition_functions)
1666 {
1667  using Constant::pow2;
1668  using Constant::pow3;
1669 
1670  // Number of points to compute
1671  const Index nf = f_grid.nelem();
1672 
1673  // Set to zero
1674  Vector absorption(nf, 0);
1675 
1676  constexpr Numeric sq_ln2 = Constant::sqrt_ln_2;
1677  constexpr Numeric sq_ln2pi = sq_ln2 / Constant::sqrt_pi;
1678  const Numeric dens = vmrs[2] * number_density(P, T);
1679  constexpr Numeric u_pi = Constant::inv_pi;
1680  constexpr Numeric u_sqln2pi = 1 / sq_ln2pi;
1681 
1682  for (Index iband=0; iband<bands.nelem(); iband++) {
1683  auto tp = convtp(vmrs, hitran, bands[iband], T, P, partition_functions.getParamType(bands[iband].QuantumIdentity()), partition_functions.getParam(bands[iband].QuantumIdentity()));
1684  const Numeric GD_div_F0 = Linefunctions::DopplerConstant(T, bands[iband].SpeciesMass());
1685 
1686  const bool nolm = not bands[iband].DoLineMixing(P);
1687  const bool sdvp = bands[iband].LineShapeType() == LineShape::Type::SDVP;
1688  const bool vp = bands[iband].LineShapeType() == LineShape::Type::VP;
1689  const bool rosenkranz = bands[iband].Population() == Absorption::PopulationType::ByHITRANRosenkranzRelmat;
1690  const bool full = bands[iband].Population() == Absorption::PopulationType::ByHITRANFullRelmat;
1691 
1692  for (Index iv=0; iv<nf; iv++) {
1693  const Numeric f = f_grid[iv];
1694 
1695  Numeric a=0;
1696  for (Index iline=0; iline<bands[iband].NumLines(); iline++) {
1697  const Numeric gamd=GD_div_F0 * tp.f0[iline];
1698  const Numeric gamd_mod=GD_div_F0 * tp.eqv.zval[iline].real();
1699  const Numeric cte = sq_ln2 / gamd;
1700  const Numeric cte_mod = sq_ln2 / gamd_mod;
1701  const Numeric popudipo = tp.pop[iline] * pow2(tp.dip[iline]);
1702 
1703  if (rosenkranz and sdvp and nolm) {
1704  const Complex w = qsdv_si(tp.f0[iline], gamd, tp.hwt[iline], tp.hwt2[iline], tp.shft[iline], 0, f);
1705  a += popudipo * w.real() * u_sqln2pi;
1706  } else if (rosenkranz and sdvp) {
1707  const Complex w = qsdv_si(tp.f0[iline], gamd, tp.hwt[iline], tp.hwt2[iline], tp.shft[iline], 0, f);
1708  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
1709  } else if ((rosenkranz and vp and nolm) or (full and nolm)) {
1710  const Numeric yy = tp.hwt[iline] * cte;
1711  const Numeric xx = (tp.f0[iline]+tp.shft[iline] - f) * cte;
1712  const Complex w = Faddeeva::w(Complex(xx, yy));
1713  a += popudipo * w.real() / gamd;
1714  } else if (rosenkranz and vp) {
1715  const Numeric yy = tp.hwt[iline] * cte;
1716  const Numeric xx = (tp.f0[iline]+tp.shft[iline] - f) * cte;
1717  const Complex w = Faddeeva::w(Complex(xx, yy));
1718  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
1719  } else if (full and vp) {
1720  const Complex z = (tp.eqv.zval[iline]-f) * cte_mod;
1721  const Complex w = Faddeeva::w(z);
1722  a += (tp.eqv.zstr[iline] * w).real() / gamd_mod;
1723  } else if (full) {
1724  a += u_sqln2pi * u_pi * (tp.eqv.zstr[iline] / (f - tp.eqv.zval[iline])).imag();
1725  } else {
1726  throw std::runtime_error("Cannot understand the combination of calculations requested...");
1727  }
1728  }
1729 
1730  // Guard to not allow negative absorption inside a band
1731  if (a > 0)
1732  absorption[iv] += a;
1733  }
1734  }
1735 
1736  for (Index iv=0; iv<nf; iv++) {
1737  const Numeric f = f_grid[iv];
1738  const Numeric fact = f * (1 - stimulated_emission(T, f));
1739  absorption[iv] *= fact * dens * sq_ln2pi;
1740  }
1741 
1742  return absorption;
1743 }
1744 
1745 void detband(CommonBlock& cmn,
1746  const Numeric& sgminr,
1747  const Numeric& sgmaxr,
1748  const Numeric& stotmax,
1749  const String& basedir="data_new/")
1750 {
1751  cmn.Bands.nBand=0;
1752  std::ifstream fortranfile;
1753  fortranfile.open(basedir + "/BandInfo.dat");
1754 
1755  if (not fortranfile.is_open())
1756  throw std::runtime_error("Cannot read the file. Please make sure you point at BandInfo.dat basedir and have the right to read.");
1757 
1758  String line;
1759  getline(fortranfile, line);
1760  while (fortranfile.good()) {
1761  Numeric stot, sgmin, sgmax;
1762  Index isotr, lfr, lir, jmxp, jmxq, jmxr;
1763  char c11, c12, c21, c22, c31, c32, c41, c42, c51, c52, x;
1764  sscanf(line.c_str(),
1765  "%1ld"
1766  "%c%c" "%1ld" "%c%c"
1767  "%c%c" "%1ld" "%c%c"
1768  "%c%c"
1769  "%12lf"
1770  "%c" "%12lf"
1771  "%c" "%12lf"
1772  "%c%c%c%c%c%c%c%c"
1773  "%4ld" "%4ld" "%4ld",
1774  &isotr,
1775  &c11, &c12, &lfr, &c21, &c22,
1776  &c31, &c32, &lir, &c41, &c42,
1777  &c51, &c52,
1778  &stot,
1779  &x, &sgmin,
1780  &x, &sgmax,
1781  &x, &x, &x, &x, &x, &x, &x, &x,
1782  &jmxp, &jmxq, &jmxr);
1783 
1784  getline(fortranfile, line);
1785  if (stot < stotmax)
1786  continue;
1787 
1788  if ((sgminr < sgmax) and (sgmaxr > sgmin)) {
1789  cmn.Bands.Isot[cmn.Bands.nBand] = isotr;
1790  if (isotr == 0)
1791  cmn.Bands.Isot[cmn.Bands.nBand] = 10;
1792 
1793  cmn.Bands.li[cmn.Bands.nBand] = lir;
1794  cmn.Bands.lf[cmn.Bands.nBand] = lfr;
1795 
1796  if (jmxp < 40 or jmxr < 40) cmn.Bands.li[cmn.Bands.nBand] = 55;
1797  if (jmxq >= 0 and jmxq < 40) cmn.Bands.li[cmn.Bands.nBand] = 55;
1798  char name[15];
1799  sprintf(name, "S%ld%c%c%ld%c%c%c%c%ld%c%c%c%c", isotr, c11, c12, lfr, c21, c22, c31, c32, lir, c41, c42, c51, c52);
1800  cmn.Bands.BandFile[cmn.Bands.nBand] = name;
1801  cmn.Bands.nBand++;
1802  if (cmn.Bands.nBand > parameters::nBmx) {
1803  throw std::runtime_error("Too many bands");
1804  }
1805  }
1806  }
1807  fortranfile.close();
1808 }
1809 
1810 void readw(CommonBlock& cmn, const String& basedir="data_new/")
1811 {
1812  for (Index l=0; l<=8; l++) {
1813  for (Index ideltal=0; ideltal<=1; ideltal++) {
1814  const Index lli = l;
1815  const Index llf = l + ideltal;
1816 
1817  const String cr = std::to_string(lli) + std::to_string(llf);
1818 
1819  std::ifstream fortranfile;
1820  const String fname = basedir + String("/WTfit") + cr + String(".dat");
1821  fortranfile.open(fname.c_str());
1822 
1823  String line;
1824  getline(fortranfile, line);
1825  while (fortranfile.good()) {
1826  char sw0r[21], sb0r[21];
1827  Numeric dmaxdt, wtmax;
1828  Index jic, jfc, jipc, jfpc;
1829  sscanf(line.c_str(),
1830  "%20s" "%20s"
1831  "%14lf" "%14lf"
1832  "%4ld" "%4ld" "%4ld" "%4ld",
1833  sw0r, sb0r, &dmaxdt, &wtmax, &jic, &jfc, &jipc, &jfpc);
1834  String ssw0r = sw0r;
1835  std::replace(ssw0r.begin(), ssw0r.end(), 'D', 'E');
1836  const Numeric w0r = std::stod(ssw0r);
1837  String ssb0r = sb0r;
1838  std::replace(ssb0r.begin(), ssb0r.end(), 'D', 'E');
1839  const Numeric b0r = std::stod(ssb0r);
1840 
1841  getline(fortranfile, line);
1842 
1843  if (lli > cmn.Wfittedp.W0pp.nbooks() or
1844  llf > cmn.Wfittedp.W0pp.npages() or
1845  jic > cmn.Wfittedp.W0pp.nrows() or
1846  jipc > cmn.Wfittedp.W0pp.ncols())
1847  throw std::runtime_error("Out of bounds in reading...");
1848 
1849  if (jic > jfc and jipc > jfpc) {
1850  cmn.Wfittedp.W0pp(lli,llf,jic,jipc) = w0r;
1851  cmn.Bfittedp.B0pp(lli,llf,jic,jipc) = b0r;
1852  } else if(jic > jfc and jipc == jfpc) {
1853  cmn.Wfittedp.W0pq(lli,llf,jic,jipc) = w0r;
1854  cmn.Bfittedp.B0pq(lli,llf,jic,jipc) = b0r;
1855  } else if(jic > jfc and jipc < jfpc) {
1856  cmn.Wfittedp.W0pr(lli,llf,jic,jipc) = w0r;
1857  cmn.Bfittedp.B0pr(lli,llf,jic,jipc) = b0r;
1858  } else if(jic < jfc and jipc > jfpc) {
1859  cmn.Wfittedr.W0rp(lli,llf,jic,jipc) = w0r;
1860  cmn.Bfittedr.B0rp(lli,llf,jic,jipc) = b0r;
1861  } else if(jic < jfc and jipc == jfpc) {
1862  cmn.Wfittedr.W0rq(lli,llf,jic,jipc) = w0r;
1863  cmn.Bfittedr.B0rq(lli,llf,jic,jipc) = b0r;
1864  } else if(jic < jfc and jipc < jfpc) {
1865  cmn.Wfittedr.W0rr(lli,llf,jic,jipc) = w0r;
1866  cmn.Bfittedr.B0rr(lli,llf,jic,jipc) = b0r;
1867  } else if(jic == jfc and jipc > jfpc) {
1868  cmn.Wfittedq.W0qp(lli,llf,jic,jipc) = w0r;
1869  cmn.Bfittedq.B0qp(lli,llf,jic,jipc) = b0r;
1870  } else if(jic == jfc and jipc == jfpc) {
1871  cmn.Wfittedq.W0qq(lli,llf,jic,jipc) = w0r;
1872  cmn.Bfittedq.B0qq(lli,llf,jic,jipc) = b0r;
1873  } else if(jic == jfc and jipc < jfpc) {
1874  cmn.Wfittedq.W0qr(lli,llf,jic,jipc) = w0r;
1875  cmn.Bfittedq.B0qr(lli,llf,jic,jipc) = b0r;
1876  }
1877  }
1878  }
1879  }
1880 }
1881 
1882 Vector compute(const Numeric p, const Numeric t, const Numeric xco2, const Numeric xh2o, const ConstVectorView invcm_grid, const Numeric stotmax, const calctype type)
1883 {
1884  const Index n = invcm_grid.nelem();
1885  Vector absorption(n);
1886  if (not n)
1887  return absorption;
1888 
1889  // Common setup and main IO
1890  CommonBlock cmn;
1891  detband(cmn, invcm_grid[0], invcm_grid[n-1], stotmax);
1892  readw(cmn);
1893  readlines(cmn);
1894 
1895  Vector absv(n);
1896  Vector absy(n);
1897  Vector absw(n);
1898  if (type == calctype::FullVP or type == calctype::FullRosenkranz or type == calctype::FullW)
1899  compabs(cmn, t,p, xco2, xh2o, invcm_grid, false, true, absv, absy, absw);
1900  else if (type == calctype::SDVP or type == calctype::SDRosenkranz or type == calctype::SDW)
1901  compabs(cmn, t,p, xco2, xh2o, invcm_grid, true, false, absv, absy, absw);
1902  else if (type == calctype::NoneVP or type == calctype::NoneRosenkranz or type == calctype::NoneW)
1903  compabs(cmn, t,p, xco2, xh2o, invcm_grid, false, false, absv, absy, absw);
1904 
1905  switch(type) {
1906  case calctype::FullVP:
1907  case calctype::SDVP:
1908  case calctype::NoneVP:
1909  for (Index i=0; i<n; i++)
1910  absorption[i] = absv[i];
1911  break;
1912  case calctype::FullRosenkranz:
1913  case calctype::SDRosenkranz:
1914  case calctype::NoneRosenkranz:
1915  for (Index i=0; i<n; i++)
1916  absorption[i] = absy[i];
1917  break;
1918  case calctype::FullW:
1919  case calctype::SDW:
1920  case calctype::NoneW:
1921  for (Index i=0; i<n; i++)
1922  absorption[i] = absw[i];
1923  break;
1924  }
1925 
1926  return absorption;
1927 }
1928 
1929 Vector compute(const HitranRelaxationMatrixData& hitran,
1930  const ArrayOfAbsorptionLines& bands,
1931  const Numeric P,
1932  const Numeric T,
1933  const ConstVectorView vmrs,
1934  const ConstVectorView f_grid,
1935  const SpeciesAuxData& partition_functions)
1936 {
1937  return f_grid.nelem() ? compabs(T, P, hitran, bands, vmrs, f_grid, partition_functions) : Vector(0);
1938 }
1939 
1940 void read(HitranRelaxationMatrixData& hitran, ArrayOfAbsorptionLines& bands, const String& basedir, const Numeric linemixinglimit, const Numeric fmin, const Numeric fmax, const Numeric stot, const ModeOfLineMixing mode)
1941 {
1942  String newbase = basedir;
1943  if (newbase.nelem() == 0)
1944  newbase = ".";
1945 
1946  CommonBlock cmn;
1947  detband(cmn, Conversion::freq2kaycm(fmin), Conversion::freq2kaycm(fmax), Conversion::arts2hitran_linestrength(stot), newbase);
1948  readw(cmn, newbase);
1949  readlines(cmn, newbase);
1950 
1951  Numeric linemixinglimit_internal=0;
1952  switch(mode) {
1953  case ModeOfLineMixing::VP:
1954  case ModeOfLineMixing::SDVP: break;
1955  case ModeOfLineMixing::FullW:
1956  case ModeOfLineMixing::VP_W:
1957  case ModeOfLineMixing::VP_Y:
1958  case ModeOfLineMixing::SDVP_Y: linemixinglimit_internal=linemixinglimit; break;
1959  default: throw std::runtime_error("Bad mode input. Must update function.");
1960  }
1961 
1962  const auto lstype = typeLP(mode) ?
1963  LineShape::Type::LP :
1964  (typeVP(mode) ?
1965  LineShape::Type::VP :
1966  LineShape::Type::SDVP) ;
1967  const auto poptype = typeFull(mode) ?
1968  Absorption::PopulationType::ByHITRANFullRelmat :
1969  Absorption::PopulationType::ByHITRANRosenkranzRelmat ;
1970 
1971  auto specs=stdarrayify(
1972  SpeciesTag("CO2-626"),
1973  SpeciesTag("CO2-636"),
1974  SpeciesTag("CO2-628"),
1975  SpeciesTag("CO2-627"),
1976  SpeciesTag("CO2-638"),
1977  SpeciesTag("CO2-637"),
1978  SpeciesTag("CO2-828"),
1979  SpeciesTag("CO2-728"),
1980  SpeciesTag("CO2-727"),
1981  SpeciesTag("CO2-838"),
1982  SpeciesTag("CO2-837"));
1983 
1984  // Move data from Fortran-style common block to ARTS
1985  hitran.B0pp = std::move(cmn.Bfittedp.B0pp);
1986  hitran.B0pq = std::move(cmn.Bfittedp.B0pq);
1987  hitran.B0pr = std::move(cmn.Bfittedp.B0pr);
1988  hitran.W0pp = std::move(cmn.Wfittedp.W0pp);
1989  hitran.W0pq = std::move(cmn.Wfittedp.W0pq);
1990  hitran.W0pr = std::move(cmn.Wfittedp.W0pr);
1991  hitran.B0qp = std::move(cmn.Bfittedq.B0qp);
1992  hitran.B0qq = std::move(cmn.Bfittedq.B0qq);
1993  hitran.B0qr = std::move(cmn.Bfittedq.B0qr);
1994  hitran.W0qp = std::move(cmn.Wfittedq.W0qp);
1995  hitran.W0qq = std::move(cmn.Wfittedq.W0qq);
1996  hitran.W0qr = std::move(cmn.Wfittedq.W0qr);
1997  hitran.B0rp = std::move(cmn.Bfittedr.B0rp);
1998  hitran.B0rq = std::move(cmn.Bfittedr.B0rq);
1999  hitran.B0rr = std::move(cmn.Bfittedr.B0rr);
2000  hitran.W0rp = std::move(cmn.Wfittedr.W0rp);
2001  hitran.W0rq = std::move(cmn.Wfittedr.W0rq);
2002  hitran.W0rr = std::move(cmn.Wfittedr.W0rr);
2003 
2004  // Reshape to band count size
2005  bands.resize(cmn.Bands.nBand);
2006  for (Index i{0}; i<cmn.Bands.nBand; i++) {
2007  QuantumNumbers outer_upper;
2008  outer_upper[QuantumNumberType::l1] = Rational(cmn.Bands.li[i]);
2009  outer_upper[QuantumNumberType::v1] = cmn.UnusedBandParams.iv1[i];
2010  outer_upper[QuantumNumberType::v2] = cmn.UnusedBandParams.iv2[i];
2011  outer_upper[QuantumNumberType::l2] = cmn.UnusedBandParams.il2[i];
2012  outer_upper[QuantumNumberType::v3] = cmn.UnusedBandParams.iv3[i];
2013  outer_upper[QuantumNumberType::r] = cmn.UnusedBandParams.ir[i];
2014  QuantumNumbers outer_lower;
2015  outer_lower[QuantumNumberType::l1] = Rational(cmn.Bands.lf[i]);
2016  outer_lower[QuantumNumberType::v1] = cmn.UnusedBandParams.fv1[i];
2017  outer_lower[QuantumNumberType::v2] = cmn.UnusedBandParams.fv2[i];
2018  outer_lower[QuantumNumberType::l2] = cmn.UnusedBandParams.fl2[i];
2019  outer_lower[QuantumNumberType::v3] = cmn.UnusedBandParams.fv3[i];
2020  outer_lower[QuantumNumberType::r] = cmn.UnusedBandParams.fr[i];
2021 
2022  bands[i] = {true,
2023  true,
2024  Absorption::CutoffType::None,
2025  Absorption::MirroringType::None,
2026  poptype,
2027  Absorption::NormalizationType::None,
2028  lstype,
2029  296,
2030  -1,
2031  linemixinglimit_internal,
2032  {specs[cmn.Bands.Isot[i]-1].Species(),
2033  specs[cmn.Bands.Isot[i]-1].Isotopologue(),
2034  outer_upper, outer_lower},
2035  {QuantumNumberType::J},
2036  {SpeciesTag("N2"), SpeciesTag("H2O"), specs[cmn.Bands.Isot[i]-1]}};
2037 
2038  bands[i].AllLines().resize(cmn.Bands.nLines[i]);
2039  for (Index j{0}; j<cmn.Bands.nLines[i]; j++) {
2040  const std::vector<Rational> inner_upper{Rational(cmn.Jiln.Ji(j, i))};
2041  const std::vector<Rational> inner_lower{Rational(cmn.Jfln.Jf(j, i))};
2042 
2043  const LineShape::ModelParameters G0_sdvp_air{LineShape::TemperatureModel::T1,
2044  Conversion::hitran2arts_broadening(cmn.GamSDVT0AIR.HWSDVT0AIR(j, i)),
2045  cmn.DTGAMAIR.BHWAIR(j, i)};
2046  const LineShape::ModelParameters G0_sdvp_h2o{LineShape::TemperatureModel::T1,
2047  Conversion::hitran2arts_broadening(cmn.GamSDVT0H2O.HWSDVT0H2O(j, i)),
2048  cmn.DTGAMH2O.BHWH2O(j, i)};
2049  const LineShape::ModelParameters G0_sdvp_co2{LineShape::TemperatureModel::T1,
2050  Conversion::hitran2arts_broadening(cmn.GamSDVT0CO2.HWSDVT0SELF(j, i)),
2051  cmn.DTGAMCO2.BHWSELF(j, i)};
2052 
2053  LineShape::ModelParameters G2_sdvp_air{G0_sdvp_air}; G2_sdvp_air.X0 *= cmn.GamSDVT0AIR.rHWT0AIR(j, i);
2054  LineShape::ModelParameters G2_sdvp_h2o{G0_sdvp_h2o}; G2_sdvp_h2o.X0 *= cmn.GamSDVT0H2O.rHWT0H2O(j, i);
2055  LineShape::ModelParameters G2_sdvp_co2{G0_sdvp_co2}; G2_sdvp_co2.X0 *= cmn.GamSDVT0CO2.rHWT0SELF(j, i);
2056 
2057  const LineShape::ModelParameters D0{LineShape::TemperatureModel::T0, Conversion::hitran2arts_broadening(cmn.SHIFT0.shft0(j, i))};
2058 
2059  const LineShape::ModelParameters G0_vp_air{LineShape::TemperatureModel::T1,
2060  Conversion::hitran2arts_broadening(cmn.GamVT0AIR.HWVT0AIR(j, i)),
2061  cmn.DTGAMAIR.BHWAIR(j, i)};
2062  const LineShape::ModelParameters G0_vp_h2o{LineShape::TemperatureModel::T1,
2063  Conversion::hitran2arts_broadening(cmn.GamVT0H2O.HWVT0H2O(j, i)),
2064  cmn.DTGAMH2O.BHWH2O(j, i)};
2065  const LineShape::ModelParameters G0_vp_co2{LineShape::TemperatureModel::T1,
2066  Conversion::hitran2arts_broadening(cmn.GamVT0CO2.HWVT0SELF(j, i)),
2067  cmn.DTGAMCO2.BHWSELF(j, i)};
2068 
2069  const LineShape::SingleSpeciesModel sdvp_air{G0_sdvp_air, D0, G2_sdvp_air};
2070  const LineShape::SingleSpeciesModel sdvp_h2o{G0_sdvp_h2o, D0, G2_sdvp_h2o};
2071  const LineShape::SingleSpeciesModel sdvp_co2{G0_sdvp_co2, D0, G2_sdvp_co2};
2072  const LineShape::SingleSpeciesModel vp_air{G0_vp_air, D0};
2073  const LineShape::SingleSpeciesModel vp_h2o{G0_vp_h2o, D0};
2074  const LineShape::SingleSpeciesModel vp_co2{G0_vp_co2, D0};
2075  const auto lsmodel = typeVP(mode) ?
2076  LineShape::Model{{vp_air, vp_h2o, vp_co2}} :
2077  LineShape::Model{{sdvp_air, sdvp_h2o, sdvp_co2}};
2078 
2079  Numeric qt0_co2, gsi0;
2080  qt_co2(parameters::T0, cmn.Bands.Isot[i], gsi0, qt0_co2);
2081 
2082  // Should probably use the isotopologue ratio
2083  bands[i].Line(j) = {Conversion::kaycm2freq(cmn.LineSg.Sig(j, i)),
2084  Conversion::hitran2arts_linestrength(cmn.UnusedBandParams.intens(j, i)),
2085  Conversion::hitran2arts_energy(cmn.Energy.E(j, i)),
2086  gsi0*Numeric(cmn.Jfln.Jf(j, i) * 2 + 1),
2087  gsi0*Numeric(cmn.Jiln.Ji(j, i) * 2 + 1),
2088  cmn.UnusedBandParams.eina(j, i),
2089  ZeemanModel{},
2090  lsmodel,
2091  inner_lower,
2092  inner_upper};
2093  }
2094  }
2095 }
2096 };
2097 
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
lm_hitran_2017
Definition: linemixing_hitran.cc:35
lin_alg.h
Linear algebra functions.
linemixing_hitran.h
Namespace and functions to deal with HITRAN linemixing.
linefunctions.h
Stuff related to lineshape functions.
parameters
Parameters parameters
Holds the command line parameters.
Definition: parameters.cc:41
physics_funcs.h
This file contains declerations of functions of physical character.