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