ARTS 2.5.10 (git: 2f1c442c)
lineshape.cc
Go to the documentation of this file.
1#include "debug.h"
2#include "enums.h"
3#include "lineshapemodel.h"
4#include "partfun.h"
5
6#include <algorithm>
7#include <cmath>
8
9#include "lineshape.h"
10#include "physics_funcs.h"
11#include "species.h"
12
13#include <Faddeeva/Faddeeva.hh>
14
15using Constant::inv_pi;
16using Constant::inv_sqrt_pi;
17using Constant::pi;
18using Constant::sqrt_ln_2;
19using Constant::sqrt_pi;
20using Math::pow2;
21using Math::pow3;
22using Math::pow4;
23
24constexpr Numeric ln_16 =
25 2.772588722239781237668928485832706272302000537441021016482720037973574487879;
26
27namespace LineShape {
28Complex Doppler::operator()(Numeric f) noexcept {
29 x = (f - mF0) * invGD;
30 F = invGD * inv_sqrt_pi * std::exp(-pow2(x));
31 return F;
32}
33
34Complex Voigt::operator()(Numeric f) noexcept {
35 real_val(z) = invGD * (f - mF0);
36 F = inv_sqrt_pi * invGD * Faddeeva::w(z);
37 dF = 2 * invGD * (Complex(0, inv_pi * invGD) - z * F);
38 return F;
39}
40
41SpeedDependentVoigt::SpeedDependentVoigt(Numeric F0_noshift, const Output &ls,
42 Numeric GD_div_F0, Numeric dZ) noexcept
43 : mF0(F0_noshift + dZ + ls.D0 - 1.5 * ls.D2),
44 invGD(sqrt_ln_2 / nonstd::abs(GD_div_F0 * mF0)),
45 invc2(1.0 / Complex(ls.G2, ls.D2)), dx(Complex(ls.G0 - 1.5 * ls.G2, mF0)),
46 x(dx * invc2), sqrty(invc2 / (2 * invGD)),
47 calcs(init(Complex(ls.G2, ls.D2))) {
48 calc();
49}
50
51Complex SpeedDependentVoigt::dFdf() const noexcept {
52 switch (calcs) {
53 case CalcType::Full:
54 return invGD * invc2 * (dw1 - dw2) / (2 * sqrt_pi * sq);
55 case CalcType::Voigt:
56 return dw1 * pow2(invGD) * inv_sqrt_pi;
57 case CalcType::LowXandHighY:
58 return dw1 * pow2(invGD) * inv_sqrt_pi -
59 dw2 * invGD * invc2 / (2 * sqrt_pi * sq);
60 case CalcType::LowYandLowX:
61 return pow2(invc2) * (-dw1 * sq + 1i * w1) / (sqrt_pi * sq);
62 case CalcType::LowYandHighX:
63 return 1i * pow2(invc2) * (x - 3) / (pi * pow3(x));
64 }
65 return {};
66}
67
68Complex SpeedDependentVoigt::dFdF0() const noexcept {
69 switch (calcs) {
70 case CalcType::Full:
71 return (4 * pow2(invGD) * (-w1 + w2) * sq +
72 1i * invc2 *
73 (dw1 * (Complex(0, 2 * mF0 * pow2(invGD)) - 2 * invGD * sq +
74 invc2) -
75 dw2 * (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
76 invc2))) /
77 (4 * sqrt_pi * invGD * mF0 * sq);
78 case CalcType::Voigt:
79 return -invGD * (Complex(0, invGD) * dw1 * (dx - Complex(0, mF0)) + w1) /
80 (sqrt_pi * mF0);
81 case CalcType::LowXandHighY:
82 return (4 * pow2(invGD) * (-w1 + w2) * sq -
83 1i * (4 * dw1 * pow3(invGD) * (dx - Complex(0, mF0)) * sq +
84 dw2 * invc2 *
85 (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
86 invc2))) /
87 (4 * sqrt_pi * invGD * mF0 * sq);
88 case CalcType::LowYandLowX:
89 return pow2(invc2) * (dw1 * sq - 1i * w1) / (sqrt_pi * sq);
90 case CalcType::LowYandHighX:
91 return 1i * pow2(invc2) * (3 - x) / (pi * pow3(x));
92 }
93 return {};
94}
95
96Complex SpeedDependentVoigt::dFdD0(Numeric dD0dD0) const noexcept {
97 switch (calcs) {
98 case CalcType::Full:
99 return -dD0dD0 *
100 (4 * pow2(invGD) * (w1 - w2) * sq +
101 1i * invc2 *
102 (-dw1 * (Complex(0, 2 * mF0 * pow2(invGD)) - 2 * invGD * sq +
103 invc2) +
104 dw2 * (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
105 invc2))) /
106 (4 * sqrt_pi * invGD * mF0 * sq);
107 case CalcType::Voigt:
108 return -dD0dD0 * invGD *
109 (Complex(0, invGD) * dw1 * (dx - Complex(0, mF0)) + w1) /
110 (sqrt_pi * mF0);
111 case CalcType::LowXandHighY:
112 return -dD0dD0 *
113 (4 * pow2(invGD) * (w1 - w2) * sq +
114 1i * (4 * dw1 * pow3(invGD) * (dx - Complex(0, mF0)) * sq +
115 dw2 * invc2 *
116 (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
117 invc2))) /
118 (4 * sqrt_pi * invGD * mF0 * sq);
119 case CalcType::LowYandLowX:
120 return dD0dD0 * pow2(invc2) * (dw1 * sq - 1i * w1) / (sqrt_pi * sq);
121 case CalcType::LowYandHighX:
122 return -Complex(0, dD0dD0) * pow2(invc2) * (x - 3) / (pi * pow3(x));
123 }
124 return {};
125}
126
127Complex SpeedDependentVoigt::dFdG0(Numeric dG0dG0) const noexcept {
128 switch (calcs) {
129 case CalcType::Full:
130 return Complex(0, dG0dG0) * invGD * invc2 * (dw1 - dw2) /
131 (2 * sqrt_pi * sq);
132 case CalcType::Voigt:
133 return Complex(0, dG0dG0) * dw1 * pow2(invGD) * inv_sqrt_pi;
134 case CalcType::LowXandHighY:
135 return Complex(0, dG0dG0) * invGD * (2 * dw1 * invGD * sq - dw2 * invc2) /
136 (2 * sqrt_pi * sq);
137 case CalcType::LowYandLowX:
138 return -dG0dG0 * pow2(invc2) * (1i * dw1 * sq + w1) / (sqrt_pi * sq);
139 case CalcType::LowYandHighX:
140 return -dG0dG0 * pow2(invc2) * (x - 3) / (pi * pow3(x));
141 }
142 return {};
143}
144
145Complex SpeedDependentVoigt::dFdD2(Numeric dD2dD2) const noexcept {
146 switch (calcs) {
147 case CalcType::Full:
148 return dD2dD2 *
149 (12 * pow2(invGD) * (w1 - w2) * sq +
150 1i * invc2 *
151 (dw1 * (-Complex(0, 2 * mF0 * pow2(invGD)) *
152 (2 * dx * invc2 + 3) +
153 4 * Complex(0, invGD) * invc2 * mF0 * sq +
154 6 * invGD * sq - Complex(0, 2 * mF0) * pow2(invc2) -
155 3 * invc2) +
156 dw2 *
157 (Complex(0, 2 * mF0 * pow2(invGD)) * (2 * dx * invc2 + 3) +
158 4 * Complex(0, invGD) * invc2 * mF0 * sq +
159 6 * invGD * sq + Complex(0, 2 * mF0) * pow2(invc2) +
160 3 * invc2))) /
161 (8 * sqrt_pi * invGD * mF0 * sq);
162 case CalcType::Voigt:
163 return 3 * dD2dD2 * invGD *
164 (Complex(0, invGD) * dw1 * (dx - Complex(0, mF0)) + w1) /
165 (2 * sqrt_pi * mF0);
166 case CalcType::LowXandHighY:
167 return dD2dD2 *
168 (12 * pow2(invGD) * (w1 - w2) * sq +
169 1i *
170 (12 * dw1 * pow3(invGD) * (dx - Complex(0, mF0)) * sq +
171 dw2 * invc2 *
172 (Complex(0, 2 * mF0 * pow2(invGD)) * (2 * dx * invc2 + 3) +
173 4 * Complex(0, invGD) * invc2 * mF0 * sq +
174 6 * invGD * sq + Complex(0, 2 * mF0) * pow2(invc2) +
175 3 * invc2))) /
176 (8 * sqrt_pi * invGD * mF0 * sq);
177 case CalcType::LowYandLowX:
178 return dD2dD2 * pow2(invc2) *
179 (4 * 1i * sq * (sqrt_pi * w1 * sq - 1) -
180 sqrt_pi * (dw1 * sq - 1i * w1) * (2 * dx * invc2 + 3)) /
181 (2 * pi * sq);
182 case CalcType::LowYandHighX:
183 return Complex(0, dD2dD2) * pow2(invc2) *
184 (-x * (2 * x - 3) + (x - 3) * (2 * dx * invc2 + 3)) /
185 (2 * pi * pow3(x));
186 }
187 return {};
188}
189
190Complex SpeedDependentVoigt::dFdG2(Numeric dG2dG2) const noexcept {
191 switch (calcs) {
192 case CalcType::Full:
193 return Complex(0, dG2dG2) * invc2 *
194 (dw1 * (-pow2(invGD) * (2 * dx * invc2 + 3) +
195 2 * invGD * invc2 * sq - pow2(invc2)) +
196 dw2 * (pow2(invGD) * (2 * dx * invc2 + 3) + 2 * invGD * invc2 * sq +
197 pow2(invc2))) /
198 (4 * sqrt_pi * invGD * sq);
199 case CalcType::Voigt:
200 return -3 * Complex(0, dG2dG2) * dw1 * pow2(invGD) / (2 * sqrt_pi);
201 case CalcType::LowXandHighY:
202 return Complex(0, dG2dG2) *
203 (-6 * dw1 * pow3(invGD) * sq +
204 dw2 * invc2 *
205 (pow2(invGD) * (2 * dx * invc2 + 3) + 2 * invGD * invc2 * sq +
206 pow2(invc2))) /
207 (4 * sqrt_pi * invGD * sq);
208 case CalcType::LowYandLowX:
209 return dG2dG2 * pow2(invc2) *
210 (4 * sq * (sqrt_pi * w1 * sq - 1) +
211 sqrt_pi * (2 * dx * invc2 + 3) * (1i * dw1 * sq + w1)) /
212 (2 * pi * sq);
213 case CalcType::LowYandHighX:
214 return dG2dG2 * pow2(invc2) *
215 (-x * (2 * x - 3) + (x - 3) * (2 * dx * invc2 + 3)) /
216 (2 * pi * pow3(x));
217 }
218 return {};
219}
220
221Complex SpeedDependentVoigt::dFdH(Numeric dZ) const noexcept {
222 switch (calcs) {
223 case CalcType::Full:
224 return -dZ *
225 (4 * pow2(invGD) * (w1 - w2) * sq +
226 1i * invc2 *
227 (-dw1 * (Complex(0, 2 * mF0 * pow2(invGD)) - 2 * invGD * sq +
228 invc2) +
229 dw2 * (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
230 invc2))) /
231 (4 * sqrt_pi * invGD * mF0 * sq);
232 case CalcType::Voigt:
233 return -dZ * invGD *
234 (Complex(0, invGD) * dw1 * (dx - Complex(0, mF0)) + w1) /
235 (sqrt_pi * mF0);
236 case CalcType::LowXandHighY:
237 return -dZ *
238 (4 * pow2(invGD) * (w1 - w2) * sq +
239 1i * (4 * dw1 * pow3(invGD) * (dx - Complex(0, mF0)) * sq +
240 dw2 * invc2 *
241 (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
242 invc2))) /
243 (4 * sqrt_pi * invGD * mF0 * sq);
244 case CalcType::LowYandLowX:
245 return dZ * pow2(invc2) * (dw1 * sq - 1i * w1) / (sqrt_pi * sq);
246 case CalcType::LowYandHighX:
247 return -Complex(0, dZ) * pow2(invc2) * (x - 3) / (pi * pow3(x));
248 }
249 return {};
250}
251
252Complex SpeedDependentVoigt::dFdVMR(const Output &d) const noexcept {
253 switch (calcs) {
254 case CalcType::Full:
255 return (-4 * pow2(invGD) * (2 * d.D0 - 3 * d.D2) * (w1 - w2) * sq +
256 1i * invc2 *
257 (dw1 * (-2 * pow2(invGD) * mF0 *
258 (Complex(3 * d.G2 - 2 * d.G0, 3 * d.D2 - 2 * d.D0) +
259 2 * dx * invc2 * Complex(d.G2, d.D2)) +
260 4 * invGD * invc2 * mF0 * sq * Complex(d.G2, d.D2) -
261 2 * invGD * (2 * d.D0 - 3 * d.D2) * sq -
262 2 * pow2(invc2) * mF0 * Complex(d.G2, d.D2) +
263 invc2 * (2 * d.D0 - 3 * d.D2)) -
264 dw2 * (-2 * pow2(invGD) * mF0 *
265 (Complex(3 * d.G2 - 2 * d.G0, 3 * d.D2 - 2 * d.D0) +
266 2 * dx * invc2 * Complex(d.G2, d.D2)) -
267 4 * invGD * invc2 * mF0 * sq * Complex(d.G2, d.D2) +
268 2 * invGD * (2 * d.D0 - 3 * d.D2) * sq -
269 2 * pow2(invc2) * mF0 * Complex(d.G2, d.D2) +
270 invc2 * (2 * d.D0 - 3 * d.D2)))) /
271 (8 * sqrt_pi * invGD * mF0 * sq);
272 case CalcType::Voigt:
273 return -invGD *
274 (Complex(0, invGD) * dw1 *
275 (dx * (2 * d.D0 - 3 * d.D2) -
276 mF0 * Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2)) +
277 w1 * (2 * d.D0 - 3 * d.D2)) /
278 (2 * sqrt_pi * mF0);
279 case CalcType::LowXandHighY:
280 return -(4 * pow2(invGD) * (2 * d.D0 - 3 * d.D2) * (w1 - w2) * sq +
281 1i * (4 * dw1 * pow3(invGD) * sq *
282 (dx * (2 * d.D0 - 3 * d.D2) -
283 mF0 *
284 Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2)) +
285 dw2 * invc2 *
286 (2 * pow2(invGD) * mF0 *
287 (Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) -
288 2 * dx * invc2 * Complex(d.G2, d.D2)) -
289 4 * invGD * invc2 * mF0 * sq * Complex(d.G2, d.D2) +
290 2 * invGD * (2 * d.D0 - 3 * d.D2) * sq -
291 2 * pow2(invc2) * mF0 * Complex(d.G2, d.D2) +
292 invc2 * (2 * d.D0 - 3 * d.D2)))) /
293 (8 * sqrt_pi * invGD * mF0 * sq);
294 case CalcType::LowYandLowX:
295 return pow2(invc2) *
296 (4 * sq * Complex(d.G2, d.D2) * (sqrt_pi * w1 * sq - 1) -
297 sqrt_pi * (1i * dw1 * sq + w1) *
298 (Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) -
299 2 * dx * invc2 * Complex(d.G2, d.D2))) /
300 (2 * pi * sq);
301 case CalcType::LowYandHighX:
302 return -pow2(invc2) *
303 (x * (2 * x - 3) * Complex(d.G2, d.D2) +
304 (x - 3) * (Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) -
305 2 * dx * invc2 * Complex(d.G2, d.D2))) /
306 (2 * pi * pow3(x));
307 }
308 return {};
309}
310
311Complex SpeedDependentVoigt::dFdT(const Output &d, Numeric T) const noexcept {
312 switch (calcs) {
313 case CalcType::Full:
314 return (-pow2(invGD) * (w1 - w2) * sq *
315 (T * (2 * d.D0 - 3 * d.D2) * sqrt_ln_2 - pow3(invGD) * mF0) *
316 ln_16 +
317 1i * invc2 *
318 (dw1 *
319 (T * invGD * invc2 * mF0 * sq * Complex(d.G2, d.D2) *
320 ln_16 -
321 2 * invGD * sq *
322 (T * (2 * d.D0 - 3 * d.D2) * sqrt_ln_2 -
323 pow3(invGD) * mF0) *
324 sqrt_ln_2 -
325 (2 * T * pow2(invGD) * mF0 *
326 (Complex(3 * d.G2 - 2 * d.G0, 3 * d.D2 - 2 * d.D0) +
327 2 * dx * invc2 * Complex(d.G2, d.D2)) *
328 sqrt_ln_2 +
329 2 * T * pow2(invc2) * mF0 * Complex(d.G2, d.D2) *
330 sqrt_ln_2 +
331 invc2 * (T * (-2 * d.D0 + 3 * d.D2) * sqrt_ln_2 +
332 pow3(invGD) * mF0)) *
333 sqrt_ln_2) +
334 dw2 *
335 (T * invGD * invc2 * mF0 * sq * Complex(d.G2, d.D2) *
336 ln_16 +
337 2 * invGD * sq *
338 (T * (-2 * d.D0 + 3 * d.D2) * sqrt_ln_2 +
339 pow3(invGD) * mF0) *
340 sqrt_ln_2 +
341 (-2 * T * pow2(invGD) * mF0 *
342 (Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) -
343 2 * dx * invc2 * Complex(d.G2, d.D2)) *
344 sqrt_ln_2 +
345 2 * T * pow2(invc2) * mF0 * Complex(d.G2, d.D2) *
346 sqrt_ln_2 +
347 invc2 * (T * (-2 * d.D0 + 3 * d.D2) * sqrt_ln_2 +
348 pow3(invGD) * mF0)) *
349 sqrt_ln_2)) *
350 sqrt_ln_2) /
351 (8 * sqrt_pi * T * invGD * mF0 * sq * pow3(sqrt_ln_2));
352 case CalcType::Voigt:
353 return -invGD *
354 (-Complex(0, invGD) * dw1 *
355 (T * mF0 * Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) *
356 sqrt_ln_2 -
357 dx * (T * (2 * d.D0 - 3 * d.D2) * sqrt_ln_2 -
358 pow3(invGD) * mF0)) +
359 w1 * (T * (2 * d.D0 - 3 * d.D2) * sqrt_ln_2 - pow3(invGD) * mF0)) /
360 (2 * sqrt_pi * T * mF0 * sqrt_ln_2);
361 case CalcType::LowXandHighY:
362 return (-pow2(invGD) * (w1 - w2) * sq *
363 (T * (2 * d.D0 - 3 * d.D2) * sqrt_ln_2 - pow3(invGD) * mF0) *
364 ln_16 +
365 1i * (dw1 * pow3(invGD) * sq *
366 (T * mF0 *
367 Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) *
368 sqrt_ln_2 -
369 dx * (T * (2 * d.D0 - 3 * d.D2) * sqrt_ln_2 -
370 pow3(invGD) * mF0)) *
371 ln_16 +
372 dw2 * invc2 *
373 (T * invGD * invc2 * mF0 * sq * Complex(d.G2, d.D2) *
374 ln_16 +
375 2 * invGD * sq *
376 (T * (-2 * d.D0 + 3 * d.D2) * sqrt_ln_2 +
377 pow3(invGD) * mF0) *
378 sqrt_ln_2 +
379 (-2 * T * pow2(invGD) * mF0 *
380 (Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) -
381 2 * dx * invc2 * Complex(d.G2, d.D2)) *
382 sqrt_ln_2 +
383 2 * T * pow2(invc2) * mF0 * Complex(d.G2, d.D2) *
384 sqrt_ln_2 +
385 invc2 * (T * (-2 * d.D0 + 3 * d.D2) * sqrt_ln_2 +
386 pow3(invGD) * mF0)) *
387 sqrt_ln_2) *
388 sqrt_ln_2)) /
389 (8 * sqrt_pi * T * invGD * mF0 * sq * pow3(sqrt_ln_2));
390 case CalcType::LowYandLowX:
391 return pow2(invc2) *
392 (4 * sq * Complex(d.G2, d.D2) * (sqrt_pi * w1 * sq - 1) -
393 sqrt_pi * (1i * dw1 * sq + w1) *
394 (Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) -
395 2 * dx * invc2 * Complex(d.G2, d.D2))) /
396 (2 * pi * sq);
397 case CalcType::LowYandHighX:
398 return -pow2(invc2) *
399 (x * (2 * x - 3) * Complex(d.G2, d.D2) +
400 (x - 3) * (Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) -
401 2 * dx * invc2 * Complex(d.G2, d.D2))) /
402 (2 * pi * pow3(x));
403 }
404 return {};
405}
406
407Complex SpeedDependentVoigt::operator()(Numeric f) noexcept {
408 imag_val(dx) = mF0 - f;
409 x = dx * invc2;
410 update_calcs();
411 calc();
412 return F;
413}
414
415constexpr Numeric abs_squared(Complex z) noexcept {
416 return pow2(z.real()) + pow2(z.imag());
417}
418
419SpeedDependentVoigt::CalcType
420SpeedDependentVoigt::init(const Complex c2) const noexcept {
421 if (abs_squared(c2) == 0)
422 return CalcType::Voigt;
423 if (abs_squared(x) <= 9e-16 * abs_squared(sqrty * sqrty))
424 return CalcType::LowXandHighY;
425 if ((abs_squared(sqrty * sqrty) <= 1.e-30 * abs_squared(x)) and
426 abs_squared(std::sqrt(x)) <= 16.e6)
427 return CalcType::LowYandLowX; // Weird case, untested
428 if ((abs_squared(sqrty * sqrty) <= 1.e-30 * abs_squared(x)))
429 return CalcType::LowYandHighX;
430 return CalcType::Full;
431}
432
433void SpeedDependentVoigt::update_calcs() noexcept {
434 if (calcs not_eq CalcType::Voigt)
435 calcs = init(Complex(1, 1));
436}
437
438void SpeedDependentVoigt::calc() noexcept {
439 switch (calcs) {
440 case CalcType::Full:
441 sq = std::sqrt(x + sqrty * sqrty);
442 w1 = Faddeeva::w(1i * (sq - sqrty));
443 w2 = Faddeeva::w(1i * (sq + sqrty));
444 F = inv_sqrt_pi * invGD * (w1 - w2);
445 dw1 = 2i * (inv_sqrt_pi - (sq - sqrty) * w1);
446 dw2 = 2i * (inv_sqrt_pi - (sq + sqrty) * w2);
447 break;
448 case CalcType::Voigt:
449 w1 = Faddeeva::w(1i * dx * invGD);
450 F = inv_sqrt_pi * invGD * w1;
451 dw1 = 2i * (inv_sqrt_pi - dx * invGD * w1);
452 break;
453 case CalcType::LowXandHighY:
454 sq = std::sqrt(x + sqrty * sqrty);
455 w1 = Faddeeva::w(1i * dx * invGD);
456 w2 = Faddeeva::w(1i * (sq + sqrty));
457 F = inv_sqrt_pi * invGD * (w1 - w2);
458 dw1 = 2i * (inv_sqrt_pi - dx * invGD * w1);
459 dw2 = 2i * (inv_sqrt_pi - (sq + sqrty) * w2);
460 break;
461 case CalcType::LowYandLowX:
462 sq = std::sqrt(x);
463 w1 = Faddeeva::w(1i * sq);
464 F = 2 * inv_pi * invc2 * (1 - sqrt_pi * sq * w1);
465 dw1 = 2i * (inv_sqrt_pi - sq * w1);
466 break;
467 case CalcType::LowYandHighX:
468 F = inv_pi * invc2 * (1 / x - 1.5 / pow2(x));
469 break;
470 }
471}
472
473HartmannTran::HartmannTran(Numeric F0_noshift, const Output &ls,
474 Numeric GD_div_F0, Numeric dZ) noexcept
475 : G0(ls.G0), D0(ls.D0), G2(ls.G2), D2(ls.D2), FVC(ls.FVC), ETA(ls.ETA),
476 mF0(F0_noshift + dZ + (1 - ls.ETA) * (ls.D0 - 1.5 * ls.D2)),
477 invGD(sqrt_ln_2 / nonstd::abs(GD_div_F0 * mF0)),
478 deltax(ls.FVC + (1 - ls.ETA) * (ls.G0 - 3 * ls.G2 / 2), mF0),
479 sqrty(1 / (2 * (1 - ls.ETA) * Complex(ls.G2, ls.D2) * invGD)) {
480 calc();
481}
482
483Complex HartmannTran::dFdf() const noexcept {
484 constexpr Complex ddeltax = -1i;
485 const Complex dx = -ddeltax / ((ETA - 1) * Complex(G2, D2));
486 const Complex dsqrtxy = dx / (2 * sqrtxy);
487
488 switch (calcs) {
489 case CalcType::Full: {
490 const Complex dz1 = dsqrtxy;
491 const Complex dz2 = dsqrtxy;
492 const Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
493 const Complex dB =
494 sqrt_pi *
495 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
496 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) /
497 (2 * sqrty * (ETA - 1) * Complex(G2, D2));
498 const Complex dK = ETA * Complex(G2, D2) * dB +
499 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
500 return inv_pi * (-A * dK + K * dA) / pow2(K);
501 }
502 case CalcType::Noc2tLowZ: {
503 const Complex dz1 = invGD * ddeltax;
504 const Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
505 const Complex dB =
506 -invGD *
507 (sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) - dz1);
508 const Complex dK = ETA * Complex(G2, D2) * dB +
509 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
510 return inv_pi * (-A * dK + K * dA) / pow2(K);
511 }
512 case CalcType::Noc2tHighZ: {
513 const Complex dz1 = invGD * ddeltax;
514 const Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
515 const Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
516 invGD * dz1 / (2 * pow2(z1)) +
517 9 * invGD * dz1 / (4 * pow4(z1));
518 const Complex dK = ETA * Complex(G2, D2) * dB +
519 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
520 return inv_pi * (-A * dK + K * dA) / pow2(K);
521 }
522 case CalcType::LowXandHighY: {
523 const Complex dz1 = invGD * ddeltax;
524 const Complex dz2 = dsqrtxy;
525 const Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
526 const Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
527 invGD * dz1 / (2 * pow2(z1)) +
528 9 * invGD * dz1 / (4 * pow4(z1));
529 const Complex dK = ETA * Complex(G2, D2) * dB +
530 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
531 return inv_pi * (-A * dK + K * dA) / pow2(K);
532 }
533 case CalcType::LowYandLowX: {
534 const Complex dz1 = dsqrtxy;
535 const Complex dz2 = dx / (2 * sqrtx);
536 const Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
537 ((ETA - 1) * Complex(G2, D2));
538 const Complex dB =
539 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
540 (2 * pow2(sqrty) + x - 1) +
541 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
542 2 * (sqrt_pi * w2 * z2 - 1) * dx) /
543 ((ETA - 1) * Complex(G2, D2));
544 const Complex dK = ETA * Complex(G2, D2) * dB +
545 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
546 return inv_pi * (-A * dK + K * dA) / pow2(K);
547 }
548 case CalcType::LowYandHighX: {
549 const Complex dz1 = dsqrtxy;
550 const Complex dA = (x - 3) * dx / ((ETA - 1) * Complex(G2, D2) * pow3(x));
551 const Complex dB =
552 (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) -
553 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx + (2 * x - 3) * x * dx / 2) /
554 ((ETA - 1) * Complex(G2, D2) * pow3(x));
555 const Complex dK = ETA * Complex(G2, D2) * dB +
556 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
557 return inv_pi * (-A * dK + K * dA) / pow2(K);
558 }
559 }
560 return {};
561}
562
563Complex HartmannTran::dFdF0() const noexcept {
564 const Numeric dGD = (1 / (invGD * mF0));
565 const Numeric dinvGD = -dGD * pow2(invGD);
566 const Complex dsqrty =
567 dinvGD / (2 * (ETA - 1) * Complex(G2, D2) * pow2(invGD));
568 constexpr Complex ddeltax = 1i;
569 const Complex dx = -ddeltax / ((ETA - 1) * Complex(G2, D2));
570 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
571
572 switch (calcs) {
573 case CalcType::Full: {
574 const Complex dz1 = dsqrtxy - dsqrty;
575 const Complex dz2 = dsqrtxy + dsqrty;
576 const Complex dA =
577 sqrt_pi *
578 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
579 const Complex dB =
580 sqrt_pi *
581 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
582 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
583 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) *
584 sqrty) /
585 (2 * (ETA - 1) * Complex(G2, D2) * pow2(sqrty));
586 const Complex dK = ETA * Complex(G2, D2) * dB +
587 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
588 return inv_pi * (-A * dK + K * dA) / pow2(K);
589 }
590 case CalcType::Noc2tLowZ: {
591 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
592 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
593 const Complex dB =
594 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
595 dz1) *
596 invGD -
597 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
598 const Complex dK = ETA * Complex(G2, D2) * dB +
599 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
600 return inv_pi * (-A * dK + K * dA) / pow2(K);
601 }
602 case CalcType::Noc2tHighZ: {
603 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
604 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
605 const Complex dB =
606 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
607 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
608 9 * dz1) *
609 invGD) /
610 (4 * pow4(z1));
611 const Complex dK = ETA * Complex(G2, D2) * dB +
612 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
613 return inv_pi * (-A * dK + K * dA) / pow2(K);
614 }
615 case CalcType::LowXandHighY: {
616 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
617 const Complex dz2 = dsqrtxy + dsqrty;
618 const Complex dA =
619 sqrt_pi *
620 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
621 const Complex dB =
622 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
623 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
624 9 * dz1) *
625 invGD) /
626 (4 * pow4(z1));
627 const Complex dK = ETA * Complex(G2, D2) * dB +
628 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
629 return inv_pi * (-A * dK + K * dA) / pow2(K);
630 }
631 case CalcType::LowYandLowX: {
632 const Complex dz1 = dsqrtxy;
633 const Complex dz2 = dx / (2 * sqrtx);
634 const Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
635 ((ETA - 1) * Complex(G2, D2));
636 const Complex dB =
637 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
638 (2 * pow2(sqrty) + x - 1) +
639 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
640 2 * (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) /
641 ((ETA - 1) * Complex(G2, D2));
642 const Complex dK = ETA * Complex(G2, D2) * dB +
643 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
644 return inv_pi * (-A * dK + K * dA) / pow2(K);
645 }
646 case CalcType::LowYandHighX: {
647 const Complex dz1 = dsqrtxy;
648 const Complex dA = (x - 3) * dx / ((ETA - 1) * Complex(G2, D2) * pow3(x));
649 const Complex dB =
650 (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
651 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x / 2 -
652 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) /
653 ((ETA - 1) * Complex(G2, D2) * pow3(x));
654 const Complex dK = ETA * Complex(G2, D2) * dB +
655 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
656 return inv_pi * (-A * dK + K * dA) / pow2(K);
657 }
658 }
659 return {};
660}
661
662Complex HartmannTran::dFdD0(Numeric dD0) const noexcept {
663 const Numeric dmF0 = (1 - ETA) * dD0;
664 const Numeric dGD = (dmF0 / (invGD * mF0));
665 const Numeric dinvGD = -dGD * pow2(invGD);
666 const Complex dsqrty =
667 dinvGD / (2 * (ETA - 1) * Complex(G2, D2) * pow2(invGD));
668 const Complex ddeltax = Complex(0, 1 - ETA) * dD0;
669 const Complex dx = -ddeltax / ((ETA - 1) * Complex(G2, D2));
670 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
671
672 switch (calcs) {
673 case CalcType::Full: {
674 const Complex dz1 = dsqrtxy - dsqrty;
675 const Complex dz2 = dsqrtxy + dsqrty;
676 const Complex dA =
677 sqrt_pi *
678 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
679 const Complex dB =
680 sqrt_pi *
681 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
682 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
683 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) *
684 sqrty) /
685 (2 * (ETA - 1) * Complex(G2, D2) * pow2(sqrty));
686 const Complex dK = ETA * Complex(G2, D2) * dB + Complex(0, ETA * dD0) * A +
687 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
688 return inv_pi * (-A * dK + K * dA) / pow2(K);
689 }
690 case CalcType::Noc2tLowZ: {
691 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
692 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
693 const Complex dB =
694 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
695 dz1) *
696 invGD -
697 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
698 const Complex dK = ETA * Complex(G2, D2) * dB + Complex(0, ETA * dD0) * A +
699 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
700 return inv_pi * (-A * dK + K * dA) / pow2(K);
701 }
702 case CalcType::Noc2tHighZ: {
703 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
704 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
705 const Complex dB =
706 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
707 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
708 9 * dz1) *
709 invGD) /
710 (4 * pow4(z1));
711 const Complex dK = ETA * Complex(G2, D2) * dB + Complex(0, ETA * dD0) * A +
712 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
713 return inv_pi * (-A * dK + K * dA) / pow2(K);
714 }
715 case CalcType::LowXandHighY: {
716 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
717 const Complex dz2 = dsqrtxy + dsqrty;
718 const Complex dA =
719 sqrt_pi *
720 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
721 const Complex dB =
722 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
723 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
724 9 * dz1) *
725 invGD) /
726 (4 * pow4(z1));
727 const Complex dK = ETA * Complex(G2, D2) * dB + Complex(0, ETA * dD0) * A +
728 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
729 return inv_pi * (-A * dK + K * dA) / pow2(K);
730 }
731 case CalcType::LowYandLowX: {
732 const Complex dz1 = dsqrtxy;
733 const Complex dz2 = dx / (2 * sqrtx);
734 const Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
735 ((ETA - 1) * Complex(G2, D2));
736 const Complex dB =
737 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
738 (2 * pow2(sqrty) + x - 1) +
739 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
740 2 * (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) /
741 ((ETA - 1) * Complex(G2, D2));
742 const Complex dK = ETA * Complex(G2, D2) * dB + Complex(0, ETA * dD0) * A +
743 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
744 return inv_pi * (-A * dK + K * dA) / pow2(K);
745 }
746 case CalcType::LowYandHighX: {
747 const Complex dz1 = dsqrtxy;
748 const Complex dA = (x - 3) * dx / ((ETA - 1) * Complex(G2, D2) * pow3(x));
749 const Complex dB =
750 (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
751 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x / 2 -
752 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) /
753 ((ETA - 1) * Complex(G2, D2) * pow3(x));
754 const Complex dK = ETA * Complex(G2, D2) * dB + Complex(0, ETA * dD0) * A +
755 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
756 return inv_pi * (-A * dK + K * dA) / pow2(K);
757 }
758 }
759 return {};
760}
761
762Complex HartmannTran::dFdG0(Numeric dG0) const noexcept {
763 const Numeric ddeltax = (1 - ETA) * dG0;
764 const Complex dx = -ddeltax / ((ETA - 1) * Complex(G2, D2));
765 const Complex dsqrtxy = dx / (2 * sqrtxy);
766
767 switch (calcs) {
768 case CalcType::Full: {
769 const Complex dz1 = dsqrtxy;
770 const Complex dz2 = dsqrtxy;
771 const Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
772 const Complex dB =
773 sqrt_pi *
774 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
775 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) /
776 (2 * sqrty * (ETA - 1) * Complex(G2, D2));
777 const Complex dK = ETA * Complex(G2, D2) * dB + ETA * A * dG0 +
778 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
779 return inv_pi * (-A * dK + K * dA) / pow2(K);
780 }
781 case CalcType::Noc2tLowZ: {
782 const Complex dz1 = invGD * ddeltax;
783 const Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
784 const Complex dB =
785 -invGD *
786 (sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) - dz1);
787 const Complex dK = ETA * Complex(G2, D2) * dB + ETA * A * dG0 +
788 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
789 return inv_pi * (-A * dK + K * dA) / pow2(K);
790 }
791 case CalcType::Noc2tHighZ: {
792 const Complex dz1 = invGD * ddeltax;
793 const Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
794 const Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
795 invGD * dz1 / (2 * pow2(z1)) +
796 9 * invGD * dz1 / (4 * pow4(z1));
797 const Complex dK = ETA * Complex(G2, D2) * dB + ETA * A * dG0 +
798 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
799 return inv_pi * (-A * dK + K * dA) / pow2(K);
800 }
801 case CalcType::LowXandHighY: {
802 const Complex dz1 = invGD * ddeltax;
803 const Complex dz2 = dsqrtxy;
804 const Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
805 const Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
806 invGD * dz1 / (2 * pow2(z1)) +
807 9 * invGD * dz1 / (4 * pow4(z1));
808 const Complex dK = ETA * Complex(G2, D2) * dB + ETA * A * dG0 +
809 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
810 return inv_pi * (-A * dK + K * dA) / pow2(K);
811 }
812 case CalcType::LowYandLowX: {
813 const Complex dz1 = dsqrtxy;
814 const Complex dz2 = dx / (2 * sqrtx);
815 const Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
816 ((ETA - 1) * Complex(G2, D2));
817 const Complex dB =
818 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
819 (2 * pow2(sqrty) + x - 1) +
820 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
821 2 * (sqrt_pi * w2 * z2 - 1) * dx) /
822 ((ETA - 1) * Complex(G2, D2));
823 const Complex dK = ETA * Complex(G2, D2) * dB + ETA * A * dG0 +
824 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
825 return inv_pi * (-A * dK + K * dA) / pow2(K);
826 }
827 case CalcType::LowYandHighX: {
828 const Complex dz1 = dsqrtxy;
829 const Complex dA = (x - 3) * dx / ((ETA - 1) * Complex(G2, D2) * pow3(x));
830 const Complex dB =
831 (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) -
832 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx + (2 * x - 3) * x * dx / 2) /
833 ((ETA - 1) * Complex(G2, D2) * pow3(x));
834 const Complex dK = ETA * Complex(G2, D2) * dB + ETA * A * dG0 +
835 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
836 return inv_pi * (-A * dK + K * dA) / pow2(K);
837 }
838 }
839 return {};
840}
841
842Complex HartmannTran::dFdD2(Numeric dD2) const noexcept {
843 const Numeric dmF0 = -3 * (1 - ETA) * dD2 / 2;
844 const Numeric dGD = (dmF0 / (invGD * mF0));
845 const Numeric dinvGD = -dGD * pow2(invGD);
846 const Complex dsqrty = (Complex(G2, D2) * dinvGD + Complex(0, invGD) * dD2) /
847 (2 * (ETA - 1) * pow2(Complex(G2, D2)) * pow2(invGD));
848 const Complex ddeltax = 1.5 * Complex(0, ETA - 1) * dD2;
849 const Complex dx = (-Complex(G2, D2) * ddeltax + Complex(0, dD2) * deltax) /
850 ((ETA - 1) * pow2(Complex(G2, D2)));
851 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
852
853 switch (calcs) {
854 case CalcType::Full: {
855 const Complex dz1 = dsqrtxy - dsqrty;
856 const Complex dz2 = dsqrtxy + dsqrty;
857 const Complex dA =
858 sqrt_pi *
859 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
860 const Complex dB =
861 (sqrt_pi * Complex(G2, D2) *
862 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
863 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
864 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
865 2 * w2 * z2 * dz2) *
866 sqrty) -
867 1i *
868 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
869 2 * sqrty) *
870 sqrty * dD2) /
871 (2 * (ETA - 1) * pow2(Complex(G2, D2)) * pow2(sqrty));
872 const Complex dK = ETA * Complex(G2, D2) * dB -
873 Complex(0, 1.5 * dD2 * ETA) * A +
874 Complex(0, dD2 * ETA) * B +
875 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
876 return inv_pi * (-A * dK + K * dA) / pow2(K);
877 }
878 case CalcType::Noc2tLowZ: {
879 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
880 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
881 const Complex dB =
882 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
883 dz1) *
884 invGD -
885 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
886 const Complex dK = ETA * Complex(G2, D2) * dB -
887 Complex(0, 1.5 * dD2 * ETA) * A +
888 Complex(0, dD2 * ETA) * B +
889 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
890 return inv_pi * (-A * dK + K * dA) / pow2(K);
891 }
892 case CalcType::Noc2tHighZ: {
893 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
894 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
895 const Complex dB =
896 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
897 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
898 9 * dz1) *
899 invGD) /
900 (4 * pow4(z1));
901 const Complex dK = ETA * Complex(G2, D2) * dB -
902 Complex(0, 1.5 * dD2 * ETA) * A +
903 Complex(0, dD2 * ETA) * B +
904 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
905 return inv_pi * (-A * dK + K * dA) / pow2(K);
906 }
907 case CalcType::LowXandHighY: {
908 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
909 const Complex dz2 = dsqrtxy + dsqrty;
910 const Complex dA =
911 sqrt_pi *
912 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
913 const Complex dB =
914 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
915 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
916 9 * dz1) *
917 invGD) /
918 (4 * pow4(z1));
919 const Complex dK = ETA * Complex(G2, D2) * dB -
920 Complex(0, 1.5 * dD2 * ETA) * A +
921 Complex(0, dD2 * ETA) * B +
922 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
923 return inv_pi * (-A * dK + K * dA) / pow2(K);
924 }
925 case CalcType::LowYandLowX: {
926 const Complex dz1 = dsqrtxy;
927 const Complex dz2 = dx / (2 * sqrtx);
928 const Complex dA =
929 2 *
930 (sqrt_pi * Complex(G2, D2) * (w2 * dz2 + z2 * 1i * dw2 * dz2) -
931 1i * (sqrt_pi * w2 * z2 - 1) * dD2) /
932 ((ETA - 1) * pow2(Complex(G2, D2)));
933 const Complex dB =
934 (-2 * Complex(G2, D2) *
935 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
936 (2 * pow2(sqrty) + x - 1) +
937 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
938 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
939 1i *
940 (2 * sqrt_pi * w1 * z1 +
941 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
942 dD2) /
943 ((ETA - 1) * pow2(Complex(G2, D2)));
944 const Complex dK = ETA * Complex(G2, D2) * dB -
945 Complex(0, 1.5 * dD2 * ETA) * A +
946 Complex(0, dD2 * ETA) * B +
947 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
948 return inv_pi * (-A * dK + K * dA) / pow2(K);
949 }
950 case CalcType::LowYandHighX: {
951 const Complex dz1 = dsqrtxy;
952 const Complex dA =
953 (Complex(G2, D2) * (x - 3) * dx + 1i * (2 * x - 3) * x * dD2 / 2) /
954 ((ETA - 1) * pow2(Complex(G2, D2)) * pow3(x));
955 const Complex dB =
956 (Complex(G2, D2) *
957 (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
958 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
959 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) -
960 1i *
961 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
962 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
963 x * dD2) /
964 (2 * (ETA - 1) * pow2(Complex(G2, D2)) * pow3(x));
965 const Complex dK = ETA * Complex(G2, D2) * dB -
966 Complex(0, 1.5 * dD2 * ETA) * A +
967 Complex(0, dD2 * ETA) * B +
968 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
969 return inv_pi * (-A * dK + K * dA) / pow2(K);
970 }
971 }
972 return {};
973}
974
975Complex HartmannTran::dFdG2(Numeric dG2) const noexcept {
976 const Complex dsqrty = dG2 / (2 * invGD * (ETA - 1) * pow2(Complex(G2, D2)));
977 const Numeric ddeltax = 3 * (ETA - 1) * dG2 / 2;
978 const Complex dx = (-Complex(G2, D2) * ddeltax + deltax * dG2) /
979 ((ETA - 1) * pow2(Complex(G2, D2)));
980 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
981
982 switch (calcs) {
983 case CalcType::Full: {
984 const Complex dz1 = dsqrtxy - dsqrty;
985 const Complex dz2 = dsqrtxy + dsqrty;
986 const Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
987 const Complex dB =
988 (sqrt_pi * Complex(G2, D2) *
989 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
990 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
991 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
992 2 * w2 * z2 * dz2) *
993 sqrty) -
994 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
995 2 * sqrty) *
996 sqrty * dG2) /
997 (2 * (ETA - 1) * pow2(Complex(G2, D2)) * pow2(sqrty));
998 const Complex dK = ETA * Complex(G2, D2) * dB - 3 * ETA * A * dG2 / 2 +
999 ETA * B * dG2 +
1000 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1001 return inv_pi * (-A * dK + K * dA) / pow2(K);
1002 }
1003 case CalcType::Noc2tLowZ: {
1004 const Complex dz1 = invGD * ddeltax;
1005 const Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
1006 const Complex dB =
1007 -invGD *
1008 (sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) - dz1);
1009 const Complex dK = ETA * Complex(G2, D2) * dB - 3 * ETA * A * dG2 / 2 +
1010 ETA * B * dG2 +
1011 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1012 return inv_pi * (-A * dK + K * dA) / pow2(K);
1013 }
1014 case CalcType::Noc2tHighZ: {
1015 const Complex dz1 = invGD * ddeltax;
1016 const Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
1017 const Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
1018 invGD * dz1 / (2 * pow2(z1)) +
1019 9 * invGD * dz1 / (4 * pow4(z1));
1020 const Complex dK = ETA * Complex(G2, D2) * dB - 3 * ETA * A * dG2 / 2 +
1021 ETA * B * dG2 +
1022 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1023 return inv_pi * (-A * dK + K * dA) / pow2(K);
1024 }
1025 case CalcType::LowXandHighY: {
1026 const Complex dz1 = invGD * ddeltax;
1027 const Complex dz2 = dsqrtxy + dsqrty;
1028 const Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
1029 const Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
1030 invGD * dz1 / (2 * pow2(z1)) +
1031 9 * invGD * dz1 / (4 * pow4(z1));
1032 const Complex dK = ETA * Complex(G2, D2) * dB - 3 * ETA * A * dG2 / 2 +
1033 ETA * B * dG2 +
1034 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1035 return inv_pi * (-A * dK + K * dA) / pow2(K);
1036 }
1037 case CalcType::LowYandLowX: {
1038 const Complex dz1 = dsqrtxy;
1039 const Complex dz2 = dx / (2 * sqrtx);
1040 const Complex dA =
1041 2 *
1042 (sqrt_pi * Complex(G2, D2) * (w2 * dz2 + z2 * 1i * dw2 * dz2) -
1043 (sqrt_pi * w2 * z2 - 1) * dG2) /
1044 ((ETA - 1) * pow2(Complex(G2, D2)));
1045 const Complex dB =
1046 (-2 * Complex(G2, D2) *
1047 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1048 (2 * pow2(sqrty) + x - 1) +
1049 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
1050 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
1051 (2 * sqrt_pi * w1 * z1 +
1052 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
1053 dG2) /
1054 ((ETA - 1) * pow2(Complex(G2, D2)));
1055 const Complex dK = ETA * Complex(G2, D2) * dB - 3 * ETA * A * dG2 / 2 +
1056 ETA * B * dG2 +
1057 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1058 return inv_pi * (-A * dK + K * dA) / pow2(K);
1059 }
1060 case CalcType::LowYandHighX: {
1061 const Complex dz1 = dsqrtxy;
1062 const Complex dA =
1063 (Complex(G2, D2) * (x - 3) * dx + (2 * x - 3) * x * dG2 / 2) /
1064 ((ETA - 1) * pow2(Complex(G2, D2)) * pow3(x));
1065 const Complex dB =
1066 (Complex(G2, D2) *
1067 (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
1068 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
1069 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) -
1070 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1071 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1072 x * dG2) /
1073 (2 * (ETA - 1) * pow2(Complex(G2, D2)) * pow3(x));
1074 const Complex dK = ETA * Complex(G2, D2) * dB - 3 * ETA * A * dG2 / 2 +
1075 ETA * B * dG2 +
1076 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1077 return inv_pi * (-A * dK + K * dA) / pow2(K);
1078 }
1079 }
1080 return {};
1081}
1082
1083Complex HartmannTran::dFdFVC(Numeric dFVC) const noexcept {
1084 const Numeric ddeltax = dFVC;
1085 const Complex dx = -ddeltax / ((ETA - 1) * Complex(G2, D2));
1086 const Complex dsqrtxy = dx / (2 * sqrtxy);
1087
1088 switch (calcs) {
1089 case CalcType::Full: {
1090 const Complex dz1 = dsqrtxy;
1091 const Complex dz2 = dsqrtxy;
1092 const Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
1093 const Complex dB =
1094 sqrt_pi *
1095 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
1096 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) /
1097 (2 * sqrty * (ETA - 1) * Complex(G2, D2));
1098 const Complex dK =
1099 ETA * Complex(G2, D2) * dB +
1100 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA - A * dFVC;
1101 return inv_pi * (-A * dK + K * dA) / pow2(K);
1102 }
1103 case CalcType::Noc2tLowZ: {
1104 const Complex dz1 = invGD * ddeltax;
1105 const Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
1106 const Complex dB =
1107 -invGD *
1108 (sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) - dz1);
1109 const Complex dK =
1110 ETA * Complex(G2, D2) * dB +
1111 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA - A * dFVC;
1112 return inv_pi * (-A * dK + K * dA) / pow2(K);
1113 }
1114 case CalcType::Noc2tHighZ: {
1115 const Complex dz1 = invGD * ddeltax;
1116 const Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
1117 const Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
1118 invGD * dz1 / (2 * pow2(z1)) +
1119 9 * invGD * dz1 / (4 * pow4(z1));
1120 const Complex dK =
1121 ETA * Complex(G2, D2) * dB +
1122 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA - A * dFVC;
1123 return inv_pi * (-A * dK + K * dA) / pow2(K);
1124 }
1125 case CalcType::LowXandHighY: {
1126 const Complex dz1 = invGD * ddeltax;
1127 const Complex dz2 = dsqrtxy;
1128 const Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
1129 const Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
1130 invGD * dz1 / (2 * pow2(z1)) +
1131 9 * invGD * dz1 / (4 * pow4(z1));
1132 const Complex dK =
1133 ETA * Complex(G2, D2) * dB +
1134 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA - A * dFVC;
1135 return inv_pi * (-A * dK + K * dA) / pow2(K);
1136 }
1137 case CalcType::LowYandLowX: {
1138 const Complex dz1 = dsqrtxy;
1139 const Complex dz2 = dx / (2 * sqrtx);
1140 const Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
1141 ((ETA - 1) * Complex(G2, D2));
1142 const Complex dB =
1143 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1144 (2 * pow2(sqrty) + x - 1) +
1145 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
1146 2 * (sqrt_pi * w2 * z2 - 1) * dx) /
1147 ((ETA - 1) * Complex(G2, D2));
1148 const Complex dK =
1149 ETA * Complex(G2, D2) * dB +
1150 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA - A * dFVC;
1151 return inv_pi * (-A * dK + K * dA) / pow2(K);
1152 }
1153 case CalcType::LowYandHighX: {
1154 const Complex dz1 = dsqrtxy;
1155 const Complex dA = (x - 3) * dx / ((ETA - 1) * Complex(G2, D2) * pow3(x));
1156 const Complex dB =
1157 (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) -
1158 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx + (2 * x - 3) * x * dx / 2) /
1159 ((ETA - 1) * Complex(G2, D2) * pow3(x));
1160 const Complex dK =
1161 ETA * Complex(G2, D2) * dB +
1162 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA - A * dFVC;
1163 return inv_pi * (-A * dK + K * dA) / pow2(K);
1164 }
1165 }
1166 return {};
1167}
1168
1169Complex HartmannTran::dFdETA(Numeric dETA) const noexcept {
1170 const Numeric dmF0 = -(D0 - 3 * D2 / 2) * dETA;
1171 const Numeric dGD = (dmF0 / (invGD * mF0));
1172 const Numeric dinvGD = -dGD * pow2(invGD);
1173 const Complex dsqrty = ((ETA - 1) * dinvGD + invGD * dETA) /
1174 (2 * Complex(G2, D2) * pow2(ETA - 1) * pow2(invGD));
1175 const Complex ddeltax = -dETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2);
1176 const Complex dx = (-(ETA - 1) * ddeltax + deltax * dETA) /
1177 (Complex(G2, D2) * pow2(ETA - 1));
1178 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
1179
1180 switch (calcs) {
1181 case CalcType::Full: {
1182 const Complex dz1 = dsqrtxy - dsqrty;
1183 const Complex dz2 = dsqrtxy + dsqrty;
1184 const Complex dA =
1185 sqrt_pi *
1186 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1187 const Complex dB =
1188 (sqrt_pi *
1189 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
1190 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
1191 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
1192 2 * w2 * z2 * dz2) *
1193 sqrty) *
1194 (ETA - 1) -
1195 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
1196 2 * sqrty) *
1197 sqrty * dETA) /
1198 (2 * Complex(G2, D2) * pow2(ETA - 1) * pow2(sqrty));
1199 const Complex dK =
1200 (-FVC + Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA) * dA +
1201 Complex(G2, D2) * B * dETA + Complex(G2, D2) * ETA * dB -
1202 Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * A * dETA;
1203 return inv_pi * (-A * dK + K * dA) / pow2(K);
1204 }
1205 case CalcType::Noc2tLowZ: {
1206 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1207 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1208 const Complex dB =
1209 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
1210 dz1) *
1211 invGD -
1212 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
1213 const Complex dK =
1214 (-FVC + Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA) * dA +
1215 Complex(G2, D2) * B * dETA + Complex(G2, D2) * ETA * dB -
1216 Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * A * dETA;
1217 return inv_pi * (-A * dK + K * dA) / pow2(K);
1218 }
1219 case CalcType::Noc2tHighZ: {
1220 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1221 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1222 const Complex dB =
1223 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1224 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1225 9 * dz1) *
1226 invGD) /
1227 (4 * pow4(z1));
1228 const Complex dK =
1229 (-FVC + Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA) * dA +
1230 Complex(G2, D2) * B * dETA + Complex(G2, D2) * ETA * dB -
1231 Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * A * dETA;
1232 return inv_pi * (-A * dK + K * dA) / pow2(K);
1233 }
1234 case CalcType::LowXandHighY: {
1235 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1236 const Complex dz2 = dsqrtxy + dsqrty;
1237 const Complex dA =
1238 sqrt_pi *
1239 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1240 const Complex dB =
1241 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1242 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1243 9 * dz1) *
1244 invGD) /
1245 (4 * pow4(z1));
1246 const Complex dK =
1247 (-FVC + Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA) * dA +
1248 Complex(G2, D2) * B * dETA + Complex(G2, D2) * ETA * dB -
1249 Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * A * dETA;
1250 return inv_pi * (-A * dK + K * dA) / pow2(K);
1251 }
1252 case CalcType::LowYandLowX: {
1253 const Complex dz1 = dsqrtxy;
1254 const Complex dz2 = dx / (2 * sqrtx);
1255 const Complex dA = 2 *
1256 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) * (ETA - 1) -
1257 (sqrt_pi * w2 * z2 - 1) * dETA) /
1258 (Complex(G2, D2) * pow2(ETA - 1));
1259 const Complex dB =
1260 (-2 * (ETA - 1) *
1261 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1262 (2 * pow2(sqrty) + x - 1) +
1263 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
1264 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
1265 (2 * sqrt_pi * w1 * z1 +
1266 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
1267 dETA) /
1268 (Complex(G2, D2) * pow2(ETA - 1));
1269 const Complex dK =
1270 (-FVC + Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA) * dA +
1271 Complex(G2, D2) * B * dETA + Complex(G2, D2) * ETA * dB -
1272 Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * A * dETA;
1273 return inv_pi * (-A * dK + K * dA) / pow2(K);
1274 }
1275 case CalcType::LowYandHighX: {
1276 const Complex dz1 = dsqrtxy;
1277 const Complex dA = ((ETA - 1) * (x - 3) * dx + (2 * x - 3) * x * dETA / 2) /
1278 (Complex(G2, D2) * pow2(ETA - 1) * pow3(x));
1279 const Complex dB =
1280 (-(2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1281 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1282 x * dETA +
1283 (ETA - 1) *
1284 (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
1285 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
1286 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx)) /
1287 (2 * Complex(G2, D2) * pow2(ETA - 1) * pow3(x));
1288 const Complex dK =
1289 (-FVC + Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA) * dA +
1290 Complex(G2, D2) * B * dETA + Complex(G2, D2) * ETA * dB -
1291 Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * A * dETA;
1292 return inv_pi * (-A * dK + K * dA) / pow2(K);
1293 }
1294 }
1295 return {};
1296}
1297
1298Complex HartmannTran::dFdH(Numeric dZ) const noexcept {
1299 const Numeric dmF0 = dZ;
1300 const Numeric dGD = (dmF0 / (invGD * mF0));
1301 const Numeric dinvGD = -dGD * pow2(invGD);
1302 const Complex dsqrty =
1303 dinvGD / (2 * (ETA - 1) * Complex(G2, D2) * pow2(invGD));
1304 const Complex ddeltax = Complex(0, dZ);
1305 const Complex dx = -ddeltax / ((ETA - 1) * Complex(G2, D2));
1306 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
1307
1308 switch (calcs) {
1309 case CalcType::Full: {
1310 const Complex dz1 = dsqrtxy - dsqrty;
1311 const Complex dz2 = dsqrtxy + dsqrty;
1312 const Complex dA =
1313 sqrt_pi *
1314 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1315 const Complex dB =
1316 sqrt_pi *
1317 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
1318 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
1319 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) *
1320 sqrty) /
1321 (2 * (ETA - 1) * Complex(G2, D2) * pow2(sqrty));
1322 const Complex dK = ETA * Complex(G2, D2) * dB +
1323 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1324 return inv_pi * (-A * dK + K * dA) / pow2(K);
1325 }
1326 case CalcType::Noc2tLowZ: {
1327 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1328 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1329 const Complex dB =
1330 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
1331 dz1) *
1332 invGD -
1333 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
1334 const Complex dK = ETA * Complex(G2, D2) * dB +
1335 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1336 return inv_pi * (-A * dK + K * dA) / pow2(K);
1337 }
1338 case CalcType::Noc2tHighZ: {
1339 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1340 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1341 const Complex dB =
1342 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1343 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1344 9 * dz1) *
1345 invGD) /
1346 (4 * pow4(z1));
1347 const Complex dK = ETA * Complex(G2, D2) * dB +
1348 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1349 return inv_pi * (-A * dK + K * dA) / pow2(K);
1350 }
1351 case CalcType::LowXandHighY: {
1352 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1353 const Complex dz2 = dsqrtxy + dsqrty;
1354 const Complex dA =
1355 sqrt_pi *
1356 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1357 const Complex dB =
1358 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1359 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1360 9 * dz1) *
1361 invGD) /
1362 (4 * pow4(z1));
1363 const Complex dK = ETA * Complex(G2, D2) * dB +
1364 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1365 return inv_pi * (-A * dK + K * dA) / pow2(K);
1366 }
1367 case CalcType::LowYandLowX: {
1368 const Complex dz1 = dsqrtxy;
1369 const Complex dz2 = dx / (2 * sqrtx);
1370 const Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
1371 ((ETA - 1) * Complex(G2, D2));
1372 const Complex dB =
1373 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1374 (2 * pow2(sqrty) + x - 1) +
1375 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
1376 2 * (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) /
1377 ((ETA - 1) * Complex(G2, D2));
1378 const Complex dK = ETA * Complex(G2, D2) * dB +
1379 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1380 return inv_pi * (-A * dK + K * dA) / pow2(K);
1381 }
1382 case CalcType::LowYandHighX: {
1383 const Complex dz1 = dsqrtxy;
1384 const Complex dA = (x - 3) * dx / ((ETA - 1) * Complex(G2, D2) * pow3(x));
1385 const Complex dB =
1386 (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
1387 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x / 2 -
1388 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) /
1389 ((ETA - 1) * Complex(G2, D2) * pow3(x));
1390 const Complex dK = ETA * Complex(G2, D2) * dB +
1391 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1392 return inv_pi * (-A * dK + K * dA) / pow2(K);
1393 }
1394 }
1395 return {};
1396}
1397
1398Complex HartmannTran::dFdVMR(const Output &d) const noexcept {
1399 const Numeric dmF0 =
1400 (1 - ETA) * (d.D0 - 3 * d.D2 / 2) - (D0 - 3 * D2 / 2) * d.ETA;
1401 const Numeric dGD = (dmF0 / (invGD * mF0));
1402 const Numeric dinvGD = -dGD * pow2(invGD);
1403 const Complex dsqrty =
1404 (Complex(G2, D2) * (ETA - 1) * dinvGD + Complex(G2, D2) * invGD * d.ETA +
1405 Complex(d.G2, d.D2) * (ETA - 1) * invGD) /
1406 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow2(invGD));
1407 const Complex ddeltax =
1408 -(ETA - 1) * Complex(d.G0 - 1.5 * d.G2, d.D0 - 1.5 * d.D2) -
1409 Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * d.ETA + d.FVC;
1410 const Complex dx = (-Complex(G2, D2) * (ETA - 1) * ddeltax +
1411 Complex(G2, D2) * deltax * d.ETA +
1412 Complex(d.G2, d.D2) * (ETA - 1) * deltax) /
1413 (pow2(Complex(G2, D2)) * pow2(ETA - 1));
1414 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
1415
1416 switch (calcs) {
1417 case CalcType::Full: {
1418 const Complex dz1 = dsqrtxy - dsqrty;
1419 const Complex dz2 = dsqrtxy + dsqrty;
1420 const Complex dA =
1421 sqrt_pi *
1422 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1423 const Complex dB =
1424 (sqrt_pi * Complex(G2, D2) *
1425 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
1426 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
1427 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
1428 2 * w2 * z2 * dz2) *
1429 sqrty) *
1430 (ETA - 1) -
1431 Complex(G2, D2) *
1432 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
1433 2 * sqrty) *
1434 sqrty * d.ETA -
1435 Complex(d.G2, d.D2) * (ETA - 1) *
1436 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
1437 2 * sqrty) *
1438 sqrty) /
1439 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow2(sqrty));
1440 const Complex dK =
1441 Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1442 Complex(d.G2, d.D2) * B * ETA +
1443 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1444 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1445 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1446 A;
1447 return inv_pi * (-A * dK + K * dA) / pow2(K);
1448 }
1449 case CalcType::Noc2tLowZ: {
1450 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1451 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1452 const Complex dB =
1453 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
1454 dz1) *
1455 invGD -
1456 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
1457 const Complex dK =
1458 Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1459 Complex(d.G2, d.D2) * B * ETA +
1460 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1461 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1462 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1463 A;
1464 return inv_pi * (-A * dK + K * dA) / pow2(K);
1465 }
1466 case CalcType::Noc2tHighZ: {
1467 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1468 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1469 const Complex dB =
1470 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1471 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1472 9 * dz1) *
1473 invGD) /
1474 (4 * pow4(z1));
1475 const Complex dK =
1476 Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1477 Complex(d.G2, d.D2) * B * ETA +
1478 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1479 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1480 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1481 A;
1482 return inv_pi * (-A * dK + K * dA) / pow2(K);
1483 }
1484 case CalcType::LowXandHighY: {
1485 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1486 const Complex dz2 = dsqrtxy + dsqrty;
1487 const Complex dA =
1488 sqrt_pi *
1489 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1490 const Complex dB =
1491 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1492 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1493 9 * dz1) *
1494 invGD) /
1495 (4 * pow4(z1));
1496 const Complex dK =
1497 Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1498 Complex(d.G2, d.D2) * B * ETA +
1499 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1500 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1501 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1502 A;
1503 return inv_pi * (-A * dK + K * dA) / pow2(K);
1504 }
1505 case CalcType::LowYandLowX: {
1506 const Complex dz1 = dsqrtxy;
1507 const Complex dz2 = dx / (2 * sqrtx);
1508 const Complex dA =
1509 2 *
1510 (sqrt_pi * Complex(G2, D2) * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1511 (ETA - 1) -
1512 Complex(G2, D2) * (sqrt_pi * w2 * z2 - 1) * d.ETA -
1513 Complex(d.G2, d.D2) * (sqrt_pi * w2 * z2 - 1) * (ETA - 1)) /
1514 (pow2(Complex(G2, D2)) * pow2(ETA - 1));
1515 const Complex dB =
1516 (-2 * Complex(G2, D2) * (ETA - 1) *
1517 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1518 (2 * pow2(sqrty) + x - 1) +
1519 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
1520 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
1521 Complex(G2, D2) *
1522 (2 * sqrt_pi * w1 * z1 +
1523 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
1524 d.ETA +
1525 Complex(d.G2, d.D2) * (ETA - 1) *
1526 (2 * sqrt_pi * w1 * z1 +
1527 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1)) /
1528 (pow2(Complex(G2, D2)) * pow2(ETA - 1));
1529 const Complex dK =
1530 Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1531 Complex(d.G2, d.D2) * B * ETA +
1532 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1533 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1534 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1535 A;
1536 return inv_pi * (-A * dK + K * dA) / pow2(K);
1537 }
1538 case CalcType::LowYandHighX: {
1539 const Complex dz1 = dsqrtxy;
1540 const Complex dA = (2 * Complex(G2, D2) * (ETA - 1) * (x - 3) * dx +
1541 Complex(G2, D2) * (2 * x - 3) * x * d.ETA +
1542 Complex(d.G2, d.D2) * (ETA - 1) * (2 * x - 3) * x) /
1543 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow3(x));
1544 const Complex dB =
1545 (-Complex(G2, D2) *
1546 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1547 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1548 x * d.ETA +
1549 Complex(G2, D2) * (ETA - 1) *
1550 (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
1551 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
1552 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) -
1553 Complex(d.G2, d.D2) *
1554 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1555 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1556 (ETA - 1) * x) /
1557 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow3(x));
1558 const Complex dK =
1559 Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1560 Complex(d.G2, d.D2) * B * ETA +
1561 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1562 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1563 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1564 A;
1565 return inv_pi * (-A * dK + K * dA) / pow2(K);
1566 }
1567 }
1568 return {};
1569}
1570
1571Complex HartmannTran::dFdT(const Output &d, Numeric T) const noexcept {
1572 const Numeric dmF0 =
1573 (1 - ETA) * (d.D0 - 3 * d.D2 / 2) - (D0 - 3 * D2 / 2) * d.ETA;
1574 const Numeric dGD =
1575 (dmF0 / (invGD * mF0)) - invGD * invGD / (2 * T * sqrt_ln_2);
1576 const Numeric dinvGD = -dGD * pow2(invGD);
1577 const Complex dsqrty =
1578 (Complex(G2, D2) * (ETA - 1) * dinvGD + Complex(G2, D2) * invGD * d.ETA +
1579 Complex(d.G2, d.D2) * (ETA - 1) * invGD) /
1580 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow2(invGD));
1581 const Complex ddeltax =
1582 -(ETA - 1) * Complex(d.G0 - 1.5 * d.G2, d.D0 - 1.5 * d.D2) -
1583 Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * d.ETA + d.FVC;
1584 const Complex dx = (-Complex(G2, D2) * (ETA - 1) * ddeltax +
1585 Complex(G2, D2) * deltax * d.ETA +
1586 Complex(d.G2, d.D2) * (ETA - 1) * deltax) /
1587 (pow2(Complex(G2, D2)) * pow2(ETA - 1));
1588 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
1589
1590 switch (calcs) {
1591 case CalcType::Full: {
1592 const Complex dz1 = dsqrtxy - dsqrty;
1593 const Complex dz2 = dsqrtxy + dsqrty;
1594 const Complex dA =
1595 sqrt_pi *
1596 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1597 const Complex dB =
1598 (sqrt_pi * Complex(G2, D2) *
1599 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
1600 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
1601 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
1602 2 * w2 * z2 * dz2) *
1603 sqrty) *
1604 (ETA - 1) -
1605 Complex(G2, D2) *
1606 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
1607 2 * sqrty) *
1608 sqrty * d.ETA -
1609 Complex(d.G2, d.D2) * (ETA - 1) *
1610 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
1611 2 * sqrty) *
1612 sqrty) /
1613 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow2(sqrty));
1614 const Complex dK =
1615 Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1616 Complex(d.G2, d.D2) * B * ETA +
1617 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1618 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1619 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1620 A;
1621 return inv_pi * (-A * dK + K * dA) / pow2(K);
1622 }
1623 case CalcType::Noc2tLowZ: {
1624 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1625 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1626 const Complex dB =
1627 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
1628 dz1) *
1629 invGD -
1630 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
1631 const Complex dK =
1632 Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1633 Complex(d.G2, d.D2) * B * ETA +
1634 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1635 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1636 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1637 A;
1638 return inv_pi * (-A * dK + K * dA) / pow2(K);
1639 }
1640 case CalcType::Noc2tHighZ: {
1641 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1642 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1643 const Complex dB =
1644 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1645 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1646 9 * dz1) *
1647 invGD) /
1648 (4 * pow4(z1));
1649 const Complex dK =
1650 Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1651 Complex(d.G2, d.D2) * B * ETA +
1652 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1653 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1654 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1655 A;
1656 return inv_pi * (-A * dK + K * dA) / pow2(K);
1657 }
1658 case CalcType::LowXandHighY: {
1659 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1660 const Complex dz2 = dsqrtxy + dsqrty;
1661 const Complex dA =
1662 sqrt_pi *
1663 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1664 const Complex dB =
1665 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1666 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1667 9 * dz1) *
1668 invGD) /
1669 (4 * pow4(z1));
1670 const Complex dK =
1671 Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1672 Complex(d.G2, d.D2) * B * ETA +
1673 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1674 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1675 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1676 A;
1677 return inv_pi * (-A * dK + K * dA) / pow2(K);
1678 }
1679 case CalcType::LowYandLowX: {
1680 const Complex dz1 = dsqrtxy;
1681 const Complex dz2 = dx / (2 * sqrtx);
1682 const Complex dA =
1683 2 *
1684 (sqrt_pi * Complex(G2, D2) * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1685 (ETA - 1) -
1686 Complex(G2, D2) * (sqrt_pi * w2 * z2 - 1) * d.ETA -
1687 Complex(d.G2, d.D2) * (sqrt_pi * w2 * z2 - 1) * (ETA - 1)) /
1688 (pow2(Complex(G2, D2)) * pow2(ETA - 1));
1689 const Complex dB =
1690 (-2 * Complex(G2, D2) * (ETA - 1) *
1691 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1692 (2 * pow2(sqrty) + x - 1) +
1693 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
1694 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
1695 Complex(G2, D2) *
1696 (2 * sqrt_pi * w1 * z1 +
1697 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
1698 d.ETA +
1699 Complex(d.G2, d.D2) * (ETA - 1) *
1700 (2 * sqrt_pi * w1 * z1 +
1701 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1)) /
1702 (pow2(Complex(G2, D2)) * pow2(ETA - 1));
1703 const Complex dK =
1704 Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1705 Complex(d.G2, d.D2) * B * ETA +
1706 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1707 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1708 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1709 A;
1710 return inv_pi * (-A * dK + K * dA) / pow2(K);
1711 }
1712 case CalcType::LowYandHighX: {
1713 const Complex dz1 = dsqrtxy;
1714 const Complex dA = (2 * Complex(G2, D2) * (ETA - 1) * (x - 3) * dx +
1715 Complex(G2, D2) * (2 * x - 3) * x * d.ETA +
1716 Complex(d.G2, d.D2) * (ETA - 1) * (2 * x - 3) * x) /
1717 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow3(x));
1718 const Complex dB =
1719 (-Complex(G2, D2) *
1720 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1721 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1722 x * d.ETA +
1723 Complex(G2, D2) * (ETA - 1) *
1724 (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
1725 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
1726 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) -
1727 Complex(d.G2, d.D2) *
1728 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1729 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1730 (ETA - 1) * x) /
1731 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow3(x));
1732 const Complex dK =
1733 Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1734 Complex(d.G2, d.D2) * B * ETA +
1735 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1736 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1737 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1738 A;
1739 return inv_pi * (-A * dK + K * dA) / pow2(K);
1740 }
1741 }
1742 return {};
1743}
1744
1745Complex HartmannTran::operator()(Numeric f) noexcept {
1746 imag_val(deltax) = mF0 - f;
1747 x = deltax / ((1 - ETA) * Complex(G2, D2));
1748 sqrtxy = std::sqrt(x + sqrty * sqrty);
1749 update_calcs();
1750 calc();
1751 return F;
1752}
1753
1754HartmannTran::CalcType HartmannTran::init(const Complex c2t) const noexcept {
1755 if (abs_squared(c2t) == 0)
1756 return CalcType::Noc2tHighZ; // nb. Value of high/low changes elsewhere
1757 if (abs_squared(x) <= 9e-16 * abs_squared(sqrty * sqrty))
1758 return CalcType::LowXandHighY;
1759 if ((abs_squared(sqrty * sqrty) <= 1.e-30 * abs_squared(x)) and
1760 abs_squared(std::sqrt(x)) <= 16.e6)
1761 return CalcType::LowYandLowX; // Weird case, untested
1762 if ((abs_squared(sqrty * sqrty) <= 1.e-30 * abs_squared(x)))
1763 return CalcType::LowYandHighX;
1764 return CalcType::Full;
1765}
1766
1767void HartmannTran::update_calcs() noexcept {
1768 calcs = init((1 - ETA) * Complex(G2, D2));
1769}
1770
1771void HartmannTran::calc() noexcept {
1772 switch (calcs) {
1773 case CalcType::Full:
1774 z1 = sqrtxy - sqrty;
1775 z2 = sqrtxy + sqrty;
1776 w1 = Faddeeva::w(1i * z1);
1777 w2 = Faddeeva::w(1i * z2);
1778 A = sqrt_pi * invGD * (w1 - w2);
1779 B = (-1 + sqrt_pi / (2 * sqrty) * (1 - pow2(z1)) * w1 -
1780 sqrt_pi / (2 * sqrty) * (1 - pow2(z2)) * w2) /
1781 ((1 - ETA) * Complex(G2, D2));
1782 break;
1783 case CalcType::Noc2tLowZ:
1784 case CalcType::Noc2tHighZ:
1785 z1 = deltax * invGD;
1786 w1 = Faddeeva::w(1i * z1);
1787 A = sqrt_pi * invGD * w1;
1788 if (abs_squared(z1) < 16e6) {
1789 calcs = CalcType::Noc2tLowZ;
1790 B = sqrt_pi * invGD * ((1 - pow2(z1)) * w1 + z1 / sqrt_pi);
1791 } else {
1792 calcs = CalcType::Noc2tHighZ;
1793 B = invGD * (sqrt_pi * w1 + 1 / z1 / 2 - 3 / pow3(z1) / 4);
1794 }
1795 break;
1796 case CalcType::LowXandHighY:
1797 z1 = deltax * invGD;
1798 z2 = sqrtxy + sqrty;
1799 w1 = Faddeeva::w(1i * z1);
1800 w2 = Faddeeva::w(1i * z2);
1801 A = sqrt_pi * invGD * (w1 - w2);
1802 B = invGD * (sqrt_pi * w1 + 1 / z1 / 2 - 3 / pow3(z1) / 4);
1803 break;
1804 case CalcType::LowYandLowX:
1805 sqrtx = std::sqrt(x);
1806 z1 = sqrtxy;
1807 z2 = sqrtx;
1808 w1 = Faddeeva::w(1i * z1);
1809 w2 = Faddeeva::w(1i * z2);
1810 A = (2 * sqrt_pi / ((1 - ETA) * Complex(G2, D2))) * (inv_sqrt_pi - z2 * w2);
1811 B = (1 / ((1 - ETA) * Complex(G2, D2))) *
1812 (-1 +
1813 2 * sqrt_pi * (1 - x - 2 * sqrty * sqrty) * (1 / sqrt_pi - z2 * w2) +
1814 2 * sqrt_pi * z1 * w1);
1815 break;
1816 case CalcType::LowYandHighX:
1817 z1 = sqrtxy;
1818 w1 = Faddeeva::w(1i * z1);
1819 A = (1 / ((1 - ETA) * Complex(G2, D2))) * (1 / x - 3 / pow2(x) / 2);
1820 B = (1 / ((1 - ETA) * Complex(G2, D2))) *
1821 (-1 + (1 - x - 2 * sqrty * sqrty) * (1 / x - 3 / pow2(x) / 2) +
1822 2 * sqrt_pi * z1 * w1);
1823 break;
1824 }
1825
1826 dw1 = 2i * (inv_sqrt_pi - z1 * w1);
1827 dw2 = 2i * (inv_sqrt_pi - z2 * w2);
1828 K = 1 - (FVC - ETA * (Complex(G0, D0) - 3 * Complex(G2, D2) / 2)) * A +
1829 ETA * Complex(G2, D2) * B;
1830 F = inv_pi * A / K;
1831}
1832
1833VanVleckHuber::VanVleckHuber(Numeric F0, Numeric T) noexcept
1834 : c1(Constant::h / (2.0 * Constant::k * T)), tanh_c1f0(std::tanh(c1 * F0)),
1835 inv_denom(1.0 / (F0 * tanh_c1f0)) {}
1836
1837Numeric VanVleckHuber::dNdT(Numeric T, Numeric f) const noexcept {
1838 const Numeric part1 =
1839 (1 - pow2(tanh_c1f0)) * c1 * f * tanh_c1f / (T * pow2(tanh_c1f0));
1840 const Numeric part2 = (pow2(tanh_c1f) - 1) * c1 * pow2(f) * inv_denom / T;
1841 return part1 + part2;
1842}
1843
1844Numeric VanVleckHuber::dNdf(Numeric f) const noexcept {
1845 const Numeric part1 = tanh_c1f * inv_denom;
1846 const Numeric part2 = c1 * f * inv_denom / (1 - pow2(tanh_c1f));
1847 return part1 + part2;
1848}
1849
1850Numeric VanVleckHuber::dNdF0() const noexcept {
1851 const Numeric part1 = -N * tanh_c1f0 * inv_denom;
1852 const Numeric part2 = c1 * N * (pow2(tanh_c1f0) - 1) / tanh_c1f0;
1853 return part1 + part2;
1854}
1855
1856Numeric VanVleckHuber::operator()(Numeric f) noexcept {
1857 tanh_c1f = std::tanh(c1 * f);
1858 N = f * tanh_c1f * inv_denom;
1859 return N;
1860}
1861
1862RosenkranzQuadratic::RosenkranzQuadratic(Numeric F0, Numeric T) noexcept
1863 : fac((Constant::h) / (2.0 * Constant::k * T) /
1864 std::sinh((Constant::h * F0) / (2.0 * Constant::k * T)) * (1.0 / F0)),
1865 dfacdT((fac / T) *
1866 (fac * F0 * F0 *
1867 std::cosh((Constant::h * F0) / (2.0 * Constant::k * T)) -
1868 1)),
1869 dfacdF0(-fac * (1 / F0 + fac * F0 *
1870 std::cosh((Constant::h * F0) /
1871 (2.0 * Constant::k * T)))) {}
1872
1873Numeric RosenkranzQuadratic::dNdT(Numeric, Numeric f) const noexcept {
1874 return f * f * dfacdT;
1875}
1876
1877Numeric RosenkranzQuadratic::dNdf(Numeric f) const noexcept {
1878 return 2.0 * f * fac;
1879}
1880
1881Numeric RosenkranzQuadratic::dNdF0() const noexcept {
1882 return N * dfacdF0 / fac;
1883}
1884
1885Numeric RosenkranzQuadratic::operator()(Numeric f) noexcept {
1886 N = fac * f * f;
1887 return N;
1888}
1889
1890Numeric SimpleFrequencyScaling::dNdT(Numeric t_ [[maybe_unused]],
1891 Numeric f) const ARTS_NOEXCEPT {
1892 ARTS_ASSERT(t_ == T,
1893 "Temperature is stored internally in SimpleFrequencyScaling\n"
1894 "but for interface reasons this function nevertheless takes "
1895 "temprature as an input\n"
1896 "For some reason, the two temperatures disagree, so regardless, "
1897 "you have encountered\n"
1898 "a path of the code that should never be encountered. The two "
1899 "temperatures are: ",
1900 T, " and ", t_, " K")
1901
1902 return -N * Constant::h * F0 * expF0 / (Constant::k * t_ * t_ * expm1F0) +
1903 Constant::h * f * f *
1904 std::exp(-(Constant::h * f) / (Constant::k * t_)) /
1905 (F0 * Constant::k * t_ * t_ * expm1F0);
1906}
1907
1908Numeric SimpleFrequencyScaling::dNdf(Numeric f) const noexcept {
1909 return N / f - N * Constant::h *
1910 std::exp(-(Constant::h * f) / (Constant::k * T)) /
1911 (Constant::k * T *
1912 std::expm1(-(Constant::h * f) / (Constant::k * T)));
1913}
1914
1915Numeric SimpleFrequencyScaling::operator()(Numeric f) noexcept {
1916 N = (f / F0) * (std::expm1(-Constant::h * f / (Constant::k * T)) / expm1F0);
1917 return N;
1918}
1919
1920LocalThermodynamicEquilibrium::LocalThermodynamicEquilibrium(
1921 Numeric I0, Numeric T0, Numeric T, Numeric F0, Numeric E0, Numeric QT,
1922 Numeric QT0, Numeric dQTdT, Numeric r, Numeric drdSELFVMR,
1923 const Numeric drdT) noexcept
1924 : LocalThermodynamicEquilibrium(
1925 I0, r, drdSELFVMR, drdT, QT0, QT, dQTdT, boltzman_ratio(T, T0, E0),
1926 dboltzman_ratio_dT_div_boltzmann_ratio(T, E0),
1927 stimulated_relative_emission(F0, T0, T),
1928 dstimulated_relative_emission_dT(F0, T0, T),
1929 dstimulated_relative_emission_dF0(F0, T0, T)) {}
1930
1931struct FullNonLocalThermodynamicEquilibriumInitialization {
1932 Numeric k, dkdF0, dkdr1, dkdr2, e, dedF0, dedr2, B, dBdT, dBdF0;
1933
1934 FullNonLocalThermodynamicEquilibriumInitialization(Numeric F0, Numeric A21,
1935 Numeric T, Numeric r1,
1936 Numeric r2, Numeric c2,
1937 Numeric c3,
1938 Numeric x) noexcept
1939 : k(c3 * (r1 * x - r2) * (A21 / c2)), dkdF0(-2.0 * k / F0),
1940 dkdr1(c3 * x * (A21 / c2)), dkdr2(-c3 * (A21 / c2)), e(c3 * r2 * A21),
1941 dedF0(e / F0), dedr2(c3 * A21),
1942 B(2 * Constant::h / Math::pow2(Constant::c) * Math::pow3(F0) /
1943 std::expm1((Constant::h / Constant::k * F0) / T)),
1944 dBdT(Math::pow2(B) * Math::pow2(Constant::c) *
1945 std::exp((Constant::h / Constant::k * F0) / T) /
1946 (2 * Constant::k * Math::pow2(F0 * T))),
1947 dBdF0(3 * B / F0 - Math::pow2(B) * Math::pow2(Constant::c) *
1948 std::exp((Constant::h / Constant::k * F0) / T) /
1949 (2 * Constant::k * T * Math::pow3(F0))) {}
1950
1951 constexpr FullNonLocalThermodynamicEquilibrium
1952 operator()(Numeric r, Numeric drdSELFVMR, Numeric drdT) &&noexcept {
1953 return {r, drdSELFVMR, drdT, k, dkdF0, dkdr1, dkdr2,
1954 e, dedF0, dedr2, B, dBdT, dBdF0};
1955 }
1956};
1957
1958FullNonLocalThermodynamicEquilibrium::FullNonLocalThermodynamicEquilibrium(
1959 Numeric F0, Numeric A21, Numeric T, Numeric g1, Numeric g2, Numeric r1,
1960 Numeric r2, Numeric r, Numeric drdSELFVMR, Numeric drdT) noexcept
1961 : FullNonLocalThermodynamicEquilibrium(
1962 FullNonLocalThermodynamicEquilibriumInitialization(
1963 F0, A21, T, r1, r2, c0 * F0 * F0 * F0, c1 * F0,
1964 g2 / g1)(r, drdSELFVMR, drdT)) {}
1965
1966struct VibrationalTemperaturesNonLocalThermodynamicEquilibriumInitializer {
1967 Numeric K1, dK1dT, K2, dK2dT, dK2dF0, K3, dK3dT, dK3dF0, dK3dTl, dK3dTu, K4,
1968 dK4dT, dK4dTu, B, dBdT, dBdF0;
1969
1970 VibrationalTemperaturesNonLocalThermodynamicEquilibriumInitializer(
1971 Numeric T, Numeric T0, Numeric F0, Numeric E0, Numeric Tl, Numeric Evl,
1972 Numeric Tu, Numeric Evu, Numeric gamma, Numeric gamma_ref, Numeric r_low,
1973 Numeric r_upp) noexcept
1974 : K1(boltzman_ratio(T, T0, E0)), dK1dT(dboltzman_ratio_dT(K1, T, E0)),
1975 K2(stimulated_relative_emission(gamma, gamma_ref)),
1976 dK2dT(dstimulated_relative_emission_dT(gamma, gamma_ref, F0, T)),
1977 dK2dF0(dstimulated_relative_emission_dF0(gamma, gamma_ref, T, T0)),
1978 K3(absorption_nlte_ratio(gamma, r_upp, r_low)),
1979 dK3dT(dabsorption_nlte_rate_dT(gamma, T, F0, Evl, Evu, K4, r_low)),
1980 dK3dF0(dabsorption_nlte_rate_dF0(gamma, T, K4, r_low)),
1981 dK3dTl(dabsorption_nlte_rate_dTl(gamma, T, Tl, Evl, r_low)),
1982 dK3dTu(dabsorption_nlte_rate_dTu(gamma, T, Tu, Evu, K4)),
1983 K4(boltzman_ratio(Tu, T, Evu)), dK4dT(dboltzman_ratio_dT(K4, T, Evu)),
1984 dK4dTu(dboltzman_ratio_dT(K4, Tu, Evu)),
1985 B(2 * Constant::h / Math::pow2(Constant::c) * Math::pow3(F0) /
1986 std::expm1((Constant::h / Constant::k * F0) / T)),
1987 dBdT(Math::pow2(B) * Math::pow2(Constant::c) *
1988 std::exp((Constant::h / Constant::k * F0) / T) /
1989 (2 * Constant::k * Math::pow2(F0 * T))),
1990 dBdF0(3 * B / F0 - Math::pow2(B) * Math::pow2(Constant::c) *
1991 std::exp((Constant::h / Constant::k * F0) / T) /
1992 (2 * Constant::k * T * Math::pow3(F0))) {}
1993
1994 constexpr VibrationalTemperaturesNonLocalThermodynamicEquilibrium
1995 operator()(Numeric I0, Numeric QT0, Numeric QT, Numeric dQTdT, Numeric r,
1996 Numeric drdSELFVMR, Numeric drdT) &&noexcept {
1997 return {I0, QT0, QT, dQTdT, r, drdSELFVMR, drdT, K1,
1998 dK1dT, K2, dK2dT, dK2dF0, K3, dK3dT, dK3dF0, dK3dTl,
1999 dK3dTu, K4, dK4dT, dK4dTu, B, dBdT, dBdF0};
2000 }
2001};
2002
2003VibrationalTemperaturesNonLocalThermodynamicEquilibrium::
2004 VibrationalTemperaturesNonLocalThermodynamicEquilibrium(
2005 Numeric I0, Numeric T0, Numeric T, Numeric Tl, Numeric Tu, Numeric F0,
2006 Numeric E0, Numeric Evl, Numeric Evu, Numeric QT, Numeric QT0,
2007 Numeric dQTdT, Numeric r, Numeric drdSELFVMR, Numeric drdT) noexcept
2008 : VibrationalTemperaturesNonLocalThermodynamicEquilibrium(
2009 VibrationalTemperaturesNonLocalThermodynamicEquilibriumInitializer(
2010 T, T0, F0, E0, Tl, Evl, Tu, Evu, stimulated_emission(T, F0),
2011 stimulated_emission(T0, F0), boltzman_ratio(Tl, T, Evl),
2012 boltzman_ratio(Tu, T, Evu))(I0, QT0, QT, dQTdT, r, drdSELFVMR,
2013 drdT)) {}
2014
2016struct CutoffRange {
2017 Index start, size;
2018};
2019
2031CutoffRange limited_range(const Numeric fl, const Numeric fu,
2032 const Vector &f_grid) ARTS_NOEXCEPT {
2033 ARTS_ASSERT(fu > fl);
2034 const Numeric *it0 = f_grid.get_c_array();
2035 const Numeric *itn = it0 + f_grid.size();
2036 const Numeric *itl = std::lower_bound(it0, itn, fl);
2037 return CutoffRange{std::distance(it0, itl),
2038 std::distance(itl, std::upper_bound(itl, itn, fu))};
2039}
2040
2042struct SparseLimitRange {
2043 Index start, size;
2044 Index start_low, size_low;
2045 Index start_upp, size_upp;
2046};
2047
2048SparseLimitRange linear_sparse_limited_range(
2049 const Numeric flc, const Numeric fuc, const Numeric fls, const Numeric fus,
2050 const Vector &f_grid, const Vector &sparse_f_grid) ARTS_NOEXCEPT {
2051 ARTS_ASSERT(fls > flc);
2052 ARTS_ASSERT(fus > fls);
2053 ARTS_ASSERT(fuc > fus);
2054
2055 const Index nvs = sparse_f_grid.size();
2056
2057 // Find bounds in sparse
2058 const Numeric *it0s = sparse_f_grid.get_c_array();
2059 const Numeric *itns = it0s + nvs;
2060 const Numeric *itlc =
2061 std::lower_bound(it0s, itns, std::nextafter(flc, fuc)); // lower cutoff
2062 const Numeric *ituc =
2063 std::upper_bound(itlc, itns, std::nextafter(fuc, flc)); // upper cutoff
2064 const Numeric *itls = std::upper_bound(itlc, ituc, fls); // lower sparse
2065 const Numeric *itus = std::lower_bound(itls, ituc, fus); // upper sparse
2066
2067 /* Start and size of sparse adjusted to the 2-grid
2068 *
2069 * The interface to the dense grid is altered slightly so that
2070 * a complete set-of-2 quadratic interopolation points becomes
2071 * a guarantee. Their start/end points ignore this restriction
2072 * but have been defined to not contain anything extra
2073 */
2074 Index const beg_lr =
2075 std::distance(it0s, itlc); /*while (beg_lr % 2) --beg_lr;*/
2076 Index end_lr = std::distance(it0s, itls);
2077 while (end_lr % 2)
2078 --end_lr;
2079 Index beg_ur = std::distance(it0s, itus);
2080 while (beg_ur % 2)
2081 ++beg_ur;
2082 Index const end_ur =
2083 std::distance(it0s, ituc); /*while (end_ur % 2) ++end_ur;*/
2084
2085 // Find new limits
2086 const Numeric fl =
2087 (end_lr <= 0 or end_lr >= nvs) ? flc : sparse_f_grid[end_lr];
2088 const Numeric fu =
2089 (beg_ur <= 0 or beg_ur >= nvs) ? fuc : sparse_f_grid[beg_ur];
2090
2091 // Find bounds in dense
2092 const Numeric *it0 = f_grid.get_c_array();
2093 const Numeric *itn = it0 + f_grid.size();
2094 const Numeric *itl = std::lower_bound(it0, itn, fl); // include boundary
2095 const Numeric *itu =
2096 std::upper_bound(itl, itn, std::nextafter(fu, fl)); // dismiss boundary
2097
2098 return {std::distance(it0, itl),
2099 std::distance(itl, itu),
2100 beg_lr,
2101 end_lr - beg_lr,
2102 beg_ur,
2103 end_ur - beg_ur};
2104}
2105
2106SparseLimitRange quad_sparse_limited_range(
2107 const Numeric flc, const Numeric fuc, const Numeric fls, const Numeric fus,
2108 const Vector &f_grid, const Vector &sparse_f_grid) ARTS_NOEXCEPT {
2109 ARTS_ASSERT(fls > flc);
2110 ARTS_ASSERT(fus > fls);
2111 ARTS_ASSERT(fuc > fus);
2112
2113 const Index nvs = sparse_f_grid.size();
2114 const Index nv = f_grid.size();
2115
2116 // Find bounds in sparse
2117 const Numeric *const it0s = sparse_f_grid.get_c_array();
2118 const Numeric *const itns = it0s + nvs;
2119 const Numeric *itlc =
2120 std::lower_bound(it0s, itns, std::nextafter(flc, fuc)); // lower cutoff
2121 const Numeric *ituc =
2122 std::upper_bound(itlc, itns, std::nextafter(fuc, flc)); // upper cutoff
2123 const Numeric *itls =
2124 std::upper_bound(itlc, ituc, std::nextafter(fls, flc)); // lower sparse
2125 const Numeric *itus =
2126 std::lower_bound(itls, ituc, std::nextafter(fus, fuc)); // upper sparse
2127
2128 /* Start and size of sparse adjusted to the 3-grid
2129 *
2130 * The interface to the dense grid is altered slightly so that
2131 * a complete set-of-3 quadratic interpolation points becomes
2132 * a guarantee. Their end points are modified to contain
2133 * only set-of-3. If the range is not modified correctly,
2134 * you risk having negative interpolation out of your range.
2135 *
2136 * To avoid negative absorption entirely, the range between
2137 * the true cutoff and the outer-most sparse grid point must
2138 * have at least two values. If they have only one value,
2139 * the quadratic interpolation will lead to negative values
2140 * in the outer range. So there's a small range that has to
2141 * be ignored
2142 */
2143 while (std::distance(it0s, itls) % 3)
2144 --itls;
2145 while (std::distance(itlc, itls) % 3 == 1)
2146 ++itlc; // skip some cutoff
2147 while (std::distance(it0s, itus) % 3)
2148 ++itus;
2149 while (std::distance(itus, ituc) % 3 == 1)
2150 --ituc; // skip some cutoff
2151
2152 // Find bounds in dense
2153 const Numeric *const it0 = f_grid.get_c_array();
2154 const Numeric *const itn = it0 + nv;
2155 const Numeric *itl;
2156 const Numeric *itu;
2157
2158 if (itls not_eq itns) {
2159 itl = std::lower_bound(it0, itn, *itls); // include boundary
2160 } else {
2161 itl = std::lower_bound(it0, itn, flc); // include boundary
2162 }
2163
2164 if (itus not_eq itns and itl not_eq itn) {
2165 itu = std::upper_bound(itl, itn,
2166 std::nextafter(*itus, *itl)); // dismiss boundary
2167 } else {
2168 itu = std::lower_bound(itl, itn,
2169 std::nextafter(fuc, flc)); // include boundary
2170 }
2171
2172 return {std::distance(it0, itl), std::distance(itl, itu),
2173 std::distance(it0s, itlc), std::distance(itlc, itls),
2174 std::distance(it0s, itus), std::distance(itus, ituc)};
2175}
2176
2182struct Derivatives {
2183 Quantum::Number::StateMatch target;
2184 union Values {
2185 Output o;
2186 Numeric n;
2187 constexpr Values() noexcept : n() {}
2188 } value;
2189 Index jac_pos{};
2190 const RetrievalQuantity *deriv{nullptr};
2191 Derivatives() noexcept : target(), value() {}
2192};
2193
2195using ArrayOfDerivatives = Array<Derivatives>;
2196
2198Index active_nelem(const ArrayOfDerivatives &derivs) noexcept {
2199 return std::distance(
2200 derivs.cbegin(),
2201 std::find_if(derivs.cbegin(), derivs.cend(),
2202 [](auto &deriv) { return deriv.deriv == nullptr; }));
2203}
2204
2205struct ComputeValues {
2206 Complex *const F;
2207 Complex *const dF;
2208 Complex *const N;
2209 Complex *const dN;
2210 const Numeric *const f;
2211 const Index size;
2212 const ArrayOfDerivatives &derivs;
2213 const Index jac_size;
2214 const Index max_jac_size;
2215 const bool do_nlte;
2216
2217 [[nodiscard]] Index jac_pos(Index iv, Index ij) const noexcept {
2218 return jac_size * iv + ij;
2219 }
2220
2221 ComputeValues(ComplexVector &F_, ComplexMatrix &dF_, ComplexVector &N_,
2222 ComplexMatrix &dN_, const Vector &f_grid, const Index start,
2223 const Index nv, const ArrayOfDerivatives &derivs_,
2224 const bool do_nlte_) noexcept
2225 : F(F_.get_c_array() + start),
2226 dF(dF_.get_c_array() + start * derivs_.size()),
2227 N(N_.get_c_array() + start),
2228 dN(dN_.get_c_array() + start * derivs_.size()),
2229 f(f_grid.get_c_array() + start), size(nv), derivs(derivs_),
2230 jac_size(derivs_.size()), max_jac_size(active_nelem(derivs)),
2231 do_nlte(do_nlte_) {}
2232
2233 ComputeValues(Complex &F_, std::vector<Complex> &dF_, Complex &N_,
2234 std::vector<Complex> &dN_, const Numeric &f_lim,
2235 const ArrayOfDerivatives &derivs_, const bool do_nlte_) noexcept
2236 : F(&F_), dF(dF_.data()), N(&N_), dN(dN_.data()), f(&f_lim), size(1),
2237 derivs(derivs_), jac_size(derivs.nelem()),
2238 max_jac_size(active_nelem(derivs)), do_nlte(do_nlte_) {}
2239
2240 ComputeValues &operator-=(const ComputeValues &cut) ARTS_NOEXCEPT {
2241 ARTS_ASSERT(cut.size == 1, "Not a cutoff limit")
2242 ARTS_ASSERT(cut.jac_size == jac_size, "Not from the same Jacobian type")
2243
2244 for (Index iv = 0; iv < size; iv++) {
2245 F[iv] -= *cut.F;
2246 for (Index ij = 0; ij < max_jac_size; ij++) {
2247 dF[jac_pos(iv, derivs[ij].jac_pos)] -= cut.dF[derivs[ij].jac_pos];
2248 }
2249 }
2250 if (do_nlte) {
2251 for (Index iv = 0; iv < size; iv++) {
2252 N[iv] -= *cut.N;
2253 for (Index ij = 0; ij < max_jac_size; ij++) {
2254 dN[jac_pos(iv, derivs[ij].jac_pos)] -= cut.dN[derivs[ij].jac_pos];
2255 }
2256 }
2257 }
2258 return *this;
2259 }
2260};
2261
2262Complex Calculator::dFdT(const Output &dXdT, Numeric T) const noexcept {
2263 return std::visit([&](auto &&LS) { return LS.dFdT(dXdT, T); }, ls);
2264}
2265
2266Complex Calculator::dFdf() const noexcept {
2267 return std::visit([](auto &&LS) { return LS.dFdf(); }, ls);
2268}
2269
2270Complex Calculator::dFdF0() const noexcept {
2271 return std::visit([](auto &&LS) { return LS.dFdF0(); }, ls);
2272}
2273
2274Complex Calculator::dFdH(Numeric dfdH) const noexcept {
2275 return std::visit([dfdH](auto &&LS) { return LS.dFdH(dfdH); }, ls);
2276}
2277
2278Complex Calculator::dFdVMR(const Output &dXdVMR) const noexcept {
2279 return std::visit([&](auto &&LS) { return LS.dFdVMR(dXdVMR); }, ls);
2280}
2281
2282Complex Calculator::dFdFVC(Numeric d) const noexcept {
2283 return std::visit([d](auto &&LS) { return LS.dFdFVC(d); }, ls);
2284}
2285
2286Complex Calculator::dFdETA(Numeric d) const noexcept {
2287 return std::visit([d](auto &&LS) { return LS.dFdETA(d); }, ls);
2288}
2289
2290Complex Calculator::dFdDV(Numeric d) const noexcept {
2291 return std::visit([d](auto &&LS) { return LS.dFdDV(d); }, ls);
2292}
2293
2294Complex Calculator::dFdD0(Numeric d) const noexcept {
2295 return std::visit([d](auto &&LS) { return LS.dFdD0(d); }, ls);
2296}
2297
2298Complex Calculator::dFdG0(Numeric d) const noexcept {
2299 return std::visit([d](auto &&LS) { return LS.dFdG0(d); }, ls);
2300}
2301
2302Complex Calculator::dFdD2(Numeric d) const noexcept {
2303 return std::visit([d](auto &&LS) { return LS.dFdD2(d); }, ls);
2304}
2305
2306Complex Calculator::dFdG2(Numeric d) const noexcept {
2307 return std::visit([d](auto &&LS) { return LS.dFdG2(d); }, ls);
2308}
2309
2310Complex Calculator::F() const noexcept {
2311 return std::visit([](auto &&LS) { return LS.F; }, ls);
2312}
2313
2314Complex Calculator::operator()(Numeric f) noexcept {
2315 return std::visit([f](auto &&LS) { return LS(f); }, ls);
2316}
2317
2318Calculator::Calculator(const Type type, const Numeric F0, const Output &X,
2319 const Numeric DC, const Numeric DZ,
2320 bool manually_mirrored) noexcept
2321 : ls(Noshape{}) {
2322 if (not manually_mirrored) {
2323 switch (type) {
2324 case Type::DP:
2325 ls = Doppler(F0, DC, DZ);
2326 break;
2327 case Type::LP: [[fallthrough]];
2328 case Type::SplitLP:
2329 ls = Lorentz(F0, X);
2330 break;
2331 case Type::VP: [[fallthrough]];
2332 case Type::SplitVP:
2333 ls = Voigt(F0, X, DC, DZ);
2334 break;
2335 case Type::SDVP: [[fallthrough]];
2336 case Type::SplitSDVP:
2337 ls = SpeedDependentVoigt(F0, X, DC, DZ);
2338 break;
2339 case Type::HTP: [[fallthrough]];
2340 case Type::SplitHTP:
2341 ls = HartmannTran(F0, X, DC, DZ);
2342 break;
2343 case Type::FINAL: { /*leave last*/
2344 }
2345 }
2346 }
2347}
2348
2349Calculator::Calculator(const Absorption::MirroringType mirror, const Type type,
2350 const Numeric F0, const Output &X, const Numeric DC,
2351 const Numeric DZ)
2352 : ls(Noshape{}) {
2353 switch (mirror) {
2354 case Absorption::MirroringType::Lorentz:
2355 ls = Lorentz(-F0, mirroredOutput(X));
2356 break;
2357 case Absorption::MirroringType::SameAsLineShape:
2358 *this = {type, -F0, mirroredOutput(X), -DC, -DZ, false};
2359 break;
2360 case Absorption::MirroringType::Manual:
2361 *this = {type, F0, mirroredOutput(X), -DC, -DZ, false};
2362 break;
2363 case Absorption::MirroringType::None:
2364 break;
2365 case Absorption::MirroringType::FINAL: { /*leave last*/
2366 }
2367 }
2368}
2369
2370Numeric Normalizer::dNdT(Numeric T, Numeric f) const noexcept {
2371 return std::visit([T, f](auto &&LSN) { return LSN.dNdT(T, f); }, ls_norm);
2372}
2373
2374Numeric Normalizer::dNdf(Numeric f) const noexcept {
2375 return std::visit([f](auto &&LS) { return LS.dNdf(f); }, ls_norm);
2376}
2377
2378Numeric Normalizer::dNdF0() const noexcept {
2379 return std::visit([](auto &&LS) { return LS.dNdF0(); }, ls_norm);
2380}
2381
2382Numeric Normalizer::operator()(Numeric f) noexcept {
2383 return std::visit([f](auto &&LSN) { return LSN(f); }, ls_norm);
2384}
2385
2386Normalizer::Normalizer(const Absorption::NormalizationType type,
2387 const Numeric F0, const Numeric T) noexcept
2388 : ls_norm(Nonorm{}) {
2389 switch (type) {
2390 case Absorption::NormalizationType::None:
2391 break;
2392 case Absorption::NormalizationType::RQ:
2393 ls_norm = RosenkranzQuadratic(F0, T);
2394 break;
2395 case Absorption::NormalizationType::VVH:
2396 ls_norm = VanVleckHuber(F0, T);
2397 break;
2398 case Absorption::NormalizationType::VVW:
2399 ls_norm = VanVleckWeisskopf(F0);
2400 break;
2401 case Absorption::NormalizationType::SFS:
2402 ls_norm = SimpleFrequencyScaling(F0, T);
2403 break;
2404 case Absorption::NormalizationType::FINAL: { /*leave last*/
2405 }
2406 }
2407}
2408
2409Numeric IntensityCalculator::S() const noexcept {
2410 return scale * std::visit([](auto &&S) { return S.S; }, ls_str);
2411}
2412
2413Numeric IntensityCalculator::dSdT() const noexcept {
2414 return scale * std::visit([](auto &&LSN) { return LSN.dSdT(); }, ls_str);
2415}
2416
2417Numeric IntensityCalculator::dSdI0() const noexcept {
2418 return scale * std::visit([](auto &&LS) { return LS.dSdI0(); }, ls_str);
2419}
2420
2421Numeric IntensityCalculator::dSdF0() const noexcept {
2422 return scale * std::visit([](auto &&LS) { return LS.dSdF0(); }, ls_str);
2423}
2424
2425Numeric IntensityCalculator::dSdNLTEu() const noexcept {
2426 return scale * std::visit([](auto &&LS) { return LS.dSdNLTEu(); }, ls_str);
2427}
2428
2429Numeric IntensityCalculator::dSdNLTEl() const noexcept {
2430 return scale * std::visit([](auto &&LS) { return LS.dSdNLTEl(); }, ls_str);
2431}
2432
2433Numeric IntensityCalculator::dSdSELFVMR() const noexcept {
2434 return (self_species == scaling_species
2435 ? std::visit([](auto &&S) { return S.S; }, ls_str)
2436 : 0) +
2437 scale * std::visit([](auto &&LS) { return LS.dSdSELFVMR(); }, ls_str);
2438}
2439
2440Numeric IntensityCalculator::dSdOTHERVMR_if() const noexcept {
2441 return self_species == scaling_species
2442 ? 0
2443 : std::visit([](auto &&S) { return S.S; }, ls_str);
2444}
2445
2446Numeric IntensityCalculator::N() const noexcept {
2447 return scale * std::visit([](auto &&S) { return S.N; }, ls_str);
2448}
2449
2450Numeric IntensityCalculator::dNdT() const noexcept {
2451 return scale * std::visit([](auto &&LSN) { return LSN.dNdT(); }, ls_str);
2452}
2453
2454Numeric IntensityCalculator::dNdI0() const noexcept {
2455 return scale * std::visit([](auto &&LS) { return LS.dNdI0(); }, ls_str);
2456}
2457
2458Numeric IntensityCalculator::dNdF0() const noexcept {
2459 return scale * std::visit([](auto &&LS) { return LS.dNdF0(); }, ls_str);
2460}
2461
2462Numeric IntensityCalculator::dNdNLTEu() const noexcept {
2463 return scale * std::visit([](auto &&LS) { return LS.dNdNLTEu(); }, ls_str);
2464}
2465
2466Numeric IntensityCalculator::dNdNLTEl() const noexcept {
2467 return scale * std::visit([](auto &&LS) { return LS.dNdNLTEl(); }, ls_str);
2468}
2469
2470Numeric IntensityCalculator::dNdSELFVMR() const noexcept {
2471 return (self_species == scaling_species
2472 ? std::visit([](auto &&S) { return S.N; }, ls_str)
2473 : 0) +
2474 scale * std::visit([](auto &&LS) { return LS.dNdSELFVMR(); }, ls_str);
2475}
2476
2477Numeric IntensityCalculator::dNdOTHERVMR_if() const noexcept {
2478 return self_species == scaling_species
2479 ? 0
2480 : std::visit([](auto &&S) { return S.N; }, ls_str);
2481}
2482
2483IntensityCalculator &
2484IntensityCalculator::adaptive_scaling(Numeric x, Species::Species self,
2485 Species::Species other) noexcept {
2486 scale = x;
2487 self_species = self;
2488 scaling_species = other;
2489 return *this;
2490}
2491
2492IntensityCalculator::IntensityCalculator(
2493 const Numeric T, const Numeric QT, const Numeric QT0, const Numeric dQTdT,
2494 const Numeric r, const Numeric drdSELFVMR, const Numeric drdT,
2495 const EnergyLevelMap &nlte, const Absorption::Lines &band,
2496 const Index line_index) noexcept
2497 : ls_str(Nostrength{}) {
2498 const auto &line = band.lines[line_index];
2499 switch (band.population) {
2500 case Absorption::PopulationType::ByHITRANFullRelmat:
2501 case Absorption::PopulationType::ByHITRANRosenkranzRelmat:
2502 case Absorption::PopulationType::ByMakarovFullRelmat:
2503 case Absorption::PopulationType::ByRovibLinearDipoleLineMixing:
2504 case Absorption::PopulationType::LTE:
2505 ls_str =
2506 LocalThermodynamicEquilibrium(line.I0, band.T0, T, line.F0, line.E0, QT,
2507 QT0, dQTdT, r, drdSELFVMR, drdT);
2508 break;
2509 case Absorption::PopulationType::NLTE: {
2510 const auto [r_low, r_upp] = nlte.get_ratio_params(band, line_index);
2511 ls_str = FullNonLocalThermodynamicEquilibrium(line.F0, line.A, T, line.glow,
2512 line.gupp, r_low, r_upp, r,
2513 drdSELFVMR, drdT);
2514 } break;
2515 case Absorption::PopulationType::VibTemps: {
2516 const auto [E_low, E_upp, T_low, T_upp] = nlte.get_vibtemp_params(band, T);
2517 ls_str = VibrationalTemperaturesNonLocalThermodynamicEquilibrium(
2518 line.I0, band.T0, T, T_low, T_upp, line.F0, line.E0, E_low, E_upp, QT,
2519 QT0, dQTdT, r, drdSELFVMR, drdT);
2520 } break;
2521 case Absorption::PopulationType::FINAL: { /*leave last*/
2522 }
2523 }
2524}
2525
2526#define CutInternalDerivativesImpl(X, Y) \
2527 else if (deriv == Jacobian::Line::Shape##X##Y) { \
2528 const Numeric d = value.n; \
2529 const Complex dFm = std::conj(ls_mirr.dFd##X(d) - ls_mirr_cut.dFd##X(d)); \
2530 const Complex dFls = ls.dFd##X(d) - ls_cut.dFd##X(d) + dFm; \
2531 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls; \
2532 if (do_nlte) { \
2533 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls; \
2534 } \
2535 }
2536
2537#define CutInternalDerivatives(X) \
2538 CutInternalDerivativesImpl(X, X0) CutInternalDerivativesImpl(X, X1) \
2539 CutInternalDerivativesImpl(X, X2) CutInternalDerivativesImpl(X, X3)
2540
2541#define InternalDerivativesImpl(X, Y) \
2542 else if (deriv == Jacobian::Line::Shape##X##Y) { \
2543 const Numeric d = value.n; \
2544 const Complex dFm = std::conj(ls_mirr.dFd##X(d)); \
2545 const Complex dFls = ls.dFd##X(d) + dFm; \
2546 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls; \
2547 if (do_nlte) { \
2548 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls; \
2549 } \
2550 }
2551
2552#define InternalDerivatives(X) \
2553 InternalDerivativesImpl(X, X0) InternalDerivativesImpl(X, X1) \
2554 InternalDerivativesImpl(X, X2) InternalDerivativesImpl(X, X3)
2555
2556#define InternalDerivativesGImpl(X) \
2557 else if (deriv == Jacobian::Line::ShapeG##X) { \
2558 const Numeric dLM = value.n; \
2559 com.dF[com.jac_pos(iv, ij)] += dLM * S * Fls; \
2560 if (do_nlte) { \
2561 com.dN[com.jac_pos(iv, ij)] += dLM * DS * Fls; \
2562 } \
2563 }
2564
2565#define InternalDerivativesG \
2566 InternalDerivativesGImpl(X0) InternalDerivativesGImpl(X1) \
2567 InternalDerivativesGImpl(X2) InternalDerivativesGImpl(X3)
2568
2569#define InternalDerivativesYImpl(X) \
2570 else if (deriv == Jacobian::Line::ShapeY##X) { \
2571 const Complex dLM = Complex(0, -value.n); \
2572 com.dF[com.jac_pos(iv, ij)] += dLM * S * Fls; \
2573 if (do_nlte) { \
2574 com.dN[com.jac_pos(iv, ij)] += dLM * DS * Fls; \
2575 } \
2576 }
2577
2578#define InternalDerivativesY \
2579 InternalDerivativesYImpl(X0) InternalDerivativesYImpl(X1) \
2580 InternalDerivativesYImpl(X2) InternalDerivativesYImpl(X3)
2581
2582#define InternalDerivativesSetupImpl(X, Y, ...) \
2583 else if (deriv == Jacobian::Line::Shape##X##Y) { \
2584 const Index pos = \
2585 band.BroadeningSpeciesPosition(deriv.Target().species_id); \
2586 if (pos >= 0 __VA_OPT__(and pos == __VA_ARGS__)) { \
2587 derivs[ij].value.n = (__VA_OPT__(1 ? 1 :) vmrs[pos]) * \
2588 band.lines[i].lineshape[pos].dX( \
2589 T, band.T0, P, Jacobian::Line::Shape##X##Y); \
2590 } else { \
2591 derivs[ij].value.n = 0; \
2592 } \
2593 }
2594
2595#define InternalDerivativesSetup(X, ...) \
2596 InternalDerivativesSetupImpl(X, X0 __VA_OPT__(, __VA_ARGS__)) \
2597 InternalDerivativesSetupImpl(X, X1 __VA_OPT__(, __VA_ARGS__)) \
2598 InternalDerivativesSetupImpl(X, X2 __VA_OPT__(, __VA_ARGS__)) \
2599 InternalDerivativesSetupImpl(X, X3 __VA_OPT__(, __VA_ARGS__))
2600
2663void cutoff_frequency_loop(ComputeValues &com, Calculator &ls,
2664 Calculator &ls_mirr, Normalizer &ls_norm,
2665 const Calculator &ls_cut,
2666 const Calculator &ls_mirr_cut,
2667 const IntensityCalculator &ls_str,
2668 const ArrayOfDerivatives &derivs, const Complex LM,
2669 const Numeric &T, const Numeric &dfdH,
2670 const Numeric &Sz,
2671 const Species::Species self_species) ARTS_NOEXCEPT {
2672 const Index nv = com.size;
2673 const bool do_nlte = com.do_nlte;
2674
2675 const Numeric Si = ls_str.S();
2676 const Numeric DNi = ls_str.N();
2677
2678 for (Index iv = 0; iv < nv; iv++) {
2679 const Numeric f = com.f[iv];
2680
2681 const Numeric Sn = ls_norm(f);
2682 const Numeric S = Sz * Sn * Si;
2683 const Numeric DS = Sz * Sn * DNi;
2684 const Complex Fm = std::conj(ls_mirr(f) - ls_mirr_cut.F());
2685 const Complex Fls = ls(f) - ls_cut.F() + Fm;
2686 com.F[iv] += S * LM * Fls;
2687 if (do_nlte) {
2688 com.N[iv] += DS * LM * Fls;
2689 }
2690
2691 for (const auto &[lt, value, ij, deriv_ptr] : derivs) {
2692 if (not deriv_ptr)
2693 break;
2694 const auto &deriv = *deriv_ptr;
2695
2696 if (deriv == Jacobian::Atm::Temperature) {
2697 const auto &dXdT = value.o;
2698 const Numeric dSn = ls_norm.dNdT(T, f);
2699 const Numeric dSi = ls_str.dSdT();
2700 const Complex dLM(dXdT.G, -dXdT.Y);
2701 const Complex dFm =
2702 std::conj(ls_mirr.dFdT(mirroredOutput(dXdT), T) -
2703 ls_mirr_cut.dFdT(mirroredOutput(dXdT), T));
2704 const Complex dFls = ls.dFdT(dXdT, T) - ls_cut.dFdT(dXdT, T) + dFm;
2705 com.dF[com.jac_pos(iv, ij)] +=
2706 S * (LM * dFls + dLM * Fls) + Sz * (dSn * Si + Sn * dSi) * LM * Fls;
2707 if (do_nlte) {
2708 const Numeric dDSi = ls_str.dNdT();
2709 com.dN[com.jac_pos(iv, ij)] += DS * (LM * dFls + dLM * Fls) +
2710 Sz * (dSn * Si + Sn * dDSi) * LM * Fls;
2711 }
2712 } else if (deriv.is_wind()) {
2713 const Complex dFm = std::conj(ls_mirr.dFdf() - ls_mirr_cut.dFdf());
2714 const Complex dFls = ls.dFdf() - ls_cut.dFdf() + dFm;
2715 const Numeric dS = Sz * ls_norm.dNdf(f) * Si;
2716 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dS * LM * Fls;
2717 if (do_nlte) {
2718 const Numeric dDS = Sz * ls_norm.dNdf(f) * DNi;
2719 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + dDS * LM * Fls;
2720 }
2721 } else if (deriv.is_mag()) {
2722 const Complex dFm =
2723 std::conj(ls_mirr.dFdH(-dfdH) - ls_mirr_cut.dFdH(-dfdH));
2724 const Complex dFls = ls.dFdH(dfdH) - ls_cut.dFdH(dfdH) + dFm;
2725 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls;
2726 if (do_nlte) {
2727 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls;
2728 }
2729 } else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR) {
2730 com.dF[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dSdSELFVMR() * LM * Fls;
2731 if (do_nlte) {
2732 com.dN[com.jac_pos(iv, ij)] +=
2733 Sz * Sn * ls_str.dNdSELFVMR() * LM * Fls;
2734 }
2735 } else if (deriv.Target().needQuantumIdentity()) {
2736 if (deriv == Jacobian::Line::VMR) {
2737 const auto &dXdVMR = value.o;
2738 const Complex dFm =
2739 std::conj(ls_mirr.dFdVMR(mirroredOutput(dXdVMR)) -
2740 ls_mirr_cut.dFdVMR(mirroredOutput(dXdVMR)));
2741 const Complex dLM(dXdVMR.G, -dXdVMR.Y);
2742 const Complex dFls = ls.dFdVMR(dXdVMR) - ls_cut.dFdVMR(dXdVMR) + dFm;
2743 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dLM * S * Fls;
2744
2745 if (self_species == deriv.Target().species_id) {
2746 com.dF[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dSdSELFVMR() * LM * Fls;
2747 }
2748
2749 if (ls_str.scaler() == deriv.Target().species_id) {
2750 com.dF[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dSdOTHERVMR_if() * LM * Fls;
2751 }
2752
2753 if (ls_str.scaler() == Species::Species::Bath) {
2754 com.dF[com.jac_pos(iv, ij)] -= Sz * Sn * ls_str.dSdOTHERVMR_if() * LM * Fls;
2755 }
2756
2757 if (do_nlte) {
2758 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + dLM * DS * Fls;
2759
2760 if (self_species == deriv.QuantumIdentity().Species()) {
2761 com.dN[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dNdSELFVMR() * LM * Fls;
2762 }
2763
2764 if (ls_str.scaler() == deriv.QuantumIdentity().Species()) {
2765 com.dN[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dNdOTHERVMR_if() * LM * Fls;
2766 }
2767
2768 if (ls_str.scaler() == Species::Species::Bath) {
2769 com.dN[com.jac_pos(iv, ij)] -= Sz * Sn * ls_str.dNdOTHERVMR_if() * LM * Fls;
2770 }
2771 }
2772 } else {
2773 if (lt == Quantum::Number::StateMatchType::Full) {
2774 if (deriv == Jacobian::Line::Center) {
2775 const Numeric d = ls_norm.dNdF0() * Si + ls_str.dSdF0() * Sn;
2776 const Complex dFm =
2777 std::conj(ls_mirr.dFdF0() - ls_mirr_cut.dFdF0());
2778 const Complex dFls = ls.dFdF0() - ls_cut.dFdF0() + dFm;
2779 const Numeric dS = Sz * d;
2780 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dS * LM * Fls;
2781 if (do_nlte) {
2782 const Numeric Dd = ls_norm.dNdF0() * DNi + ls_str.dNdF0() * Sn;
2783 const Numeric DdS = Sz * Dd;
2784 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + DdS * LM * Fls;
2785 }
2786 } else if (deriv == Jacobian::Line::Strength) {
2787 const Numeric dS = Sz * Sn * ls_str.dSdI0();
2788 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
2789 if (do_nlte) {
2790 const Numeric DdS = Sz * Sn * ls_str.dNdI0();
2791 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
2792 }
2793 }
2794 // All line shape derivatives
2795 CutInternalDerivatives(G0) CutInternalDerivatives(D0)
2796 CutInternalDerivatives(G2) CutInternalDerivatives(D2)
2797 CutInternalDerivatives(ETA) CutInternalDerivatives(FVC)
2798 InternalDerivativesY InternalDerivativesG
2799 CutInternalDerivatives(DV)
2800 } else if (lt == Quantum::Number::StateMatchType::Level) {
2801 if (deriv == Jacobian::Line::NLTE) {
2802 if (lt.upp) {
2803 const Numeric dS = Sz * Sn * ls_str.dSdNLTEu();
2804 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
2805
2806 if (do_nlte) {
2807 const Numeric DdS = Sz * Sn * ls_str.dNdNLTEu();
2808 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
2809 }
2810 }
2811
2812 if (lt.low) {
2813 const Numeric dS = Sz * Sn * ls_str.dSdNLTEl();
2814 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
2815
2816 if (do_nlte) {
2817 const Numeric DdS = Sz * Sn * ls_str.dNdNLTEl();
2818 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
2819 }
2820 }
2821 }
2822 }
2823 }
2824 }
2825 }
2826 }
2827}
2828
2890void frequency_loop(ComputeValues &com, Calculator &ls, Calculator &ls_mirr,
2891 Normalizer &ls_norm, const IntensityCalculator &ls_str,
2892 const ArrayOfDerivatives &derivs, const Complex LM,
2893 const Numeric &T, const Numeric &dfdH, const Numeric &Sz,
2894 const Species::Species self_species) ARTS_NOEXCEPT {
2895 const Index nv = com.size;
2896 const bool do_nlte = com.do_nlte;
2897
2898 const Numeric Si = ls_str.S();
2899 const Numeric DNi = ls_str.N();
2900
2901 for (Index iv = 0; iv < nv; iv++) {
2902 const Numeric f = com.f[iv];
2903
2904 const Numeric Sn = ls_norm(f);
2905 const Numeric S = Sz * Sn * Si;
2906 const Numeric DS = Sz * Sn * DNi;
2907 const Complex Fm = std::conj(ls_mirr(f));
2908 const Complex Fls = ls(f) + Fm;
2909 com.F[iv] += S * LM * Fls;
2910 if (do_nlte) {
2911 com.N[iv] += DS * LM * Fls;
2912 }
2913
2914 for (const auto &[lt, value, ij, deriv_ptr] : derivs) {
2915 if (not deriv_ptr)
2916 break;
2917 const auto &deriv = *deriv_ptr;
2918
2919 if (deriv == Jacobian::Atm::Temperature) {
2920 const auto &dXdT = value.o;
2921 const Numeric dSn = ls_norm.dNdT(T, f);
2922 const Numeric dSi = ls_str.dSdT();
2923 const Complex dLM(dXdT.G, -dXdT.Y);
2924 const Complex dFm = std::conj(ls_mirr.dFdT(mirroredOutput(dXdT), T));
2925 const Complex dFls = ls.dFdT(dXdT, T) + dFm;
2926 com.dF[com.jac_pos(iv, ij)] +=
2927 S * (LM * dFls + dLM * Fls) + Sz * (dSn * Si + Sn * dSi) * LM * Fls;
2928 if (do_nlte) {
2929 const Numeric dDSi = ls_str.dNdT();
2930 com.dN[com.jac_pos(iv, ij)] += DS * (LM * dFls + dLM * Fls) +
2931 Sz * (dSn * Si + Sn * dDSi) * LM * Fls;
2932 }
2933 } else if (deriv.is_wind()) {
2934 const Complex dFm = std::conj(ls_mirr.dFdf());
2935 const Complex dFls = ls.dFdf() + dFm;
2936 const Numeric dS = Sz * ls_norm.dNdf(f) * Si;
2937 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dS * LM * Fls;
2938 if (do_nlte) {
2939 const Numeric dDS = Sz * ls_norm.dNdf(f) * DNi;
2940 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + dDS * LM * Fls;
2941 }
2942 } else if (deriv.is_mag()) {
2943 const Complex dFm = std::conj(ls_mirr.dFdH(-dfdH));
2944 const Complex dFls = ls.dFdH(dfdH) + dFm;
2945 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls;
2946 if (do_nlte) {
2947 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls;
2948 }
2949 } else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR) {
2950 com.dF[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dSdSELFVMR() * LM * Fls;
2951 if (do_nlte) {
2952 com.dN[com.jac_pos(iv, ij)] +=
2953 Sz * Sn * ls_str.dNdSELFVMR() * LM * Fls;
2954 }
2955 } else if (deriv.Target().needQuantumIdentity()) {
2956 if (deriv == Jacobian::Line::VMR) {
2957 const auto &dXdVMR = value.o;
2958 const Complex dFm = std::conj(ls_mirr.dFdVMR(mirroredOutput(dXdVMR)));
2959 const Complex dLM(dXdVMR.G, -dXdVMR.Y);
2960 const Complex dFls = ls.dFdVMR(dXdVMR) + dFm;
2961 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dLM * S * Fls;
2962
2963 if (self_species == deriv.Target().species_id) {
2964 com.dF[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dSdSELFVMR() * LM * Fls;
2965 }
2966
2967 if (ls_str.scaler() == deriv.Target().species_id) {
2968 com.dF[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dSdOTHERVMR_if() * LM * Fls;
2969 }
2970
2971 if (ls_str.scaler() == Species::Species::Bath) {
2972 com.dF[com.jac_pos(iv, ij)] -= Sz * Sn * ls_str.dSdOTHERVMR_if() * LM * Fls;
2973 }
2974
2975 if (do_nlte) {
2976 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + dLM * DS * Fls;
2977
2978 if (self_species == deriv.QuantumIdentity().Species()) {
2979 com.dN[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dNdSELFVMR() * LM * Fls;
2980 }
2981
2982 if (ls_str.scaler() == deriv.QuantumIdentity().Species()) {
2983 com.dN[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dNdOTHERVMR_if() * LM * Fls;
2984 }
2985
2986 if (ls_str.scaler() == Species::Species::Bath) {
2987 com.dN[com.jac_pos(iv, ij)] -= Sz * Sn * ls_str.dNdOTHERVMR_if() * LM * Fls;
2988 }
2989 }
2990 } else {
2991 if (lt == Quantum::Number::StateMatchType::Full) {
2992 if (deriv == Jacobian::Line::Center) {
2993 const Numeric d = ls_norm.dNdF0() * Si + ls_str.dSdF0() * Sn;
2994 const Complex dFm = std::conj(ls_mirr.dFdF0());
2995 const Complex dFls = ls.dFdF0() + dFm;
2996 const Numeric dS = Sz * d;
2997 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dS * LM * Fls;
2998 if (do_nlte) {
2999 const Numeric Dd = ls_norm.dNdF0() * DNi + ls_str.dNdF0() * Sn;
3000 const Numeric DdS = Sz * Dd;
3001 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + DdS * LM * Fls;
3002 }
3003 } else if (deriv == Jacobian::Line::Strength) {
3004 const Numeric dS = Sz * Sn * ls_str.dSdI0();
3005 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
3006 if (do_nlte) {
3007 const Numeric DdS = Sz * Sn * ls_str.dNdI0();
3008 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
3009 }
3010 }
3011 // All line shape derivatives
3012 InternalDerivatives(G0) InternalDerivatives(D0)
3013 InternalDerivatives(G2) InternalDerivatives(D2)
3014 InternalDerivatives(ETA) InternalDerivatives(FVC)
3015 InternalDerivativesY InternalDerivativesG
3016 InternalDerivatives(DV)
3017 } else if (lt == Quantum::Number::StateMatchType::Level) {
3018 if (deriv == Jacobian::Line::NLTE) {
3019 if (lt.upp) {
3020 const Numeric dS = Sz * Sn * ls_str.dSdNLTEu();
3021 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
3022
3023 if (do_nlte) {
3024 const Numeric DdS = Sz * Sn * ls_str.dNdNLTEu();
3025 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
3026 }
3027 }
3028
3029 if (lt.low) {
3030 const Numeric dS = Sz * Sn * ls_str.dSdNLTEl();
3031 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
3032
3033 if (do_nlte) {
3034 const Numeric DdS = Sz * Sn * ls_str.dNdNLTEl();
3035 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
3036 }
3037 }
3038 }
3039 }
3040 }
3041 }
3042 }
3043 }
3044}
3045
3084void cutoff_loop_sparse_linear(
3085 ComputeData &com, ComputeData &sparse_com, Normalizer ls_norm,
3086 const IntensityCalculator ls_str, const AbsorptionLines &band,
3087 const ArrayOfDerivatives &derivs, const Output X, const Numeric &T,
3088 const Numeric &H, const Numeric &sparse_lim, const Numeric &DC,
3089 const Index i,
3090 const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT {
3091 // Basic settings
3092 const bool do_nlte = ls_str.do_nlte();
3093 const bool do_cutoff = band.cutoff not_eq Absorption::CutoffType::None;
3094 const Numeric fu = band.CutoffFreq(i, X.D0);
3095 const Numeric fl = band.CutoffFreqMinus(i, X.D0);
3096 const Numeric fus = band.lines[i].F0 + sparse_lim;
3097 const Numeric fls = band.lines[i].F0 - sparse_lim;
3098
3099 // Find sparse and dense ranges
3100 const auto [dense_start, dense_size, sparse_low_start, sparse_low_size,
3101 sparse_upp_start, sparse_upp_size] =
3102 linear_sparse_limited_range(fl, fu, fls, fus, com.f_grid,
3103 sparse_com.f_grid);
3104 if ((dense_size + sparse_low_size + sparse_upp_size) == 0)
3105 return;
3106
3107 // Get the compute data view
3108 ComputeValues comval(com.F, com.dF, com.N, com.dN, com.f_grid, dense_start,
3109 dense_size, derivs, do_nlte);
3110
3111 // Get views of the sparse data
3112 ComputeValues sparse_low_range(
3113 sparse_com.F, sparse_com.dF, sparse_com.N, sparse_com.dN,
3114 sparse_com.f_grid, sparse_low_start, sparse_low_size, derivs, do_nlte);
3115 ComputeValues sparse_upp_range(
3116 sparse_com.F, sparse_com.dF, sparse_com.N, sparse_com.dN,
3117 sparse_com.f_grid, sparse_upp_start, sparse_upp_size, derivs, do_nlte);
3118
3119 const Index nz = band.ZeemanCount(i, zeeman_polarization);
3120 for (Index iz = 0; iz < nz; iz++) {
3121 const Numeric dfdH = band.ZeemanSplitting(i, zeeman_polarization, iz);
3122 const Numeric Sz = band.ZeemanStrength(i, zeeman_polarization, iz);
3123 const Complex LM = Complex(1 + X.G, -X.Y);
3124 Calculator ls(band.lineshapetype, band.lines[i].F0, X, DC, dfdH * H,
3125 band.mirroring == Absorption::MirroringType::Manual);
3126 Calculator ls_mirr(band.mirroring, band.lineshapetype, band.lines[i].F0, X,
3127 DC, dfdH * H);
3128
3129 if (do_cutoff) {
3130 // Initialize and set the cutoff values
3131 Calculator ls_cut = ls;
3132 Calculator ls_mirr_cut = ls_mirr;
3133 ls_cut(fu);
3134 ls_mirr_cut(fu);
3135
3136 cutoff_frequency_loop(comval, ls, ls_mirr, ls_norm, ls_cut, ls_mirr_cut,
3137 ls_str, derivs, LM, T, dfdH, Sz, band.Species());
3138
3139 if (sparse_low_size) {
3140 cutoff_frequency_loop(sparse_low_range, ls, ls_mirr, ls_norm, ls_cut,
3141 ls_mirr_cut, ls_str, derivs, LM, T, dfdH, Sz,
3142 band.Species());
3143 }
3144
3145 if (sparse_upp_size) {
3146 cutoff_frequency_loop(sparse_upp_range, ls, ls_mirr, ls_norm, ls_cut,
3147 ls_mirr_cut, ls_str, derivs, LM, T, dfdH, Sz,
3148 band.Species());
3149 }
3150
3151 } else {
3152 frequency_loop(comval, ls, ls_mirr, ls_norm, ls_str, derivs, LM, T, dfdH,
3153 Sz, band.Species());
3154
3155 if (sparse_low_size) {
3156 frequency_loop(sparse_low_range, ls, ls_mirr, ls_norm, ls_str, derivs,
3157 LM, T, dfdH, Sz, band.Species());
3158 }
3159
3160 if (sparse_upp_size) {
3161 frequency_loop(sparse_upp_range, ls, ls_mirr, ls_norm, ls_str, derivs,
3162 LM, T, dfdH, Sz, band.Species());
3163 }
3164 }
3165 }
3166}
3167
3168void cutoff_loop_sparse_triple(
3169 ComputeData &com, ComputeData &sparse_com, Normalizer ls_norm,
3170 const IntensityCalculator ls_str, const AbsorptionLines &band,
3171 const ArrayOfDerivatives &derivs, const Output X, const Numeric &T,
3172 const Numeric &H, const Numeric &sparse_lim, const Numeric &DC,
3173 const Index i,
3174 const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT {
3175 // Basic settings
3176 const bool do_nlte = ls_str.do_nlte();
3177 const bool do_cutoff = band.cutoff not_eq Absorption::CutoffType::None;
3178 const Numeric fu = band.CutoffFreq(i, X.D0);
3179 const Numeric fl = band.CutoffFreqMinus(i, X.D0);
3180 const Numeric fus = band.lines[i].F0 + sparse_lim;
3181 const Numeric fls = band.lines[i].F0 - sparse_lim;
3182
3183 // Find sparse and dense ranges
3184 const auto [dense_start, dense_size, sparse_low_start, sparse_low_size,
3185 sparse_upp_start, sparse_upp_size] =
3186 quad_sparse_limited_range(fl, fu, fls, fus, com.f_grid,
3187 sparse_com.f_grid);
3188 if ((dense_size + sparse_low_size + sparse_upp_size) == 0)
3189 return;
3190
3191 // Get the compute data view
3192 ComputeValues comval(com.F, com.dF, com.N, com.dN, com.f_grid, dense_start,
3193 dense_size, derivs, do_nlte);
3194
3195 // Get views of the sparse data
3196 ComputeValues sparse_low_range(
3197 sparse_com.F, sparse_com.dF, sparse_com.N, sparse_com.dN,
3198 sparse_com.f_grid, sparse_low_start, sparse_low_size, derivs, do_nlte);
3199 ComputeValues sparse_upp_range(
3200 sparse_com.F, sparse_com.dF, sparse_com.N, sparse_com.dN,
3201 sparse_com.f_grid, sparse_upp_start, sparse_upp_size, derivs, do_nlte);
3202
3203 const Index nz = band.ZeemanCount(i, zeeman_polarization);
3204 for (Index iz = 0; iz < nz; iz++) {
3205 const Numeric dfdH = band.ZeemanSplitting(i, zeeman_polarization, iz);
3206 const Numeric Sz = band.ZeemanStrength(i, zeeman_polarization, iz);
3207 const Complex LM = Complex(1 + X.G, -X.Y);
3208 Calculator ls(band.lineshapetype, band.lines[i].F0, X, DC, dfdH * H,
3209 band.mirroring == Absorption::MirroringType::Manual);
3210 Calculator ls_mirr(band.mirroring, band.lineshapetype, band.lines[i].F0, X,
3211 DC, dfdH * H);
3212
3213 if (do_cutoff) {
3214 // Initialize and set the cutoff values
3215 Calculator ls_cut = ls;
3216 Calculator ls_mirr_cut = ls_mirr;
3217 ls_cut(fu);
3218 ls_mirr_cut(fu);
3219
3220 cutoff_frequency_loop(comval, ls, ls_mirr, ls_norm, ls_cut, ls_mirr_cut,
3221 ls_str, derivs, LM, T, dfdH, Sz, band.Species());
3222
3223 if (sparse_low_size) {
3224 cutoff_frequency_loop(sparse_low_range, ls, ls_mirr, ls_norm, ls_cut,
3225 ls_mirr_cut, ls_str, derivs, LM, T, dfdH, Sz,
3226 band.Species());
3227 }
3228
3229 if (sparse_upp_size) {
3230 cutoff_frequency_loop(sparse_upp_range, ls, ls_mirr, ls_norm, ls_cut,
3231 ls_mirr_cut, ls_str, derivs, LM, T, dfdH, Sz,
3232 band.Species());
3233 }
3234
3235 } else {
3236 frequency_loop(comval, ls, ls_mirr, ls_norm, ls_str, derivs, LM, T, dfdH,
3237 Sz, band.Species());
3238
3239 if (sparse_low_size) {
3240 frequency_loop(sparse_low_range, ls, ls_mirr, ls_norm, ls_str, derivs,
3241 LM, T, dfdH, Sz, band.Species());
3242 }
3243
3244 if (sparse_upp_size) {
3245 frequency_loop(sparse_upp_range, ls, ls_mirr, ls_norm, ls_str, derivs,
3246 LM, T, dfdH, Sz, band.Species());
3247 }
3248 }
3249 }
3250}
3251
3290void cutoff_loop(ComputeData &com, Normalizer ls_norm,
3291 const IntensityCalculator ls_str, const AbsorptionLines &band,
3292 const ArrayOfDerivatives &derivs, const Output X,
3293 const Numeric &T, const Numeric &H, const Numeric &DC,
3294 const Index i,
3295 const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT {
3296 // Basic settings
3297 const bool do_nlte = ls_str.do_nlte();
3298 const bool do_cutoff = band.cutoff not_eq Absorption::CutoffType::None;
3299 const Numeric fu = band.CutoffFreq(i, X.D0);
3300 const Numeric fl = band.CutoffFreqMinus(i, X.D0);
3301
3302 // Only for the cutoff-range
3303 const auto [cutstart, cutsize] = limited_range(fl, fu, com.f_grid);
3304 if (not cutsize)
3305 return;
3306
3307 // Get the compute data view
3308 ComputeValues comval(com.F, com.dF, com.N, com.dN, com.f_grid, cutstart,
3309 cutsize, derivs, do_nlte);
3310
3311 const Index nz = band.ZeemanCount(i, zeeman_polarization);
3312 for (Index iz = 0; iz < nz; iz++) {
3313 const Numeric dfdH = band.ZeemanSplitting(i, zeeman_polarization, iz);
3314 const Numeric Sz = band.ZeemanStrength(i, zeeman_polarization, iz);
3315 const Complex LM = Complex(1 + X.G, -X.Y);
3316 Calculator ls(band.lineshapetype, band.lines[i].F0, X, DC, dfdH * H,
3317 band.mirroring == Absorption::MirroringType::Manual);
3318 Calculator ls_mirr(band.mirroring, band.lineshapetype, band.lines[i].F0, X,
3319 DC, dfdH * H);
3320
3321 if (do_cutoff) {
3322 // Initialize and set the cutoff values
3323 Calculator ls_cut = ls;
3324 Calculator ls_mirr_cut = ls_mirr;
3325 ls_cut(fu);
3326 ls_mirr_cut(fu);
3327
3328 cutoff_frequency_loop(comval, ls, ls_mirr, ls_norm, ls_cut, ls_mirr_cut,
3329 ls_str, derivs, LM, T, dfdH, Sz, band.Species());
3330
3331 } else {
3332 frequency_loop(comval, ls, ls_mirr, ls_norm, ls_str, derivs, LM, T, dfdH,
3333 Sz, band.Species());
3334 }
3335 }
3336}
3337
3369void line_loop(ComputeData &com, ComputeData &sparse_com,
3370 const AbsorptionLines &band,
3371 const ArrayOfRetrievalQuantity &jacobian_quantities,
3372 const EnergyLevelMap &nlte, const Vector &vmrs,
3373 const ArrayOfSpeciesTag &self_tag, const Numeric &P,
3374 const Numeric &T, const Numeric &H, const Numeric &sparse_lim,
3375 const Numeric QT, const Numeric QT0, const Numeric dQTdT,
3376 const Numeric r, const Numeric drdSELFVMR, const Numeric drdT,
3377 const Zeeman::Polarization zeeman_polarization,
3378 const Options::LblSpeedup speedup_type) ARTS_NOEXCEPT {
3379 const Index nj = jacobian_quantities.nelem();
3380 const Index nl = band.NumLines();
3381
3382 // Derivatives are allocated ahead of all loops
3383 ArrayOfDerivatives derivs(nj);
3384
3385 // Doppler constant
3386 const Numeric DC = band.DopplerConstant(T);
3387
3388 if (not independent_per_broadener(band.lineshapetype)) {
3389 for (Index i = 0; i < nl; i++) {
3390 // Pre-compute the derivatives
3391 for (Index ij = 0; ij < nj; ij++) {
3392 const auto &deriv = jacobian_quantities[ij];
3393 derivs[ij].jac_pos = -1;
3394 derivs[ij].deriv = nullptr;
3395
3396 if (not deriv.propmattype())
3397 continue;
3398 derivs[ij].jac_pos = ij;
3399 derivs[ij].deriv = &deriv;
3400
3401 if (deriv == Jacobian::Atm::Temperature) {
3402 derivs[ij].value.o = band.ShapeParameters_dT(i, T, P, vmrs);
3403 } else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR) {
3404 if (not(deriv == self_tag)) { // Remove if its not good
3405 derivs[ij].jac_pos = -1;
3406 derivs[ij].deriv = nullptr;
3407 }
3408 } else if (deriv.Target().needQuantumIdentity()) {
3409 if (deriv == Jacobian::Line::VMR) {
3410 derivs[ij].value.o =
3411 band.ShapeParameters_dVMR(i, T, P, deriv.QuantumIdentity());
3412 } else {
3413 auto &lt = derivs[ij].target = {deriv.Target().qid,
3414 band.lines[i].localquanta,
3415 band.quantumidentity};
3416 if (lt == Quantum::Number::StateMatchType::Full) {
3417 if constexpr (false) { /*skip so the rest can be a else-if block*/
3418 }
3419 // All line shape derivatives
3420 InternalDerivativesSetup(G0) InternalDerivativesSetup(D0)
3421 InternalDerivativesSetup(G2) InternalDerivativesSetup(D2)
3422 InternalDerivativesSetup(ETA)
3423 InternalDerivativesSetup(FVC)
3424 InternalDerivativesSetup(Y)
3425 InternalDerivativesSetup(G)
3426 InternalDerivativesSetup(DV)
3427 }
3428 }
3429 }
3430 }
3431 std::remove_if(derivs.begin(), derivs.end(),
3432 [](Derivatives &dd) { return dd.deriv == nullptr; });
3433
3434 // Call cut off loop with or without sparsity
3435 switch (speedup_type) {
3436 case Options::LblSpeedup::None:
3437 cutoff_loop(com, Normalizer(band.normalization, band.lines[i].F0, T),
3438 IntensityCalculator(T, QT, QT0, dQTdT, r, drdSELFVMR, drdT,
3439 nlte, band, i),
3440 band, derivs, band.ShapeParameters(i, T, P, vmrs), T, H, DC,
3441 i, zeeman_polarization);
3442 break;
3443 case Options::LblSpeedup::QuadraticIndependent:
3444 cutoff_loop_sparse_triple(
3445 com, sparse_com,
3446 Normalizer(band.normalization, band.lines[i].F0, T),
3447 IntensityCalculator(T, QT, QT0, dQTdT, r, drdSELFVMR, drdT, nlte,
3448 band, i),
3449 band, derivs, band.ShapeParameters(i, T, P, vmrs), T, H, sparse_lim,
3450 DC, i, zeeman_polarization);
3451 break;
3452 case Options::LblSpeedup::LinearIndependent:
3453 cutoff_loop_sparse_linear(
3454 com, sparse_com,
3455 Normalizer(band.normalization, band.lines[i].F0, T),
3456 IntensityCalculator(T, QT, QT0, dQTdT, r, drdSELFVMR, drdT, nlte,
3457 band, i),
3458 band, derivs, band.ShapeParameters(i, T, P, vmrs), T, H, sparse_lim,
3459 DC, i, zeeman_polarization);
3460 break;
3461 case Options::LblSpeedup::FINAL: { /* Leave last */
3462 }
3463 }
3464 }
3465 } else {
3466 for (Index i = 0; i < nl; i++) {
3467 for (Index ib=0; ib<band.NumBroadeners(); ib++) {
3468 // Pre-compute the derivatives
3469 for (Index ij = 0; ij < nj; ij++) {
3470 const auto &deriv = jacobian_quantities[ij];
3471 derivs[ij].jac_pos = -1;
3472 derivs[ij].deriv = nullptr;
3473
3474 if (not deriv.propmattype())
3475 continue;
3476 derivs[ij].jac_pos = ij;
3477 derivs[ij].deriv = &deriv;
3478
3479 if (deriv == Jacobian::Atm::Temperature) {
3480 derivs[ij].value.o = band.ShapeParameters_dT(i, T, P, ib);
3481 } else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR) {
3482 if (not(deriv == self_tag)) { // Remove if its not good
3483 derivs[ij].jac_pos = -1;
3484 derivs[ij].deriv = nullptr;
3485 }
3486 } else if (deriv.Target().needQuantumIdentity()) {
3487 if (deriv == Jacobian::Line::VMR) {
3488 // The VMR goes into the line strength in this mode
3489 derivs[ij].value.o = LineShape::Output{};
3490 } else {
3491 auto &lt = derivs[ij].target = {deriv.Target().qid,
3492 band.lines[i].localquanta,
3493 band.quantumidentity};
3494 if (lt == Quantum::Number::StateMatchType::Full) {
3495 if constexpr (false) { /*skip so the rest can be a else-if block*/
3496 }
3497 // All line shape derivatives
3498 InternalDerivativesSetup(G0, ib)
3499 InternalDerivativesSetup(D0, ib)
3500 InternalDerivativesSetup(G2, ib)
3501 InternalDerivativesSetup(D2, ib)
3502 InternalDerivativesSetup(ETA, ib)
3503 InternalDerivativesSetup(FVC, ib)
3504 InternalDerivativesSetup(Y, ib)
3505 InternalDerivativesSetup(G, ib)
3506 InternalDerivativesSetup(DV, ib)
3507 }
3508 }
3509 }
3510 }
3511 std::remove_if(derivs.begin(), derivs.end(),
3512 [](Derivatives &dd) { return dd.deriv == nullptr; });
3513
3514 // The line shape strength rescaled by VMR of the broadener
3515 const auto ls_str = IntensityCalculator(T, QT, QT0, dQTdT, r,
3516 drdSELFVMR, drdT, nlte, band, i)
3517 .adaptive_scaling(vmrs[ib], band.Species(),
3518 band.broadeningspecies[ib]);
3519
3520 // Call cut off loop with or without sparsity
3521 switch (speedup_type) {
3522 case Options::LblSpeedup::None:
3523 cutoff_loop(com, Normalizer(band.normalization, band.lines[i].F0, T),
3524 ls_str, band, derivs, band.ShapeParameters(i, T, P, ib),
3525 T, H, DC, i, zeeman_polarization);
3526 break;
3527 case Options::LblSpeedup::QuadraticIndependent:
3528 cutoff_loop_sparse_triple(
3529 com, sparse_com,
3530 Normalizer(band.normalization, band.lines[i].F0, T), ls_str, band,
3531 derivs, band.ShapeParameters(i, T, P, ib), T, H, sparse_lim, DC,
3532 i, zeeman_polarization);
3533 break;
3534 case Options::LblSpeedup::LinearIndependent:
3535 cutoff_loop_sparse_linear(
3536 com, sparse_com,
3537 Normalizer(band.normalization, band.lines[i].F0, T), ls_str, band,
3538 derivs, band.ShapeParameters(i, T, P, ib), T, H, sparse_lim, DC,
3539 i, zeeman_polarization);
3540 break;
3541 case Options::LblSpeedup::FINAL: { /* Leave last */
3542 }
3543 }
3544 }
3545 }
3546 }
3547}
3548
3549void compute(ComputeData &com, ComputeData &sparse_com,
3550 const AbsorptionLines &band,
3551 const ArrayOfRetrievalQuantity &jacobian_quantities,
3552 const EnergyLevelMap &nlte, const Vector &vmrs,
3553 const ArrayOfSpeciesTag &self_tag, const Numeric &self_vmr,
3554 const Numeric &isot_ratio, const Numeric &P, const Numeric &T,
3555 const Numeric &H, const Numeric &sparse_lim,
3556 const Zeeman::Polarization zeeman_polarization,
3557 const Options::LblSpeedup speedup_type,
3558 const bool robust) ARTS_NOEXCEPT {
3559 [[maybe_unused]] const Index nj = jacobian_quantities.nelem();
3560 const Index nl = band.NumLines();
3561 const Index nv = com.f_grid.nelem();
3562
3563 // Tests that must be true while calling this function
3564 ARTS_ASSERT(H >= 0, "Only for positive H. You provided: ", H)
3565 ARTS_ASSERT(P > 0, "Only for abs positive P. You provided: ", P)
3566 ARTS_ASSERT(T > 0, "Only for abs positive T. You provided: ", T)
3567 ARTS_ASSERT(band.OK(), "Band is poorly constructed. You need to use "
3568 "a detailed debugger to find out why.")
3569 ARTS_ASSERT(com.F.size() == nv, "F is wrong size. Size is (", com.F.size(),
3570 ") but should be: (", nv, ')')
3571 ARTS_ASSERT(not com.do_nlte or com.N.size() == nv,
3572 "N is wrong size. Size is (", com.N.size(), ") but should be (",
3573 nv, ')')
3574 ARTS_ASSERT(nj == 0 or (com.dF.nrows() == nv and com.dF.ncols() == nj),
3575 "dF is wrong size. Size is (", com.dF.nrows(), " x ",
3576 com.dF.ncols(), ") but should be: (", nv, " x ", nj, ")")
3577 ARTS_ASSERT(nj == 0 or not com.do_nlte or
3578 (com.dN.nrows() == nv and com.dN.ncols() == nj),
3579 "dN is wrong size. Size is (", com.dN.nrows(), " x ",
3580 com.dN.ncols(), ") but should be: (", nv, " x ", nj, ")")
3581 ARTS_ASSERT((sparse_lim > 0 and sparse_com.f_grid.size() > 1) or
3582 (sparse_lim == 0),
3583 "Sparse limit is either 0, or the sparse frequency grid has to "
3584 "have upper and lower values")
3585
3586 // Early return test
3587 if (nv == 0 or nl == 0 or
3588 (Absorption::relaxationtype_relmat(band.population) and
3589 band.DoLineMixing(P))) {
3590 return; // No line-by-line computations required/wanted
3591 }
3592
3593 const Numeric dnumdensdVMR = isot_ratio * number_density(P, T);
3594
3595 if (robust and band.DoLineMixing(P) and band.AnyLinemixing()) {
3596 ComputeData com_safe(com.f_grid, jacobian_quantities, com.do_nlte);
3597 ComputeData sparse_com_safe(sparse_com.f_grid, jacobian_quantities,
3598 sparse_com.do_nlte);
3599
3600 line_loop(com_safe, sparse_com_safe, band, jacobian_quantities, nlte, vmrs,
3601 self_tag, P, T, H, sparse_lim,
3602 single_partition_function(T, band.Isotopologue()),
3603 single_partition_function(band.T0, band.Isotopologue()),
3604 dsingle_partition_function_dT(T, band.Isotopologue()),
3605 self_vmr * dnumdensdVMR, dnumdensdVMR,
3606 self_vmr * isot_ratio * dnumber_density_dt(P, T),
3607 zeeman_polarization, speedup_type);
3608
3609 com_safe.enforce_positive_absorption();
3610 sparse_com_safe.enforce_positive_absorption();
3611
3612 com += com_safe;
3613 sparse_com += sparse_com_safe;
3614 } else {
3615 line_loop(com, sparse_com, band, jacobian_quantities, nlte, vmrs, self_tag,
3616 P, T, H, sparse_lim,
3617 single_partition_function(T, band.Isotopologue()),
3618 single_partition_function(band.T0, band.Isotopologue()),
3619 dsingle_partition_function_dT(T, band.Isotopologue()),
3620 self_vmr * dnumdensdVMR, dnumdensdVMR,
3621 self_vmr * isot_ratio * dnumber_density_dt(P, T),
3622 zeeman_polarization, speedup_type);
3623 }
3624}
3625
3626#undef InternalDerivatives
3627#undef InternalDerivativesG
3628#undef InternalDerivativesY
3629
3630Index sparse_f_grid_red(const Vector &f_grid,
3631 const Numeric &sparse_df) noexcept {
3632 if (f_grid.nelem())
3633 return f_grid.nelem() /
3634 Index(1 +
3635 std::abs(f_grid[f_grid.nelem() - 1] - f_grid[0]) / sparse_df);
3636 return 0;
3637}
3638
3639Vector linear_sparse_f_grid(const Vector &f_grid,
3640 const Numeric &sparse_df) ARTS_NOEXCEPT {
3641 const Index n = sparse_f_grid_red(f_grid, sparse_df);
3642 const Index nv = f_grid.nelem();
3643
3644 if (nv and n) {
3645 std::vector<Numeric> sparse_f_grid;
3646 for (Index iv = 0; iv < nv - n; iv += n) {
3647 sparse_f_grid.emplace_back(f_grid[iv]);
3648 sparse_f_grid.emplace_back(f_grid[iv + n]);
3649 }
3650
3651 const Numeric f0 = sparse_f_grid.back();
3652 if (f0 not_eq f_grid[nv - 1]) {
3653 sparse_f_grid.emplace_back(f0);
3654 sparse_f_grid.emplace_back(f_grid[nv - 1]);
3655 }
3656
3657 return sparse_f_grid;
3658 }
3659 return Vector(0);
3660}
3661
3662bool good_linear_sparse_f_grid(const Vector &f_grid_dense,
3663 const Vector &f_grid_sparse) noexcept {
3664 const Index nf_sparse = f_grid_sparse.nelem();
3665 const Index nf_dense = f_grid_dense.nelem();
3666
3667 if (nf_sparse == 1)
3668 return false;
3669
3670 if (nf_sparse and nf_dense)
3671 return f_grid_dense[0] >= f_grid_sparse[0] and
3672 f_grid_dense[nf_dense - 1] <= f_grid_sparse[nf_sparse - 1];
3673
3674 return true;
3675}
3676
3677Vector triple_sparse_f_grid(const Vector &f_grid,
3678 const Numeric &sparse_df) noexcept {
3679 const Index n = sparse_f_grid_red(f_grid, sparse_df);
3680 const Index nv = f_grid.nelem();
3681
3682 if (nv and n > 2) {
3683 std::vector<Numeric> sparse_f_grid;
3684 for (Index iv = 0; iv < nv - n; iv += n) {
3685 sparse_f_grid.emplace_back(f_grid[iv]);
3686 sparse_f_grid.emplace_back(f_grid[iv] +
3687 0.5 * (f_grid[iv + n] - f_grid[iv]));
3688 sparse_f_grid.emplace_back(f_grid[iv + n]);
3689 }
3690
3691 const Numeric f0 = sparse_f_grid.back();
3692 if (f0 not_eq f_grid[nv - 1]) {
3693 sparse_f_grid.emplace_back(f0);
3694 sparse_f_grid.emplace_back(f0 + 0.5 * (f_grid[nv - 1] - f0));
3695 sparse_f_grid.emplace_back(f_grid[nv - 1]);
3696 }
3697
3698 return sparse_f_grid;
3699 }
3700 return Vector(0);
3701}
3702
3703void ComputeData::interp_add_even(const ComputeData &sparse) ARTS_NOEXCEPT {
3704 const Index nv = f_grid.nelem();
3705 const Index sparse_nv = sparse.f_grid.nelem();
3706 const Index nj = dF.ncols();
3707
3708 ARTS_ASSERT(do_nlte == sparse.do_nlte, "Must have the same NLTE status")
3709 ARTS_ASSERT(sparse_nv > 1, "Must have at least two sparse grid-points")
3710 ARTS_ASSERT(nv == 0 or (f_grid[0] == sparse.f_grid[0] and
3711 f_grid[nv - 1] >= sparse.f_grid[sparse_nv - 1]),
3712 "If there are any dense frequency points, then the sparse "
3713 "frequency points must fully cover them")
3714 ARTS_ASSERT(not(sparse_nv % 2), "Must be multiple of to")
3715
3716 Index sparse_iv = 0;
3717 Numeric f1 = sparse.f_grid[sparse_iv + 1];
3718 Numeric f0 = sparse.f_grid[sparse_iv];
3719 Numeric inv = 1.0 / (f1 - f0);
3720 for (Index iv = 0; iv < nv; iv++) {
3721 if (sparse_iv < (sparse_nv - 2) and f1 == f_grid[iv]) {
3722 sparse_iv += 2;
3723 f1 = sparse.f_grid[sparse_iv + 1];
3724 f0 = sparse.f_grid[sparse_iv];
3725 inv = 1.0 / (f1 - f0);
3726 }
3727
3728 const Numeric xm0 = f_grid[iv] - f0;
3729 const Numeric xm1 = f_grid[iv] - f1;
3730 const Numeric l0 = -xm1 * inv;
3731 const Numeric l1 = xm0 * inv;
3732
3733 F[iv] += l0 * sparse.F[sparse_iv] + l1 * sparse.F[sparse_iv + 1];
3734 for (Index ij = 0; ij < nj; ij++) {
3735 dF(iv, ij) +=
3736 l0 * sparse.dF(sparse_iv, ij) + l1 * sparse.dF(sparse_iv + 1, ij);
3737 }
3738 if (do_nlte) {
3739 N[iv] += l0 * sparse.N[sparse_iv] + l1 * sparse.N[sparse_iv + 1];
3740 for (Index ij = 0; ij < nj; ij++) {
3741 dN(iv, ij) +=
3742 l0 * sparse.dN(sparse_iv, ij) + l1 * sparse.dN(sparse_iv + 1, ij);
3743 }
3744 }
3745 }
3746}
3747
3748void ComputeData::interp_add_triplequad(const ComputeData &sparse)
3749 ARTS_NOEXCEPT {
3750 const Index nv = f_grid.nelem();
3751 const Index sparse_nv = sparse.f_grid.nelem();
3752 const Index nj = dF.ncols();
3753
3754 ARTS_ASSERT(do_nlte == sparse.do_nlte, "Must have the same NLTE status")
3755 ARTS_ASSERT(sparse_nv > 2, "Must have at least three sparse grid-points")
3756 ARTS_ASSERT(nv == 0 or (f_grid[0] == sparse.f_grid[0] and
3757 f_grid[nv - 1] >= sparse.f_grid[sparse_nv - 1]),
3758 "If there are any dense frequency points, then the sparse "
3759 "frequency points must fully cover them")
3760 ARTS_ASSERT(not(sparse_nv % 3), "Must be multiple of three")
3761
3762 Index sparse_iv = 0;
3763 Numeric f2 = sparse.f_grid[sparse_iv + 2];
3764 Numeric f1 = sparse.f_grid[sparse_iv + 1];
3765 Numeric f0 = sparse.f_grid[sparse_iv];
3766 Numeric inv = 1.0 / Math::pow2(f1 - f0);
3767 for (Index iv = 0; iv < nv; iv++) {
3768 if (sparse_iv < (sparse_nv - 3) and f2 == f_grid[iv]) {
3769 sparse_iv += 3;
3770 f2 = sparse.f_grid[sparse_iv + 2];
3771 f1 = sparse.f_grid[sparse_iv + 1];
3772 f0 = sparse.f_grid[sparse_iv];
3773 inv = 1.0 / Math::pow2(f1 - f0);
3774 }
3775 ARTS_ASSERT(
3776 f_grid[iv] >= f0 and (f_grid[iv] < f2 or (f2 == f_grid[iv] and
3777 sparse_iv == sparse_nv - 3)),
3778 "Out of range frequency grid. Must be caught earlier.\n"
3779 "The sparse range is from: ",
3780 f0, " to ", f2, " with ", f1,
3781 " as the half-way grid point.\n"
3782 "The dense frequency is ",
3783 f_grid[iv], " and the indices are: sparse_iv=", sparse_iv, "; iv=", iv)
3784
3785 const Numeric xm0 = f_grid[iv] - f0;
3786 const Numeric xm1 = f_grid[iv] - f1;
3787 const Numeric xm2 = f_grid[iv] - f2;
3788 const Numeric l0 = 0.5 * xm1 * xm2 * inv; // --
3789 const Numeric l1 = -xm0 * xm2 * inv; // +-
3790 const Numeric l2 = 0.5 * xm0 * xm1 * inv; // ++
3791
3792 F[iv] += l0 * sparse.F[sparse_iv] + l1 * sparse.F[sparse_iv + 1] +
3793 l2 * sparse.F[sparse_iv + 2];
3794 for (Index ij = 0; ij < nj; ij++) {
3795 dF(iv, ij) += l0 * sparse.dF(sparse_iv, ij) +
3796 l1 * sparse.dF(sparse_iv + 1, ij) +
3797 l2 * sparse.dF(sparse_iv + 2, ij);
3798 }
3799 if (do_nlte) {
3800 N[iv] += l0 * sparse.N[sparse_iv] + l1 * sparse.N[sparse_iv + 1] +
3801 l2 * sparse.N[sparse_iv + 2];
3802 for (Index ij = 0; ij < nj; ij++) {
3803 dN(iv, ij) += l0 * sparse.dN(sparse_iv, ij) +
3804 l1 * sparse.dN(sparse_iv + 1, ij) +
3805 l2 * sparse.dN(sparse_iv + 2, ij);
3806 }
3807 }
3808 }
3809}
3810} // namespace LineShape
Array
This can be used to make arrays out of anything.
Definition: array.h:48
Array::nelem
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
ComplexMatrix
The ComplexMatrix class.
Definition: matpack_complex.h:915
ComplexVector
The ComplexVector class.
Definition: matpack_complex.h:593
LineShape::Calculator
Line shape calculator.
Definition: lineshape.h:661
LineShape::Calculator::dFdH
Complex dFdH(Numeric dfdH) const noexcept
Definition: lineshape.cc:2274
LineShape::Calculator::dFdF0
Complex dFdF0() const noexcept
Definition: lineshape.cc:2270
LineShape::Calculator::F
Complex F() const noexcept
Definition: lineshape.cc:2310
LineShape::Calculator::dFdT
Complex dFdT(const Output &dXdT, Numeric T) const noexcept
Definition: lineshape.cc:2262
LineShape::Calculator::dFdD0
Complex dFdD0(Numeric d) const noexcept
Definition: lineshape.cc:2294
LineShape::Calculator::dFdVMR
Complex dFdVMR(const Output &dXdVMR) const noexcept
Definition: lineshape.cc:2278
LineShape::Calculator::dFdG0
Complex dFdG0(Numeric d) const noexcept
Definition: lineshape.cc:2298
LineShape::Calculator::dFdG2
Complex dFdG2(Numeric d) const noexcept
Definition: lineshape.cc:2306
LineShape::Calculator::dFdFVC
Complex dFdFVC(Numeric d) const noexcept
Definition: lineshape.cc:2282
LineShape::Calculator::Calculator
Calculator(const Type type, const Numeric F0, const Output &X, const Numeric DC, const Numeric DZ, bool manually_mirrored) noexcept
Definition: lineshape.cc:2318
LineShape::Calculator::dFdD2
Complex dFdD2(Numeric d) const noexcept
Definition: lineshape.cc:2302
LineShape::Calculator::operator()
Complex operator()(Numeric f) noexcept
Call operator on frequency. Must call this before any of the derivatives.
Definition: lineshape.cc:2314
LineShape::Calculator::dFdETA
Complex dFdETA(Numeric d) const noexcept
Definition: lineshape.cc:2286
LineShape::Calculator::dFdf
Complex dFdf() const noexcept
Definition: lineshape.cc:2266
LineShape::Calculator::dFdDV
Complex dFdDV(Numeric d) const noexcept
Definition: lineshape.cc:2290
LineShape::IntensityCalculator
Class encapsulating all supported types of intensity calculations of individual absorption lines.
Definition: lineshape.h:743
LineShape::IntensityCalculator::N
Numeric N() const noexcept
The line strength source offset.
Definition: lineshape.cc:2446
LineShape::IntensityCalculator::adaptive_scaling
IntensityCalculator & adaptive_scaling(Numeric x, Species::Species self, Species::Species other) noexcept
Rescale the line strength parameters by x.
Definition: lineshape.cc:2484
LineShape::IntensityCalculator::dSdI0
Numeric dSdI0() const noexcept
The line strength absorption derivative wrt the reference line strength.
Definition: lineshape.cc:2417
LineShape::IntensityCalculator::dSdNLTEl
Numeric dSdNLTEl() const noexcept
The line strength absorption derivative wrt either the lower state number density distribution or its...
Definition: lineshape.cc:2429
LineShape::IntensityCalculator::dNdF0
Numeric dNdF0() const noexcept
The line strength source offset derivative wrt the line center.
Definition: lineshape.cc:2458
LineShape::IntensityCalculator::dSdF0
Numeric dSdF0() const noexcept
The line strength absorption derivative wrt the line center.
Definition: lineshape.cc:2421
LineShape::IntensityCalculator::dSdNLTEu
Numeric dSdNLTEu() const noexcept
The line strength absorption derivative wrt either the upper state number density distribution or its...
Definition: lineshape.cc:2425
LineShape::IntensityCalculator::dNdOTHERVMR_if
Numeric dNdOTHERVMR_if() const noexcept
The line source offset derivative wrt the VMR of the rescaling species from adaptive_scaling (if self...
Definition: lineshape.cc:2477
LineShape::IntensityCalculator::S
Numeric S() const noexcept
The line strength absorption.
Definition: lineshape.cc:2409
LineShape::IntensityCalculator::dSdOTHERVMR_if
Numeric dSdOTHERVMR_if() const noexcept
The line strength derivative wrt the VMR of the rescaling species from adaptive_scaling (if self not_...
Definition: lineshape.cc:2440
LineShape::IntensityCalculator::dNdSELFVMR
Numeric dNdSELFVMR() const noexcept
The line source offset derivative wrt the VMR of the band's species.
Definition: lineshape.cc:2470
LineShape::IntensityCalculator::self_species
Species::Species self_species
Definition: lineshape.h:752
LineShape::IntensityCalculator::scaling_species
Species::Species scaling_species
Definition: lineshape.h:753
LineShape::IntensityCalculator::dNdT
Numeric dNdT() const noexcept
The line strength source offset derivative wrt temperature.
Definition: lineshape.cc:2450
LineShape::IntensityCalculator::dNdI0
Numeric dNdI0() const noexcept
The line strength source offset derivative wrt the reference line strength.
Definition: lineshape.cc:2454
LineShape::IntensityCalculator::dNdNLTEu
Numeric dNdNLTEu() const noexcept
The line strength source offset derivative wrt either the upper state number density distribution or ...
Definition: lineshape.cc:2462
LineShape::IntensityCalculator::dSdSELFVMR
Numeric dSdSELFVMR() const noexcept
The line strength derivative wrt the VMR of the band's species.
Definition: lineshape.cc:2433
LineShape::IntensityCalculator::dNdNLTEl
Numeric dNdNLTEl() const noexcept
The line strength source offset derivative wrt either the lower state number density distribution or ...
Definition: lineshape.cc:2466
LineShape::IntensityCalculator::dSdT
Numeric dSdT() const noexcept
The line strength absorption derivative wrt temperature.
Definition: lineshape.cc:2413
LineShape::IntensityCalculator::IntensityCalculator
IntensityCalculator(const Numeric T, const Numeric QT, const Numeric QT0, const Numeric dQTdT, const Numeric r, const Numeric drdSELFVMR, const Numeric drdT, const EnergyLevelMap &nlte, const Absorption::Lines &band, const Index line_index) noexcept
Definition: lineshape.cc:2492
LineShape::Normalizer::dNdF0
Numeric dNdF0() const noexcept
Definition: lineshape.cc:2378
LineShape::Normalizer::operator()
Numeric operator()(Numeric f) noexcept
Definition: lineshape.cc:2382
LineShape::Normalizer::dNdf
Numeric dNdf(Numeric f) const noexcept
Definition: lineshape.cc:2374
LineShape::Normalizer::dNdT
Numeric dNdT(Numeric T, Numeric f) const noexcept
Definition: lineshape.cc:2370
LineShape::Normalizer::Normalizer
Normalizer(const Absorption::NormalizationType type, const Numeric F0, const Numeric T) noexcept
Definition: lineshape.cc:2386
RetrievalQuantity
Deals with internal derivatives, Jacobian definition, and OEM calculations.
Definition: jacobian.h:325
Vector
The Vector class.
Definition: matpackI.h:910
debug.h
Helper macros for debugging.
ARTS_NOEXCEPT
#define ARTS_NOEXCEPT
Definition: debug.h:99
ARTS_ASSERT
#define ARTS_ASSERT(condition,...)
Definition: debug.h:102
inv
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:167
stimulated_emission
Numeric stimulated_emission(Numeric T, Numeric F0)
Computes exp(-hf/kT)
Definition: linescaling.cc:33
dabsorption_nlte_rate_dT
Numeric dabsorption_nlte_rate_dT(const Numeric &gamma, const Numeric &T, const Numeric &F0, const Numeric &El, const Numeric &Eu, const Numeric &r_upp, const Numeric &r_low)
Computes temperature derivatives of (r_low - r_upp * gamma) / (1 - gamma)
Definition: linescaling.cc:135
dabsorption_nlte_rate_dF0
Numeric dabsorption_nlte_rate_dF0(const Numeric &gamma, const Numeric &T, const Numeric &r_upp, const Numeric &r_low)
Computes frequency derivative of (r_low - r_upp * gamma) / (1 - gamma)
Definition: linescaling.cc:157
dsingle_partition_function_dT
Numeric dsingle_partition_function_dT(const Numeric &T, const Species::IsotopeRecord &ir)
Computes the partition function temperature derivative.
Definition: linescaling.cc:28
dstimulated_relative_emission_dF0
Numeric dstimulated_relative_emission_dF0(const Numeric F0, const Numeric T0, const Numeric T) noexcept
Computes.
Definition: linescaling.cc:61
boltzman_ratio
Numeric boltzman_ratio(const Numeric &T, const Numeric &T0, const Numeric &E0)
Computes exp(E0/c (T - T0) / (T * T0))
Definition: linescaling.cc:94
dabsorption_nlte_rate_dTl
Numeric dabsorption_nlte_rate_dTl(const Numeric &gamma, const Numeric &T, const Numeric &Tl, const Numeric &El, const Numeric &r_low)
Computes lower state temperature derivative of (r_low - r_upp * gamma) / (1 - gamma)
Definition: linescaling.cc:166
single_partition_function
Numeric single_partition_function(const Numeric &T, const Species::IsotopeRecord &ir)
Computes the partition function at one temperature.
Definition: linescaling.cc:23
dboltzman_ratio_dT
Numeric dboltzman_ratio_dT(const Numeric &boltzmann_ratio, const Numeric &T, const Numeric &E0)
Computes temperature derivatives exp(E0/k (T - T0) / (T * T0))
Definition: linescaling.cc:100
stimulated_relative_emission
Numeric stimulated_relative_emission(const Numeric F0, const Numeric T0, const Numeric T) noexcept
Computes.
Definition: linescaling.cc:52
absorption_nlte_ratio
Numeric absorption_nlte_ratio(const Numeric &gamma, const Numeric &r_upp, const Numeric &r_low) noexcept
Computes (r_low - r_upp * gamma) / (1 - gamma)
Definition: linescaling.cc:129
dabsorption_nlte_rate_dTu
Numeric dabsorption_nlte_rate_dTu(const Numeric &gamma, const Numeric &T, const Numeric &Tu, const Numeric &Eu, const Numeric &r_upp)
Computes upper state temperature derivative of (r_low - r_upp * gamma) / (1 - gamma)
Definition: linescaling.cc:177
dstimulated_relative_emission_dT
Numeric dstimulated_relative_emission_dT(const Numeric F0, const Numeric T0, const Numeric T) noexcept
Computes.
Definition: linescaling.cc:56
dboltzman_ratio_dT_div_boltzmann_ratio
constexpr Numeric dboltzman_ratio_dT_div_boltzmann_ratio(Numeric T, Numeric E0)
Computes temperature derivatives exp(E0/k (T - T0) / (T * T0)) / exp(E0/c (T - T0) / (T * T0))
Definition: linescaling.h:191
CutInternalDerivatives
#define CutInternalDerivatives(X)
Definition: lineshape.cc:2537
InternalDerivatives
#define InternalDerivatives(X)
Definition: lineshape.cc:2552
InternalDerivativesSetup
#define InternalDerivativesSetup(X,...)
Definition: lineshape.cc:2595
InternalDerivativesY
#define InternalDerivativesY
Definition: lineshape.cc:2578
ln_16
constexpr Numeric ln_16
Definition: lineshape.cc:24
InternalDerivativesG
#define InternalDerivativesG
Definition: lineshape.cc:2565
lineshapemodel.h
Contains the line shape namespace.
fac
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:65
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
real_val
Numeric & real_val(Complex &c) noexcept
Return a reference to the real value of c.
Definition: matpack_complex.h:40
imag_val
Numeric & imag_val(Complex &c) noexcept
Return a reference to the imaginary value of c.
Definition: matpack_complex.h:43
Complex
std::complex< Numeric > Complex
Definition: matpack_complex.h:36
Absorption::relaxationtype_relmat
constexpr bool relaxationtype_relmat(PopulationType in) noexcept
Definition: absorptionlines.h:145
Constant::inv_sqrt_pi
constexpr Numeric inv_sqrt_pi
Inverse of the square root of pi.
Definition: arts_constants.h:90
Constant::sqrt_ln_2
constexpr Numeric sqrt_ln_2
Square root of natural logarithm of 2.
Definition: arts_constants.h:117
Constant::pi
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
Definition: arts_constants.h:75
Constant::k
constexpr Numeric k
Boltzmann constant convenience name [J/K].
Definition: arts_constants.h:170
Constant::sqrt_pi
constexpr Numeric sqrt_pi
Square root of pi.
Definition: arts_constants.h:87
Constant::c
constexpr Numeric c
Speed of light convenience name [m/s].
Definition: arts_constants.h:137
Constant::inv_pi
constexpr Numeric inv_pi
Inverse of pi.
Definition: arts_constants.h:78
Constant::h
constexpr Numeric h
Planck constant convenience name [J s].
Definition: arts_constants.h:146
LineShape
Computations of line shape derived parameters.
Definition: lineshape.cc:27
LineShape::frequency_loop
void frequency_loop(ComputeValues &com, Calculator &ls, Calculator &ls_mirr, Normalizer &ls_norm, const IntensityCalculator &ls_str, const ArrayOfDerivatives &derivs, const Complex LM, const Numeric &T, const Numeric &dfdH, const Numeric &Sz, const Species::Species self_species) ARTS_NOEXCEPT
Frequency loop of the line shape call.
Definition: lineshape.cc:2890
LineShape::sparse_f_grid_red
Index sparse_f_grid_red(const Vector &f_grid, const Numeric &sparse_df) noexcept
Definition: lineshape.cc:3630
LineShape::cutoff_frequency_loop
void cutoff_frequency_loop(ComputeValues &com, Calculator &ls, Calculator &ls_mirr, Normalizer &ls_norm, const Calculator &ls_cut, const Calculator &ls_mirr_cut, const IntensityCalculator &ls_str, const ArrayOfDerivatives &derivs, const Complex LM, const Numeric &T, const Numeric &dfdH, const Numeric &Sz, const Species::Species self_species) ARTS_NOEXCEPT
Cutoff frequency loop of the line shape call.
Definition: lineshape.cc:2663
LineShape::active_nelem
Index active_nelem(const ArrayOfDerivatives &derivs) noexcept
Helper function to find the last relevant derivative.
Definition: lineshape.cc:2198
LineShape::abs_squared
constexpr Numeric abs_squared(Complex z) noexcept
Definition: lineshape.cc:415
LineShape::mirroredOutput
constexpr Output mirroredOutput(Output x) noexcept
Output to be used by mirroring calls.
Definition: lineshapemodel.h:564
LineShape::compute
void compute(ComputeData &com, ComputeData &sparse_com, const AbsorptionLines &band, const ArrayOfRetrievalQuantity &jacobian_quantities, const EnergyLevelMap &nlte, const Vector &vmrs, const ArrayOfSpeciesTag &self_tag, const Numeric &self_vmr, const Numeric &isot_ratio, const Numeric &P, const Numeric &T, const Numeric &H, const Numeric &sparse_lim, const Zeeman::Polarization zeeman_polarization, const Options::LblSpeedup speedup_type, const bool robust) ARTS_NOEXCEPT
Compute the absorption of an absorption band.
Definition: lineshape.cc:3549
LineShape::cutoff_loop_sparse_linear
void cutoff_loop_sparse_linear(ComputeData &com, ComputeData &sparse_com, Normalizer ls_norm, const IntensityCalculator ls_str, const AbsorptionLines &band, const ArrayOfDerivatives &derivs, const Output X, const Numeric &T, const Numeric &H, const Numeric &sparse_lim, const Numeric &DC, const Index i, const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT
Cutoff considerations of the line shape.
Definition: lineshape.cc:3084
LineShape::line_loop
void line_loop(ComputeData &com, ComputeData &sparse_com, const AbsorptionLines &band, const ArrayOfRetrievalQuantity &jacobian_quantities, const EnergyLevelMap &nlte, const Vector &vmrs, const ArrayOfSpeciesTag &self_tag, const Numeric &P, const Numeric &T, const Numeric &H, const Numeric &sparse_lim, const Numeric QT, const Numeric QT0, const Numeric dQTdT, const Numeric r, const Numeric drdSELFVMR, const Numeric drdT, const Zeeman::Polarization zeeman_polarization, const Options::LblSpeedup speedup_type) ARTS_NOEXCEPT
Loop all the lines of the band.
Definition: lineshape.cc:3369
LineShape::cutoff_loop_sparse_triple
void cutoff_loop_sparse_triple(ComputeData &com, ComputeData &sparse_com, Normalizer ls_norm, const IntensityCalculator ls_str, const AbsorptionLines &band, const ArrayOfDerivatives &derivs, const Output X, const Numeric &T, const Numeric &H, const Numeric &sparse_lim, const Numeric &DC, const Index i, const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT
Definition: lineshape.cc:3168
LineShape::good_linear_sparse_f_grid
bool good_linear_sparse_f_grid(const Vector &f_grid_dense, const Vector &f_grid_sparse) noexcept
Definition: lineshape.cc:3662
LineShape::quad_sparse_limited_range
SparseLimitRange quad_sparse_limited_range(const Numeric flc, const Numeric fuc, const Numeric fls, const Numeric fus, const Vector &f_grid, const Vector &sparse_f_grid) ARTS_NOEXCEPT
Definition: lineshape.cc:2106
LineShape::cutoff_loop
void cutoff_loop(ComputeData &com, Normalizer ls_norm, const IntensityCalculator ls_str, const AbsorptionLines &band, const ArrayOfDerivatives &derivs, const Output X, const Numeric &T, const Numeric &H, const Numeric &DC, const Index i, const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT
Cutoff considerations of the line shape.
Definition: lineshape.cc:3290
LineShape::vmrs
Vector vmrs(const ConstVectorView &atmospheric_vmrs, const ArrayOfArrayOfSpeciesTag &atmospheric_species, const ArrayOfSpecies &lineshape_species) ARTS_NOEXCEPT
Returns a VMR vector for this model's main calculations.
Definition: lineshapemodel.cc:458
LineShape::triple_sparse_f_grid
Vector triple_sparse_f_grid(const Vector &f_grid, const Numeric &sparse_df) noexcept
Definition: lineshape.cc:3677
LineShape::linear_sparse_limited_range
SparseLimitRange linear_sparse_limited_range(const Numeric flc, const Numeric fuc, const Numeric fls, const Numeric fus, const Vector &f_grid, const Vector &sparse_f_grid) ARTS_NOEXCEPT
Definition: lineshape.cc:2048
LineShape::linear_sparse_f_grid
Vector linear_sparse_f_grid(const Vector &f_grid, const Numeric &sparse_df) ARTS_NOEXCEPT
Definition: lineshape.cc:3639
LineShape::limited_range
CutoffRange limited_range(const Numeric fl, const Numeric fu, const Vector &f_grid) ARTS_NOEXCEPT
Gets the start and size of a range such that.
Definition: lineshape.cc:2031
Math::pow3
constexpr auto pow3(auto x) noexcept
power of three
Definition: arts_constexpr_math.h:42
Math::pow4
constexpr auto pow4(auto x) noexcept
power of four
Definition: arts_constexpr_math.h:45
Math::pow2
constexpr auto pow2(auto x) noexcept
power of two
Definition: arts_constexpr_math.h:39
Zeeman::Polarization
Polarization
Zeeman polarization selection.
Definition: zeemandata.h:44
nonstd::abs
constexpr T abs(T x) noexcept
Definition: nonstd.h:13
physics_funcs.h
This file contains declerations of functions of physical character.
dnumber_density_dt
constexpr Numeric dnumber_density_dt(Numeric p, Numeric t) noexcept
dnumber_density_dT
Definition: physics_funcs.h:84
number_density
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Definition: physics_funcs.h:70
N
#define N
Definition: rng.cc:164
LineShape::ComputeData
Main computational data for the line shape and strength calculations.
Definition: lineshape.h:840
LineShape::ComputeData::interp_add_even
void interp_add_even(const ComputeData &sparse) ARTS_NOEXCEPT
Add a sparse grid to this grid via linear interpolation.
Definition: lineshape.cc:3703
LineShape::ComputeData::interp_add_triplequad
void interp_add_triplequad(const ComputeData &sparse) ARTS_NOEXCEPT
Add a sparse grid to this grid via square interpolation.
Definition: lineshape.cc:3748
LineShape::ComputeData::enforce_positive_absorption
void enforce_positive_absorption() noexcept
All four fields are set to zero at i if F[i].real() < 0.
Definition: lineshape.h:876
LineShape::ComputeValues::derivs
const ArrayOfDerivatives & derivs
Definition: lineshape.cc:2212
LineShape::ComputeValues::ComputeValues
ComputeValues(ComplexVector &F_, ComplexMatrix &dF_, ComplexVector &N_, ComplexMatrix &dN_, const Vector &f_grid, const Index start, const Index nv, const ArrayOfDerivatives &derivs_, const bool do_nlte_) noexcept
Definition: lineshape.cc:2221
LineShape::ComputeValues::operator-=
ComputeValues & operator-=(const ComputeValues &cut) ARTS_NOEXCEPT
Definition: lineshape.cc:2240
LineShape::ComputeValues::f
const Numeric *const f
Definition: lineshape.cc:2210
LineShape::ComputeValues::ComputeValues
ComputeValues(Complex &F_, std::vector< Complex > &dF_, Complex &N_, std::vector< Complex > &dN_, const Numeric &f_lim, const ArrayOfDerivatives &derivs_, const bool do_nlte_) noexcept
Definition: lineshape.cc:2233
LineShape::ComputeValues::jac_pos
Index jac_pos(Index iv, Index ij) const noexcept
Definition: lineshape.cc:2217
LineShape::CutoffRange
Struct to keep the cutoff limited range values.
Definition: lineshape.cc:2016
LineShape::Derivatives
Data struct for keeping derivative keys and values.
Definition: lineshape.cc:2182
LineShape::Derivatives::deriv
const RetrievalQuantity * deriv
Definition: lineshape.cc:2190
LineShape::Derivatives::value
union LineShape::Derivatives::Values value
LineShape::Derivatives::target
Quantum::Number::StateMatch target
Definition: lineshape.cc:2183
LineShape::Doppler::operator()
Complex operator()(Numeric f) noexcept
Definition: lineshape.cc:28
LineShape::FullNonLocalThermodynamicEquilibriumInitialization::operator()
constexpr FullNonLocalThermodynamicEquilibrium operator()(Numeric r, Numeric drdSELFVMR, Numeric drdT) &&noexcept
Definition: lineshape.cc:1952
LineShape::FullNonLocalThermodynamicEquilibriumInitialization::FullNonLocalThermodynamicEquilibriumInitialization
FullNonLocalThermodynamicEquilibriumInitialization(Numeric F0, Numeric A21, Numeric T, Numeric r1, Numeric r2, Numeric c2, Numeric c3, Numeric x) noexcept
Definition: lineshape.cc:1934
LineShape::FullNonLocalThermodynamicEquilibrium::FullNonLocalThermodynamicEquilibrium
constexpr FullNonLocalThermodynamicEquilibrium(Numeric r, Numeric drdSELFVMR, Numeric drdt, Numeric k, Numeric dkdF0, Numeric dkdr1, Numeric dkdr2, Numeric e, Numeric dedF0, Numeric dedr2, Numeric B, Numeric dBdT, Numeric dBdF0) noexcept
Definition: lineshape.h:488
LineShape::HartmannTran::dFdG0
Complex dFdG0(Numeric dG0) const noexcept
Definition: lineshape.cc:762
LineShape::HartmannTran::dFdH
Complex dFdH(Numeric dZ) const noexcept
Definition: lineshape.cc:1298
LineShape::HartmannTran::dFdT
Complex dFdT(const Output &d, Numeric T) const noexcept
Definition: lineshape.cc:1571
LineShape::HartmannTran::init
CalcType init(const Complex c2t) const noexcept
Definition: lineshape.cc:1754
LineShape::HartmannTran::HartmannTran
HartmannTran(Numeric F0_noshift, const Output &ls, Numeric GD_div_F0, Numeric dZ) noexcept
Definition: lineshape.cc:473