ARTS  2.3.1285(git:92a29ea9-dirty)
predefined_absorption_models.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 "predefined_absorption_models.h"
28 #include "lin_alg.h"
29 #include "linescaling.h"
30 #include "wigner_functions.h"
31 
32 
33 constexpr auto necs2020 = 38;
34 constexpr auto nlines_mpm2020 = necs2020 + 6;
35 
36 
37 constexpr LineShape::SingleSpeciesModel init_mpm2020_slsm(Numeric g00,
38  Numeric y0,
39  Numeric y1,
40  Numeric g0,
41  Numeric g1,
42  Numeric dv0,
43  Numeric dv1,
44  Numeric x)
45 {
46  return LineShape::SingleSpeciesModel(
47  {LineShape::TemperatureModel::T1, g00, x, NAN, NAN},
48  {LineShape::TemperatureModel::None, NAN, NAN, NAN, NAN},
49  {LineShape::TemperatureModel::None, NAN, NAN, NAN, NAN},
50  {LineShape::TemperatureModel::None, NAN, NAN, NAN, NAN},
51  {LineShape::TemperatureModel::None, NAN, NAN, NAN, NAN},
52  {LineShape::TemperatureModel::None, NAN, NAN, NAN, NAN},
53  {LineShape::TemperatureModel::T4, y0, y1, x, NAN},
54  {LineShape::TemperatureModel::T4, g0, g1, 2*x, NAN},
55  {LineShape::TemperatureModel::T4, dv0, dv1, 2*x, NAN});
56 }
57 
58 
59 constexpr std::array<LineShape::SingleSpeciesModel, nlines_mpm2020> init_mpm2020_lsm()
60 {
61  // Pressure broadening [1/Pa] at reference temperature
62  constexpr std::array<Numeric, nlines_mpm2020> g00 =
63  {1.685E+4, 1.703E+4, 1.513E+4, 1.495E+4, 1.433E+4,
64  1.408E+4, 1.353E+4, 1.353E+4, 1.303E+4, 1.319E+4,
65  1.262E+4, 1.265E+4, 1.238E+4, 1.217E+4, 1.207E+4,
66  1.207E+4, 1.137E+4, 1.137E+4, 1.101E+4, 1.101E+4,
67  1.037E+4, 1.038E+4, 9.96E+3, 9.96E+3, 9.55E+3,
68  9.55E+3, 9.06E+3, 9.06E+3, 8.58E+3, 8.58E+3,
69  8.11E+3, 8.11E+3, 7.64E+3, 7.64E+3, 7.17E+3,
70  7.17E+3, 6.69E+3, 6.69E+3, 1.64E+4, 1.64E+4,
71  1.60E+4, 1.60E+4, 1.62E+4, 1.47E+4,};
72 
73  // First order line mixing first coefficient [1/Pa] at reference temperature
74  constexpr std::array<Numeric, nlines_mpm2020> y0 =
75  {-4.1E-7, 0.00000277, -0.00000372, 0.00000559, -0.00000573,
76  0.00000618, -0.00000366, 0.00000278, -8.9E-7, -2.1E-7,
77  6.0E-7, -0.00000152, 0.00000216, -0.00000293, 0.00000373,
78  -0.00000436, 0.00000491, -0.00000542, 0.00000571, -0.00000613,
79  0.00000636, -0.00000670, 0.00000690, -0.00000718, 0.00000740,
80  -0.00000763, 0.00000788, -0.00000807, 0.00000834, -0.00000849,
81  0.00000876, -0.00000887, 0.00000915, -0.00000922, 0.00000950,
82  -0.00000955, 0.00000987, -0.00000988, 0.00000, 0.00000,
83  0.00000, 0.00000, 0.00000, 0.00000,};
84 
85  // First order line mixing second coefficient [1/Pa] at reference temperature
86  constexpr std::array<Numeric, nlines_mpm2020> y1 =
87  {0.00000, 0.00000124, -2E-8, 8E-8, 4.5E-7,
88  -9.3E-7, 0.00000264, -0.00000351, 0.00000359, -0.00000416,
89  0.00000326, -0.00000353, 0.00000484, -0.00000503, 0.00000579,
90  -0.00000590, 0.00000616, -0.00000619, 0.00000611, -0.00000609,
91  0.00000574, -0.00000568, 0.00000574, -0.00000566, 0.0000060,
92  -0.0000059, 0.0000063, -0.0000062, 0.0000064, -0.0000063,
93  0.0000065, -0.0000064, 0.0000065, -0.0000064, 0.0000065,
94  -0.0000064, 0.0000064, -0.0000062, 0.00000, 0.00000,
95  0.00000, 0.00000, 0.00000, 0.00000, };
96 
97  // Second order line mixing strength adjustment first coefficient [1/Pa^2] at reference temperature
98  constexpr std::array<Numeric, nlines_mpm2020> g0 =
99  {-6.95E-14, -9.0E-12, -1.03E-11, -2.39E-11, -1.72E-11,
100  -1.71E-11, 2.8E-12, 1.50E-11, 1.32E-11, 1.70E-11,
101  8.7E-12, 6.9E-12, 8.3E-12, 6.7E-12, 7E-13,
102  1.6E-12, -2.1E-12, -6.6E-12, -9.5E-12, -1.15E-11,
103  -1.18E-11, -1.40E-11, -1.73E-11, -1.86E-11, -2.17E-11,
104  -2.27E-11, -2.34E-11, -2.42E-11, -2.66E-11, -2.72E-11,
105  -3.01E-11, -3.04E-11, -3.34E-11, -3.33E-11, -3.61E-11,
106  -3.58E-11, -3.48E-11, -3.44E-11, 0E-10, 0E-10,
107  0E-10, 0E-10, 0E-10, 0E-10,};
108 
109  // Second order line mixing strength adjustment second coefficient [1/Pa^2] at reference temperature
110  constexpr std::array<Numeric, nlines_mpm2020> g1 =
111  {0E-10, -4.5E-12, 7E-13, 3.3E-12, 8.1E-12,
112  1.62E-11, 1.79E-11, 2.25E-11, 5.4E-12, 3E-13,
113  4E-14, -4.7E-12, -3.4E-12, -7.1E-12, -1.80E-11,
114  -2.10E-11, -2.85E-11, -3.23E-11, -3.63E-11, -3.80E-11,
115  -3.78E-11, -3.87E-11, -3.92E-11, -3.94E-11, -4.24E-11,
116  -4.22E-11, -4.65E-11, -4.6E-11, -5.1E-11, -5.0E-11,
117  -5.5E-11, -5.4E-11, -5.8E-11, -5.6E-11, -6.2E-11,
118  -5.9E-11, -6.8E-11, -6.5E-11, 0E-10, 0E-10,
119  0E-10, 0E-10, 0E-10, 0E-10, };
120 
121  // Second order line mixing frequency adjustment first coefficient [Hz/Pa^2] at reference temperature
122  constexpr std::array<Numeric, nlines_mpm2020> dv0 =
123  {-0.000028, 0.000597, -0.00195, 0.0032, -0.00475,
124  0.00541, -0.00232, 0.00154, 0.00007, -0.00084,
125  -0.00025, -0.00014, -0.00004, -0.00020, 0.0005,
126  -0.00066, 0.00072, -0.0008, 0.00064, -0.00070,
127  0.00056, -0.00060, 0.00047, -0.00049, 0.00040,
128  -0.00041, 0.00036, -0.00037, 0.00033, -0.00034,
129  0.00032, -0.00032, 0.00030, -0.00030, 0.00028,
130  -0.00029, 0.00029, -0.00029, 0.0, 0.0,
131  0.0, 0.0, 0.0, 0.0, };
132 
133  // Second order line mixing frequency adjustment second coefficient [Hz/Pa^2] at reference temperature
134  constexpr std::array<Numeric, nlines_mpm2020> dv1 =
135  {-0.000039, 0.0009, -0.0012, 0.0016, -0.0027,
136  0.0029, 0.0006, -0.0015, 0.0010, -0.0014,
137  -0.0013, 0.0013, 0.0004, -0.0005, 0.0010,
138  -0.0010, 0.0010, -0.0011, 0.0008, -0.0009,
139  0.0003, -0.0003, 0.00009, -0.00009, 0.00017,
140  -0.00016, 0.00024, -0.00023, 0.00024, -0.00024,
141  0.00024, -0.00020, 0.00017, -0.00016, 0.00013,
142  -0.00012, 0.00005, -0.00004, 0.0, 0.0,
143  0.0, 0.0, 0.0, 0.0, };
144 
145  // Temperature scaling exponent [-]
146  constexpr Numeric x = 0.754;
147 
148  return {init_mpm2020_slsm(g00[0], y0[0], y1[0], g0[0], g1[0], dv0[0], dv1[0], x),
149  init_mpm2020_slsm(g00[1], y0[1], y1[1], g0[1], g1[1], dv0[1], dv1[1], x),
150  init_mpm2020_slsm(g00[2], y0[2], y1[2], g0[2], g1[2], dv0[2], dv1[2], x),
151  init_mpm2020_slsm(g00[3], y0[3], y1[3], g0[3], g1[3], dv0[3], dv1[3], x),
152  init_mpm2020_slsm(g00[4], y0[4], y1[4], g0[4], g1[4], dv0[4], dv1[4], x),
153  init_mpm2020_slsm(g00[5], y0[5], y1[5], g0[5], g1[5], dv0[5], dv1[5], x),
154  init_mpm2020_slsm(g00[6], y0[6], y1[6], g0[6], g1[6], dv0[6], dv1[6], x),
155  init_mpm2020_slsm(g00[7], y0[7], y1[7], g0[7], g1[7], dv0[7], dv1[7], x),
156  init_mpm2020_slsm(g00[8], y0[8], y1[8], g0[8], g1[8], dv0[8], dv1[8], x),
157  init_mpm2020_slsm(g00[9], y0[9], y1[9], g0[9], g1[9], dv0[9], dv1[9], x),
158  init_mpm2020_slsm(g00[10], y0[10], y1[10], g0[10], g1[10], dv0[10], dv1[10], x),
159  init_mpm2020_slsm(g00[11], y0[11], y1[11], g0[11], g1[11], dv0[11], dv1[11], x),
160  init_mpm2020_slsm(g00[12], y0[12], y1[12], g0[12], g1[12], dv0[12], dv1[12], x),
161  init_mpm2020_slsm(g00[13], y0[13], y1[13], g0[13], g1[13], dv0[13], dv1[13], x),
162  init_mpm2020_slsm(g00[14], y0[14], y1[14], g0[14], g1[14], dv0[14], dv1[14], x),
163  init_mpm2020_slsm(g00[15], y0[15], y1[15], g0[15], g1[15], dv0[15], dv1[15], x),
164  init_mpm2020_slsm(g00[16], y0[16], y1[16], g0[16], g1[16], dv0[16], dv1[16], x),
165  init_mpm2020_slsm(g00[17], y0[17], y1[17], g0[17], g1[17], dv0[17], dv1[17], x),
166  init_mpm2020_slsm(g00[18], y0[18], y1[18], g0[18], g1[18], dv0[18], dv1[18], x),
167  init_mpm2020_slsm(g00[19], y0[19], y1[19], g0[19], g1[19], dv0[19], dv1[19], x),
168  init_mpm2020_slsm(g00[20], y0[20], y1[20], g0[20], g1[20], dv0[20], dv1[20], x),
169  init_mpm2020_slsm(g00[21], y0[21], y1[21], g0[21], g1[21], dv0[21], dv1[21], x),
170  init_mpm2020_slsm(g00[22], y0[22], y1[22], g0[22], g1[22], dv0[22], dv1[22], x),
171  init_mpm2020_slsm(g00[23], y0[23], y1[23], g0[23], g1[23], dv0[23], dv1[23], x),
172  init_mpm2020_slsm(g00[24], y0[24], y1[24], g0[24], g1[24], dv0[24], dv1[24], x),
173  init_mpm2020_slsm(g00[25], y0[25], y1[25], g0[25], g1[25], dv0[25], dv1[25], x),
174  init_mpm2020_slsm(g00[26], y0[26], y1[26], g0[26], g1[26], dv0[26], dv1[26], x),
175  init_mpm2020_slsm(g00[27], y0[27], y1[27], g0[27], g1[27], dv0[27], dv1[27], x),
176  init_mpm2020_slsm(g00[28], y0[28], y1[28], g0[28], g1[28], dv0[28], dv1[28], x),
177  init_mpm2020_slsm(g00[29], y0[29], y1[29], g0[29], g1[29], dv0[29], dv1[29], x),
178  init_mpm2020_slsm(g00[30], y0[30], y1[30], g0[30], g1[30], dv0[30], dv1[30], x),
179  init_mpm2020_slsm(g00[31], y0[31], y1[31], g0[31], g1[31], dv0[31], dv1[31], x),
180  init_mpm2020_slsm(g00[32], y0[32], y1[32], g0[32], g1[32], dv0[32], dv1[32], x),
181  init_mpm2020_slsm(g00[33], y0[33], y1[33], g0[33], g1[33], dv0[33], dv1[33], x),
182  init_mpm2020_slsm(g00[34], y0[34], y1[34], g0[34], g1[34], dv0[34], dv1[34], x),
183  init_mpm2020_slsm(g00[35], y0[35], y1[35], g0[35], g1[35], dv0[35], dv1[35], x),
184  init_mpm2020_slsm(g00[36], y0[36], y1[36], g0[36], g1[36], dv0[36], dv1[36], x),
185  init_mpm2020_slsm(g00[37], y0[37], y1[37], g0[37], g1[37], dv0[37], dv1[37], x),
186  init_mpm2020_slsm(g00[38], y0[38], y1[38], g0[38], g1[38], dv0[38], dv1[38], x),
187  init_mpm2020_slsm(g00[39], y0[39], y1[39], g0[39], g1[39], dv0[39], dv1[39], x),
188  init_mpm2020_slsm(g00[40], y0[40], y1[40], g0[40], g1[40], dv0[40], dv1[40], x),
189  init_mpm2020_slsm(g00[41], y0[41], y1[41], g0[41], g1[41], dv0[41], dv1[41], x),
190  init_mpm2020_slsm(g00[42], y0[42], y1[42], g0[42], g1[42], dv0[42], dv1[42], x),
191  init_mpm2020_slsm(g00[43], y0[43], y1[43], g0[43], g1[43], dv0[43], dv1[43], x),};
192 }
193 
194 
195 constexpr QuantumIdentifier init_mpm2020_qid(Index species, Index isot, Rational Jup, Rational Jlow, Rational Nup, Rational Nlow)
196 {
197  return QuantumIdentifier(species, isot, QuantumNumbers(Jup, Nup, 0_rat), QuantumNumbers(Jlow, Nlow, 0_rat));
198 }
199 
200 
201 constexpr std::array<QuantumIdentifier, nlines_mpm2020> init_mpm2020_qids(const Index& species, const Index& isot)
202 {
203  // N of upper level
204  constexpr std::array<Index, nlines_mpm2020> Np = {
205  1, 1, 3, 3, 5, 5, 7, 7, 9, 9,
206  11, 11, 13, 13, 15, 15, 17, 17,
207  19, 19, 21, 21, 23, 23, 25, 25,
208  27, 27, 29, 29, 31, 31, 33, 33,
209  35, 35, 37, 37, 1, 1, 1, 3, 3, 3};
210 
211  // N of lower level
212  constexpr std::array<Index, nlines_mpm2020> Npp = {
213  1, 1, 3, 3, 5, 5, 7, 7, 9, 9,
214  11, 11, 13, 13, 15, 15, 17, 17,
215  19, 19, 21, 21, 23, 23, 25, 25,
216  27, 27, 29, 29, 31, 31, 33, 33,
217  35, 35, 37, 37, 3, 3, 3, 5, 5, 5};
218 
219  // J of upper level
220  constexpr std::array<Index, nlines_mpm2020> Jp = {
221  1, 1, 3, 3, 5, 5, 7, 7, 9, 9,
222  11, 11, 13, 13, 15, 15, 17, 17,
223  19, 19, 21, 21, 23, 23, 25, 25,
224  27, 27, 29, 29, 31, 31, 33, 33,
225  35, 35, 37, 37, 1, 2, 2, 3, 4, 4};
226 
227  // J of lower level
228  constexpr std::array<Index, nlines_mpm2020> Jpp = {
229  0, 2, 2, 4, 4, 6, 6, 8, 8, 10,
230  10, 12, 12, 14, 14, 16, 16, 18,
231  18, 20, 20, 22, 22, 24, 24, 26,
232  26, 28, 28, 30, 30, 32, 32, 34,
233  34, 36, 36, 38, 2, 2, 3, 4, 4, 5};
234 
235  return {QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[0]), Rational(Np[0]), 0_rat), QuantumNumbers(Rational(Jpp[0]), Rational(Npp[0]), 0_rat)),
236  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[1]), Rational(Np[1]), 0_rat), QuantumNumbers(Rational(Jpp[1]), Rational(Npp[1]), 0_rat)),
237  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[2]), Rational(Np[2]), 0_rat), QuantumNumbers(Rational(Jpp[2]), Rational(Npp[2]), 0_rat)),
238  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[3]), Rational(Np[3]), 0_rat), QuantumNumbers(Rational(Jpp[3]), Rational(Npp[3]), 0_rat)),
239  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[4]), Rational(Np[4]), 0_rat), QuantumNumbers(Rational(Jpp[4]), Rational(Npp[4]), 0_rat)),
240  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[5]), Rational(Np[5]), 0_rat), QuantumNumbers(Rational(Jpp[5]), Rational(Npp[5]), 0_rat)),
241  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[6]), Rational(Np[6]), 0_rat), QuantumNumbers(Rational(Jpp[6]), Rational(Npp[6]), 0_rat)),
242  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[7]), Rational(Np[7]), 0_rat), QuantumNumbers(Rational(Jpp[7]), Rational(Npp[7]), 0_rat)),
243  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[8]), Rational(Np[8]), 0_rat), QuantumNumbers(Rational(Jpp[8]), Rational(Npp[8]), 0_rat)),
244  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[9]), Rational(Np[9]), 0_rat), QuantumNumbers(Rational(Jpp[9]), Rational(Npp[9]), 0_rat)),
245  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[10]), Rational(Np[10]), 0_rat), QuantumNumbers(Rational(Jpp[10]), Rational(Npp[10]), 0_rat)),
246  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[11]), Rational(Np[11]), 0_rat), QuantumNumbers(Rational(Jpp[11]), Rational(Npp[11]), 0_rat)),
247  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[12]), Rational(Np[12]), 0_rat), QuantumNumbers(Rational(Jpp[12]), Rational(Npp[12]), 0_rat)),
248  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[13]), Rational(Np[13]), 0_rat), QuantumNumbers(Rational(Jpp[13]), Rational(Npp[13]), 0_rat)),
249  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[14]), Rational(Np[14]), 0_rat), QuantumNumbers(Rational(Jpp[14]), Rational(Npp[14]), 0_rat)),
250  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[15]), Rational(Np[15]), 0_rat), QuantumNumbers(Rational(Jpp[15]), Rational(Npp[15]), 0_rat)),
251  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[16]), Rational(Np[16]), 0_rat), QuantumNumbers(Rational(Jpp[16]), Rational(Npp[16]), 0_rat)),
252  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[17]), Rational(Np[17]), 0_rat), QuantumNumbers(Rational(Jpp[17]), Rational(Npp[17]), 0_rat)),
253  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[18]), Rational(Np[18]), 0_rat), QuantumNumbers(Rational(Jpp[18]), Rational(Npp[18]), 0_rat)),
254  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[19]), Rational(Np[19]), 0_rat), QuantumNumbers(Rational(Jpp[19]), Rational(Npp[19]), 0_rat)),
255  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[20]), Rational(Np[20]), 0_rat), QuantumNumbers(Rational(Jpp[20]), Rational(Npp[20]), 0_rat)),
256  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[21]), Rational(Np[21]), 0_rat), QuantumNumbers(Rational(Jpp[21]), Rational(Npp[21]), 0_rat)),
257  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[22]), Rational(Np[22]), 0_rat), QuantumNumbers(Rational(Jpp[22]), Rational(Npp[22]), 0_rat)),
258  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[23]), Rational(Np[23]), 0_rat), QuantumNumbers(Rational(Jpp[23]), Rational(Npp[23]), 0_rat)),
259  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[24]), Rational(Np[24]), 0_rat), QuantumNumbers(Rational(Jpp[24]), Rational(Npp[24]), 0_rat)),
260  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[25]), Rational(Np[25]), 0_rat), QuantumNumbers(Rational(Jpp[25]), Rational(Npp[25]), 0_rat)),
261  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[26]), Rational(Np[26]), 0_rat), QuantumNumbers(Rational(Jpp[26]), Rational(Npp[26]), 0_rat)),
262  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[27]), Rational(Np[27]), 0_rat), QuantumNumbers(Rational(Jpp[27]), Rational(Npp[27]), 0_rat)),
263  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[28]), Rational(Np[28]), 0_rat), QuantumNumbers(Rational(Jpp[28]), Rational(Npp[28]), 0_rat)),
264  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[29]), Rational(Np[29]), 0_rat), QuantumNumbers(Rational(Jpp[29]), Rational(Npp[29]), 0_rat)),
265  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[30]), Rational(Np[30]), 0_rat), QuantumNumbers(Rational(Jpp[30]), Rational(Npp[30]), 0_rat)),
266  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[31]), Rational(Np[31]), 0_rat), QuantumNumbers(Rational(Jpp[31]), Rational(Npp[31]), 0_rat)),
267  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[32]), Rational(Np[32]), 0_rat), QuantumNumbers(Rational(Jpp[32]), Rational(Npp[32]), 0_rat)),
268  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[33]), Rational(Np[33]), 0_rat), QuantumNumbers(Rational(Jpp[33]), Rational(Npp[33]), 0_rat)),
269  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[34]), Rational(Np[34]), 0_rat), QuantumNumbers(Rational(Jpp[34]), Rational(Npp[34]), 0_rat)),
270  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[35]), Rational(Np[35]), 0_rat), QuantumNumbers(Rational(Jpp[35]), Rational(Npp[35]), 0_rat)),
271  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[36]), Rational(Np[36]), 0_rat), QuantumNumbers(Rational(Jpp[36]), Rational(Npp[36]), 0_rat)),
272  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[37]), Rational(Np[37]), 0_rat), QuantumNumbers(Rational(Jpp[37]), Rational(Npp[37]), 0_rat)),
273  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[38]), Rational(Np[38]), 0_rat), QuantumNumbers(Rational(Jpp[38]), Rational(Npp[38]), 0_rat)),
274  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[39]), Rational(Np[39]), 0_rat), QuantumNumbers(Rational(Jpp[39]), Rational(Npp[39]), 0_rat)),
275  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[40]), Rational(Np[40]), 0_rat), QuantumNumbers(Rational(Jpp[40]), Rational(Npp[40]), 0_rat)),
276  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[41]), Rational(Np[41]), 0_rat), QuantumNumbers(Rational(Jpp[41]), Rational(Npp[41]), 0_rat)),
277  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[42]), Rational(Np[42]), 0_rat), QuantumNumbers(Rational(Jpp[42]), Rational(Npp[42]), 0_rat)),
278  QuantumIdentifier(species, isot, QuantumNumbers(Rational(Jp[43]), Rational(Np[43]), 0_rat), QuantumNumbers(Rational(Jpp[43]), Rational(Npp[43]), 0_rat)),};
279 }
280 
281 
282 void Absorption::PredefinedModel::makarov2020_o2_lines_mpm(Matrix& xsec,
283  ArrayOfMatrix& dxsec,
284  const Vector& f,
285  const Vector& p,
286  const Vector& t,
287  const Vector& water_vmr,
288  const ArrayOfRetrievalQuantity& jacs,
289  const ArrayOfIndex& jacs_pos)
290 {
291  using Constant::pi;
292  using Constant::sqrt_pi;
293  using Constant::inv_sqrt_pi;
294  using Constant::pow2;
295  using Constant::pow3;
296 
297  // Central frequency [Hz]
298  constexpr std::array<Numeric, nlines_mpm2020> f0 = {
299  1.18750334E+11, 5.6264774E+10, 6.2486253E+10, 5.8446588E+10, 6.0306056E+10,
300  5.9590983E+10, 5.9164204E+10, 6.0434778E+10, 5.8323877E+10, 6.1150562E+10,
301  5.7612486E+10, 6.1800158E+10, 5.6968211E+10, 6.2411220E+10, 5.6363399E+10,
302  6.2997984E+10, 5.5783815E+10, 6.3568526E+10, 5.5221384E+10, 6.4127775E+10,
303  5.4671180E+10, 6.4678910E+10, 5.4130025E+10, 6.5224078E+10, 5.3595775E+10,
304  6.5764779E+10, 5.3066934E+10, 6.6302096E+10, 5.2542418E+10, 6.6836834E+10,
305  5.2021429E+10, 6.7369601E+10, 5.1503360E+10, 6.7900868E+10, 5.0987745E+10,
306  6.8431006E+10, 5.0474214E+10, 6.8960312E+10, 3.68498246E+11, 4.24763020E+11,
307  4.87249273E+11, 7.15392902E+11, 7.73839490E+11, 8.34145546E+11, };
308 
309  // Intensity [1 / Pa] (rounded to 10 digits because at most 9 digits exist in f0)
310  constexpr std::array<Numeric, nlines_mpm2020> intens = {
311  1.591521878E-21, 1.941172240E-21, 4.834543970E-21, 4.959264029E-21, 7.010386457E-21,
312  7.051673348E-21, 8.085012578E-21, 8.108262250E-21, 8.145673278E-21, 8.149757320E-21,
313  7.396406085E-21, 7.401923754E-21, 6.162286575E-21, 6.168475265E-21, 4.749226167E-21,
314  4.754435107E-21, 3.405982896E-21, 3.408455562E-21, 2.282498656E-21, 2.283934341E-21,
315  1.432692459E-21, 1.433513473E-21, 8.439995690E-22, 8.443521837E-22, 4.672706507E-22,
316  4.676049313E-22, 2.435008301E-22, 2.437304596E-22, 1.195038747E-22, 1.196873412E-22,
317  5.532759045E-23, 5.537261239E-23, 2.416832398E-23, 2.418989865E-23, 9.969285671E-24,
318  9.977543709E-24, 3.882541154E-24, 3.888101811E-24, 3.676253816E-23, 3.017524005E-22,
319  9.792882227E-23, 2.756166168E-23, 1.486462215E-22, 4.411918954E-23, };
320 
321  // Temperature intensity modifier
322  constexpr std::array<Numeric, nlines_mpm2020> a2 = {
323  0.01, 0.014, 0.083, 0.083, 0.207, 0.207, 0.387, 0.386,
324  0.621, 0.621, 0.910, 0.910, 1.255, 1.255, 1.654, 1.654,
325  2.109, 2.108, 2.618, 2.617, 3.182, 3.181, 3.800, 3.800,
326  4.474, 4.473, 5.201, 5.200, 5.983, 5.982, 6.819, 6.818,
327  7.709, 7.708, 8.653, 8.652, 9.651, 9.650, 0.048, 0.044,
328  0.049, 0.145, 0.141, 0.145};
329 
330  // Line shape model in SI units
331  constexpr auto lsm = init_mpm2020_lsm();
332 
333  // Reference temperature [K]
334  constexpr Numeric t0 = 300;
335 
336  // FIXME: Can be constexpr when SpeciesTag is constexpr
337  auto species = SpeciesTag("O2-66");
338  const std::array<QuantumIdentifier, nlines_mpm2020> qids =
339  init_mpm2020_qids(species.Species(), species.Isotopologue());
340 
341  // Model setting
342  const bool do_temp_deriv = do_temperature_jacobian(jacs);
343 
344  // Per pressure level
345  #pragma omp parallel for if (not arts_omp_in_parallel() and p.nelem() >= arts_omp_get_max_threads()) schedule(guided)
346  for (Index ip=0; ip<p.nelem(); ip++) {
347  const Numeric theta = t0 / t[ip];
348  const Numeric theta_m1 = theta - 1;
349  const Numeric theta_3 = pow3(theta);
350  const Numeric GD_div_F0 = Linefunctions::DopplerConstant(t[ip], species.SpeciesMass());
351 
352  for (Index i=0; i<nlines_mpm2020; i++) {
353  const Numeric invGD = 1 / (GD_div_F0 * f0[i]);
354  const Numeric fac = sqrt_pi * invGD;
355  const Numeric ST = theta_3 * p[ip] * intens[i] * std::exp(-a2[i] * theta_m1);
356  const Numeric G0 = (1 + 0.1*water_vmr[ip]) * p[ip] * lsm[i].compute(t[ip], t0, LineShape::Variable::G0);
357  const Numeric Y = p[ip] * lsm[i].compute(t[ip], t0, LineShape::Variable::Y);
358  const Numeric G = pow2( p[ip]) * lsm[i].compute(t[ip], t0, LineShape::Variable::G);
359  const Numeric DV = pow2(p[ip]) * lsm[i].compute(t[ip], t0, LineShape::Variable::DV);
360 
361  const Numeric dinvGD_dT = do_temp_deriv ? - invGD * Linefunctions::dDopplerConstant_dT(t[ip], GD_div_F0) : 0;
362  const Numeric dST_dT = do_temp_deriv ? (a2[i]*t0 - 3*t[ip]) / pow2(t[ip]) * ST : 0;
363  const Numeric dG0_dT = do_temp_deriv ? (1 + 0.1*water_vmr[ip]) * p[ip] * lsm[i].compute_dT(t[ip], t0, LineShape::Variable::G0) : 0;
364  const Numeric dY_dT = do_temp_deriv ? p[ip] * lsm[i].compute_dT(t[ip], t0, LineShape::Variable::Y) : 0;
365  const Numeric dG_dT = do_temp_deriv ? pow2(p[ip]) * lsm[i].compute_dT(t[ip], t0, LineShape::Variable::G) : 0;
366  const Numeric dDV_dT = do_temp_deriv ? pow2(p[ip]) * lsm[i].compute_dT(t[ip], t0, LineShape::Variable::DV) : 0;
367 
368  for (Index j=0; j<f.nelem(); j++) {
369  const Complex z = Complex(f0[i] + DV - f[j], G0) * invGD;
370  const Complex Fv = fac * Faddeeva::w(z);
371  const Complex Flm = 1 / Complex(G0, f[j] + f0[i] + DV);
372 
373  const Complex abs = std::real(
374  /* around line center */
375  Complex(1 + G, Y) * Fv +
376  /* mirrored line far from line center */
377  Complex(1 + G, -Y) * Flm);
378 
379  xsec(j, ip) += ST * pow2(f[j]) * abs.real();
380 
381  if (jacs_pos.nelem()) {
382  const Complex dw = 2 * (Complex(0, fac * inv_sqrt_pi) - z * Fv);
383  const Complex dm = - pi * pow2(Flm);
384 
385  for (Index iq=0; iq<jacs_pos.nelem(); iq++) {
386  const auto& deriv = jacs[jacs_pos[iq]];
387 
388  if (deriv == JacPropMatType::Temperature) {
389  const Complex dFv = dw * (invGD * Complex(dDV_dT, dG0_dT) - dinvGD_dT) + Fv * dinvGD_dT;
390  const Complex dFlm = dm * Complex(dG0_dT, dDV_dT);
391  dxsec[iq](j, ip) += pow2(f[j]) * (ST * std::real(
392  /* around line center */
393  Complex(1 + G, Y) * dFv + Complex(dG_dT, dY_dT) * Fv +
394  /* mirrored line far from line center */
395  Complex(1 + G, -Y) * dFlm + Complex(G, -dY_dT) * Flm) + abs.real() * dST_dT);
396  } else if (is_frequency_parameter(deriv)) {
397  const Complex dFv = - dw * invGD;
398  const Complex dFlm = Complex(0, 1) * dm;
399  dxsec[iq](j, ip) += ST * (pow2(f[j]) * std::real(
400  /* around line center */
401  Complex(1 + G, Y) * dFv +
402  /* mirrored line far from line center */
403  Complex(1 + G, -Y) * dFlm) + 2 * abs.real() * f[j]);
404  } else if (deriv.QuantumIdentity().In(qids[i])) {
405  if (deriv == JacPropMatType::LineShapeG0X0) {
406  dxsec[iq](j, ip) += ST * pow2(f[j]) * std::real(
407  /* around line center */
408  Complex(1 + G, Y) * Complex(0, 1) * dw * invGD +
409  /* mirrored line far from line center */
410  Complex(1 + G, -Y) * dm) *
411  lsm[i].compute_dX0(t[ip], t0, LineShape::Variable::G0);
412  } else if (deriv == JacPropMatType::LineShapeG0X1) {
413  dxsec[iq](j, ip) += ST * pow2(f[j]) * std::real(
414  /* around line center */
415  Complex(1 + G, Y) * Complex(0, 1) * dw * invGD +
416  /* mirrored line far from line center */
417  Complex(1 + G, -Y) * dm) *
418  lsm[i].compute_dX1(t[ip], t0, LineShape::Variable::DV);
419  } else if (deriv == JacPropMatType::LineShapeDVX0) {
420  const Complex dFv = dw * invGD;
421  const Complex dFlm = Complex(0, 1) * dm;
422  dxsec[iq](j, ip) += ST * pow2(f[j]) * std::real(
423  /* around line center */
424  Complex(1 + G, Y) * dFv +
425  /* mirrored line far from line center */
426  Complex(1 + G, -Y) * dFlm) *
427  lsm[i].compute_dX0(t[ip], t0, LineShape::Variable::DV);;
428  } else if (deriv == JacPropMatType::LineShapeDVX1) {
429  const Complex dFv = dw * invGD;
430  const Complex dFlm = Complex(0, 1) * dm;
431  dxsec[iq](j, ip) += ST * pow2(f[j]) * std::real(
432  /* around line center */
433  Complex(1 + G, Y) * dFv +
434  /* mirrored line far from line center */
435  Complex(1 + G, -Y) * dFlm) *
436  lsm[i].compute_dX1(t[ip], t0, LineShape::Variable::DV);;
437  } else if (deriv == JacPropMatType::LineShapeDVX2) {
438  const Complex dFv = dw * invGD;
439  const Complex dFlm = Complex(0, 1) * dm;
440  dxsec[iq](j, ip) += ST * pow2(f[j]) * std::real(
441  /* around line center */
442  Complex(1 + G, Y) * dFv +
443  /* mirrored line far from line center */
444  Complex(1 + G, -Y) * dFlm) *
445  lsm[i].compute_dX2(t[ip], t0, LineShape::Variable::DV);
446  } else if (deriv == JacPropMatType::LineShapeGX0) {
447  dxsec[iq](j, ip) += ST * pow2(f[j]) * std::real(
448  /* around line center */
449  Fv +
450  /* mirrored line far from line center */
451  Flm) *
452  lsm[i].compute_dX0(t[ip], t0, LineShape::Variable::G);
453  } else if (deriv == JacPropMatType::LineShapeYX0) {
454  dxsec[iq](j, ip) += ST * pow2(f[j]) * std::real(
455  /* around line center */
456  Fv +
457  /* mirrored line far from line center */
458  Flm) *
459  lsm[i].compute_dX0(t[ip], t0, LineShape::Variable::Y);
460  } else if (deriv == JacPropMatType::LineShapeGX1) {
461  dxsec[iq](j, ip) += ST * pow2(f[j]) * std::real(
462  /* around line center */
463  Fv -
464  /* mirrored line far from line center */
465  Flm) *
466  lsm[i].compute_dX1(t[ip], t0, LineShape::Variable::G);
467  } else if (deriv == JacPropMatType::LineShapeYX1) {
468  dxsec[iq](j, ip) += ST * pow2(f[j]) * std::real(
469  /* around line center */
470  Fv +
471  /* mirrored line far from line center */
472  Flm) *
473  lsm[i].compute_dX1(t[ip], t0, LineShape::Variable::Y);
474  } else if (deriv == JacPropMatType::LineShapeGX2) {
475  dxsec[iq](j, ip) += ST * pow2(f[j]) * std::real(
476  /* around line center */
477  Fv -
478  /* mirrored line far from line center */
479  Flm) *
480  lsm[i].compute_dX2(t[ip], t0, LineShape::Variable::G);
481  } else if (deriv == JacPropMatType::LineShapeYX2) {
482  dxsec[iq](j, ip) += ST * pow2(f[j]) * std::real(
483  /* around line center */
484  Fv -
485  /* mirrored line far from line center */
486  Flm) *
487  lsm[i].compute_dX2(t[ip], t0, LineShape::Variable::Y);
488  } else if (deriv == JacPropMatType::LineCenter) {
489  const Complex dFv = Fv / f0[i] - dw * invGD + dw * z / f0[i];
490  const Complex dFlm = Complex(0, 1) * dm;
491  dxsec[iq](j, ip) += ST * pow2(f[j]) * std::real(
492  /* around line center */
493  Complex(1 + G, Y) * dFv +
494  /* mirrored line far from line center */
495  Complex(1 + G, -Y) * dFlm);
496  } else if (deriv == JacPropMatType::LineStrength) {
497  dxsec[iq](j, ip) += theta_3 * p[ip] * std::exp(-a2[i] * theta_m1) * pow2(f[j]) * abs.real();
498  } else if (deriv == JacPropMatType::LineSpecialParameter1) {
499  dxsec[iq](j, ip) += -theta_m1 * ST * pow2(f[j]) * abs.real();
500  }
501  }
502  }
503  }
504  }
505  }
506  }
507 }
508 
509 
510 void normalize_relaxation_matrix(Eigen::Ref<Eigen::MatrixXcd> W,
511  const Eigen::Ref<Eigen::ArrayXd> rho,
512  const Eigen::Ref<Eigen::ArrayXd> d)
513 {
514  // FIXME:
515  // rho[i] / rho[j] or rho[j] / rho[i] in the next for-loop. [See {Prop. 1}]
516  // The current order is the same as for the second-to-next for-loop. [See {Prop. 2}]
517  // If the first loop changes, the renormalization in the second loop
518  // must change...
519 
520  // Population density balanced parameters
521  for (Index i=0; i<necs2020; i++) {
522  for (Index j=0; j<i; j++) {
523  if (i not_eq j) {
524  W(i, j) = W(j, i) * (rho[i] / rho[j]); // {Prop. 1} --- division-order should change?
525  }
526  }
527  }
528 
529  // Balance by the reduced dipole
530  for (Index i=0; i<necs2020; i++) {
531  auto norm = - std::imag(d[i] * W(i, i)) / std::imag(W.row(i) * d.abs().matrix() - std::abs(d[i]) * W(i, i)); // {Prop. 2} --- dot-product axis should to column?
532  for (Index j=0; j<necs2020; j++) {
533  if (i not_eq j) {
534  W(i, j) *= norm; // {Prop. 2} --- If yes, then change W(i, j) to W(j, i)
535  }
536  }
537  }
538 }
539 
540 
541 void Absorption::PredefinedModel::makarov2020_o2_lines_ecs(ComplexVector& I, const Vector& f, Numeric P, Numeric T, Numeric water_vmr)
542 {
543  using Constant::h;
544  using Constant::k;
545  using Constant::pi;
546  using Constant::pow3;
547 
548  constexpr auto qids = init_mpm2020_qids(0, 0);
549  constexpr auto lsm = init_mpm2020_lsm();
550  constexpr auto T0 = 300;
551 
552  // Central frequency [Hz]
553  constexpr std::array<Numeric, nlines_mpm2020> f0 = {
554  1.18750334E+11, 5.6264774E+10, 6.2486253E+10, 5.8446588E+10, 6.0306056E+10,
555  5.9590983E+10, 5.9164204E+10, 6.0434778E+10, 5.8323877E+10, 6.1150562E+10,
556  5.7612486E+10, 6.1800158E+10, 5.6968211E+10, 6.2411220E+10, 5.6363399E+10,
557  6.2997984E+10, 5.5783815E+10, 6.3568526E+10, 5.5221384E+10, 6.4127775E+10,
558  5.4671180E+10, 6.4678910E+10, 5.4130025E+10, 6.5224078E+10, 5.3595775E+10,
559  6.5764779E+10, 5.3066934E+10, 6.6302096E+10, 5.2542418E+10, 6.6836834E+10,
560  5.2021429E+10, 6.7369601E+10, 5.1503360E+10, 6.7900868E+10, 5.0987745E+10,
561  6.8431006E+10, 5.0474214E+10, 6.8960312E+10, 3.68498246E+11, 4.24763020E+11,
562  4.87249273E+11, 7.15392902E+11, 7.73839490E+11, 8.34145546E+11, };
563 
564  // Intensity [1 / Pa] (rounded to 10 digits because at most 9 digits exist in f0)
565  constexpr std::array<Numeric, nlines_mpm2020> intens = {
566  1.591521878E-21, 1.941172240E-21, 4.834543970E-21, 4.959264029E-21, 7.010386457E-21,
567  7.051673348E-21, 8.085012578E-21, 8.108262250E-21, 8.145673278E-21, 8.149757320E-21,
568  7.396406085E-21, 7.401923754E-21, 6.162286575E-21, 6.168475265E-21, 4.749226167E-21,
569  4.754435107E-21, 3.405982896E-21, 3.408455562E-21, 2.282498656E-21, 2.283934341E-21,
570  1.432692459E-21, 1.433513473E-21, 8.439995690E-22, 8.443521837E-22, 4.672706507E-22,
571  4.676049313E-22, 2.435008301E-22, 2.437304596E-22, 1.195038747E-22, 1.196873412E-22,
572  5.532759045E-23, 5.537261239E-23, 2.416832398E-23, 2.418989865E-23, 9.969285671E-24,
573  9.977543709E-24, 3.882541154E-24, 3.888101811E-24, 3.676253816E-23, 3.017524005E-22,
574  9.792882227E-23, 2.756166168E-23, 1.486462215E-22, 4.411918954E-23, };
575 
576  // Temperature intensity modifier
577  constexpr std::array<Numeric, nlines_mpm2020> a2 = {
578  0.01, 0.014, 0.083, 0.083, 0.207, 0.207, 0.387, 0.386,
579  0.621, 0.621, 0.910, 0.910, 1.255, 1.255, 1.654, 1.654,
580  2.109, 2.108, 2.618, 2.617, 3.182, 3.181, 3.800, 3.800,
581  4.474, 4.473, 5.201, 5.200, 5.983, 5.982, 6.819, 6.818,
582  7.709, 7.708, 8.653, 8.652, 9.651, 9.650, 0.048, 0.044,
583  0.049, 0.145, 0.141, 0.145};
584 
585  // Computed arrays values
586  Eigen::Array<Numeric, necs2020, 1> d;
587  Eigen::Array<Numeric, necs2020, 1> rho;
588  Eigen::Matrix<Complex, necs2020, necs2020> W;
589 
590  // Diagonal and upper triangular values
591  for (Index i=0; i<necs2020; i++) {
592  auto Ji = qids[i].UpperQuantumNumber(QuantumNumberType::J);
593  auto Jf = qids[i].LowerQuantumNumber(QuantumNumberType::J);
594  auto Ni = qids[i].UpperQuantumNumber(QuantumNumberType::N);
595  auto Nf = qids[i].LowerQuantumNumber(QuantumNumberType::N);
596  for (Index j=i; j<necs2020; j++) {
597  auto Ji_p = qids[j].UpperQuantumNumber(QuantumNumberType::J);
598  auto Jf_p = qids[j].LowerQuantumNumber(QuantumNumberType::J);
599  auto Ni_p = qids[j].UpperQuantumNumber(QuantumNumberType::N);
600  auto Nf_p = qids[j].LowerQuantumNumber(QuantumNumberType::N);
601  if (i not_eq j) {
602  W(i, j) = 1i * o2_ecs_wigner_symbol_tran(Ji, Jf, Ni, Nf, 1_rat, 1_rat, Ji_p, Jf_p, Ni_p, Nf_p, 1_rat, T);
603  } else {
604  W(i, j) = 1i * (1 + 0.1*water_vmr) * P * lsm[i].compute(T, T0, LineShape::Variable::G0);
605  }
606  }
607  }
608 
609  std::transform(qids.cbegin(), qids.cbegin()+necs2020, d.data(), [](auto& qns){return
610  o2_makarov2013_reduced_dipole (qns.UpperQuantumNumber(QuantumNumberType::J),
611  qns.LowerQuantumNumber(QuantumNumberType::J),
612  qns.UpperQuantumNumber(QuantumNumberType::N));});
613 // std::cout << "Ji = np.array([";
614 // for (Index i=0; i<necs2020; i++) {
615 // std::cout << qids[i].UpperQuantumNumber(QuantumNumberType::J) << ",\n";
616 // }
617 // std::cout<<"])\n";
618 // std::cout << "f0 = np.array([";
619 // for (Index i=0; i<necs2020; i++) {
620 // std::cout << f0[i] << ",\n";
621 // }
622 // std::cout<<"])\n";
623 // std::cout << "Jf = np.array([";
624 // for (Index i=0; i<necs2020; i++) {
625 // std::cout << qids[i].LowerQuantumNumber(QuantumNumberType::J) << ",\n";
626 // }
627 // std::cout<<"])\n";
628 // std::cout << "N = np.array([";
629 // for (Index i=0; i<necs2020; i++) {
630 // std::cout << qids[i].LowerQuantumNumber(QuantumNumberType::N) << ",\n";
631 // }
632 // std::cout<<"])\n";
633 //
634 // std::cout << "d = np.array([";
635 // for (Index i=0; i<necs2020; i++) {
636 // std::cout << d[i] << ",\n";
637 // }
638 // std::cout<<"])\n";
639 
640  std::transform(qids.cbegin(), qids.cbegin()+necs2020, rho.data(), [T](auto& qns){return
641  boltzman_factor(T, o2_ecs_erot_jn_same(qns.LowerQuantumNumber(QuantumNumberType::J)));});
642 // std::cout << "rho = np.array([";
643 // for (Index i=0; i<necs2020; i++) {
644 // std::cout << rho[i] << ",\n";
645 // }
646 // std::cout<<"])\n";
647 //
648 // std::cout << "Wprenorm = np.array([";
649 // for (Index i=0; i<necs2020; i++) {
650 // std::cout << "[";
651 // for (Index j=0; j<necs2020; j++) {
652 // std::cout << W(i, j).imag() << ",\n";
653 // }
654 // std::cout<<"],\n";
655 // }
656 // std::cout<<"])\n";
657  normalize_relaxation_matrix(W, rho, d);
658 // std::cout << "Wnorm = np.array([";
659 // for (Index i=0; i<necs2020; i++) {
660 // std::cout << "[";
661 // for (Index j=0; j<necs2020; j++) {
662 // std::cout << W(i, j).imag() << ",\n";
663 // }
664 // std::cout<<"],\n";
665 // }
666 // std::cout<<"])\n";
667 
668  // Rescale dipole to be proportional to line strength
669  const Numeric theta = T0 / T;
670  const Numeric theta_m1 = theta - 1;
671  const Numeric theta_3 = pow3(theta);
672  for (Index i=0; i<necs2020; i++) {
673  const Numeric ST = pi * theta_3 * P * intens[i] * std::exp(-a2[i] * theta_m1);
674  const Numeric sgn = std::abs(d[i]) / d[i];
675  d[i] = sgn * f0[i] * std::sqrt(ST / rho[i]);
676  }
677  //std::cout << "ST = np.array([";
678  //for (Index i=0; i<necs2020; i++) {
679  // std::cout << d[i] << ",\n";
680  //}
681  //std::cout<<"])\n";
682 
683  /*if (not full)*/
684  for (Index i=0; i<necs2020; i++) {
685  W(i, i) += f0[i];
686  }
687 
688  const Eigen::ComplexEigenSolver<Eigen::Matrix<Complex, necs2020, necs2020>> eV(W, true);
689  const auto& D = eV.eigenvalues();
690  const auto& V = eV.eigenvectors();
691  const auto& Vinv = W = eV.eigenvectors().inverse(); // Reuse W memory but with different &name
692 
693  Eigen::Array<Complex, necs2020, 1> B; B *= 0;
694  for (Index m=0; m<necs2020; m++) {
695  for (Index i=0; i<necs2020; i++) {
696  for (Index j=0; j<necs2020; j++) {
697  B[m] += rho[i] * d[i] * d[j] * V(j, m) * Vinv(m, i);
698  }
699  }
700  }
701 
702 // // Lorentz profile!
703 // const Numeric div = Constant::inv_pi / rho.cwiseProduct(d.cwiseAbs2()).sum();
704 // for (Index i=0; i<f.nelem(); i++) {
705 // I[i] = div * (B.array() / (f[i] - D.array())).sum();
706 
707  auto GD0 = Linefunctions::DopplerConstant(T, SpeciesTag("O2-66").SpeciesMass());
708  for (Index i=0; i<necs2020; i++) {
709  for (Index iv=0; iv<f.nelem(); iv++) {
710  I[iv] += (Constant::inv_sqrt_pi / (GD0 * D[i].real())) * Faddeeva::w((f[iv] - std::conj(D[i])) / (GD0 * D[i].real())) * B[i];
711  }
712  }
713 }
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
o2_makarov2013_reduced_dipole
Numeric o2_makarov2013_reduced_dipole(const Rational &Jup, const Rational &Jlo, const Rational &N)
Returns the reduced dipole moment following Makarov etal 2013.
Definition: wigner_functions.cc:267
init_mpm2020_slsm
constexpr LineShape::SingleSpeciesModel init_mpm2020_slsm(Numeric g00, Numeric y0, Numeric y1, Numeric g0, Numeric g1, Numeric dv0, Numeric dv1, Numeric x)
Definition: predefined_absorption_models.cc:37
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
Vector
The Vector class.
Definition: matpackI.h:860
abs
#define abs(x)
Definition: legacy_continua.cc:20626
fac
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:63
normalize_relaxation_matrix
void normalize_relaxation_matrix(Eigen::Ref< Eigen::MatrixXcd > W, const Eigen::Ref< Eigen::ArrayXd > rho, const Eigen::Ref< Eigen::ArrayXd > d)
Definition: predefined_absorption_models.cc:510
lin_alg.h
Linear algebra functions.
boltzman_factor
Numeric boltzman_factor(Numeric T, Numeric E0)
Computes exp(- E0/kT)
Definition: linescaling.cc:176
conj
constexpr Complex conj(Complex c)
the conjugate of c
Definition: complex.h:65
w
cmplx FADDEEVA() w(cmplx z, double relerr)
Definition: Faddeeva.cc:680
Constant::pow2
constexpr T pow2(T x)
power of two
Definition: constants.h:64
wigner_functions.h
Wigner symbol interactions.
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
is_frequency_parameter
bool is_frequency_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a frequency parameter in propagation matrix calculations.
Definition: jacobian.cc:1120
init_mpm2020_lsm
constexpr std::array< LineShape::SingleSpeciesModel, nlines_mpm2020 > init_mpm2020_lsm()
Definition: predefined_absorption_models.cc:59
init_mpm2020_qids
constexpr std::array< QuantumIdentifier, nlines_mpm2020 > init_mpm2020_qids(const Index &species, const Index &isot)
Definition: predefined_absorption_models.cc:201
Absorption::PredefinedModel::makarov2020_o2_lines_mpm
void makarov2020_o2_lines_mpm(Matrix &xsec, ArrayOfMatrix &dxsec, const Vector &f, const Vector &p, const Vector &t, const Vector &water_vmr, const ArrayOfRetrievalQuantity &jacs, const ArrayOfIndex &jacs_pos)
Adds Makarov MPM2020 O2 absorption lines to the absorption matrix.
Definition: predefined_absorption_models.cc:282
Complex
std::complex< Numeric > Complex
Definition: complex.h:33
Rational
Implements rational numbers to work with other ARTS types.
Definition: rational.h:54
SpeciesTag
A tag group can consist of the sum of several of these.
Definition: abs_species_tags.h:44
QuantumIdentifier
Class to identify and match lines by their quantum numbers.
Definition: quantum.h:390
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
dw
constexpr Complex dw(Complex z, Complex w) noexcept
The Faddeeva function partial derivative.
Definition: linefunctions.cc:45
Matrix
The Matrix class.
Definition: matpackI.h:1193
o2_ecs_erot_jn_same
Numeric o2_ecs_erot_jn_same(Rational J)
Energy of the J=N line at J.
Definition: wigner_functions.cc:159
init_mpm2020_qid
constexpr QuantumIdentifier init_mpm2020_qid(Index species, Index isot, Rational Jup, Rational Jlow, Rational Nup, Rational Nlow)
Definition: predefined_absorption_models.cc:195
Constant::pow3
constexpr T pow3(T x)
power of three
Definition: constants.h:70
Linefunctions::dDopplerConstant_dT
Numeric dDopplerConstant_dT(const Numeric &T, const Numeric &dc)
Returns the temperature derivative of the frequency-independent part of the Doppler broadening...
Definition: linefunctions.cc:811
do_temperature_jacobian
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition: jacobian.cc:1279
ComplexVector
The ComplexVector class.
Definition: complex.h:573
lm_hitran_2017::compute
Vector compute(const Numeric p, const Numeric t, const Numeric xco2, const Numeric xh2o, const ConstVectorView invcm_grid, const Numeric stotmax, const calctype type)
Definition: linemixing_hitran.cc:1882
QuantumNumbers
Container class for Quantum Numbers.
Definition: quantum.h:222
Array
This can be used to make arrays out of anything.
Definition: array.h:40
necs2020
constexpr auto necs2020
Definition: predefined_absorption_models.cc:33
linescaling.h
Constains various line scaling functions.
o2_ecs_wigner_symbol_tran
Numeric o2_ecs_wigner_symbol_tran(const Rational &Ji, const Rational &Jf, const Rational &Ni, const Rational &Nf, const Rational &Si, const Rational &Sf, const Rational &Ji_p, const Rational &Jf_p, const Rational &Ni_p, const Rational &Nf_p, const Rational &n, const Numeric &T)
Returns the wigner symbol used in Tran etal 2006.
Definition: wigner_functions.cc:218
nlines_mpm2020
constexpr auto nlines_mpm2020
Definition: predefined_absorption_models.cc:34
Linefunctions::DopplerConstant
Numeric DopplerConstant(Numeric T, Numeric mass)
Returns the frequency-independent part of the Doppler broadening.
Definition: linefunctions.cc:807
a2
#define a2
Definition: complex.h:58
Absorption::PredefinedModel::makarov2020_o2_lines_ecs
void makarov2020_o2_lines_ecs(ComplexVector &I, const Vector &f, Numeric P, Numeric T, Numeric water_vmr)
Definition: predefined_absorption_models.cc:541
transform
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1476
LineShape::SingleSpeciesModel
Compute the line shape parameters for a single broadening species.
Definition: linefunctiondata.h:204
predefined_absorption_models.h
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620