ARTS 2.5.0 (git: 9ee3ac6c)
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 Constant::pow2;
15using Constant::pow3;
16using Constant::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, const Output &ls,
38 Numeric GD_div_F0, Numeric dZ) noexcept
39 : mF0(F0_noshift + dZ + ls.D0 - 1.5 * ls.D2),
40 invGD(sqrt_ln_2 / nonstd::abs(GD_div_F0 * mF0)),
41 invc2(1.0 / Complex(ls.G2, ls.D2)), dx(Complex(ls.G0 - 1.5 * ls.G2, mF0)),
42 x(dx * invc2), sqrty(invc2 / (2 * invGD)),
43 calcs(init(Complex(ls.G2, ls.D2))) {
44 calc();
45}
46
47Complex SpeedDependentVoigt::dFdf() const noexcept {
48 switch (calcs) {
49 case CalcType::Full:
50 return invGD * invc2 * (dw1 - dw2) / (2 * sqrt_pi * sq);
51 case CalcType::Voigt:
52 return dw1 * pow2(invGD) * inv_sqrt_pi;
53 case CalcType::LowXandHighY:
54 return dw1 * pow2(invGD) * inv_sqrt_pi -
55 dw2 * invGD * invc2 / (2 * sqrt_pi * sq);
56 case CalcType::LowYandLowX:
57 return pow2(invc2) * (-dw1 * sq + 1i * w1) / (sqrt_pi * sq);
58 case CalcType::LowYandHighX:
59 return 1i * pow2(invc2) * (x - 3) / (pi * pow3(x));
60 }
61 return {};
62}
63
64Complex SpeedDependentVoigt::dFdF0() const noexcept {
65 switch (calcs) {
66 case CalcType::Full:
67 return (4 * pow2(invGD) * (-w1 + w2) * sq +
68 1i * invc2 *
69 (dw1 * (Complex(0, 2 * mF0 * pow2(invGD)) - 2 * invGD * sq +
70 invc2) -
71 dw2 * (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
72 invc2))) /
73 (4 * sqrt_pi * invGD * mF0 * sq);
74 case CalcType::Voigt:
75 return -invGD * (Complex(0, invGD) * dw1 * (dx - Complex(0, mF0)) + w1) /
76 (sqrt_pi * mF0);
77 case CalcType::LowXandHighY:
78 return (4 * pow2(invGD) * (-w1 + w2) * sq -
79 1i * (4 * dw1 * pow3(invGD) * (dx - Complex(0, mF0)) * sq +
80 dw2 * invc2 *
81 (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
82 invc2))) /
83 (4 * sqrt_pi * invGD * mF0 * sq);
84 case CalcType::LowYandLowX:
85 return pow2(invc2) * (dw1 * sq - 1i * w1) / (sqrt_pi * sq);
86 case CalcType::LowYandHighX:
87 return 1i * pow2(invc2) * (3 - x) / (pi * pow3(x));
88 }
89 return {};
90}
91
92Complex SpeedDependentVoigt::dFdD0(Numeric dD0dD0) const noexcept {
93 switch (calcs) {
94 case CalcType::Full:
95 return -dD0dD0 *
96 (4 * pow2(invGD) * (w1 - w2) * sq +
97 1i * invc2 *
98 (-dw1 * (Complex(0, 2 * mF0 * pow2(invGD)) - 2 * invGD * sq +
99 invc2) +
100 dw2 * (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
101 invc2))) /
102 (4 * sqrt_pi * invGD * mF0 * sq);
103 case CalcType::Voigt:
104 return -dD0dD0 * invGD *
105 (Complex(0, invGD) * dw1 * (dx - Complex(0, mF0)) + w1) /
106 (sqrt_pi * mF0);
107 case CalcType::LowXandHighY:
108 return -dD0dD0 *
109 (4 * pow2(invGD) * (w1 - w2) * sq +
110 1i * (4 * dw1 * pow3(invGD) * (dx - Complex(0, mF0)) * sq +
111 dw2 * invc2 *
112 (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
113 invc2))) /
114 (4 * sqrt_pi * invGD * mF0 * sq);
115 case CalcType::LowYandLowX:
116 return dD0dD0 * pow2(invc2) * (dw1 * sq - 1i * w1) / (sqrt_pi * sq);
117 case CalcType::LowYandHighX:
118 return -Complex(0, dD0dD0) * pow2(invc2) * (x - 3) / (pi * pow3(x));
119 }
120 return {};
121}
122
123Complex SpeedDependentVoigt::dFdG0(Numeric dG0dG0) const noexcept {
124 switch (calcs) {
125 case CalcType::Full:
126 return Complex(0, dG0dG0) * invGD * invc2 * (dw1 - dw2) /
127 (2 * sqrt_pi * sq);
128 case CalcType::Voigt:
129 return Complex(0, dG0dG0) * dw1 * pow2(invGD) * inv_sqrt_pi;
130 case CalcType::LowXandHighY:
131 return Complex(0, dG0dG0) * invGD * (2 * dw1 * invGD * sq - dw2 * invc2) /
132 (2 * sqrt_pi * sq);
133 case CalcType::LowYandLowX:
134 return -dG0dG0 * pow2(invc2) * (1i * dw1 * sq + w1) / (sqrt_pi * sq);
135 case CalcType::LowYandHighX:
136 return -dG0dG0 * pow2(invc2) * (x - 3) / (pi * pow3(x));
137 }
138 return {};
139}
140
141Complex SpeedDependentVoigt::dFdD2(Numeric dD2dD2) const noexcept {
142 switch (calcs) {
143 case CalcType::Full:
144 return dD2dD2 *
145 (12 * pow2(invGD) * (w1 - w2) * sq +
146 1i * invc2 *
147 (dw1 * (-Complex(0, 2 * mF0 * pow2(invGD)) *
148 (2 * dx * invc2 + 3) +
149 4 * Complex(0, invGD) * invc2 * mF0 * sq +
150 6 * invGD * sq - Complex(0, 2 * mF0) * pow2(invc2) -
151 3 * invc2) +
152 dw2 *
153 (Complex(0, 2 * mF0 * pow2(invGD)) * (2 * dx * invc2 + 3) +
154 4 * Complex(0, invGD) * invc2 * mF0 * sq +
155 6 * invGD * sq + Complex(0, 2 * mF0) * pow2(invc2) +
156 3 * invc2))) /
157 (8 * sqrt_pi * invGD * mF0 * sq);
158 case CalcType::Voigt:
159 return 3 * dD2dD2 * invGD *
160 (Complex(0, invGD) * dw1 * (dx - Complex(0, mF0)) + w1) /
161 (2 * sqrt_pi * mF0);
162 case CalcType::LowXandHighY:
163 return dD2dD2 *
164 (12 * pow2(invGD) * (w1 - w2) * sq +
165 1i *
166 (12 * dw1 * pow3(invGD) * (dx - Complex(0, mF0)) * sq +
167 dw2 * invc2 *
168 (Complex(0, 2 * mF0 * pow2(invGD)) * (2 * dx * invc2 + 3) +
169 4 * Complex(0, invGD) * invc2 * mF0 * sq +
170 6 * invGD * sq + Complex(0, 2 * mF0) * pow2(invc2) +
171 3 * invc2))) /
172 (8 * sqrt_pi * invGD * mF0 * sq);
173 case CalcType::LowYandLowX:
174 return dD2dD2 * pow2(invc2) *
175 (4 * 1i * sq * (sqrt_pi * w1 * sq - 1) -
176 sqrt_pi * (dw1 * sq - 1i * w1) * (2 * dx * invc2 + 3)) /
177 (2 * pi * sq);
178 case CalcType::LowYandHighX:
179 return Complex(0, dD2dD2) * pow2(invc2) *
180 (-x * (2 * x - 3) + (x - 3) * (2 * dx * invc2 + 3)) /
181 (2 * pi * pow3(x));
182 }
183 return {};
184}
185
186Complex SpeedDependentVoigt::dFdG2(Numeric dG2dG2) const noexcept {
187 switch (calcs) {
188 case CalcType::Full:
189 return Complex(0, dG2dG2) * invc2 *
190 (dw1 * (-pow2(invGD) * (2 * dx * invc2 + 3) +
191 2 * invGD * invc2 * sq - pow2(invc2)) +
192 dw2 * (pow2(invGD) * (2 * dx * invc2 + 3) + 2 * invGD * invc2 * sq +
193 pow2(invc2))) /
194 (4 * sqrt_pi * invGD * sq);
195 case CalcType::Voigt:
196 return -3 * Complex(0, dG2dG2) * dw1 * pow2(invGD) / (2 * sqrt_pi);
197 case CalcType::LowXandHighY:
198 return Complex(0, dG2dG2) *
199 (-6 * dw1 * pow3(invGD) * sq +
200 dw2 * invc2 *
201 (pow2(invGD) * (2 * dx * invc2 + 3) + 2 * invGD * invc2 * sq +
202 pow2(invc2))) /
203 (4 * sqrt_pi * invGD * sq);
204 case CalcType::LowYandLowX:
205 return dG2dG2 * pow2(invc2) *
206 (4 * sq * (sqrt_pi * w1 * sq - 1) +
207 sqrt_pi * (2 * dx * invc2 + 3) * (1i * dw1 * sq + w1)) /
208 (2 * pi * sq);
209 case CalcType::LowYandHighX:
210 return dG2dG2 * pow2(invc2) *
211 (-x * (2 * x - 3) + (x - 3) * (2 * dx * invc2 + 3)) /
212 (2 * pi * pow3(x));
213 }
214 return {};
215}
216
217Complex SpeedDependentVoigt::dFdH(Numeric dZ) const noexcept {
218 switch (calcs) {
219 case CalcType::Full:
220 return -dZ *
221 (4 * pow2(invGD) * (w1 - w2) * sq +
222 1i * invc2 *
223 (-dw1 * (Complex(0, 2 * mF0 * pow2(invGD)) - 2 * invGD * sq +
224 invc2) +
225 dw2 * (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
226 invc2))) /
227 (4 * sqrt_pi * invGD * mF0 * sq);
228 case CalcType::Voigt:
229 return -dZ * invGD *
230 (Complex(0, invGD) * dw1 * (dx - Complex(0, mF0)) + w1) /
231 (sqrt_pi * mF0);
232 case CalcType::LowXandHighY:
233 return -dZ *
234 (4 * pow2(invGD) * (w1 - w2) * sq +
235 1i * (4 * dw1 * pow3(invGD) * (dx - Complex(0, mF0)) * sq +
236 dw2 * invc2 *
237 (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
238 invc2))) /
239 (4 * sqrt_pi * invGD * mF0 * sq);
240 case CalcType::LowYandLowX:
241 return dZ * pow2(invc2) * (dw1 * sq - 1i * w1) / (sqrt_pi * sq);
242 case CalcType::LowYandHighX:
243 return -Complex(0, dZ) * pow2(invc2) * (x - 3) / (pi * pow3(x));
244 }
245 return {};
246}
247
248Complex SpeedDependentVoigt::dFdVMR(const Output &d) const noexcept {
249 switch (calcs) {
250 case CalcType::Full:
251 return (-4 * pow2(invGD) * (2 * d.D0 - 3 * d.D2) * (w1 - w2) * sq +
252 1i * invc2 *
253 (dw1 * (-2 * pow2(invGD) * mF0 *
254 (Complex(3 * d.G2 - 2 * d.G0, 3 * d.D2 - 2 * d.D0) +
255 2 * dx * invc2 * Complex(d.G2, d.D2)) +
256 4 * invGD * invc2 * mF0 * sq * Complex(d.G2, d.D2) -
257 2 * invGD * (2 * d.D0 - 3 * d.D2) * sq -
258 2 * pow2(invc2) * mF0 * Complex(d.G2, d.D2) +
259 invc2 * (2 * d.D0 - 3 * d.D2)) -
260 dw2 * (-2 * pow2(invGD) * mF0 *
261 (Complex(3 * d.G2 - 2 * d.G0, 3 * d.D2 - 2 * d.D0) +
262 2 * dx * invc2 * Complex(d.G2, d.D2)) -
263 4 * invGD * invc2 * mF0 * sq * Complex(d.G2, d.D2) +
264 2 * invGD * (2 * d.D0 - 3 * d.D2) * sq -
265 2 * pow2(invc2) * mF0 * Complex(d.G2, d.D2) +
266 invc2 * (2 * d.D0 - 3 * d.D2)))) /
267 (8 * sqrt_pi * invGD * mF0 * sq);
268 case CalcType::Voigt:
269 return -invGD *
270 (Complex(0, invGD) * dw1 *
271 (dx * (2 * d.D0 - 3 * d.D2) -
272 mF0 * Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2)) +
273 w1 * (2 * d.D0 - 3 * d.D2)) /
274 (2 * sqrt_pi * mF0);
275 case CalcType::LowXandHighY:
276 return -(4 * pow2(invGD) * (2 * d.D0 - 3 * d.D2) * (w1 - w2) * sq +
277 1i * (4 * dw1 * pow3(invGD) * sq *
278 (dx * (2 * d.D0 - 3 * d.D2) -
279 mF0 *
280 Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2)) +
281 dw2 * invc2 *
282 (2 * pow2(invGD) * mF0 *
283 (Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) -
284 2 * dx * invc2 * Complex(d.G2, d.D2)) -
285 4 * invGD * invc2 * mF0 * sq * Complex(d.G2, d.D2) +
286 2 * invGD * (2 * d.D0 - 3 * d.D2) * sq -
287 2 * pow2(invc2) * mF0 * Complex(d.G2, d.D2) +
288 invc2 * (2 * d.D0 - 3 * d.D2)))) /
289 (8 * sqrt_pi * invGD * mF0 * sq);
290 case CalcType::LowYandLowX:
291 return pow2(invc2) *
292 (4 * sq * Complex(d.G2, d.D2) * (sqrt_pi * w1 * sq - 1) -
293 sqrt_pi * (1i * dw1 * sq + w1) *
294 (Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) -
295 2 * dx * invc2 * Complex(d.G2, d.D2))) /
296 (2 * pi * sq);
297 case CalcType::LowYandHighX:
298 return -pow2(invc2) *
299 (x * (2 * x - 3) * Complex(d.G2, d.D2) +
300 (x - 3) * (Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) -
301 2 * dx * invc2 * Complex(d.G2, d.D2))) /
302 (2 * pi * pow3(x));
303 }
304 return {};
305}
306
307Complex SpeedDependentVoigt::dFdT(const Output &d, Numeric T) const noexcept {
308 switch (calcs) {
309 case CalcType::Full:
310 return (-pow2(invGD) * (w1 - w2) * sq *
311 (T * (2 * d.D0 - 3 * d.D2) * sqrt_ln_2 - pow3(invGD) * mF0) *
312 ln_16 +
313 1i * invc2 *
314 (dw1 *
315 (T * invGD * invc2 * mF0 * sq * Complex(d.G2, d.D2) *
316 ln_16 -
317 2 * invGD * sq *
318 (T * (2 * d.D0 - 3 * d.D2) * sqrt_ln_2 -
319 pow3(invGD) * mF0) *
320 sqrt_ln_2 -
321 (2 * T * pow2(invGD) * mF0 *
322 (Complex(3 * d.G2 - 2 * d.G0, 3 * d.D2 - 2 * d.D0) +
323 2 * dx * invc2 * Complex(d.G2, d.D2)) *
324 sqrt_ln_2 +
325 2 * T * pow2(invc2) * mF0 * Complex(d.G2, d.D2) *
326 sqrt_ln_2 +
327 invc2 * (T * (-2 * d.D0 + 3 * d.D2) * sqrt_ln_2 +
328 pow3(invGD) * mF0)) *
329 sqrt_ln_2) +
330 dw2 *
331 (T * invGD * invc2 * mF0 * sq * Complex(d.G2, d.D2) *
332 ln_16 +
333 2 * invGD * sq *
334 (T * (-2 * d.D0 + 3 * d.D2) * sqrt_ln_2 +
335 pow3(invGD) * mF0) *
336 sqrt_ln_2 +
337 (-2 * T * pow2(invGD) * mF0 *
338 (Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) -
339 2 * dx * invc2 * Complex(d.G2, d.D2)) *
340 sqrt_ln_2 +
341 2 * T * pow2(invc2) * mF0 * Complex(d.G2, d.D2) *
342 sqrt_ln_2 +
343 invc2 * (T * (-2 * d.D0 + 3 * d.D2) * sqrt_ln_2 +
344 pow3(invGD) * mF0)) *
345 sqrt_ln_2)) *
346 sqrt_ln_2) /
347 (8 * sqrt_pi * T * invGD * mF0 * sq * pow3(sqrt_ln_2));
348 case CalcType::Voigt:
349 return -invGD *
350 (-Complex(0, invGD) * dw1 *
351 (T * mF0 * Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) *
352 sqrt_ln_2 -
353 dx * (T * (2 * d.D0 - 3 * d.D2) * sqrt_ln_2 -
354 pow3(invGD) * mF0)) +
355 w1 * (T * (2 * d.D0 - 3 * d.D2) * sqrt_ln_2 - pow3(invGD) * mF0)) /
356 (2 * sqrt_pi * T * mF0 * sqrt_ln_2);
357 case CalcType::LowXandHighY:
358 return (-pow2(invGD) * (w1 - w2) * sq *
359 (T * (2 * d.D0 - 3 * d.D2) * sqrt_ln_2 - pow3(invGD) * mF0) *
360 ln_16 +
361 1i * (dw1 * pow3(invGD) * sq *
362 (T * mF0 *
363 Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) *
364 sqrt_ln_2 -
365 dx * (T * (2 * d.D0 - 3 * d.D2) * sqrt_ln_2 -
366 pow3(invGD) * mF0)) *
367 ln_16 +
368 dw2 * invc2 *
369 (T * invGD * invc2 * mF0 * sq * Complex(d.G2, d.D2) *
370 ln_16 +
371 2 * invGD * sq *
372 (T * (-2 * d.D0 + 3 * d.D2) * sqrt_ln_2 +
373 pow3(invGD) * mF0) *
374 sqrt_ln_2 +
375 (-2 * T * pow2(invGD) * mF0 *
376 (Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) -
377 2 * dx * invc2 * Complex(d.G2, d.D2)) *
378 sqrt_ln_2 +
379 2 * T * pow2(invc2) * mF0 * Complex(d.G2, d.D2) *
380 sqrt_ln_2 +
381 invc2 * (T * (-2 * d.D0 + 3 * d.D2) * sqrt_ln_2 +
382 pow3(invGD) * mF0)) *
383 sqrt_ln_2) *
384 sqrt_ln_2)) /
385 (8 * sqrt_pi * T * invGD * mF0 * sq * pow3(sqrt_ln_2));
386 case CalcType::LowYandLowX:
387 return pow2(invc2) *
388 (4 * sq * Complex(d.G2, d.D2) * (sqrt_pi * w1 * sq - 1) -
389 sqrt_pi * (1i * dw1 * sq + w1) *
390 (Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) -
391 2 * dx * invc2 * Complex(d.G2, d.D2))) /
392 (2 * pi * sq);
393 case CalcType::LowYandHighX:
394 return -pow2(invc2) *
395 (x * (2 * x - 3) * Complex(d.G2, d.D2) +
396 (x - 3) * (Complex(2 * d.G0 - 3 * d.G2, 2 * d.D0 - 3 * d.D2) -
397 2 * dx * invc2 * Complex(d.G2, d.D2))) /
398 (2 * pi * pow3(x));
399 }
400 return {};
401}
402
403Complex SpeedDependentVoigt::operator()(Numeric f) noexcept {
404 imag_val(dx) = mF0 - f;
405 x = dx * invc2;
406 update_calcs();
407 calc();
408 return F;
409}
410
411constexpr Numeric abs_squared(Complex z) noexcept {
412 return pow2(z.real()) + pow2(z.imag());
413}
414
415SpeedDependentVoigt::CalcType
416SpeedDependentVoigt::init(const Complex c2) const noexcept {
417 if (abs_squared(c2) == 0)
418 return CalcType::Voigt;
419 if (abs_squared(x) <= 9e-16 * abs_squared(sqrty * sqrty))
420 return CalcType::LowXandHighY;
421 if ((abs_squared(sqrty * sqrty) <= 1.e-30 * abs_squared(x)) and abs_squared(std::sqrt(x)) <= 16.e6)
422 return CalcType::LowYandLowX; // Weird case, untested
423 if ((abs_squared(sqrty * sqrty) <= 1.e-30 * abs_squared(x)))
424 return CalcType::LowYandHighX;
425 return CalcType::Full;
426}
427
428void SpeedDependentVoigt::update_calcs() noexcept {
429 if (calcs not_eq CalcType::Voigt)
430 calcs = init(Complex(1, 1));
431}
432
433void SpeedDependentVoigt::calc() noexcept {
434 switch (calcs) {
435 case CalcType::Full:
436 sq = std::sqrt(x + sqrty * sqrty);
437 w1 = Faddeeva::w(1i * (sq - sqrty));
438 w2 = Faddeeva::w(1i * (sq + sqrty));
439 F = inv_sqrt_pi * invGD * (w1 - w2);
440 dw1 = 2i * (inv_sqrt_pi - (sq - sqrty) * w1);
441 dw2 = 2i * (inv_sqrt_pi - (sq + sqrty) * w2);
442 break;
443 case CalcType::Voigt:
444 w1 = Faddeeva::w(1i * dx * invGD);
445 F = inv_sqrt_pi * invGD * w1;
446 dw1 = 2i * (inv_sqrt_pi - dx * invGD * w1);
447 break;
448 case CalcType::LowXandHighY:
449 sq = std::sqrt(x + sqrty * sqrty);
450 w1 = Faddeeva::w(1i * dx * invGD);
451 w2 = Faddeeva::w(1i * (sq + sqrty));
452 F = inv_sqrt_pi * invGD * (w1 - w2);
453 dw1 = 2i * (inv_sqrt_pi - dx * invGD * w1);
454 dw2 = 2i * (inv_sqrt_pi - (sq + sqrty) * w2);
455 break;
456 case CalcType::LowYandLowX:
457 sq = std::sqrt(x);
458 w1 = Faddeeva::w(1i * sq);
459 F = 2 * inv_pi * invc2 * (1 - sqrt_pi * sq * w1);
460 dw1 = 2i * (inv_sqrt_pi - sq * w1);
461 break;
462 case CalcType::LowYandHighX:
463 F = inv_pi * invc2 * (1 / x - 1.5 / pow2(x));
464 break;
465 }
466}
467
468HartmannTran::HartmannTran(Numeric F0_noshift, const Output &ls,
469 Numeric GD_div_F0, Numeric dZ) noexcept
470 : G0(ls.G0), D0(ls.D0), G2(ls.G2), D2(ls.D2), FVC(ls.FVC), ETA(ls.ETA),
471 mF0(F0_noshift + dZ + (1 - ls.ETA) * (ls.D0 - 1.5 * ls.D2)),
472 invGD(sqrt_ln_2 / nonstd::abs(GD_div_F0 * mF0)),
473 deltax(ls.FVC + (1 - ls.ETA) * (ls.G0 - 3 * ls.G2 / 2), mF0),
474 sqrty(1 / (2 * (1 - ls.ETA) * Complex(ls.G2, ls.D2) * invGD)) {
475 calc();
476}
477
478Complex HartmannTran::dFdf() const noexcept {
479 constexpr Complex ddeltax = -1i;
480 Complex dx = -ddeltax / ((ETA - 1) * Complex(G2, D2));
481 Complex dsqrtxy = dx / (2 * sqrtxy);
482
483 switch (calcs) {
484 case CalcType::Full: {
485 Complex dz1 = dsqrtxy;
486 Complex dz2 = dsqrtxy;
487 Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
488 Complex dB =
489 sqrt_pi *
490 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
491 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) /
492 (2 * sqrty * (ETA - 1) * Complex(G2, D2));
493 Complex dK = ETA * Complex(G2, D2) * dB +
494 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
495 return inv_pi * (-A * dK + K * dA) / pow2(K);
496 }
497 case CalcType::Noc2tLowZ: {
498 Complex dz1 = invGD * ddeltax;
499 Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
500 Complex dB =
501 -invGD *
502 (sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) - dz1);
503 Complex dK = ETA * Complex(G2, D2) * dB +
504 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
505 return inv_pi * (-A * dK + K * dA) / pow2(K);
506 }
507 case CalcType::Noc2tHighZ: {
508 Complex dz1 = invGD * ddeltax;
509 Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
510 Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
511 invGD * dz1 / (2 * pow2(z1)) +
512 9 * invGD * dz1 / (4 * pow4(z1));
513 Complex dK = ETA * Complex(G2, D2) * dB +
514 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
515 return inv_pi * (-A * dK + K * dA) / pow2(K);
516 }
517 case CalcType::LowXandHighY: {
518 Complex dz1 = invGD * ddeltax;
519 Complex dz2 = dsqrtxy;
520 Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
521 Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
522 invGD * dz1 / (2 * pow2(z1)) +
523 9 * invGD * dz1 / (4 * pow4(z1));
524 Complex dK = ETA * Complex(G2, D2) * dB +
525 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
526 return inv_pi * (-A * dK + K * dA) / pow2(K);
527 }
528 case CalcType::LowYandLowX: {
529 Complex dz1 = dsqrtxy;
530 Complex dz2 = dx / (2 * sqrtx);
531 Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
532 ((ETA - 1) * Complex(G2, D2));
533 Complex dB =
534 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
535 (2 * pow2(sqrty) + x - 1) +
536 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
537 2 * (sqrt_pi * w2 * z2 - 1) * dx) /
538 ((ETA - 1) * Complex(G2, D2));
539 Complex dK = ETA * Complex(G2, D2) * dB +
540 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
541 return inv_pi * (-A * dK + K * dA) / pow2(K);
542 }
543 case CalcType::LowYandHighX: {
544 Complex dz1 = dsqrtxy;
545 Complex dA = (x - 3) * dx / ((ETA - 1) * Complex(G2, D2) * pow3(x));
546 Complex dB =
547 (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) -
548 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx + (2 * x - 3) * x * dx / 2) /
549 ((ETA - 1) * Complex(G2, D2) * pow3(x));
550 Complex dK = ETA * Complex(G2, D2) * dB +
551 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
552 return inv_pi * (-A * dK + K * dA) / pow2(K);
553 }
554 }
555 return {};
556}
557
558Complex HartmannTran::dFdF0() const noexcept {
559 Numeric dGD = (1 / (invGD * mF0));
560 Numeric dinvGD = -dGD * pow2(invGD);
561 Complex dsqrty = dinvGD / (2 * (ETA - 1) * Complex(G2, D2) * pow2(invGD));
562 constexpr Complex ddeltax = 1i;
563 Complex dx = -ddeltax / ((ETA - 1) * Complex(G2, D2));
564 Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
565
566 switch (calcs) {
567 case CalcType::Full: {
568 Complex dz1 = dsqrtxy - dsqrty;
569 Complex dz2 = dsqrtxy + dsqrty;
570 Complex dA = sqrt_pi * ((w1 - w2) * dinvGD +
571 (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
572 Complex dB =
573 sqrt_pi *
574 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
575 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
576 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) *
577 sqrty) /
578 (2 * (ETA - 1) * Complex(G2, D2) * pow2(sqrty));
579 Complex dK = ETA * Complex(G2, D2) * dB +
580 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
581 return inv_pi * (-A * dK + K * dA) / pow2(K);
582 }
583 case CalcType::Noc2tLowZ: {
584 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
585 Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
586 Complex dB =
587 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
588 dz1) *
589 invGD -
590 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
591 Complex dK = ETA * Complex(G2, D2) * dB +
592 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
593 return inv_pi * (-A * dK + K * dA) / pow2(K);
594 }
595 case CalcType::Noc2tHighZ: {
596 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
597 Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
598 Complex dB =
599 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
600 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
601 9 * dz1) *
602 invGD) /
603 (4 * pow4(z1));
604 Complex dK = ETA * Complex(G2, D2) * dB +
605 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
606 return inv_pi * (-A * dK + K * dA) / pow2(K);
607 }
608 case CalcType::LowXandHighY: {
609 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
610 Complex dz2 = dsqrtxy + dsqrty;
611 Complex dA = sqrt_pi * ((w1 - w2) * dinvGD +
612 (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
613 Complex dB =
614 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
615 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
616 9 * dz1) *
617 invGD) /
618 (4 * pow4(z1));
619 Complex dK = ETA * Complex(G2, D2) * dB +
620 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
621 return inv_pi * (-A * dK + K * dA) / pow2(K);
622 }
623 case CalcType::LowYandLowX: {
624 Complex dz1 = dsqrtxy;
625 Complex dz2 = dx / (2 * sqrtx);
626 Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
627 ((ETA - 1) * Complex(G2, D2));
628 Complex dB =
629 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
630 (2 * pow2(sqrty) + x - 1) +
631 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
632 2 * (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) /
633 ((ETA - 1) * Complex(G2, D2));
634 Complex dK = ETA * Complex(G2, D2) * dB +
635 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
636 return inv_pi * (-A * dK + K * dA) / pow2(K);
637 }
638 case CalcType::LowYandHighX: {
639 Complex dz1 = dsqrtxy;
640 Complex dA = (x - 3) * dx / ((ETA - 1) * Complex(G2, D2) * pow3(x));
641 Complex dB = (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
642 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x / 2 -
643 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) /
644 ((ETA - 1) * Complex(G2, D2) * pow3(x));
645 Complex dK = ETA * Complex(G2, D2) * dB +
646 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
647 return inv_pi * (-A * dK + K * dA) / pow2(K);
648 }
649 }
650 return {};
651}
652
653Complex HartmannTran::dFdD0(Numeric dD0) const noexcept {
654 Numeric dmF0 = (1 - ETA) * dD0;
655 Numeric dGD = (dmF0 / (invGD * mF0));
656 Numeric dinvGD = -dGD * pow2(invGD);
657 Complex dsqrty = dinvGD / (2 * (ETA - 1) * Complex(G2, D2) * pow2(invGD));
658 Complex ddeltax = Complex(0, 1 - ETA) * dD0;
659 Complex dx = -ddeltax / ((ETA - 1) * Complex(G2, D2));
660 Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
661
662 switch (calcs) {
663 case CalcType::Full: {
664 Complex dz1 = dsqrtxy - dsqrty;
665 Complex dz2 = dsqrtxy + dsqrty;
666 Complex dA = sqrt_pi * ((w1 - w2) * dinvGD +
667 (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
668 Complex dB =
669 sqrt_pi *
670 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
671 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
672 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) *
673 sqrty) /
674 (2 * (ETA - 1) * Complex(G2, D2) * pow2(sqrty));
675 Complex dK = ETA * Complex(G2, D2) * dB + Complex(0, ETA * dD0) * A +
676 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
677 return inv_pi * (-A * dK + K * dA) / pow2(K);
678 }
679 case CalcType::Noc2tLowZ: {
680 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
681 Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
682 Complex dB =
683 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
684 dz1) *
685 invGD -
686 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
687 Complex dK = ETA * Complex(G2, D2) * dB + Complex(0, ETA * dD0) * A +
688 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
689 return inv_pi * (-A * dK + K * dA) / pow2(K);
690 }
691 case CalcType::Noc2tHighZ: {
692 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
693 Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
694 Complex dB =
695 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
696 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
697 9 * dz1) *
698 invGD) /
699 (4 * pow4(z1));
700 Complex dK = ETA * Complex(G2, D2) * dB + Complex(0, ETA * dD0) * A +
701 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
702 return inv_pi * (-A * dK + K * dA) / pow2(K);
703 }
704 case CalcType::LowXandHighY: {
705 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
706 Complex dz2 = dsqrtxy + dsqrty;
707 Complex dA = sqrt_pi * ((w1 - w2) * dinvGD +
708 (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
709 Complex dB =
710 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
711 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
712 9 * dz1) *
713 invGD) /
714 (4 * pow4(z1));
715 Complex dK = ETA * Complex(G2, D2) * dB + Complex(0, ETA * dD0) * A +
716 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
717 return inv_pi * (-A * dK + K * dA) / pow2(K);
718 }
719 case CalcType::LowYandLowX: {
720 Complex dz1 = dsqrtxy;
721 Complex dz2 = dx / (2 * sqrtx);
722 Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
723 ((ETA - 1) * Complex(G2, D2));
724 Complex dB =
725 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
726 (2 * pow2(sqrty) + x - 1) +
727 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
728 2 * (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) /
729 ((ETA - 1) * Complex(G2, D2));
730 Complex dK = ETA * Complex(G2, D2) * dB + Complex(0, ETA * dD0) * A +
731 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
732 return inv_pi * (-A * dK + K * dA) / pow2(K);
733 }
734 case CalcType::LowYandHighX: {
735 Complex dz1 = dsqrtxy;
736 Complex dA = (x - 3) * dx / ((ETA - 1) * Complex(G2, D2) * pow3(x));
737 Complex dB = (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
738 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x / 2 -
739 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) /
740 ((ETA - 1) * Complex(G2, D2) * pow3(x));
741 Complex dK = ETA * Complex(G2, D2) * dB + Complex(0, ETA * dD0) * A +
742 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
743 return inv_pi * (-A * dK + K * dA) / pow2(K);
744 }
745 }
746 return {};
747}
748
749Complex HartmannTran::dFdG0(Numeric dG0) const noexcept {
750 Numeric ddeltax = (1 - ETA) * dG0;
751 Complex dx = -ddeltax / ((ETA - 1) * Complex(G2, D2));
752 Complex dsqrtxy = dx / (2 * sqrtxy);
753
754 switch (calcs) {
755 case CalcType::Full: {
756 Complex dz1 = dsqrtxy;
757 Complex dz2 = dsqrtxy;
758 Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
759 Complex dB =
760 sqrt_pi *
761 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
762 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) /
763 (2 * sqrty * (ETA - 1) * Complex(G2, D2));
764 Complex dK = ETA * Complex(G2, D2) * dB + ETA * A * dG0 +
765 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
766 return inv_pi * (-A * dK + K * dA) / pow2(K);
767 }
768 case CalcType::Noc2tLowZ: {
769 Complex dz1 = invGD * ddeltax;
770 Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
771 Complex dB =
772 -invGD *
773 (sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) - dz1);
774 Complex dK = ETA * Complex(G2, D2) * dB + ETA * A * dG0 +
775 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
776 return inv_pi * (-A * dK + K * dA) / pow2(K);
777 }
778 case CalcType::Noc2tHighZ: {
779 Complex dz1 = invGD * ddeltax;
780 Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
781 Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
782 invGD * dz1 / (2 * pow2(z1)) +
783 9 * invGD * dz1 / (4 * pow4(z1));
784 Complex dK = ETA * Complex(G2, D2) * dB + ETA * A * dG0 +
785 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
786 return inv_pi * (-A * dK + K * dA) / pow2(K);
787 }
788 case CalcType::LowXandHighY: {
789 Complex dz1 = invGD * ddeltax;
790 Complex dz2 = dsqrtxy;
791 Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
792 Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
793 invGD * dz1 / (2 * pow2(z1)) +
794 9 * invGD * dz1 / (4 * pow4(z1));
795 Complex dK = ETA * Complex(G2, D2) * dB + ETA * A * dG0 +
796 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
797 return inv_pi * (-A * dK + K * dA) / pow2(K);
798 }
799 case CalcType::LowYandLowX: {
800 Complex dz1 = dsqrtxy;
801 Complex dz2 = dx / (2 * sqrtx);
802 Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
803 ((ETA - 1) * Complex(G2, D2));
804 Complex dB =
805 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
806 (2 * pow2(sqrty) + x - 1) +
807 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
808 2 * (sqrt_pi * w2 * z2 - 1) * dx) /
809 ((ETA - 1) * Complex(G2, D2));
810 Complex dK = ETA * Complex(G2, D2) * dB + ETA * A * dG0 +
811 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
812 return inv_pi * (-A * dK + K * dA) / pow2(K);
813 }
814 case CalcType::LowYandHighX: {
815 Complex dz1 = dsqrtxy;
816 Complex dA = (x - 3) * dx / ((ETA - 1) * Complex(G2, D2) * pow3(x));
817 Complex dB =
818 (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) -
819 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx + (2 * x - 3) * x * dx / 2) /
820 ((ETA - 1) * Complex(G2, D2) * pow3(x));
821 Complex dK = ETA * Complex(G2, D2) * dB + ETA * A * dG0 +
822 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
823 return inv_pi * (-A * dK + K * dA) / pow2(K);
824 }
825 }
826 return {};
827}
828
829Complex HartmannTran::dFdD2(Numeric dD2) const noexcept {
830 Numeric dmF0 = -3 * (1 - ETA) * dD2 / 2;
831 Numeric dGD = (dmF0 / (invGD * mF0));
832 Numeric dinvGD = -dGD * pow2(invGD);
833 Complex dsqrty = (Complex(G2, D2) * dinvGD + Complex(0, invGD) * dD2) /
834 (2 * (ETA - 1) * pow2(Complex(G2, D2)) * pow2(invGD));
835 Complex ddeltax = 1.5 * Complex(0, ETA - 1) * dD2;
836 Complex dx = (-Complex(G2, D2) * ddeltax + Complex(0, dD2) * deltax) /
837 ((ETA - 1) * pow2(Complex(G2, D2)));
838 Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
839
840 switch (calcs) {
841 case CalcType::Full: {
842 Complex dz1 = dsqrtxy - dsqrty;
843 Complex dz2 = dsqrtxy + dsqrty;
844 Complex dA = sqrt_pi * ((w1 - w2) * dinvGD +
845 (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
846 Complex dB = (sqrt_pi * Complex(G2, D2) *
847 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
848 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
849 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
850 2 * w2 * z2 * dz2) *
851 sqrty) -
852 1i *
853 (sqrt_pi * (pow2(z1) - 1) * w1 -
854 sqrt_pi * (pow2(z2) - 1) * w2 + 2 * sqrty) *
855 sqrty * dD2) /
856 (2 * (ETA - 1) * pow2(Complex(G2, D2)) * pow2(sqrty));
857 Complex dK = ETA * Complex(G2, D2) * dB - Complex(0, 1.5 * dD2 * ETA) * A +
858 Complex(0, dD2 * ETA) * B +
859 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
860 return inv_pi * (-A * dK + K * dA) / pow2(K);
861 }
862 case CalcType::Noc2tLowZ: {
863 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
864 Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
865 Complex dB =
866 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
867 dz1) *
868 invGD -
869 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
870 Complex dK = ETA * Complex(G2, D2) * dB - Complex(0, 1.5 * dD2 * ETA) * A +
871 Complex(0, dD2 * ETA) * B +
872 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
873 return inv_pi * (-A * dK + K * dA) / pow2(K);
874 }
875 case CalcType::Noc2tHighZ: {
876 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
877 Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
878 Complex dB =
879 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
880 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
881 9 * dz1) *
882 invGD) /
883 (4 * pow4(z1));
884 Complex dK = ETA * Complex(G2, D2) * dB - Complex(0, 1.5 * dD2 * ETA) * A +
885 Complex(0, dD2 * ETA) * B +
886 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
887 return inv_pi * (-A * dK + K * dA) / pow2(K);
888 }
889 case CalcType::LowXandHighY: {
890 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
891 Complex dz2 = dsqrtxy + dsqrty;
892 Complex dA = sqrt_pi * ((w1 - w2) * dinvGD +
893 (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
894 Complex dB =
895 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
896 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
897 9 * dz1) *
898 invGD) /
899 (4 * pow4(z1));
900 Complex dK = ETA * Complex(G2, D2) * dB - Complex(0, 1.5 * dD2 * ETA) * A +
901 Complex(0, dD2 * ETA) * B +
902 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
903 return inv_pi * (-A * dK + K * dA) / pow2(K);
904 }
905 case CalcType::LowYandLowX: {
906 Complex dz1 = dsqrtxy;
907 Complex dz2 = dx / (2 * sqrtx);
908 Complex dA = 2 *
909 (sqrt_pi * Complex(G2, D2) * (w2 * dz2 + z2 * 1i * dw2 * dz2) -
910 1i * (sqrt_pi * w2 * z2 - 1) * dD2) /
911 ((ETA - 1) * pow2(Complex(G2, D2)));
912 Complex dB =
913 (-2 * Complex(G2, D2) *
914 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
915 (2 * pow2(sqrty) + x - 1) +
916 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
917 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
918 1i *
919 (2 * sqrt_pi * w1 * z1 +
920 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
921 dD2) /
922 ((ETA - 1) * pow2(Complex(G2, D2)));
923 Complex dK = ETA * Complex(G2, D2) * dB - Complex(0, 1.5 * dD2 * ETA) * A +
924 Complex(0, dD2 * ETA) * B +
925 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
926 return inv_pi * (-A * dK + K * dA) / pow2(K);
927 }
928 case CalcType::LowYandHighX: {
929 Complex dz1 = dsqrtxy;
930 Complex dA =
931 (Complex(G2, D2) * (x - 3) * dx + 1i * (2 * x - 3) * x * dD2 / 2) /
932 ((ETA - 1) * pow2(Complex(G2, D2)) * pow3(x));
933 Complex dB =
934 (Complex(G2, D2) *
935 (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
936 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
937 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) -
938 1i *
939 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
940 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
941 x * dD2) /
942 (2 * (ETA - 1) * pow2(Complex(G2, D2)) * pow3(x));
943 Complex dK = ETA * Complex(G2, D2) * dB - Complex(0, 1.5 * dD2 * ETA) * A +
944 Complex(0, dD2 * ETA) * B +
945 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
946 return inv_pi * (-A * dK + K * dA) / pow2(K);
947 }
948 }
949 return {};
950}
951
952Complex HartmannTran::dFdG2(Numeric dG2) const noexcept {
953 Complex dsqrty = dG2 / (2 * invGD * (ETA - 1) * pow2(Complex(G2, D2)));
954 Numeric ddeltax = 3 * (ETA - 1) * dG2 / 2;
955 Complex dx = (-Complex(G2, D2) * ddeltax + deltax * dG2) /
956 ((ETA - 1) * pow2(Complex(G2, D2)));
957 Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
958
959 switch (calcs) {
960 case CalcType::Full: {
961 Complex dz1 = dsqrtxy - dsqrty;
962 Complex dz2 = dsqrtxy + dsqrty;
963 Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
964 Complex dB = (sqrt_pi * Complex(G2, D2) *
965 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
966 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
967 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
968 2 * w2 * z2 * dz2) *
969 sqrty) -
970 (sqrt_pi * (pow2(z1) - 1) * w1 -
971 sqrt_pi * (pow2(z2) - 1) * w2 + 2 * sqrty) *
972 sqrty * dG2) /
973 (2 * (ETA - 1) * pow2(Complex(G2, D2)) * pow2(sqrty));
974 Complex dK = ETA * Complex(G2, D2) * dB - 3 * ETA * A * dG2 / 2 +
975 ETA * B * dG2 +
976 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
977 return inv_pi * (-A * dK + K * dA) / pow2(K);
978 }
979 case CalcType::Noc2tLowZ: {
980 Complex dz1 = invGD * ddeltax;
981 Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
982 Complex dB =
983 -invGD *
984 (sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) - dz1);
985 Complex dK = ETA * Complex(G2, D2) * dB - 3 * ETA * A * dG2 / 2 +
986 ETA * B * dG2 +
987 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
988 return inv_pi * (-A * dK + K * dA) / pow2(K);
989 }
990 case CalcType::Noc2tHighZ: {
991 Complex dz1 = invGD * ddeltax;
992 Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
993 Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
994 invGD * dz1 / (2 * pow2(z1)) +
995 9 * invGD * dz1 / (4 * pow4(z1));
996 Complex dK = ETA * Complex(G2, D2) * dB - 3 * ETA * A * dG2 / 2 +
997 ETA * B * dG2 +
998 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
999 return inv_pi * (-A * dK + K * dA) / pow2(K);
1000 }
1001 case CalcType::LowXandHighY: {
1002 Complex dz1 = invGD * ddeltax;
1003 Complex dz2 = dsqrtxy + dsqrty;
1004 Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
1005 Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
1006 invGD * dz1 / (2 * pow2(z1)) +
1007 9 * invGD * dz1 / (4 * pow4(z1));
1008 Complex dK = ETA * Complex(G2, D2) * dB - 3 * ETA * A * dG2 / 2 +
1009 ETA * B * dG2 +
1010 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1011 return inv_pi * (-A * dK + K * dA) / pow2(K);
1012 }
1013 case CalcType::LowYandLowX: {
1014 Complex dz1 = dsqrtxy;
1015 Complex dz2 = dx / (2 * sqrtx);
1016 Complex dA = 2 *
1017 (sqrt_pi * Complex(G2, D2) * (w2 * dz2 + z2 * 1i * dw2 * dz2) -
1018 (sqrt_pi * w2 * z2 - 1) * dG2) /
1019 ((ETA - 1) * pow2(Complex(G2, D2)));
1020 Complex dB =
1021 (-2 * Complex(G2, D2) *
1022 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1023 (2 * pow2(sqrty) + x - 1) +
1024 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
1025 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
1026 (2 * sqrt_pi * w1 * z1 +
1027 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
1028 dG2) /
1029 ((ETA - 1) * pow2(Complex(G2, D2)));
1030 Complex dK = ETA * Complex(G2, D2) * dB - 3 * ETA * A * dG2 / 2 +
1031 ETA * B * dG2 +
1032 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1033 return inv_pi * (-A * dK + K * dA) / pow2(K);
1034 }
1035 case CalcType::LowYandHighX: {
1036 Complex dz1 = dsqrtxy;
1037 Complex dA = (Complex(G2, D2) * (x - 3) * dx + (2 * x - 3) * x * dG2 / 2) /
1038 ((ETA - 1) * pow2(Complex(G2, D2)) * pow3(x));
1039 Complex dB =
1040 (Complex(G2, D2) *
1041 (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
1042 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
1043 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) -
1044 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1045 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1046 x * dG2) /
1047 (2 * (ETA - 1) * pow2(Complex(G2, D2)) * pow3(x));
1048 Complex dK = ETA * Complex(G2, D2) * dB - 3 * ETA * A * dG2 / 2 +
1049 ETA * B * dG2 +
1050 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1051 return inv_pi * (-A * dK + K * dA) / pow2(K);
1052 }
1053 }
1054 return {};
1055}
1056
1057Complex HartmannTran::dFdFVC(Numeric dFVC) const noexcept {
1058 Numeric ddeltax = dFVC;
1059 Complex dx = -ddeltax / ((ETA - 1) * Complex(G2, D2));
1060 Complex dsqrtxy = dx / (2 * sqrtxy);
1061
1062 switch (calcs) {
1063 case CalcType::Full: {
1064 Complex dz1 = dsqrtxy;
1065 Complex dz2 = dsqrtxy;
1066 Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
1067 Complex dB =
1068 sqrt_pi *
1069 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
1070 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) /
1071 (2 * sqrty * (ETA - 1) * Complex(G2, D2));
1072 Complex dK = ETA * Complex(G2, D2) * dB +
1073 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA -
1074 A * dFVC;
1075 return inv_pi * (-A * dK + K * dA) / pow2(K);
1076 }
1077 case CalcType::Noc2tLowZ: {
1078 Complex dz1 = invGD * ddeltax;
1079 Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
1080 Complex dB =
1081 -invGD *
1082 (sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) - dz1);
1083 Complex dK = ETA * Complex(G2, D2) * dB +
1084 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA -
1085 A * dFVC;
1086 return inv_pi * (-A * dK + K * dA) / pow2(K);
1087 }
1088 case CalcType::Noc2tHighZ: {
1089 Complex dz1 = invGD * ddeltax;
1090 Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
1091 Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
1092 invGD * dz1 / (2 * pow2(z1)) +
1093 9 * invGD * dz1 / (4 * pow4(z1));
1094 Complex dK = ETA * Complex(G2, D2) * dB +
1095 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA -
1096 A * dFVC;
1097 return inv_pi * (-A * dK + K * dA) / pow2(K);
1098 }
1099 case CalcType::LowXandHighY: {
1100 Complex dz1 = invGD * ddeltax;
1101 Complex dz2 = dsqrtxy;
1102 Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
1103 Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
1104 invGD * dz1 / (2 * pow2(z1)) +
1105 9 * invGD * dz1 / (4 * pow4(z1));
1106 Complex dK = ETA * Complex(G2, D2) * dB +
1107 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA -
1108 A * dFVC;
1109 return inv_pi * (-A * dK + K * dA) / pow2(K);
1110 }
1111 case CalcType::LowYandLowX: {
1112 Complex dz1 = dsqrtxy;
1113 Complex dz2 = dx / (2 * sqrtx);
1114 Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
1115 ((ETA - 1) * Complex(G2, D2));
1116 Complex dB =
1117 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1118 (2 * pow2(sqrty) + x - 1) +
1119 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
1120 2 * (sqrt_pi * w2 * z2 - 1) * dx) /
1121 ((ETA - 1) * Complex(G2, D2));
1122 Complex dK = ETA * Complex(G2, D2) * dB +
1123 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA -
1124 A * dFVC;
1125 return inv_pi * (-A * dK + K * dA) / pow2(K);
1126 }
1127 case CalcType::LowYandHighX: {
1128 Complex dz1 = dsqrtxy;
1129 Complex dA = (x - 3) * dx / ((ETA - 1) * Complex(G2, D2) * pow3(x));
1130 Complex dB =
1131 (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) -
1132 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx + (2 * x - 3) * x * dx / 2) /
1133 ((ETA - 1) * Complex(G2, D2) * pow3(x));
1134 Complex dK = ETA * Complex(G2, D2) * dB +
1135 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA -
1136 A * dFVC;
1137 return inv_pi * (-A * dK + K * dA) / pow2(K);
1138 }
1139 }
1140 return {};
1141}
1142
1143Complex HartmannTran::dFdETA(Numeric dETA) const noexcept {
1144 Numeric dmF0 = -(D0 - 3 * D2 / 2) * dETA;
1145 Numeric dGD = (dmF0 / (invGD * mF0));
1146 Numeric dinvGD = -dGD * pow2(invGD);
1147 Complex dsqrty = ((ETA - 1) * dinvGD + invGD * dETA) /
1148 (2 * Complex(G2, D2) * pow2(ETA - 1) * pow2(invGD));
1149 Complex ddeltax = -dETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2);
1150 Complex dx = (-(ETA - 1) * ddeltax + deltax * dETA) /
1151 (Complex(G2, D2) * pow2(ETA - 1));
1152 Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
1153
1154 switch (calcs) {
1155 case CalcType::Full: {
1156 Complex dz1 = dsqrtxy - dsqrty;
1157 Complex dz2 = dsqrtxy + dsqrty;
1158 Complex dA = sqrt_pi * ((w1 - w2) * dinvGD +
1159 (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1160 Complex dB = (sqrt_pi *
1161 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
1162 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
1163 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
1164 2 * w2 * z2 * dz2) *
1165 sqrty) *
1166 (ETA - 1) -
1167 (sqrt_pi * (pow2(z1) - 1) * w1 -
1168 sqrt_pi * (pow2(z2) - 1) * w2 + 2 * sqrty) *
1169 sqrty * dETA) /
1170 (2 * Complex(G2, D2) * pow2(ETA - 1) * pow2(sqrty));
1171 Complex dK = (-FVC + Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA) * dA +
1172 Complex(G2, D2) * B * dETA + Complex(G2, D2) * ETA * dB -
1173 Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * A * dETA;
1174 return inv_pi * (-A * dK + K * dA) / pow2(K);
1175 }
1176 case CalcType::Noc2tLowZ: {
1177 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1178 Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1179 Complex dB =
1180 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
1181 dz1) *
1182 invGD -
1183 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
1184 Complex dK = (-FVC + Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA) * dA +
1185 Complex(G2, D2) * B * dETA + Complex(G2, D2) * ETA * dB -
1186 Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * A * dETA;
1187 return inv_pi * (-A * dK + K * dA) / pow2(K);
1188 }
1189 case CalcType::Noc2tHighZ: {
1190 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1191 Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1192 Complex dB =
1193 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1194 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1195 9 * dz1) *
1196 invGD) /
1197 (4 * pow4(z1));
1198 Complex dK = (-FVC + Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA) * dA +
1199 Complex(G2, D2) * B * dETA + Complex(G2, D2) * ETA * dB -
1200 Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * A * dETA;
1201 return inv_pi * (-A * dK + K * dA) / pow2(K);
1202 }
1203 case CalcType::LowXandHighY: {
1204 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1205 Complex dz2 = dsqrtxy + dsqrty;
1206 Complex dA = sqrt_pi * ((w1 - w2) * dinvGD +
1207 (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1208 Complex dB =
1209 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1210 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1211 9 * dz1) *
1212 invGD) /
1213 (4 * pow4(z1));
1214 Complex dK = (-FVC + Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA) * dA +
1215 Complex(G2, D2) * B * dETA + Complex(G2, D2) * ETA * dB -
1216 Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * A * dETA;
1217 return inv_pi * (-A * dK + K * dA) / pow2(K);
1218 }
1219 case CalcType::LowYandLowX: {
1220 Complex dz1 = dsqrtxy;
1221 Complex dz2 = dx / (2 * sqrtx);
1222 Complex dA = 2 *
1223 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) * (ETA - 1) -
1224 (sqrt_pi * w2 * z2 - 1) * dETA) /
1225 (Complex(G2, D2) * pow2(ETA - 1));
1226 Complex dB =
1227 (-2 * (ETA - 1) *
1228 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1229 (2 * pow2(sqrty) + x - 1) +
1230 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
1231 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
1232 (2 * sqrt_pi * w1 * z1 +
1233 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
1234 dETA) /
1235 (Complex(G2, D2) * pow2(ETA - 1));
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::LowYandHighX: {
1242 Complex dz1 = dsqrtxy;
1243 Complex dA = ((ETA - 1) * (x - 3) * dx + (2 * x - 3) * x * dETA / 2) /
1244 (Complex(G2, D2) * pow2(ETA - 1) * pow3(x));
1245 Complex dB = (-(2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1246 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1247 x * dETA +
1248 (ETA - 1) * (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) *
1249 pow3(x) +
1250 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
1251 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx)) /
1252 (2 * Complex(G2, D2) * pow2(ETA - 1) * pow3(x));
1253 Complex dK = (-FVC + Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA) * dA +
1254 Complex(G2, D2) * B * dETA + Complex(G2, D2) * ETA * dB -
1255 Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * A * dETA;
1256 return inv_pi * (-A * dK + K * dA) / pow2(K);
1257 }
1258 }
1259 return {};
1260}
1261
1262Complex HartmannTran::dFdH(Numeric dZ) const noexcept {
1263 Numeric dmF0 = dZ;
1264 Numeric dGD = (dmF0 / (invGD * mF0));
1265 Numeric dinvGD = -dGD * pow2(invGD);
1266 Complex dsqrty = dinvGD / (2 * (ETA - 1) * Complex(G2, D2) * pow2(invGD));
1267 Complex ddeltax = Complex(0, dZ);
1268 Complex dx = -ddeltax / ((ETA - 1) * Complex(G2, D2));
1269 Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
1270
1271 switch (calcs) {
1272 case CalcType::Full: {
1273 Complex dz1 = dsqrtxy - dsqrty;
1274 Complex dz2 = dsqrtxy + dsqrty;
1275 Complex dA = sqrt_pi * ((w1 - w2) * dinvGD +
1276 (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1277 Complex dB =
1278 sqrt_pi *
1279 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
1280 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
1281 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) *
1282 sqrty) /
1283 (2 * (ETA - 1) * Complex(G2, D2) * pow2(sqrty));
1284 Complex dK = ETA * Complex(G2, D2) * dB +
1285 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1286 return inv_pi * (-A * dK + K * dA) / pow2(K);
1287 }
1288 case CalcType::Noc2tLowZ: {
1289 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1290 Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1291 Complex dB =
1292 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
1293 dz1) *
1294 invGD -
1295 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
1296 Complex dK = ETA * Complex(G2, D2) * dB +
1297 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1298 return inv_pi * (-A * dK + K * dA) / pow2(K);
1299 }
1300 case CalcType::Noc2tHighZ: {
1301 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1302 Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1303 Complex dB =
1304 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1305 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1306 9 * dz1) *
1307 invGD) /
1308 (4 * pow4(z1));
1309 Complex dK = ETA * Complex(G2, D2) * dB +
1310 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1311 return inv_pi * (-A * dK + K * dA) / pow2(K);
1312 }
1313 case CalcType::LowXandHighY: {
1314 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1315 Complex dz2 = dsqrtxy + dsqrty;
1316 Complex dA = sqrt_pi * ((w1 - w2) * dinvGD +
1317 (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1318 Complex dB =
1319 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1320 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1321 9 * dz1) *
1322 invGD) /
1323 (4 * pow4(z1));
1324 Complex dK = ETA * Complex(G2, D2) * dB +
1325 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1326 return inv_pi * (-A * dK + K * dA) / pow2(K);
1327 }
1328 case CalcType::LowYandLowX: {
1329 Complex dz1 = dsqrtxy;
1330 Complex dz2 = dx / (2 * sqrtx);
1331 Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
1332 ((ETA - 1) * Complex(G2, D2));
1333 Complex dB =
1334 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1335 (2 * pow2(sqrty) + x - 1) +
1336 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
1337 2 * (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) /
1338 ((ETA - 1) * Complex(G2, D2));
1339 Complex dK = ETA * Complex(G2, D2) * dB +
1340 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1341 return inv_pi * (-A * dK + K * dA) / pow2(K);
1342 }
1343 case CalcType::LowYandHighX: {
1344 Complex dz1 = dsqrtxy;
1345 Complex dA = (x - 3) * dx / ((ETA - 1) * Complex(G2, D2) * pow3(x));
1346 Complex dB = (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
1347 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x / 2 -
1348 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) /
1349 ((ETA - 1) * Complex(G2, D2) * pow3(x));
1350 Complex dK = ETA * Complex(G2, D2) * dB +
1351 (ETA * Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) - FVC) * dA;
1352 return inv_pi * (-A * dK + K * dA) / pow2(K);
1353 }
1354 }
1355 return {};
1356}
1357
1358Complex HartmannTran::dFdVMR(const Output &d) const noexcept {
1359 Numeric dmF0 = (1 - ETA) * (d.D0 - 3 * d.D2 / 2) - (D0 - 3 * D2 / 2) * d.ETA;
1360 Numeric dGD = (dmF0 / (invGD * mF0));
1361 Numeric dinvGD = -dGD * pow2(invGD);
1362 Complex dsqrty =
1363 (Complex(G2, D2) * (ETA - 1) * dinvGD + Complex(G2, D2) * invGD * d.ETA +
1364 Complex(d.G2, d.D2) * (ETA - 1) * invGD) /
1365 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow2(invGD));
1366 Complex ddeltax = -(ETA - 1) * Complex(d.G0 - 1.5 * d.G2, d.D0 - 1.5 * d.D2) -
1367 Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * d.ETA + d.FVC;
1368 Complex dx = (-Complex(G2, D2) * (ETA - 1) * ddeltax +
1369 Complex(G2, D2) * deltax * d.ETA +
1370 Complex(d.G2, d.D2) * (ETA - 1) * deltax) /
1371 (pow2(Complex(G2, D2)) * pow2(ETA - 1));
1372 Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
1373
1374 switch (calcs) {
1375 case CalcType::Full: {
1376 Complex dz1 = dsqrtxy - dsqrty;
1377 Complex dz2 = dsqrtxy + dsqrty;
1378 Complex dA = sqrt_pi * ((w1 - w2) * dinvGD +
1379 (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1380 Complex dB = (sqrt_pi * Complex(G2, D2) *
1381 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
1382 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
1383 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
1384 2 * w2 * z2 * dz2) *
1385 sqrty) *
1386 (ETA - 1) -
1387 Complex(G2, D2) *
1388 (sqrt_pi * (pow2(z1) - 1) * w1 -
1389 sqrt_pi * (pow2(z2) - 1) * w2 + 2 * sqrty) *
1390 sqrty * d.ETA -
1391 Complex(d.G2, d.D2) * (ETA - 1) *
1392 (sqrt_pi * (pow2(z1) - 1) * w1 -
1393 sqrt_pi * (pow2(z2) - 1) * w2 + 2 * sqrty) *
1394 sqrty) /
1395 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow2(sqrty));
1396 Complex dK = Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1397 Complex(d.G2, d.D2) * B * ETA +
1398 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1399 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1400 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1401 A;
1402 return inv_pi * (-A * dK + K * dA) / pow2(K);
1403 }
1404 case CalcType::Noc2tLowZ: {
1405 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1406 Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1407 Complex dB =
1408 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
1409 dz1) *
1410 invGD -
1411 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
1412 Complex dK = Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1413 Complex(d.G2, d.D2) * B * ETA +
1414 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1415 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1416 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1417 A;
1418 return inv_pi * (-A * dK + K * dA) / pow2(K);
1419 }
1420 case CalcType::Noc2tHighZ: {
1421 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1422 Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1423 Complex dB =
1424 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1425 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1426 9 * dz1) *
1427 invGD) /
1428 (4 * pow4(z1));
1429 Complex dK = Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1430 Complex(d.G2, d.D2) * B * ETA +
1431 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1432 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1433 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1434 A;
1435 return inv_pi * (-A * dK + K * dA) / pow2(K);
1436 }
1437 case CalcType::LowXandHighY: {
1438 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1439 Complex dz2 = dsqrtxy + dsqrty;
1440 Complex dA = sqrt_pi * ((w1 - w2) * dinvGD +
1441 (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1442 Complex dB =
1443 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1444 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1445 9 * dz1) *
1446 invGD) /
1447 (4 * pow4(z1));
1448 Complex dK = Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1449 Complex(d.G2, d.D2) * B * ETA +
1450 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1451 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1452 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1453 A;
1454 return inv_pi * (-A * dK + K * dA) / pow2(K);
1455 }
1456 case CalcType::LowYandLowX: {
1457 Complex dz1 = dsqrtxy;
1458 Complex dz2 = dx / (2 * sqrtx);
1459 Complex dA = 2 *
1460 (sqrt_pi * Complex(G2, D2) * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1461 (ETA - 1) -
1462 Complex(G2, D2) * (sqrt_pi * w2 * z2 - 1) * d.ETA -
1463 Complex(d.G2, d.D2) * (sqrt_pi * w2 * z2 - 1) * (ETA - 1)) /
1464 (pow2(Complex(G2, D2)) * pow2(ETA - 1));
1465 Complex dB =
1466 (-2 * Complex(G2, D2) * (ETA - 1) *
1467 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1468 (2 * pow2(sqrty) + x - 1) +
1469 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
1470 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
1471 Complex(G2, D2) *
1472 (2 * sqrt_pi * w1 * z1 +
1473 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
1474 d.ETA +
1475 Complex(d.G2, d.D2) * (ETA - 1) *
1476 (2 * sqrt_pi * w1 * z1 +
1477 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1)) /
1478 (pow2(Complex(G2, D2)) * pow2(ETA - 1));
1479 Complex dK = Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1480 Complex(d.G2, d.D2) * B * ETA +
1481 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1482 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1483 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1484 A;
1485 return inv_pi * (-A * dK + K * dA) / pow2(K);
1486 }
1487 case CalcType::LowYandHighX: {
1488 Complex dz1 = dsqrtxy;
1489 Complex dA = (2 * Complex(G2, D2) * (ETA - 1) * (x - 3) * dx +
1490 Complex(G2, D2) * (2 * x - 3) * x * d.ETA +
1491 Complex(d.G2, d.D2) * (ETA - 1) * (2 * x - 3) * x) /
1492 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow3(x));
1493 Complex dB =
1494 (-Complex(G2, D2) *
1495 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1496 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1497 x * d.ETA +
1498 Complex(G2, D2) * (ETA - 1) *
1499 (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
1500 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
1501 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) -
1502 Complex(d.G2, d.D2) *
1503 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1504 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1505 (ETA - 1) * x) /
1506 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow3(x));
1507 Complex dK = 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 }
1516 return {};
1517}
1518
1519Complex HartmannTran::dFdT(const Output &d, Numeric T) const noexcept {
1520 Numeric dmF0 = (1 - ETA) * (d.D0 - 3 * d.D2 / 2) - (D0 - 3 * D2 / 2) * d.ETA;
1521 Numeric dGD = (dmF0 / (invGD * mF0)) - invGD * invGD / (2 * T * sqrt_ln_2);
1522 Numeric dinvGD = -dGD * pow2(invGD);
1523 Complex dsqrty =
1524 (Complex(G2, D2) * (ETA - 1) * dinvGD + Complex(G2, D2) * invGD * d.ETA +
1525 Complex(d.G2, d.D2) * (ETA - 1) * invGD) /
1526 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow2(invGD));
1527 Complex ddeltax = -(ETA - 1) * Complex(d.G0 - 1.5 * d.G2, d.D0 - 1.5 * d.D2) -
1528 Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * d.ETA + d.FVC;
1529 Complex dx = (-Complex(G2, D2) * (ETA - 1) * ddeltax +
1530 Complex(G2, D2) * deltax * d.ETA +
1531 Complex(d.G2, d.D2) * (ETA - 1) * deltax) /
1532 (pow2(Complex(G2, D2)) * pow2(ETA - 1));
1533 Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
1534
1535 switch (calcs) {
1536 case CalcType::Full: {
1537 Complex dz1 = dsqrtxy - dsqrty;
1538 Complex dz2 = dsqrtxy + dsqrty;
1539 Complex dA = sqrt_pi * ((w1 - w2) * dinvGD +
1540 (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1541 Complex dB = (sqrt_pi * Complex(G2, D2) *
1542 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
1543 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
1544 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
1545 2 * w2 * z2 * dz2) *
1546 sqrty) *
1547 (ETA - 1) -
1548 Complex(G2, D2) *
1549 (sqrt_pi * (pow2(z1) - 1) * w1 -
1550 sqrt_pi * (pow2(z2) - 1) * w2 + 2 * sqrty) *
1551 sqrty * d.ETA -
1552 Complex(d.G2, d.D2) * (ETA - 1) *
1553 (sqrt_pi * (pow2(z1) - 1) * w1 -
1554 sqrt_pi * (pow2(z2) - 1) * w2 + 2 * sqrty) *
1555 sqrty) /
1556 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow2(sqrty));
1557 Complex dK = Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1558 Complex(d.G2, d.D2) * B * ETA +
1559 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1560 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1561 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1562 A;
1563 return inv_pi * (-A * dK + K * dA) / pow2(K);
1564 }
1565 case CalcType::Noc2tLowZ: {
1566 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1567 Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1568 Complex dB =
1569 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
1570 dz1) *
1571 invGD -
1572 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
1573 Complex dK = Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1574 Complex(d.G2, d.D2) * B * ETA +
1575 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1576 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1577 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1578 A;
1579 return inv_pi * (-A * dK + K * dA) / pow2(K);
1580 }
1581 case CalcType::Noc2tHighZ: {
1582 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1583 Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1584 Complex dB =
1585 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1586 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1587 9 * dz1) *
1588 invGD) /
1589 (4 * pow4(z1));
1590 Complex dK = Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1591 Complex(d.G2, d.D2) * B * ETA +
1592 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1593 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1594 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1595 A;
1596 return inv_pi * (-A * dK + K * dA) / pow2(K);
1597 }
1598 case CalcType::LowXandHighY: {
1599 Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1600 Complex dz2 = dsqrtxy + dsqrty;
1601 Complex dA = sqrt_pi * ((w1 - w2) * dinvGD +
1602 (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1603 Complex dB =
1604 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1605 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1606 9 * dz1) *
1607 invGD) /
1608 (4 * pow4(z1));
1609 Complex dK = Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1610 Complex(d.G2, d.D2) * B * ETA +
1611 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1612 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1613 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1614 A;
1615 return inv_pi * (-A * dK + K * dA) / pow2(K);
1616 }
1617 case CalcType::LowYandLowX: {
1618 Complex dz1 = dsqrtxy;
1619 Complex dz2 = dx / (2 * sqrtx);
1620 Complex dA = 2 *
1621 (sqrt_pi * Complex(G2, D2) * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1622 (ETA - 1) -
1623 Complex(G2, D2) * (sqrt_pi * w2 * z2 - 1) * d.ETA -
1624 Complex(d.G2, d.D2) * (sqrt_pi * w2 * z2 - 1) * (ETA - 1)) /
1625 (pow2(Complex(G2, D2)) * pow2(ETA - 1));
1626 Complex dB =
1627 (-2 * Complex(G2, D2) * (ETA - 1) *
1628 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1629 (2 * pow2(sqrty) + x - 1) +
1630 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
1631 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
1632 Complex(G2, D2) *
1633 (2 * sqrt_pi * w1 * z1 +
1634 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
1635 d.ETA +
1636 Complex(d.G2, d.D2) * (ETA - 1) *
1637 (2 * sqrt_pi * w1 * z1 +
1638 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1)) /
1639 (pow2(Complex(G2, D2)) * pow2(ETA - 1));
1640 Complex dK = Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1641 Complex(d.G2, d.D2) * B * ETA +
1642 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1643 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1644 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1645 A;
1646 return inv_pi * (-A * dK + K * dA) / pow2(K);
1647 }
1648 case CalcType::LowYandHighX: {
1649 Complex dz1 = dsqrtxy;
1650 Complex dA = (2 * Complex(G2, D2) * (ETA - 1) * (x - 3) * dx +
1651 Complex(G2, D2) * (2 * x - 3) * x * d.ETA +
1652 Complex(d.G2, d.D2) * (ETA - 1) * (2 * x - 3) * x) /
1653 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow3(x));
1654 Complex dB =
1655 (-Complex(G2, D2) *
1656 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1657 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1658 x * d.ETA +
1659 Complex(G2, D2) * (ETA - 1) *
1660 (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
1661 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
1662 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) -
1663 Complex(d.G2, d.D2) *
1664 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1665 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1666 (ETA - 1) * x) /
1667 (2 * pow2(Complex(G2, D2)) * pow2(ETA - 1) * pow3(x));
1668 Complex dK = Complex(G2, D2) * B * d.ETA + Complex(G2, D2) * ETA * dB +
1669 Complex(d.G2, d.D2) * B * ETA +
1670 (Complex(G0 - 1.5 * G2, D0 - 1.5 * D2) * ETA - FVC) * dA +
1671 (-Complex(1.5 * G2 - G0, 1.5 * D2 - D0) * d.ETA -
1672 Complex(1.5 * d.G2 - d.G0, 1.5 * d.D2 - d.D0) * ETA - d.FVC) *
1673 A;
1674 return inv_pi * (-A * dK + K * dA) / pow2(K);
1675 }
1676 }
1677 return {};
1678}
1679
1680Complex HartmannTran::operator()(Numeric f) noexcept {
1681 imag_val(deltax) = mF0 - f;
1682 x = deltax / ((1 - ETA) * Complex(G2, D2));
1683 sqrtxy = std::sqrt(x + sqrty * sqrty);
1684 update_calcs();
1685 calc();
1686 return F;
1687}
1688
1689HartmannTran::CalcType HartmannTran::init(const Complex c2t) const noexcept {
1690 if (abs_squared(c2t) == 0)
1691 return CalcType::Noc2tHighZ; // nb. Value of high/low changes elsewhere
1692 if (abs_squared(x) <= 9e-16 * abs_squared(sqrty * sqrty))
1693 return CalcType::LowXandHighY;
1694 if ((abs_squared(sqrty * sqrty) <= 1.e-30 * abs_squared(x)) and
1695 abs_squared(std::sqrt(x)) <= 16.e6)
1696 return CalcType::LowYandLowX; // Weird case, untested
1697 if ((abs_squared(sqrty * sqrty) <= 1.e-30 * abs_squared(x)))
1698 return CalcType::LowYandHighX;
1699 return CalcType::Full;
1700}
1701
1702void HartmannTran::update_calcs() noexcept {
1703 calcs = init((1 - ETA) * Complex(G2, D2));
1704}
1705
1706void HartmannTran::calc() noexcept {
1707 switch (calcs) {
1708 case CalcType::Full:
1709 z1 = sqrtxy - sqrty;
1710 z2 = sqrtxy + sqrty;
1711 w1 = Faddeeva::w(1i * z1);
1712 w2 = Faddeeva::w(1i * z2);
1713 A = sqrt_pi * invGD * (w1 - w2);
1714 B = (-1 + sqrt_pi / (2 * sqrty) * (1 - pow2(z1)) * w1 -
1715 sqrt_pi / (2 * sqrty) * (1 - pow2(z2)) * w2) /
1716 ((1 - ETA) * Complex(G2, D2));
1717 break;
1718 case CalcType::Noc2tLowZ:
1719 case CalcType::Noc2tHighZ:
1720 z1 = deltax * invGD;
1721 w1 = Faddeeva::w(1i * z1);
1722 A = sqrt_pi * invGD * w1;
1723 if (abs_squared(z1) < 16e6) {
1724 calcs = CalcType::Noc2tLowZ;
1725 B = sqrt_pi * invGD * ((1 - pow2(z1)) * w1 + z1 / sqrt_pi);
1726 } else {
1727 calcs = CalcType::Noc2tHighZ;
1728 B = invGD * (sqrt_pi * w1 + 1 / z1 / 2 - 3 / pow3(z1) / 4);
1729 }
1730 break;
1731 case CalcType::LowXandHighY:
1732 z1 = deltax * invGD;
1733 z2 = sqrtxy + sqrty;
1734 w1 = Faddeeva::w(1i * z1);
1735 w2 = Faddeeva::w(1i * z2);
1736 A = sqrt_pi * invGD * (w1 - w2);
1737 B = invGD * (sqrt_pi * w1 + 1 / z1 / 2 - 3 / pow3(z1) / 4);
1738 break;
1739 case CalcType::LowYandLowX:
1740 sqrtx = std::sqrt(x);
1741 z1 = sqrtxy;
1742 z2 = sqrtx;
1743 w1 = Faddeeva::w(1i * z1);
1744 w2 = Faddeeva::w(1i * z2);
1745 A = (2 * sqrt_pi / ((1 - ETA) * Complex(G2, D2))) * (inv_sqrt_pi - z2 * w2);
1746 B = (1 / ((1 - ETA) * Complex(G2, D2))) *
1747 (-1 +
1748 2 * sqrt_pi * (1 - x - 2 * sqrty * sqrty) * (1 / sqrt_pi - z2 * w2) +
1749 2 * sqrt_pi * z1 * w1);
1750 break;
1751 case CalcType::LowYandHighX:
1752 z1 = sqrtxy;
1753 w1 = Faddeeva::w(1i * z1);
1754 A = (1 / ((1 - ETA) * Complex(G2, D2))) * (1 / x - 3 / pow2(x) / 2);
1755 B = (1 / ((1 - ETA) * Complex(G2, D2))) *
1756 (-1 + (1 - x - 2 * sqrty * sqrty) * (1 / x - 3 / pow2(x) / 2) +
1757 2 * sqrt_pi * z1 * w1);
1758 break;
1759 }
1760
1761 dw1 = 2i * (inv_sqrt_pi - z1 * w1);
1762 dw2 = 2i * (inv_sqrt_pi - z2 * w2);
1763 K = 1 - (FVC - ETA * (Complex(G0, D0) - 3 * Complex(G2, D2) / 2)) * A +
1764 ETA * Complex(G2, D2) * B;
1765 F = inv_pi * A / K;
1766}
1767
1768VanVleckHuber::VanVleckHuber(Numeric F0, Numeric T) noexcept
1769 : c1(Constant::h / (2.0 * Constant::k * T)), tanh_c1f0(std::tanh(c1 * F0)),
1770 inv_denom(1.0 / (F0 * tanh_c1f0)) {}
1771
1772Numeric VanVleckHuber::dNdT(Numeric T, Numeric f) const noexcept {
1773 const Numeric part1 =
1774 (1 - pow2(tanh_c1f0)) * c1 * f * tanh_c1f / (T * pow2(tanh_c1f0));
1775 const Numeric part2 = (pow2(tanh_c1f) - 1) * c1 * pow2(f) * inv_denom / T;
1776 return part1 + part2;
1777}
1778
1779Numeric VanVleckHuber::dNdf(Numeric f) const noexcept {
1780 const Numeric part1 = tanh_c1f * inv_denom;
1781 const Numeric part2 = c1 * f * inv_denom / (1 - pow2(tanh_c1f));
1782 return part1 + part2;
1783}
1784
1785Numeric VanVleckHuber::dNdF0() const noexcept {
1786 const Numeric part1 = -N * tanh_c1f0 * inv_denom;
1787 const Numeric part2 = c1 * N * (pow2(tanh_c1f0) - 1) / tanh_c1f0;
1788 return part1 + part2;
1789}
1790
1791Numeric VanVleckHuber::operator()(Numeric f) noexcept {
1792 tanh_c1f = std::tanh(c1 * f);
1793 N = f * tanh_c1f * inv_denom;
1794 return N;
1795}
1796
1797RosenkranzQuadratic::RosenkranzQuadratic(Numeric F0, Numeric T) noexcept
1798 : fac((Constant::h) / (2.0 * Constant::k * T) /
1799 std::sinh((Constant::h * F0) / (2.0 * Constant::k * T)) * (1.0 / F0)),
1800 dfacdT((fac / T) *
1801 (fac * F0 * F0 *
1802 std::cosh((Constant::h * F0) / (2.0 * Constant::k * T)) -
1803 1)),
1804 dfacdF0(-fac * (1 / F0 + fac * F0 *
1805 std::cosh((Constant::h * F0) /
1806 (2.0 * Constant::k * T)))) {}
1807
1808Numeric RosenkranzQuadratic::dNdT(Numeric, Numeric f) const noexcept {
1809 return f * f * dfacdT;
1810}
1811
1812Numeric RosenkranzQuadratic::dNdf(Numeric f) const noexcept {
1813 return 2.0 * f * fac;
1814}
1815
1816Numeric RosenkranzQuadratic::dNdF0() const noexcept {
1817 return N * dfacdF0 / fac;
1818}
1819
1820Numeric RosenkranzQuadratic::operator()(Numeric f) noexcept {
1821 N = fac * f * f;
1822 return N;
1823}
1824
1825Numeric SimpleFrequencyScaling::dNdT(Numeric t_ [[maybe_unused]], Numeric f) const ARTS_NOEXCEPT {
1826 ARTS_ASSERT(t_ == T, "Temperature is stored internally in SimpleFrequencyScaling\n"
1827 "but for interface reasons this function nevertheless takes temprature as an input\n"
1828 "For some reason, the two temperatures disagree, so regardless, you have encountered\n"
1829 "a path of the code that should never be encountered. The two temperatures are: ",
1830 T, " and ", t_, " K")
1831
1832 return
1833 - N * Constant::h * F0 * expF0 / (Constant::k * t_ * t_ * expm1F0)
1834 + Constant::h * f * f * std::exp(- (Constant::h * f) / (Constant::k * t_)) /
1835 (F0 * Constant::k * t_ * t_ * expm1F0);
1836}
1837
1838Numeric SimpleFrequencyScaling::dNdf(Numeric f) const noexcept {
1839 return N / f - N * Constant::h *
1840 std::exp(- (Constant::h * f) / (Constant::k * T)) /
1841 (Constant::k * T * std::expm1(- (Constant::h * f) / (Constant::k * T)));
1842}
1843
1844Numeric SimpleFrequencyScaling::operator()(Numeric f) noexcept {
1845 N = (f / F0) * (std::expm1(- Constant::h * f / (Constant::k * T)) / expm1F0);
1846 return N;
1847}
1848
1849LocalThermodynamicEquilibrium::LocalThermodynamicEquilibrium(
1850 Numeric I0, Numeric T0, Numeric T, Numeric F0, Numeric E0, Numeric QT,
1851 Numeric QT0, Numeric dQTdT, Numeric r, Numeric drdSELFVMR, const Numeric drdT) noexcept
1852 : LocalThermodynamicEquilibrium(
1853 I0, r, drdSELFVMR, drdT, QT0, QT, dQTdT, boltzman_ratio(T, T0, E0),
1854 dboltzman_ratio_dT_div_boltzmann_ratio(T, E0),
1855 stimulated_relative_emission(F0, T0, T),
1856 dstimulated_relative_emission_dT(F0, T0, T),
1857 dstimulated_relative_emission_dF0(F0, T0, T)) {
1858 }
1859
1860struct FullNonLocalThermodynamicEquilibriumInitialization {
1861 Numeric k, dkdF0, dkdr1, dkdr2,
1862 e, dedF0, dedr2,
1863 B, dBdT, dBdF0;
1864
1865 FullNonLocalThermodynamicEquilibriumInitialization(Numeric F0, Numeric A21, Numeric T, Numeric r1, Numeric r2, Numeric c2, Numeric c3, Numeric x) noexcept :
1866 k(c3 * (r1 * x - r2) * (A21 / c2)),
1867 dkdF0(- 2.0 * k / F0),
1868 dkdr1(c3 * x * (A21 / c2)),
1869 dkdr2(- c3 * (A21 / c2)),
1870 e(c3 * r2 * A21),
1871 dedF0(e / F0),
1872 dedr2(c3 * A21),
1873 B(2 * Constant::h / Constant::pow2(Constant::c) * Constant::pow3(F0) / std::expm1((Constant::h / Constant::k * F0) / T)),
1874 dBdT(Constant::pow2(B) * Constant::pow2(Constant::c) * std::exp((Constant::h / Constant::k * F0) / T) / (2*Constant::k*Constant::pow2(F0*T))),
1875 dBdF0(3 * B / F0 - Constant::pow2(B) * Constant::pow2(Constant::c) * std::exp((Constant::h / Constant::k * F0) / T) / (2*Constant::k*T*Constant::pow3(F0)))
1876 {}
1877
1878 constexpr FullNonLocalThermodynamicEquilibrium operator()(
1879 Numeric r, Numeric drdSELFVMR, Numeric drdT
1880 ) && noexcept {
1881 return FullNonLocalThermodynamicEquilibrium(r, drdSELFVMR, drdT,
1882 k, dkdF0, dkdr1, dkdr2,
1883 e, dedF0, dedr2,
1884 B, dBdT, dBdF0);
1885 }
1886};
1887
1888FullNonLocalThermodynamicEquilibrium::FullNonLocalThermodynamicEquilibrium(
1889 Numeric F0, Numeric A21, Numeric T,
1890 Numeric g1, Numeric g2, Numeric r1,
1891 Numeric r2, Numeric r, Numeric drdSELFVMR, Numeric drdT) noexcept :
1892 FullNonLocalThermodynamicEquilibrium(FullNonLocalThermodynamicEquilibriumInitialization(F0, A21, T, r1, r2, c0 * F0 * F0 * F0, c1 * F0, g2 / g1)(r, drdSELFVMR, drdT)) {}
1893
1894struct VibrationalTemperaturesNonLocalThermodynamicEquilibriumInitializer {
1895 Numeric K1, dK1dT,
1896 K2, dK2dT, dK2dF0,
1897 K3, dK3dT, dK3dF0, dK3dTl, dK3dTu,
1898 K4, dK4dT, dK4dTu,
1899 B, dBdT, dBdF0;
1900
1901 VibrationalTemperaturesNonLocalThermodynamicEquilibriumInitializer(Numeric T, Numeric T0, Numeric F0, Numeric E0, Numeric Tl, Numeric Evl, Numeric Tu, Numeric Evu, Numeric gamma, Numeric gamma_ref, Numeric r_low, Numeric r_upp
1902 ) noexcept :
1903 K1(boltzman_ratio(T, T0, E0)),
1904 dK1dT(dboltzman_ratio_dT(K1, T, E0)),
1905 K2(stimulated_relative_emission(gamma, gamma_ref)),
1906 dK2dT(dstimulated_relative_emission_dT(gamma, gamma_ref, F0, T)),
1907 dK2dF0(dstimulated_relative_emission_dF0(gamma, gamma_ref, T, T0)),
1908 K3(absorption_nlte_ratio(gamma, r_upp, r_low)),
1909 dK3dT(dabsorption_nlte_rate_dT(gamma, T, F0, Evl, Evu, K4, r_low)),
1910 dK3dF0(dabsorption_nlte_rate_dF0(gamma, T, K4, r_low)),
1911 dK3dTl(dabsorption_nlte_rate_dTl(gamma, T, Tl, Evl, r_low)),
1912 dK3dTu(dabsorption_nlte_rate_dTu(gamma, T, Tu, Evu, K4)),
1913 K4(boltzman_ratio(Tu, T, Evu)),
1914 dK4dT(dboltzman_ratio_dT(K4, T, Evu)),
1915 dK4dTu(dboltzman_ratio_dT(K4, Tu, Evu)),
1916 B(2 * Constant::h / Constant::pow2(Constant::c) * Constant::pow3(F0) / std::expm1((Constant::h / Constant::k * F0) / T)),
1917 dBdT(Constant::pow2(B) * Constant::pow2(Constant::c) * std::exp((Constant::h / Constant::k * F0) / T) / (2*Constant::k*Constant::pow2(F0*T))),
1918 dBdF0(3 * B / F0 - Constant::pow2(B) * Constant::pow2(Constant::c) * std::exp((Constant::h / Constant::k * F0) / T) / (2*Constant::k*T*Constant::pow3(F0)))
1919 {}
1920
1921 constexpr VibrationalTemperaturesNonLocalThermodynamicEquilibrium operator()(Numeric I0,
1922 Numeric QT0, Numeric QT, Numeric dQTdT,
1923 Numeric r, Numeric drdSELFVMR, Numeric drdT) && noexcept {
1924 return VibrationalTemperaturesNonLocalThermodynamicEquilibrium(I0,
1925 QT0, QT, dQTdT,
1926 r, drdSELFVMR, drdT,
1927 K1, dK1dT,
1928 K2, dK2dT, dK2dF0,
1929 K3, dK3dT, dK3dF0, dK3dTl, dK3dTu,
1930 K4, dK4dT, dK4dTu,
1931 B, dBdT, dBdF0);
1932 }
1933};
1934
1935VibrationalTemperaturesNonLocalThermodynamicEquilibrium::
1936 VibrationalTemperaturesNonLocalThermodynamicEquilibrium(
1937 Numeric I0, Numeric T0, Numeric T, Numeric Tl, Numeric Tu, Numeric F0,
1938 Numeric E0, Numeric Evl, Numeric Evu, Numeric QT, Numeric QT0,
1939 Numeric dQTdT, Numeric r, Numeric drdSELFVMR, Numeric drdT) noexcept
1940 :
1941 VibrationalTemperaturesNonLocalThermodynamicEquilibrium(
1942 VibrationalTemperaturesNonLocalThermodynamicEquilibriumInitializer(T, T0, F0, E0, Tl, Evl, Tu, Evu,
1943 stimulated_emission(T, F0),
1944 stimulated_emission(T0, F0),
1945 boltzman_ratio(Tl, T, Evl),
1946 boltzman_ratio(Tu, T, Evu)
1947 )(I0, QT0, QT, dQTdT, r, drdSELFVMR, drdT)
1948 ) {}
1949
1950#define CutInternalDerivativesImpl(X, Y) \
1951 else if (deriv == Jacobian::Line::Shape##X##Y) { \
1952 const Numeric d = value.n; \
1953 const Complex dFm = \
1954 std::conj(ls_mirr.dFd##X(d) - ls_mirr_cut.dFd##X(d)); \
1955 const Complex dFls = ls.dFd##X(d) - ls_cut.dFd##X(d) + dFm; \
1956 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls; \
1957 if (do_nlte) { \
1958 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls; \
1959 } \
1960 }
1961
1962#define CutInternalDerivatives(X) \
1963CutInternalDerivativesImpl(X, X0) CutInternalDerivativesImpl(X, X1) \
1964 CutInternalDerivativesImpl(X, X2) CutInternalDerivativesImpl(X, X3)
1965
1966#define InternalDerivativesImpl(X, Y) \
1967 else if (deriv == Jacobian::Line::Shape##X##Y) { \
1968 const Numeric d = value.n; \
1969 const Complex dFm = std::conj(ls_mirr.dFd##X(d)); \
1970 const Complex dFls = ls.dFd##X(d) + dFm; \
1971 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls; \
1972 if (do_nlte) { \
1973 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls; \
1974 } \
1975 }
1976
1977#define InternalDerivatives(X) \
1978 InternalDerivativesImpl(X, X0) InternalDerivativesImpl(X, X1) \
1979 InternalDerivativesImpl(X, X2) InternalDerivativesImpl(X, X3)
1980
1981#define InternalDerivativesSetupImpl(X, Y) \
1982 else if (deriv == Jacobian::Line::Shape##X##Y) { \
1983 const Index pos = \
1984 band.BroadeningSpeciesPosition(deriv.Target().LineSpecies()); \
1985 if (pos >= 0) { \
1986 derivs[ij].value.n = band.Line(i).LineShape().d##X##_d##Y( \
1987 T, band.T0(), P, pos, vmrs); \
1988 } else { \
1989 derivs[ij].value.n = 0; \
1990 } \
1991 }
1992
1993#define InternalDerivativesSetup(X) \
1994 InternalDerivativesSetupImpl(X, X0) InternalDerivativesSetupImpl(X, X1) \
1995 InternalDerivativesSetupImpl(X, X2) InternalDerivativesSetupImpl(X, X3)
1996
1997#define InternalDerivativesGImpl(X) \
1998 else if (deriv == Jacobian::Line::ShapeG##X) { \
1999 const Numeric dLM = value.n; \
2000 com.dF[com.jac_pos(iv, ij)] += dLM * S * Fls; \
2001 if (do_nlte) { \
2002 com.dN[com.jac_pos(iv, ij)] += dLM * DS * Fls; \
2003 } \
2004 }
2005
2006#define InternalDerivativesG \
2007 InternalDerivativesGImpl(X0) InternalDerivativesGImpl(X1) \
2008 InternalDerivativesGImpl(X2) InternalDerivativesGImpl(X3)
2009
2010#define InternalDerivativesYImpl(X) \
2011 else if (deriv == Jacobian::Line::ShapeY##X) { \
2012 const Complex dLM = Complex(0, -value.n); \
2013 com.dF[com.jac_pos(iv, ij)] += dLM * S * Fls; \
2014 if (do_nlte) { \
2015 com.dN[com.jac_pos(iv, ij)] += dLM * DS * Fls; \
2016 } \
2017 }
2018
2019#define InternalDerivativesY \
2020 InternalDerivativesYImpl(X0) InternalDerivativesYImpl(X1) \
2021 InternalDerivativesYImpl(X2) InternalDerivativesYImpl(X3)
2022
2024struct CutoffRange {
2025 Index start, size;
2026};
2027
2039CutoffRange limited_range(const Numeric fl, const Numeric fu, const Vector &f_grid) ARTS_NOEXCEPT {
2040 ARTS_ASSERT(fu > fl);
2041 const Numeric * it0 = f_grid.get_c_array();
2042 const Numeric * itn = it0 + f_grid.size();
2043 const Numeric * itl = std::lower_bound(it0, itn, fl);
2044 return CutoffRange{std::distance(it0, itl), std::distance(itl, std::upper_bound(itl, itn, fu))};
2045}
2046
2048struct SparseLimitRange {
2049 Index start, size;
2050 Index start_low, size_low;
2051 Index start_upp, size_upp;
2052};
2053
2054SparseLimitRange linear_sparse_limited_range(const Numeric flc,
2055 const Numeric fuc,
2056 const Numeric fls,
2057 const Numeric fus,
2058 const Vector& f_grid,
2059 const Vector& sparse_f_grid) ARTS_NOEXCEPT {
2060 ARTS_ASSERT(fls > flc);
2061 ARTS_ASSERT(fus > fls);
2062 ARTS_ASSERT(fuc > fus);
2063
2064 const Index nvs = sparse_f_grid.size();
2065
2066 // Find bounds in sparse
2067 const Numeric * it0s = sparse_f_grid.get_c_array();
2068 const Numeric * itns = it0s + nvs;
2069 const Numeric * itlc = std::lower_bound(it0s, itns, std::nextafter(flc, fuc)); // lower cutoff
2070 const Numeric * ituc = std::upper_bound(itlc, itns, std::nextafter(fuc, flc)); // upper cutoff
2071 const Numeric * itls = std::upper_bound(itlc, ituc, fls); // lower sparse
2072 const Numeric * itus = std::lower_bound(itls, ituc, fus); // upper sparse
2073
2074 /* Start and size of sparse adjusted to the 2-grid
2075 *
2076 * The interface to the dense grid is altered slightly so that
2077 * a complete set-of-2 quadratic interopolation points becomes
2078 * a guarantee. Their start/end points ignore this restriction
2079 * but have been defined to not contain anything extra
2080 */
2081 Index beg_lr = std::distance(it0s, itlc); /*while (beg_lr % 2) --beg_lr;*/
2082 Index end_lr = std::distance(it0s, itls); while (end_lr % 2) --end_lr;
2083 Index beg_ur = std::distance(it0s, itus); while (beg_ur % 2) ++beg_ur;
2084 Index end_ur = std::distance(it0s, ituc); /*while (end_ur % 2) ++end_ur;*/
2085
2086 // Find new limits
2087 const Numeric fl = (end_lr <= 0 or end_lr >= nvs) ? flc : sparse_f_grid[end_lr];
2088 const Numeric fu = (beg_ur <= 0 or beg_ur >= nvs) ? fuc : sparse_f_grid[beg_ur];
2089
2090 // Find bounds in dense
2091 const Numeric * it0 = f_grid.get_c_array();
2092 const Numeric * itn = it0 + f_grid.size();
2093 const Numeric * itl = std::lower_bound(it0, itn, fl); // include boundary
2094 const Numeric * itu = std::upper_bound(itl, itn, std::nextafter(fu, fl)); // dismiss boundary
2095
2096 return {
2097 std::distance(it0, itl), std::distance(itl, itu),
2098 beg_lr, end_lr - beg_lr,
2099 beg_ur, end_ur - beg_ur
2100 };
2101}
2102
2103SparseLimitRange quad_sparse_limited_range(const Numeric flc,
2104 const Numeric fuc,
2105 const Numeric fls,
2106 const Numeric fus,
2107 const Vector& f_grid,
2108 const Vector& sparse_f_grid) ARTS_NOEXCEPT {
2109 ARTS_ASSERT(fls > flc);
2110 ARTS_ASSERT(fus > fls);
2111 ARTS_ASSERT(fuc > fus);
2112
2113 const Index nvs = sparse_f_grid.size();
2114 const Index nv = f_grid.size();
2115
2116 // Find bounds in sparse
2117 const Numeric * const it0s = sparse_f_grid.get_c_array();
2118 const Numeric * const itns = it0s + nvs;
2119 const Numeric * itlc = std::lower_bound(it0s, itns, std::nextafter(flc, fuc)); // lower cutoff
2120 const Numeric * ituc = std::upper_bound(itlc, itns, std::nextafter(fuc, flc)); // upper cutoff
2121 const Numeric * itls = std::upper_bound(itlc, ituc, std::nextafter(fls, flc)); // lower sparse
2122 const Numeric * itus = std::lower_bound(itls, ituc, std::nextafter(fus, fuc)); // upper sparse
2123
2124 /* Start and size of sparse adjusted to the 3-grid
2125 *
2126 * The interface to the dense grid is altered slightly so that
2127 * a complete set-of-3 quadratic interpolation points becomes
2128 * a guarantee. Their end points are modified to contain
2129 * only set-of-3. If the range is not modified correctly,
2130 * you risk having negative interpolation out of your range.
2131 *
2132 * To avoid negative absorption entirely, the range between
2133 * the true cutoff and the outer-most sparse grid point must
2134 * have at least two values. If they have only one value,
2135 * the quadratic interpolation will lead to negative values
2136 * in the outer range. So there's a small range that has to
2137 * be ignored
2138 */
2139 while (std::distance(it0s, itls) % 3) --itls;
2140 while (std::distance(itlc, itls) % 3 == 1) ++itlc; // skip some cutoff
2141 while (std::distance(it0s, itus) % 3) ++itus;
2142 while (std::distance(itus, ituc) % 3 == 1) --ituc; // skip some cutoff
2143
2144 // Find bounds in dense
2145 const Numeric * const it0 = f_grid.get_c_array();
2146 const Numeric * const itn = it0 + nv;
2147 const Numeric * itl;
2148 const Numeric * itu;
2149
2150 if (itls not_eq itns) {
2151 itl = std::lower_bound(it0, itn, *itls); // include boundary
2152 } else {
2153 itl = std::lower_bound(it0, itn, flc); // include boundary
2154 }
2155
2156 if (itus not_eq itns and itl not_eq itn) {
2157 itu = std::upper_bound(itl, itn, std::nextafter(*itus, *itl)); // dismiss boundary
2158 } else {
2159 itu = std::lower_bound(itl, itn, std::nextafter(fuc, flc)); // include boundary
2160 }
2161
2162 return {
2163 std::distance(it0, itl), std::distance(itl, itu),
2164 std::distance(it0s, itlc), std::distance(itlc, itls),
2165 std::distance(it0s, itus), std::distance(itus, ituc)
2166 };
2167}
2168
2174struct Derivatives {
2175 Absorption::QuantumIdentifierLineTarget target;
2176 union Values {
2177 Output o;
2178 Numeric n;
2179 constexpr Values() noexcept : n() {}
2180 } value;
2181 Index jac_pos{};
2182 const RetrievalQuantity * deriv{nullptr};
2183 Derivatives() noexcept : target(), value() {}
2184};
2185
2187using ArrayOfDerivatives = Array<Derivatives>;
2188
2190Index active_nelem(const ArrayOfDerivatives& derivs) noexcept {
2191 return std::distance(derivs.cbegin(),
2192 std::find_if(derivs.cbegin(), derivs.cend(),
2193 [](auto& deriv) {
2194 return deriv.deriv == nullptr;
2195 }));
2196}
2197
2198struct ComputeValues {
2199 Complex * const F;
2200 Complex * const dF;
2201 Complex * const N;
2202 Complex * const dN;
2203 const Numeric * const f;
2204 const Index size;
2205 const ArrayOfDerivatives& derivs;
2206 const Index jac_size;
2207 const Index max_jac_size;
2208 const bool do_nlte;
2209
2210 [[nodiscard]] Index jac_pos(Index iv, Index ij) const noexcept {return jac_size * iv + ij;}
2211
2212 ComputeValues(ComplexVector &F_,
2213 ComplexMatrix &dF_,
2214 ComplexVector &N_,
2215 ComplexMatrix &dN_,
2216 const Vector &f_grid,
2217 const Index start,
2218 const Index nv,
2219 const ArrayOfDerivatives& derivs_,
2220 const bool do_nlte_) noexcept :
2221 F(F_.get_c_array()+start),
2222 dF(dF_.get_c_array()+start*derivs_.size()),
2223 N(N_.get_c_array()+start),
2224 dN(dN_.get_c_array()+start*derivs_.size()),
2225 f(f_grid.get_c_array()+start),
2226 size(nv), derivs(derivs_), jac_size(derivs_.size()), max_jac_size(active_nelem(derivs)), do_nlte(do_nlte_) {
2227 }
2228
2229 ComputeValues(Complex &F_, std::vector<Complex> &dF_,
2230 Complex &N_, std::vector<Complex> &dN_,
2231 const Numeric &f_lim,
2232 const ArrayOfDerivatives& derivs_, const bool do_nlte_) noexcept :
2233 F(&F_), dF(dF_.data()), N(&N_), dN(dN_.data()), f(&f_lim), size(1),
2234 derivs(derivs_), jac_size(derivs.nelem()), max_jac_size(active_nelem(derivs)), do_nlte(do_nlte_) {
2235 }
2236
2237 ComputeValues& operator-=(const ComputeValues& cut) ARTS_NOEXCEPT {
2238 ARTS_ASSERT(cut.size == 1, "Not a cutoff limit")
2239 ARTS_ASSERT(cut.jac_size == jac_size, "Not from the same Jacobian type")
2240
2241 for (Index iv=0; iv<size; iv++) {
2242 F[iv] -= *cut.F;
2243 for (Index ij=0; ij<max_jac_size; ij++) {
2244 dF[jac_pos(iv, derivs[ij].jac_pos)] -= cut.dF[derivs[ij].jac_pos];
2245 }
2246 }
2247 if (do_nlte) {
2248 for (Index iv=0; iv<size; iv++) {
2249 N[iv] -= *cut.N;
2250 for (Index ij=0; ij<max_jac_size; ij++) {
2251 dN[jac_pos(iv, derivs[ij].jac_pos)] -= cut.dN[derivs[ij].jac_pos];
2252 }
2253 }
2254 }
2255 return *this;
2256 }
2257};
2258
2321void cutoff_frequency_loop(ComputeValues &com,
2322 Calculator &ls,
2323 Calculator &ls_mirr,
2324 Normalizer &ls_norm,
2325 const Calculator &ls_cut,
2326 const Calculator &ls_mirr_cut,
2327 const IntensityCalculator &ls_str,
2328 const ArrayOfDerivatives &derivs,
2329 const Complex LM,
2330 const Numeric &T,
2331 const Numeric &dfdH,
2332 const Numeric &Sz,
2333 const Species::Species self_species) ARTS_NOEXCEPT {
2334 const Index nv = com.size;
2335 const bool do_nlte = com.do_nlte;
2336
2337 const Numeric Si = ls_str.S();
2338 const Numeric DNi = ls_str.N();
2339
2340 for (Index iv = 0; iv < nv; iv++) {
2341 const Numeric f = com.f[iv];
2342
2343 const Numeric Sn = ls_norm(f);
2344 const Numeric S = Sz * Sn * Si;
2345 const Numeric DS = Sz * Sn * DNi;
2346 const Complex Fm = std::conj(ls_mirr(f) - ls_mirr_cut.F());
2347 const Complex Fls = ls(f) - ls_cut.F() + Fm;
2348 com.F[iv] += S * LM * Fls;
2349 if (do_nlte) {
2350 com.N[iv] += DS * LM * Fls;
2351 }
2352
2353 for (const auto& [lt, value, ij, deriv_ptr]: derivs) {
2354 if (not deriv_ptr) break;
2355 const auto& deriv = *deriv_ptr;
2356
2357 if (deriv == Jacobian::Atm::Temperature) {
2358 const auto &dXdT = value.o;
2359 const Numeric dSn = ls_norm.dNdT(T, f);
2360 const Numeric dSi = ls_str.dSdT();
2361 const Complex dLM(dXdT.G, -dXdT.Y);
2362 const Complex dFm = std::conj(ls_mirr.dFdT(mirroredOutput(dXdT), T) - ls_mirr_cut.dFdT(mirroredOutput(dXdT), T));
2363 const Complex dFls = ls.dFdT(dXdT, T) - ls_cut.dFdT(dXdT, T) + dFm;
2364 com.dF[com.jac_pos(iv, ij)] +=
2365 S * (LM * dFls + dLM * Fls) + Sz * (dSn * Si + Sn * dSi) * LM * Fls;
2366 if (do_nlte) {
2367 const Numeric dDSi = ls_str.dNdT();
2368 com.dN[com.jac_pos(iv, ij)] += DS * (LM * dFls + dLM * Fls) +
2369 Sz * (dSn * Si + Sn * dDSi) * LM * Fls;
2370 }
2371 } else if (deriv.is_wind()) {
2372 const Complex dFm = std::conj(ls_mirr.dFdf() - ls_mirr_cut.dFdf());
2373 const Complex dFls = ls.dFdf() - ls_cut.dFdf() + dFm;
2374 const Numeric dS = Sz * ls_norm.dNdf(f) * Si;
2375 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dS * LM * Fls;
2376 if (do_nlte) {
2377 const Numeric dDS = Sz * ls_norm.dNdf(f) * DNi;
2378 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + dDS * LM * Fls;
2379 }
2380 } else if (deriv.is_mag()) {
2381 const Complex dFm = std::conj(ls_mirr.dFdH(-dfdH) - ls_mirr_cut.dFdH(-dfdH));
2382 const Complex dFls = ls.dFdH(dfdH) - ls_cut.dFdH(dfdH) + dFm;
2383 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls;
2384 if (do_nlte) {
2385 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls;
2386 }
2387 } else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR) {
2388 com.dF[com.jac_pos(iv, ij)] += ls_str.dSdSELFVMR() * LM * Fls;
2389 if (do_nlte) {
2390 com.dN[com.jac_pos(iv, ij)] += ls_str.dNdSELFVMR() * LM * Fls;
2391 }
2392 } else if (deriv.Target().needQuantumIdentity()) {
2393 if (deriv == Jacobian::Line::VMR) {
2394 const auto &dXdVMR = value.o;
2395 const Complex dFm = std::conj(ls_mirr.dFdVMR(mirroredOutput(dXdVMR)) - ls_mirr_cut.dFdVMR(mirroredOutput(dXdVMR)));
2396 const Complex dLM(dXdVMR.G, -dXdVMR.Y);
2397 const Complex dFls = ls.dFdVMR(dXdVMR) - ls_cut.dFdVMR(dXdVMR) + dFm;
2398 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dLM * S * Fls;
2399 if (self_species == deriv.QuantumIdentity().Species()) {
2400 com.dF[com.jac_pos(iv, ij)] += ls_str.dSdSELFVMR() * LM * Fls;
2401 }
2402 if (do_nlte) {
2403 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + dLM * DS * Fls;
2404 if (self_species == deriv.QuantumIdentity().Species()) {
2405 com.dF[com.jac_pos(iv, ij)] += ls_str.dNdSELFVMR() * LM * Fls;
2406 }
2407 }
2408 } else {
2409 if (lt == Absorption::QuantumIdentifierLineTargetType::Line) {
2410 if (deriv == Jacobian::Line::Center) {
2411 const Numeric d = ls_norm.dNdF0() * Si + ls_str.dSdF0() * Sn;
2412 const Complex dFm = std::conj(ls_mirr.dFdF0() - ls_mirr_cut.dFdF0());
2413 const Complex dFls = ls.dFdF0() - ls_cut.dFdF0() + dFm;
2414 const Numeric dS = Sz * d;
2415 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dS * LM * Fls;
2416 if (do_nlte) {
2417 const Numeric Dd = ls_norm.dNdF0() * DNi + ls_str.dNdF0() * Sn;
2418 const Numeric DdS = Sz * Dd;
2419 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + DdS * LM * Fls;
2420 }
2421 } else if (deriv == Jacobian::Line::Strength) {
2422 const Numeric dS = Sz * Sn * ls_str.dSdI0();
2423 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
2424 if (do_nlte) {
2425 const Numeric DdS = Sz * Sn * ls_str.dNdI0();
2426 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
2427 }
2428 }
2429 // All line shape derivatives
2430 CutInternalDerivatives(G0)
2431 CutInternalDerivatives(D0)
2432 CutInternalDerivatives(G2)
2433 CutInternalDerivatives(D2)
2434 CutInternalDerivatives(ETA)
2435 CutInternalDerivatives(FVC)
2436 InternalDerivativesY
2437 InternalDerivativesG
2438 CutInternalDerivatives(DV)
2439 } else if (lt == Absorption::QuantumIdentifierLineTargetType::Level) {
2440 if (deriv == Jacobian::Line::NLTE) {
2441 if (lt.upper) {
2442 const Numeric dS = Sz * Sn * ls_str.dSdNLTEu();
2443 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
2444
2445 if (do_nlte) {
2446 const Numeric DdS = Sz * Sn * ls_str.dNdNLTEu();
2447 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
2448 }
2449 }
2450
2451 if (lt.lower) {
2452 const Numeric dS = Sz * Sn * ls_str.dSdNLTEl();
2453 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
2454
2455 if (do_nlte) {
2456 const Numeric DdS = Sz * Sn * ls_str.dNdNLTEl();
2457 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
2458 }
2459 }
2460 }
2461 }
2462 }
2463 }
2464 }
2465 }
2466}
2467
2529void frequency_loop(ComputeValues &com,
2530 Calculator &ls,
2531 Calculator &ls_mirr,
2532 Normalizer &ls_norm,
2533 const IntensityCalculator &ls_str,
2534 const ArrayOfDerivatives &derivs,
2535 const Complex LM,
2536 const Numeric &T,
2537 const Numeric &dfdH,
2538 const Numeric &Sz,
2539 const Species::Species self_species) ARTS_NOEXCEPT {
2540 const Index nv = com.size;
2541 const bool do_nlte = com.do_nlte;
2542
2543 const Numeric Si = ls_str.S();
2544 const Numeric DNi = ls_str.N();
2545
2546 for (Index iv = 0; iv < nv; iv++) {
2547 const Numeric f = com.f[iv];
2548
2549 const Numeric Sn = ls_norm(f);
2550 const Numeric S = Sz * Sn * Si;
2551 const Numeric DS = Sz * Sn * DNi;
2552 const Complex Fm = std::conj(ls_mirr(f));
2553 const Complex Fls = ls(f) + Fm;
2554 com.F[iv] += S * LM * Fls;
2555 if (do_nlte) {
2556 com.N[iv] += DS * LM * Fls;
2557 }
2558
2559 for (const auto& [lt, value, ij, deriv_ptr]: derivs) {
2560 if (not deriv_ptr) break;
2561 const auto& deriv = *deriv_ptr;
2562
2563 if (deriv == Jacobian::Atm::Temperature) {
2564 const auto &dXdT = value.o;
2565 const Numeric dSn = ls_norm.dNdT(T, f);
2566 const Numeric dSi = ls_str.dSdT();
2567 const Complex dLM(dXdT.G, -dXdT.Y);
2568 const Complex dFm = std::conj(ls_mirr.dFdT(mirroredOutput(dXdT), T));
2569 const Complex dFls = ls.dFdT(dXdT, T) + dFm;
2570 com.dF[com.jac_pos(iv, ij)] +=
2571 S * (LM * dFls + dLM * Fls) + Sz * (dSn * Si + Sn * dSi) * LM * Fls;
2572 if (do_nlte) {
2573 const Numeric dDSi = ls_str.dNdT();
2574 com.dN[com.jac_pos(iv, ij)] += DS * (LM * dFls + dLM * Fls) +
2575 Sz * (dSn * Si + Sn * dDSi) * LM * Fls;
2576 }
2577 } else if (deriv.is_wind()) {
2578 const Complex dFm = std::conj(ls_mirr.dFdf());
2579 const Complex dFls = ls.dFdf() + dFm;
2580 const Numeric dS = Sz * ls_norm.dNdf(f) * Si;
2581 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dS * LM * Fls;
2582 if (do_nlte) {
2583 const Numeric dDS = Sz * ls_norm.dNdf(f) * DNi;
2584 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + dDS * LM * Fls;
2585 }
2586 } else if (deriv.is_mag()) {
2587 const Complex dFm = std::conj(ls_mirr.dFdH(-dfdH));
2588 const Complex dFls = ls.dFdH(dfdH) + dFm;
2589 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls;
2590 if (do_nlte) {
2591 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls;
2592 }
2593 } else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR) {
2594 com.dF[com.jac_pos(iv, ij)] += ls_str.dSdSELFVMR() * LM * Fls;
2595 if (do_nlte) {
2596 com.dN[com.jac_pos(iv, ij)] += ls_str.dNdSELFVMR() * LM * Fls;
2597 }
2598 } else if (deriv.Target().needQuantumIdentity()) {
2599 if (deriv == Jacobian::Line::VMR) {
2600 const auto &dXdVMR = value.o;
2601 const Complex dFm = std::conj(ls_mirr.dFdVMR(mirroredOutput(dXdVMR)));
2602 const Complex dLM(dXdVMR.G, -dXdVMR.Y);
2603 const Complex dFls = ls.dFdVMR(dXdVMR) + dFm;
2604 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dLM * S * Fls;
2605 if (self_species == deriv.QuantumIdentity().Species()) {
2606 com.dF[com.jac_pos(iv, ij)] += ls_str.dSdSELFVMR() * LM * Fls;
2607 }
2608 if (do_nlte) {
2609 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + dLM * DS * Fls;
2610 if (self_species == deriv.QuantumIdentity().Species()) {
2611 com.dN[com.jac_pos(iv, ij)] += ls_str.dNdSELFVMR() * LM * Fls;
2612 }
2613 }
2614 } else {
2615 if (lt == Absorption::QuantumIdentifierLineTargetType::Line) {
2616 if (deriv == Jacobian::Line::Center) {
2617 const Numeric d = ls_norm.dNdF0() * Si + ls_str.dSdF0() * Sn;
2618 const Complex dFm = std::conj(ls_mirr.dFdF0());
2619 const Complex dFls = ls.dFdF0() + dFm;
2620 const Numeric dS = Sz * d;
2621 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dS * LM * Fls;
2622 if (do_nlte) {
2623 const Numeric Dd = ls_norm.dNdF0() * DNi + ls_str.dNdF0() * Sn;
2624 const Numeric DdS = Sz * Dd;
2625 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + DdS * LM * Fls;
2626 }
2627 } else if (deriv == Jacobian::Line::Strength) {
2628 const Numeric dS = Sz * Sn * ls_str.dSdI0();
2629 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
2630 if (do_nlte) {
2631 const Numeric DdS = Sz * Sn * ls_str.dNdI0();
2632 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
2633 }
2634 }
2635 // All line shape derivatives
2636 InternalDerivatives(G0)
2637 InternalDerivatives(D0)
2638 InternalDerivatives(G2)
2639 InternalDerivatives(D2)
2640 InternalDerivatives(ETA)
2641 InternalDerivatives(FVC)
2642 InternalDerivativesY
2643 InternalDerivativesG
2644 InternalDerivatives(DV)
2645 } else if (lt == Absorption::QuantumIdentifierLineTargetType::Level) {
2646 if (deriv == Jacobian::Line::NLTE) {
2647 if (lt.upper) {
2648 const Numeric dS = Sz * Sn * ls_str.dSdNLTEu();
2649 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
2650
2651 if (do_nlte) {
2652 const Numeric DdS = Sz * Sn * ls_str.dNdNLTEu();
2653 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
2654 }
2655 }
2656
2657 if (lt.lower) {
2658 const Numeric dS = Sz * Sn * ls_str.dSdNLTEl();
2659 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
2660
2661 if (do_nlte) {
2662 const Numeric DdS = Sz * Sn * ls_str.dNdNLTEl();
2663 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
2664 }
2665 }
2666 }
2667 }
2668 }
2669 }
2670 }
2671 }
2672}
2673
2712void cutoff_loop_sparse_linear(ComputeData &com,
2713 ComputeData &sparse_com,
2714 Normalizer ls_norm,
2715 const IntensityCalculator ls_str,
2716 const AbsorptionLines &band,
2717 const ArrayOfDerivatives &derivs,
2718 const Output X,
2719 const Numeric &T,
2720 const Numeric &H,
2721 const Numeric &sparse_lim,
2722 const Numeric &DC,
2723 const Index i,
2724 const bool do_zeeman,
2725 const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT {
2726 // Basic settings
2727 const bool do_nlte = ls_str.do_nlte();
2728 const bool do_cutoff = band.Cutoff() not_eq Absorption::CutoffType::None;
2729 const Numeric fu = band.CutoffFreq(i, X.D0);
2730 const Numeric fl = band.CutoffFreqMinus(i, X.D0);
2731 const Numeric fus = band.F0(i) + sparse_lim;
2732 const Numeric fls = band.F0(i) - sparse_lim;
2733
2734 // Find sparse and dense ranges
2735 const auto [dense_start, dense_size,
2736 sparse_low_start, sparse_low_size,
2737 sparse_upp_start, sparse_upp_size] = linear_sparse_limited_range(fl, fu, fls, fus,
2738 com.f_grid,
2739 sparse_com.f_grid);
2740 if ((dense_size + sparse_low_size + sparse_upp_size) == 0) return;
2741
2742 // Get the compute data view
2743 ComputeValues comval(com.F, com.dF, com.N, com.dN, com.f_grid, dense_start, dense_size, derivs, do_nlte);
2744
2745 // Get views of the sparse data
2746 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);
2747 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);
2748
2749 const Index nz = do_zeeman ? band.ZeemanCount(i, zeeman_polarization) : 1;
2750 for (Index iz = 0; iz < nz; iz++) {
2751 const Numeric dfdH = do_zeeman ? band.ZeemanSplitting(i, zeeman_polarization, iz) : 0;
2752 const Numeric Sz = do_zeeman ? band.ZeemanStrength(i, zeeman_polarization, iz) : 1;
2753 const Complex LM = Complex(1 + X.G, -X.Y);
2754 Calculator ls(band.LineShapeType(), band.F0(i), X, DC, dfdH * H, band.Mirroring() == Absorption::MirroringType::Manual);
2755 Calculator ls_mirr(band.Mirroring(), band.LineShapeType(), band.F0(i), X, DC, dfdH * H);
2756
2757 if (do_cutoff) {
2758 // Initialize and set the cutoff values
2759 Calculator ls_cut = ls;
2760 Calculator ls_mirr_cut = ls_mirr;
2761 ls_cut(fu);
2762 ls_mirr_cut(fu);
2763
2764 cutoff_frequency_loop(comval, ls, ls_mirr, ls_norm, ls_cut, ls_mirr_cut, ls_str,
2765 derivs, LM, T, dfdH, Sz, band.Species());
2766
2767 if (sparse_low_size) {
2768 cutoff_frequency_loop(sparse_low_range, ls, ls_mirr, ls_norm, ls_cut, ls_mirr_cut, ls_str,
2769 derivs, LM, T, dfdH, Sz, band.Species());
2770 }
2771
2772 if (sparse_upp_size) {
2773 cutoff_frequency_loop(sparse_upp_range, ls, ls_mirr, ls_norm, ls_cut, ls_mirr_cut, ls_str,
2774 derivs, LM, T, dfdH, Sz, band.Species());
2775 }
2776
2777 } else {
2778 frequency_loop(comval, ls, ls_mirr, ls_norm, ls_str,
2779 derivs, LM, T, dfdH, Sz, band.Species());
2780
2781 if (sparse_low_size) {
2782 frequency_loop(sparse_low_range, ls, ls_mirr, ls_norm, ls_str,
2783 derivs, LM, T, dfdH, Sz, band.Species());
2784 }
2785
2786 if (sparse_upp_size) {
2787 frequency_loop(sparse_upp_range, ls, ls_mirr, ls_norm, ls_str,
2788 derivs, LM, T, dfdH, Sz, band.Species());
2789 }
2790 }
2791 }
2792}
2793
2794void cutoff_loop_sparse_triple(ComputeData &com,
2795 ComputeData &sparse_com,
2796 Normalizer ls_norm,
2797 const IntensityCalculator ls_str,
2798 const AbsorptionLines &band,
2799 const ArrayOfDerivatives &derivs,
2800 const Output X,
2801 const Numeric &T,
2802 const Numeric &H,
2803 const Numeric &sparse_lim,
2804 const Numeric &DC,
2805 const Index i,
2806 const bool do_zeeman,
2807 const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT {
2808 // Basic settings
2809 const bool do_nlte = ls_str.do_nlte();
2810 const bool do_cutoff = band.Cutoff() not_eq Absorption::CutoffType::None;
2811 const Numeric fu = band.CutoffFreq(i, X.D0);
2812 const Numeric fl = band.CutoffFreqMinus(i, X.D0);
2813 const Numeric fus = band.F0(i) + sparse_lim;
2814 const Numeric fls = band.F0(i) - sparse_lim;
2815
2816 // Find sparse and dense ranges
2817 const auto [dense_start, dense_size,
2818 sparse_low_start, sparse_low_size,
2819 sparse_upp_start, sparse_upp_size] = quad_sparse_limited_range(fl, fu, fls, fus,
2820 com.f_grid,
2821 sparse_com.f_grid);
2822 if ((dense_size + sparse_low_size + sparse_upp_size) == 0) return;
2823
2824 // Get the compute data view
2825 ComputeValues comval(com.F, com.dF, com.N, com.dN, com.f_grid, dense_start, dense_size, derivs, do_nlte);
2826
2827 // Get views of the sparse data
2828 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);
2829 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);
2830
2831 const Index nz = do_zeeman ? band.ZeemanCount(i, zeeman_polarization) : 1;
2832 for (Index iz = 0; iz < nz; iz++) {
2833 const Numeric dfdH = do_zeeman ? band.ZeemanSplitting(i, zeeman_polarization, iz) : 0;
2834 const Numeric Sz = do_zeeman ? band.ZeemanStrength(i, zeeman_polarization, iz) : 1;
2835 const Complex LM = Complex(1 + X.G, -X.Y);
2836 Calculator ls(band.LineShapeType(), band.F0(i), X, DC, dfdH * H, band.Mirroring() == Absorption::MirroringType::Manual);
2837 Calculator ls_mirr(band.Mirroring(), band.LineShapeType(), band.F0(i), X, DC, dfdH * H);
2838
2839 if (do_cutoff) {
2840 // Initialize and set the cutoff values
2841 Calculator ls_cut = ls;
2842 Calculator ls_mirr_cut = ls_mirr;
2843 ls_cut(fu);
2844 ls_mirr_cut(fu);
2845
2846 cutoff_frequency_loop(comval, ls, ls_mirr, ls_norm, ls_cut, ls_mirr_cut, ls_str,
2847 derivs, LM, T, dfdH, Sz, band.Species());
2848
2849 if (sparse_low_size) {
2850 cutoff_frequency_loop(sparse_low_range, ls, ls_mirr, ls_norm, ls_cut, ls_mirr_cut, ls_str,
2851 derivs, LM, T, dfdH, Sz, band.Species());
2852 }
2853
2854 if (sparse_upp_size) {
2855 cutoff_frequency_loop(sparse_upp_range, ls, ls_mirr, ls_norm, ls_cut, ls_mirr_cut, ls_str,
2856 derivs, LM, T, dfdH, Sz, band.Species());
2857 }
2858
2859 } else {
2860 frequency_loop(comval, ls, ls_mirr, ls_norm, ls_str,
2861 derivs, LM, T, dfdH, Sz, band.Species());
2862
2863 if (sparse_low_size) {
2864 frequency_loop(sparse_low_range, ls, ls_mirr, ls_norm, ls_str,
2865 derivs, LM, T, dfdH, Sz, band.Species());
2866 }
2867
2868 if (sparse_upp_size) {
2869 frequency_loop(sparse_upp_range, ls, ls_mirr, ls_norm, ls_str,
2870 derivs, LM, T, dfdH, Sz, band.Species());
2871 }
2872 }
2873 }
2874}
2875
2914void cutoff_loop(ComputeData &com,
2915 Normalizer ls_norm,
2916 const IntensityCalculator ls_str,
2917 const AbsorptionLines &band,
2918 const ArrayOfDerivatives &derivs,
2919 const Output X,
2920 const Numeric &T,
2921 const Numeric &H,
2922 const Numeric &DC,
2923 const Index i,
2924 const bool do_zeeman,
2925 const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT {
2926 // Basic settings
2927 const bool do_nlte = ls_str.do_nlte();
2928 const bool do_cutoff = band.Cutoff() not_eq Absorption::CutoffType::None;
2929 const Numeric fu = band.CutoffFreq(i, X.D0);
2930 const Numeric fl = band.CutoffFreqMinus(i, X.D0);
2931
2932 // Only for the cutoff-range
2933 const auto [cutstart, cutsize] = limited_range(fl, fu, com.f_grid);
2934 if (not cutsize) return;
2935
2936 // Get the compute data view
2937 ComputeValues comval(com.F, com.dF, com.N, com.dN, com.f_grid, cutstart, cutsize, derivs, do_nlte);
2938
2939 const Index nz = do_zeeman ? band.ZeemanCount(i, zeeman_polarization) : 1;
2940 for (Index iz = 0; iz < nz; iz++) {
2941 const Numeric dfdH = do_zeeman ? band.ZeemanSplitting(i, zeeman_polarization, iz) : 0;
2942 const Numeric Sz = do_zeeman ? band.ZeemanStrength(i, zeeman_polarization, iz) : 1;
2943 const Complex LM = Complex(1 + X.G, -X.Y);
2944 Calculator ls(band.LineShapeType(), band.F0(i), X, DC, dfdH * H, band.Mirroring() == Absorption::MirroringType::Manual);
2945 Calculator ls_mirr(band.Mirroring(), band.LineShapeType(), band.F0(i), X, DC, dfdH * H);
2946
2947 if (do_cutoff) {
2948 // Initialize and set the cutoff values
2949 Calculator ls_cut = ls;
2950 Calculator ls_mirr_cut = ls_mirr;
2951 ls_cut(fu);
2952 ls_mirr_cut(fu);
2953
2954 cutoff_frequency_loop(comval, ls, ls_mirr, ls_norm, ls_cut, ls_mirr_cut, ls_str,
2955 derivs, LM, T, dfdH, Sz, band.Species());
2956
2957 } else {
2958 frequency_loop(comval, ls, ls_mirr, ls_norm, ls_str,
2959 derivs, LM, T, dfdH, Sz, band.Species());
2960 }
2961 }
2962}
2963
2995void line_loop(ComputeData &com,
2996 ComputeData &sparse_com,
2997 const AbsorptionLines &band,
2998 const ArrayOfRetrievalQuantity &jacobian_quantities,
2999 const EnergyLevelMap &nlte,
3000 const Vector &vmrs,
3001 const ArrayOfSpeciesTag& self_tag,
3002 const Numeric &P,
3003 const Numeric &T,
3004 const Numeric &H,
3005 const Numeric &sparse_lim,
3006 const Numeric QT,
3007 const Numeric QT0,
3008 const Numeric dQTdT,
3009 const Numeric r,
3010 const Numeric drdSELFVMR,
3011 const Numeric drdT,
3012 const bool do_zeeman,
3013 const Zeeman::Polarization zeeman_polarization,
3014 const Options::LblSpeedup speedup_type) ARTS_NOEXCEPT {
3015 const Index nj = jacobian_quantities.nelem();
3016 const Index nl = band.NumLines();
3017
3018 // Derivatives are allocated ahead of all loops
3019 ArrayOfDerivatives derivs(nj);
3020
3021 // Doppler constant
3022 const Numeric DC = band.DopplerConstant(T);
3023
3024 for (Index i = 0; i < nl; i++) {
3025 // Pre-compute the derivatives
3026 for (Index ij = 0; ij < nj; ij++) {
3027 const auto& deriv = jacobian_quantities[ij];
3028 derivs[ij].jac_pos = -1;
3029 derivs[ij].deriv = nullptr;
3030
3031 if (not deriv.propmattype()) continue;
3032 derivs[ij].jac_pos = ij;
3033 derivs[ij].deriv = &deriv;
3034
3035 if (deriv == Jacobian::Atm::Temperature) {
3036 derivs[ij].value.o = band.ShapeParameters_dT(i, T, P, vmrs);
3037 } else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR) {
3038 if (not (deriv == self_tag)) { // Remove if its not good
3039 derivs[ij].jac_pos = -1;
3040 derivs[ij].deriv = nullptr;
3041 }
3042 } else if (deriv.Target().needQuantumIdentity()) {
3043 if (deriv == Jacobian::Line::VMR) {
3044 derivs[ij].value.o = band.ShapeParameters_dVMR(i, T, P, deriv.QuantumIdentity());
3045 } else {
3046 auto &lt =
3047 derivs[ij].target = {deriv.Target().QuantumIdentity(), band, i};
3048 if (lt == Absorption::QuantumIdentifierLineTargetType::Line) {
3049 if constexpr (false) {/*skip so the rest can be a else-if block*/}
3050 // All line shape derivatives
3051 InternalDerivativesSetup(G0)
3052 InternalDerivativesSetup(D0)
3053 InternalDerivativesSetup(G2)
3054 InternalDerivativesSetup(D2)
3055 InternalDerivativesSetup(ETA)
3056 InternalDerivativesSetup(FVC)
3057 InternalDerivativesSetup(Y)
3058 InternalDerivativesSetup(G)
3059 InternalDerivativesSetup(DV)
3060 }
3061 }
3062 }
3063 }
3064 std::remove_if(derivs.begin(), derivs.end(), [](Derivatives& dd) { return dd.deriv == nullptr; });
3065
3066 // Call cut off loop with or without sparsity
3067 switch (speedup_type) {
3068 case Options::LblSpeedup::None:
3069 cutoff_loop(com,
3070 Normalizer(band.Normalization(), band.F0(i), T),
3071 IntensityCalculator(T, QT, QT0, dQTdT, r, drdSELFVMR, drdT, nlte, band, i),
3072 band, derivs, band.ShapeParameters(i, T, P, vmrs), T, H,
3073 DC, i, do_zeeman, zeeman_polarization);
3074 break;
3075 case Options::LblSpeedup::QuadraticIndependent:
3076 cutoff_loop_sparse_triple(com, sparse_com,
3077 Normalizer(band.Normalization(), band.F0(i), T),
3078 IntensityCalculator(T, QT, QT0, dQTdT, r, drdSELFVMR, drdT, nlte, band, i),
3079 band, derivs, band.ShapeParameters(i, T, P, vmrs), T, H, sparse_lim,
3080 DC, i, do_zeeman, zeeman_polarization);
3081 break;
3082 case Options::LblSpeedup::LinearIndependent:
3083 cutoff_loop_sparse_linear(com, sparse_com,
3084 Normalizer(band.Normalization(), band.F0(i), T),
3085 IntensityCalculator(T, QT, QT0, dQTdT, r, drdSELFVMR, drdT, nlte, band, i),
3086 band, derivs, band.ShapeParameters(i, T, P, vmrs), T, H, sparse_lim,
3087 DC, i, do_zeeman, zeeman_polarization);
3088 break;
3089 case Options::LblSpeedup::FINAL: { /* Leave last */ }
3090 }
3091 }
3092}
3093
3116void compute(ComputeData &com,
3117 ComputeData &sparse_com,
3118 const AbsorptionLines &band,
3119 const ArrayOfRetrievalQuantity &jacobian_quantities,
3120 const EnergyLevelMap &nlte,
3121 const Vector &vmrs,
3122 const ArrayOfSpeciesTag& self_tag,
3123 const Numeric& self_vmr, const Numeric &isot_ratio, const Numeric &P, const Numeric &T, const Numeric &H,
3124 const Numeric &sparse_lim,
3125 const bool do_zeeman, const Zeeman::Polarization zeeman_polarization,
3126 const Options::LblSpeedup speedup_type) ARTS_NOEXCEPT {
3127 [[maybe_unused]] const Index nj = jacobian_quantities.nelem();
3128 const Index nl = band.NumLines();
3129 const Index nv = com.f_grid.nelem();
3130
3131 // Tests that must be true while calling this function
3132 ARTS_ASSERT(H >= 0, "Only for positive H. You provided: ", H)
3133 ARTS_ASSERT(P > 0, "Only for abs positive P. You provided: ", P)
3134 ARTS_ASSERT(T > 0, "Only for abs positive T. You provided: ", T)
3135 ARTS_ASSERT(band.OK(), "Band is poorly constructed. You need to use "
3136 "a detailed debugger to find out why.")
3137 ARTS_ASSERT(com.F.size() == nv, "F is wrong size. Size is (", com.F.size(),
3138 ") but should be: (", nv, ')')
3139 ARTS_ASSERT(not com.do_nlte or com.N.size() == nv, "N is wrong size. Size is (",
3140 com.N.size(), ") but should be (", nv, ')')
3141 ARTS_ASSERT(nj == 0 or (com.dF.nrows() == nv and com.dF.ncols() == nj),
3142 "dF is wrong size. Size is (", com.dF.nrows(), " x ", com.dF.ncols(),
3143 ") but should be: (", nv, " x ", nj, ")")
3144 ARTS_ASSERT(nj == 0 or not com.do_nlte or (com.dN.nrows() == nv and com.dN.ncols() == nj),
3145 "dN is wrong size. Size is (", com.dN.nrows(), " x ", com.dN.ncols(),
3146 ") but should be: (", nv, " x ", nj, ")")
3147 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")
3148
3149 // Early return test
3150 if (nv == 0 or nl == 0 or (Absorption::relaxationtype_relmat(band.Population()) and band.DoLineMixing(P))) {
3151 return; // No line-by-line computations required/wanted
3152 }
3153
3154 const Numeric dnumdensdVMR = isot_ratio * number_density(P, T);
3155 line_loop(com, sparse_com, band, jacobian_quantities, nlte, vmrs, self_tag, P, T, H, sparse_lim,
3156 single_partition_function(T, band.Isotopologue()),
3157 single_partition_function(band.T0(), band.Isotopologue()),
3158 dsingle_partition_function_dT(T, band.Isotopologue()),
3159 self_vmr * dnumdensdVMR, dnumdensdVMR, self_vmr * isot_ratio * dnumber_density_dt(P, T),
3160 do_zeeman, zeeman_polarization, speedup_type);
3161}
3162
3163#undef InternalDerivatives
3164#undef InternalDerivativesG
3165#undef InternalDerivativesY
3166
3167
3168Index sparse_f_grid_red(const Vector& f_grid, const Numeric& sparse_df) noexcept {
3169 if (f_grid.nelem())
3170 return f_grid.nelem() / Index(1 + std::abs(f_grid[f_grid.nelem() - 1] - f_grid[0]) / sparse_df);
3171 return 0;
3172}
3173
3174Vector linear_sparse_f_grid(const Vector& f_grid, const Numeric& sparse_df) ARTS_NOEXCEPT {
3175 const Index n = sparse_f_grid_red(f_grid, sparse_df);
3176 const Index nv = f_grid.nelem();
3177
3178 if (nv and n) {
3179 std::vector<Numeric> sparse_f_grid;
3180 for (Index iv=0; iv<nv-n; iv+=n) {
3181 sparse_f_grid.emplace_back(f_grid[iv]);
3182 sparse_f_grid.emplace_back(f_grid[iv + n]);
3183 }
3184
3185 const Numeric f0 = sparse_f_grid.back();
3186 if (f0 not_eq f_grid[nv-1]) {
3187 sparse_f_grid.emplace_back(f0);
3188 sparse_f_grid.emplace_back(f_grid[nv-1]);
3189 }
3190
3191 return sparse_f_grid;
3192 }
3193 return Vector(0);
3194}
3195
3196bool good_linear_sparse_f_grid(const Vector& f_grid_dense, const Vector& f_grid_sparse) noexcept {
3197 const Index nf_sparse=f_grid_sparse.nelem();
3198 const Index nf_dense=f_grid_dense.nelem();
3199
3200 if (nf_sparse == 1)
3201 return false;
3202
3203 if(nf_sparse and nf_dense)
3204 return f_grid_dense[0] >= f_grid_sparse[0] and f_grid_dense[nf_dense-1] <= f_grid_sparse[nf_sparse-1];
3205
3206 return true;
3207}
3208
3209Vector triple_sparse_f_grid(const Vector& f_grid, const Numeric& sparse_df) noexcept {
3210 const Index n = sparse_f_grid_red(f_grid, sparse_df);
3211 const Index nv = f_grid.nelem();
3212
3213 if (nv and n > 2) {
3214 std::vector<Numeric> sparse_f_grid;
3215 for (Index iv=0; iv<nv-n; iv+=n) {
3216 sparse_f_grid.emplace_back(f_grid[iv]);
3217 sparse_f_grid.emplace_back(f_grid[iv] + 0.5 * (f_grid[iv+n] - f_grid[iv]));
3218 sparse_f_grid.emplace_back(f_grid[iv+n]);
3219 }
3220
3221 const Numeric f0 = sparse_f_grid.back();
3222 if (f0 not_eq f_grid[nv-1]) {
3223 sparse_f_grid.emplace_back(f0);
3224 sparse_f_grid.emplace_back(f0 + 0.5 * (f_grid[nv-1] - f0));
3225 sparse_f_grid.emplace_back(f_grid[nv-1]);
3226 }
3227
3228 return sparse_f_grid;
3229 }
3230 return Vector(0);
3231}
3232
3233void ComputeData::interp_add_even(const ComputeData& sparse) ARTS_NOEXCEPT {
3234 const Index nv = f_grid.nelem();
3235 const Index sparse_nv = sparse.f_grid.nelem();
3236 const Index nj = dF.ncols();
3237
3238 ARTS_ASSERT(do_nlte == sparse.do_nlte, "Must have the same NLTE status")
3239 ARTS_ASSERT(sparse_nv > 1, "Must have at least two sparse grid-points")
3240 ARTS_ASSERT(nv == 0 or (f_grid[0] == sparse.f_grid[0] and f_grid[nv - 1] >= sparse.f_grid[sparse_nv - 1]),
3241 "If there are any dense frequency points, then the sparse frequency points must fully cover them")
3242 ARTS_ASSERT(not (sparse_nv % 2), "Must be multiple of to")
3243
3244 Index sparse_iv=0;
3245 Numeric f1 = sparse.f_grid[sparse_iv + 1];
3246 Numeric f0 = sparse.f_grid[sparse_iv];
3247 Numeric inv = 1.0 / (f1 - f0);
3248 for (Index iv=0; iv<nv; iv++) {
3249 if (sparse_iv < (sparse_nv - 2) and f1 == f_grid[iv]) {
3250 sparse_iv += 2;
3251 f1 = sparse.f_grid[sparse_iv + 1];
3252 f0 = sparse.f_grid[sparse_iv];
3253 inv = 1.0 / (f1 - f0);
3254 }
3255
3256 const Numeric xm0 = f_grid[iv] - f0;
3257 const Numeric xm1 = f_grid[iv] - f1;
3258 const Numeric l0 = - xm1 * inv;
3259 const Numeric l1 = xm0 * inv;
3260
3261 F[iv] += l0 * sparse.F[sparse_iv] + l1 * sparse.F[sparse_iv + 1];
3262 for (Index ij=0; ij<nj; ij++) {
3263 dF(iv, ij) += l0 * sparse.dF(sparse_iv, ij) + l1 * sparse.dF(sparse_iv + 1, ij);
3264 }
3265 if (do_nlte) {
3266 N[iv] += l0 * sparse.N[sparse_iv] + l1 * sparse.N[sparse_iv + 1];
3267 for (Index ij=0; ij<nj; ij++) {
3268 dN(iv, ij) += l0 * sparse.dN(sparse_iv, ij) + l1 * sparse.dN(sparse_iv + 1, ij);
3269 }
3270 }
3271 }
3272}
3273
3274void ComputeData::interp_add_triplequad(const ComputeData& sparse) ARTS_NOEXCEPT {
3275 const Index nv = f_grid.nelem();
3276 const Index sparse_nv = sparse.f_grid.nelem();
3277 const Index nj = dF.ncols();
3278
3279 ARTS_ASSERT(do_nlte == sparse.do_nlte, "Must have the same NLTE status")
3280 ARTS_ASSERT(sparse_nv > 2, "Must have at least three sparse grid-points")
3281 ARTS_ASSERT(nv == 0 or (f_grid[0] == sparse.f_grid[0] and f_grid[nv - 1] >= sparse.f_grid[sparse_nv - 1]),
3282 "If there are any dense frequency points, then the sparse frequency points must fully cover them")
3283 ARTS_ASSERT(not (sparse_nv % 3), "Must be multiple of three")
3284
3285 Index sparse_iv=0;
3286 Numeric f2 = sparse.f_grid[sparse_iv + 2];
3287 Numeric f1 = sparse.f_grid[sparse_iv + 1];
3288 Numeric f0 = sparse.f_grid[sparse_iv];
3289 Numeric inv = 1.0 / Constant::pow2(f1 - f0);
3290 for (Index iv = 0; iv < nv; iv++) {
3291 if (sparse_iv < (sparse_nv - 3) and f2 == f_grid[iv]) {
3292 sparse_iv += 3;
3293 f2 = sparse.f_grid[sparse_iv + 2];
3294 f1 = sparse.f_grid[sparse_iv + 1];
3295 f0 = sparse.f_grid[sparse_iv];
3296 inv = 1.0 / Constant::pow2(f1 - f0);
3297 }
3298 ARTS_ASSERT (f_grid[iv] >= f0 and (f_grid[iv] < f2 or (f2 == f_grid[iv] and sparse_iv == sparse_nv - 3)),
3299 "Out of range frequency grid. Must be caught earlier.\n"
3300 "The sparse range is from: ", f0, " to ", f2, " with ", f1, " as the half-way grid point.\n"
3301 "The dense frequency is ", f_grid[iv], " and the indices are: sparse_iv=", sparse_iv, "; iv=", iv)
3302
3303 const Numeric xm0 = f_grid[iv] - f0;
3304 const Numeric xm1 = f_grid[iv] - f1;
3305 const Numeric xm2 = f_grid[iv] - f2;
3306 const Numeric l0 = 0.5 * xm1 * xm2 * inv; // --
3307 const Numeric l1 = - xm0 * xm2 * inv; // +-
3308 const Numeric l2 = 0.5 * xm0 * xm1 * inv; // ++
3309
3310 F[iv] += l0 * sparse.F[sparse_iv] + l1 * sparse.F[sparse_iv + 1] + l2 * sparse.F[sparse_iv + 2];
3311 for (Index ij=0; ij<nj; ij++) {
3312 dF(iv, ij) += l0 * sparse.dF(sparse_iv, ij) + l1 * sparse.dF(sparse_iv + 1, ij) + l2 * sparse.dF(sparse_iv + 2, ij);
3313 }
3314 if (do_nlte) {
3315 N[iv] += l0 * sparse.N[sparse_iv] + l1 * sparse.N[sparse_iv + 1] + l2 * sparse.N[sparse_iv + 2];
3316 for (Index ij=0; ij<nj; ij++) {
3317 dN(iv, ij) += l0 * sparse.dN(sparse_iv, ij) + l1 * sparse.dN(sparse_iv + 1, ij) + l2 * sparse.dN(sparse_iv + 2, ij);
3318 }
3319 }
3320 }
3321}
3322} // namespace LineShape
E0
Numeric E0
Definition: arts_api_classes.cc:271
F0
F0
Definition: arts_api_classes.cc:269
Array
This can be used to make arrays out of anything.
Definition: array.h:107
Array::nelem
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
ComplexMatrix
The ComplexMatrix class.
Definition: matpack_complex.h:884
ComplexVector
The ComplexVector class.
Definition: matpack_complex.h:592
LineShape::Calculator
Line shape calculator.
Definition: lineshape.h:535
RetrievalQuantity
Deals with internal derivatives, Jacobian definition, and OEM calculations.
Definition: jacobian.h:345
Vector
The Vector class.
Definition: matpackI.h:876
ARTS_NOEXCEPT
#define ARTS_NOEXCEPT
Definition: debug.h:80
ARTS_ASSERT
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
abs
#define abs(x)
Definition: legacy_continua.cc:20567
dx
#define dx
Definition: legacy_continua.cc:22121
inv
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:167
stimulated_emission
Numeric stimulated_emission(Numeric T, Numeric F0)
Computes exp(-hf/kT)
Definition: linescaling.cc:33
dabsorption_nlte_rate_dT
Numeric dabsorption_nlte_rate_dT(const Numeric &gamma, const Numeric &T, const Numeric &F0, const Numeric &El, const Numeric &Eu, const Numeric &r_upp, const Numeric &r_low)
Computes temperature derivatives of (r_low - r_upp * gamma) / (1 - gamma)
Definition: linescaling.cc:139
dabsorption_nlte_rate_dF0
Numeric dabsorption_nlte_rate_dF0(const Numeric &gamma, const Numeric &T, const Numeric &r_upp, const Numeric &r_low)
Computes frequency derivative of (r_low - r_upp * gamma) / (1 - gamma)
Definition: linescaling.cc:163
dsingle_partition_function_dT
Numeric dsingle_partition_function_dT(const Numeric &T, const Species::IsotopeRecord &ir)
Computes the partition function temperature derivative.
Definition: linescaling.cc:28
dstimulated_relative_emission_dF0
Numeric dstimulated_relative_emission_dF0(const Numeric F0, const Numeric T0, const Numeric T) noexcept
Computes.
Definition: linescaling.cc:60
boltzman_ratio
Numeric boltzman_ratio(const Numeric &T, const Numeric &T0, const Numeric &E0)
Computes exp(E0/c (T - T0) / (T * T0))
Definition: linescaling.cc:97
dabsorption_nlte_rate_dTl
Numeric dabsorption_nlte_rate_dTl(const Numeric &gamma, const Numeric &T, const Numeric &Tl, const Numeric &El, const Numeric &r_low)
Computes lower state temperature derivative of (r_low - r_upp * gamma) / (1 - gamma)
Definition: linescaling.cc:174
single_partition_function
Numeric single_partition_function(const Numeric &T, const Species::IsotopeRecord &ir)
Computes the partition function at one temperature.
Definition: linescaling.cc:23
dboltzman_ratio_dT
Numeric dboltzman_ratio_dT(const Numeric &boltzmann_ratio, const Numeric &T, const Numeric &E0)
Computes temperature derivatives exp(E0/k (T - T0) / (T * T0))
Definition: linescaling.cc:104
stimulated_relative_emission
Numeric stimulated_relative_emission(const Numeric F0, const Numeric T0, const Numeric T) noexcept
Computes.
Definition: linescaling.cc:51
absorption_nlte_ratio
Numeric absorption_nlte_ratio(const Numeric &gamma, const Numeric &r_upp, const Numeric &r_low) noexcept
Computes (r_low - r_upp * gamma) / (1 - gamma)
Definition: linescaling.cc:133
dabsorption_nlte_rate_dTu
Numeric dabsorption_nlte_rate_dTu(const Numeric &gamma, const Numeric &T, const Numeric &Tu, const Numeric &Eu, const Numeric &r_upp)
Computes upper state temperature derivative of (r_low - r_upp * gamma) / (1 - gamma)
Definition: linescaling.cc:187
dstimulated_relative_emission_dT
Numeric dstimulated_relative_emission_dT(const Numeric F0, const Numeric T0, const Numeric T) noexcept
Computes.
Definition: linescaling.cc:55
dboltzman_ratio_dT_div_boltzmann_ratio
constexpr Numeric dboltzman_ratio_dT_div_boltzmann_ratio(Numeric T, Numeric E0)
Computes temperature derivatives exp(E0/k (T - T0) / (T * T0)) / exp(E0/c (T - T0) / (T * T0))
Definition: linescaling.h:191
CutInternalDerivatives
#define CutInternalDerivatives(X)
Definition: lineshape.cc:1962
InternalDerivatives
#define InternalDerivatives(X)
Definition: lineshape.cc:1977
InternalDerivativesY
#define InternalDerivativesY
Definition: lineshape.cc:2019
InternalDerivativesSetup
#define InternalDerivativesSetup(X)
Definition: lineshape.cc:1993
ln_16
constexpr Numeric ln_16
Definition: lineshape.cc:20
InternalDerivativesG
#define InternalDerivativesG
Definition: lineshape.cc:2006
fac
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:63
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
real_val
Numeric & real_val(Complex &c) noexcept
Return a reference to the real value of c.
Definition: matpack_complex.h:39
conj
constexpr Complex conj(Complex c) noexcept
the conjugate of c
Definition: matpack_complex.h:80
imag_val
Numeric & imag_val(Complex &c) noexcept
Return a reference to the imaginary value of c.
Definition: matpack_complex.h:42
Complex
std::complex< Numeric > Complex
Definition: matpack_complex.h:35
ARTS::Continua::init
Workspace & init(Workspace &ws)
Definition: test_cpp_api.cc:147
Absorption::relaxationtype_relmat
constexpr bool relaxationtype_relmat(PopulationType in) noexcept
Definition: absorptionlines.h:154
Constant::pow3
constexpr auto pow3(T x) -> decltype(pow2(x) *x)
power of three
Definition: constants.h:71
Constant::pow4
constexpr auto pow4(T x) -> decltype(pow2(pow2(x)))
power of four
Definition: constants.h:77
Constant::pow2
constexpr auto pow2(T x) -> decltype(x *x)
power of two
Definition: constants.h:65
Jacobian::Temperature
Temperature
Definition: jacobian.h:57
LineShape
Computations of line shape derived parameters.
Definition: lineshape.cc:23
LineShape::frequency_loop
void frequency_loop(ComputeValues &com, Calculator &ls, Calculator &ls_mirr, Normalizer &ls_norm, const IntensityCalculator &ls_str, const ArrayOfDerivatives &derivs, const Complex LM, const Numeric &T, const Numeric &dfdH, const Numeric &Sz, const Species::Species self_species) ARTS_NOEXCEPT
Frequency loop of the line shape call.
Definition: lineshape.cc:2529
LineShape::sparse_f_grid_red
Index sparse_f_grid_red(const Vector &f_grid, const Numeric &sparse_df) noexcept
Definition: lineshape.cc:3168
LineShape::cutoff_frequency_loop
void cutoff_frequency_loop(ComputeValues &com, Calculator &ls, Calculator &ls_mirr, Normalizer &ls_norm, const Calculator &ls_cut, const Calculator &ls_mirr_cut, const IntensityCalculator &ls_str, const ArrayOfDerivatives &derivs, const Complex LM, const Numeric &T, const Numeric &dfdH, const Numeric &Sz, const Species::Species self_species) ARTS_NOEXCEPT
Cutoff frequency loop of the line shape call.
Definition: lineshape.cc:2321
LineShape::active_nelem
Index active_nelem(const ArrayOfDerivatives &derivs) noexcept
Helper function to find the last relevant derivative.
Definition: lineshape.cc:2190
LineShape::cutoff_loop_sparse_linear
void cutoff_loop_sparse_linear(ComputeData &com, ComputeData &sparse_com, Normalizer ls_norm, const IntensityCalculator ls_str, const AbsorptionLines &band, const ArrayOfDerivatives &derivs, const Output X, const Numeric &T, const Numeric &H, const Numeric &sparse_lim, const Numeric &DC, const Index i, const bool do_zeeman, const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT
Cutoff considerations of the line shape.
Definition: lineshape.cc:2712
LineShape::abs_squared
constexpr Numeric abs_squared(Complex z) noexcept
Definition: lineshape.cc:411
LineShape::mirroredOutput
constexpr Output mirroredOutput(Output x) noexcept
Output to be used by mirroring calls.
Definition: lineshapemodel.h:482
LineShape::line_loop
void line_loop(ComputeData &com, ComputeData &sparse_com, const AbsorptionLines &band, const ArrayOfRetrievalQuantity &jacobian_quantities, const EnergyLevelMap &nlte, const Vector &vmrs, const ArrayOfSpeciesTag &self_tag, const Numeric &P, const Numeric &T, const Numeric &H, const Numeric &sparse_lim, const Numeric QT, const Numeric QT0, const Numeric dQTdT, const Numeric r, const Numeric drdSELFVMR, const Numeric drdT, const bool do_zeeman, const Zeeman::Polarization zeeman_polarization, const Options::LblSpeedup speedup_type) ARTS_NOEXCEPT
Loop all the lines of the band.
Definition: lineshape.cc:2995
LineShape::cutoff_loop_sparse_triple
void cutoff_loop_sparse_triple(ComputeData &com, ComputeData &sparse_com, Normalizer ls_norm, const IntensityCalculator ls_str, const AbsorptionLines &band, const ArrayOfDerivatives &derivs, const Output X, const Numeric &T, const Numeric &H, const Numeric &sparse_lim, const Numeric &DC, const Index i, const bool do_zeeman, const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT
Definition: lineshape.cc:2794
LineShape::good_linear_sparse_f_grid
bool good_linear_sparse_f_grid(const Vector &f_grid_dense, const Vector &f_grid_sparse) noexcept
Definition: lineshape.cc:3196
LineShape::quad_sparse_limited_range
SparseLimitRange quad_sparse_limited_range(const Numeric flc, const Numeric fuc, const Numeric fls, const Numeric fus, const Vector &f_grid, const Vector &sparse_f_grid) ARTS_NOEXCEPT
Definition: lineshape.cc:2103
LineShape::cutoff_loop
void cutoff_loop(ComputeData &com, Normalizer ls_norm, const IntensityCalculator ls_str, const AbsorptionLines &band, const ArrayOfDerivatives &derivs, const Output X, const Numeric &T, const Numeric &H, const Numeric &DC, const Index i, const bool do_zeeman, const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT
Cutoff considerations of the line shape.
Definition: lineshape.cc:2914
LineShape::vmrs
Vector vmrs(const ConstVectorView &atmospheric_vmrs, const ArrayOfArrayOfSpeciesTag &atmospheric_species, const ArrayOfSpecies &lineshape_species) ARTS_NOEXCEPT
Returns a VMR vector for this model's main calculations.
Definition: lineshapemodel.cc:456
LineShape::triple_sparse_f_grid
Vector triple_sparse_f_grid(const Vector &f_grid, const Numeric &sparse_df) noexcept
Definition: lineshape.cc:3209
LineShape::linear_sparse_limited_range
SparseLimitRange linear_sparse_limited_range(const Numeric flc, const Numeric fuc, const Numeric fls, const Numeric fus, const Vector &f_grid, const Vector &sparse_f_grid) ARTS_NOEXCEPT
Definition: lineshape.cc:2054
LineShape::linear_sparse_f_grid
Vector linear_sparse_f_grid(const Vector &f_grid, const Numeric &sparse_df) ARTS_NOEXCEPT
Definition: lineshape.cc:3174
LineShape::limited_range
CutoffRange limited_range(const Numeric fl, const Numeric fu, const Vector &f_grid) ARTS_NOEXCEPT
Gets the start and size of a range such that.
Definition: lineshape.cc:2039
LineShape::compute
void compute(ComputeData &com, ComputeData &sparse_com, const AbsorptionLines &band, const ArrayOfRetrievalQuantity &jacobian_quantities, const EnergyLevelMap &nlte, const Vector &vmrs, const ArrayOfSpeciesTag &self_tag, const Numeric &self_vmr, const Numeric &isot_ratio, const Numeric &P, const Numeric &T, const Numeric &H, const Numeric &sparse_lim, const bool do_zeeman, const Zeeman::Polarization zeeman_polarization, const Options::LblSpeedup speedup_type) ARTS_NOEXCEPT
Compute the line shape in its entirety.
Definition: lineshape.cc:3116
Options::QuadraticIndependent
X3 LineCenter QuadraticIndependent
Definition: constants.h:577
Options::None
X3 LineCenter None
Definition: constants.h:576
Zeeman::start
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the lowest M for a polarization type of this transition.
Definition: zeemandata.h:78
Zeeman::Polarization
Polarization
Zeeman polarization selection.
Definition: zeemandata.h:44
nonstd::abs
constexpr T abs(T x) noexcept
Definition: nonstd.h:13
oem::LM
invlib::LevenbergMarquardt< Numeric, CovarianceMatrix, Std > LM
Levenberg-Marquardt (LM) optimization using normed ARTS QR solver.
Definition: oem.h:172
oem::Vector
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
physics_funcs.h
This file contains declerations of functions of physical character.
dnumber_density_dt
constexpr Numeric dnumber_density_dt(Numeric p, Numeric t) noexcept
dnumber_density_dT
Definition: physics_funcs.h:84
number_density
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Definition: physics_funcs.h:70
d
#define d
w
#define w
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:739
N
#define N
Definition: rng.cc:164
LineShape::ComputeData::interp_add_even
void interp_add_even(const ComputeData &sparse) ARTS_NOEXCEPT
Definition: lineshape.cc:3233
LineShape::ComputeData::interp_add_triplequad
void interp_add_triplequad(const ComputeData &sparse) ARTS_NOEXCEPT
Definition: lineshape.cc:3274
LineShape::ComputeValues::derivs
const ArrayOfDerivatives & derivs
Definition: lineshape.cc:2205
LineShape::ComputeValues::ComputeValues
ComputeValues(ComplexVector &F_, ComplexMatrix &dF_, ComplexVector &N_, ComplexMatrix &dN_, const Vector &f_grid, const Index start, const Index nv, const ArrayOfDerivatives &derivs_, const bool do_nlte_) noexcept
Definition: lineshape.cc:2212
LineShape::ComputeValues::operator-=
ComputeValues & operator-=(const ComputeValues &cut) ARTS_NOEXCEPT
Definition: lineshape.cc:2237
LineShape::ComputeValues::f
const Numeric *const f
Definition: lineshape.cc:2203
LineShape::ComputeValues::ComputeValues
ComputeValues(Complex &F_, std::vector< Complex > &dF_, Complex &N_, std::vector< Complex > &dN_, const Numeric &f_lim, const ArrayOfDerivatives &derivs_, const bool do_nlte_) noexcept
Definition: lineshape.cc:2229
LineShape::ComputeValues::jac_pos
Index jac_pos(Index iv, Index ij) const noexcept
Definition: lineshape.cc:2210
LineShape::CutoffRange
Struct to keep the cutoff limited range values.
Definition: lineshape.cc:2024
LineShape::Derivatives
Data struct for keeping derivative keys and values.
Definition: lineshape.cc:2174
LineShape::Derivatives::target
Absorption::QuantumIdentifierLineTarget target
Definition: lineshape.cc:2175
LineShape::Derivatives::deriv
const RetrievalQuantity * deriv
Definition: lineshape.cc:2182
LineShape::Derivatives::value
union LineShape::Derivatives::Values value
LineShape::Doppler::operator()
Complex operator()(Numeric f) noexcept
Definition: lineshape.cc:24
LineShape::FullNonLocalThermodynamicEquilibriumInitialization::operator()
constexpr FullNonLocalThermodynamicEquilibrium operator()(Numeric r, Numeric drdSELFVMR, Numeric drdT) &&noexcept
Definition: lineshape.cc:1878
LineShape::FullNonLocalThermodynamicEquilibriumInitialization::FullNonLocalThermodynamicEquilibriumInitialization
FullNonLocalThermodynamicEquilibriumInitialization(Numeric F0, Numeric A21, Numeric T, Numeric r1, Numeric r2, Numeric c2, Numeric c3, Numeric x) noexcept
Definition: lineshape.cc:1865
LineShape::FullNonLocalThermodynamicEquilibrium::FullNonLocalThermodynamicEquilibrium
constexpr FullNonLocalThermodynamicEquilibrium(Numeric r, Numeric drdSELFVMR, Numeric drdt, Numeric k, Numeric dkdF0, Numeric dkdr1, Numeric dkdr2, Numeric e, Numeric dedF0, Numeric dedr2, Numeric B, Numeric dBdT, Numeric dBdF0) noexcept
Definition: lineshape.h:413
LineShape::HartmannTran::dFdG0
Complex dFdG0(Numeric dG0) const noexcept
Definition: lineshape.cc:749
LineShape::HartmannTran::dFdH
Complex dFdH(Numeric dZ) const noexcept
Definition: lineshape.cc:1262
LineShape::HartmannTran::dFdT
Complex dFdT(const Output &d, Numeric T) const noexcept
Definition: lineshape.cc:1519
LineShape::HartmannTran::init
CalcType init(const Complex c2t) const noexcept
Definition: lineshape.cc:1689
LineShape::HartmannTran::HartmannTran
HartmannTran(Numeric F0_noshift, const Output &ls, Numeric GD_div_F0, Numeric dZ) noexcept
Definition: lineshape.cc:468
LineShape::HartmannTran::dFdf
Complex dFdf() const noexcept
Definition: lineshape.cc:478
LineShape::HartmannTran::dFdD2
Complex dFdD2(Numeric dD2) const noexcept
Definition: lineshape.cc:829
LineShape::HartmannTran::dFdD0
Complex dFdD0(Numeric dD0) const noexcept
Definition: lineshape.cc:653
LineShape::HartmannTran::update_calcs
void update_calcs() noexcept
Definition: lineshape.cc:1702
LineShape::HartmannTran::dFdG2
Complex dFdG2(Numeric dG2) const noexcept
Definition: lineshape.cc:952
LineShape::HartmannTran::dFdF0
Complex dFdF0() const noexcept
Definition: lineshape.cc:558
LineShape::HartmannTran::dFdETA
Complex dFdETA(Numeric dETA) const noexcept
Definition: lineshape.cc:1143
LineShape::HartmannTran::dFdFVC
Complex dFdFVC(Numeric dFVC) const noexcept
Definition: lineshape.cc:1057
LineShape::HartmannTran::calc
void calc() noexcept
Definition: lineshape.cc:1706
LineShape::HartmannTran::operator()
Complex operator()(Numeric f) noexcept
Definition: lineshape.cc:1680
LineShape::HartmannTran::dFdVMR
Complex dFdVMR(const Output &d) const noexcept
Definition: lineshape.cc:1358
LineShape::LocalThermodynamicEquilibrium::LocalThermodynamicEquilibrium
constexpr LocalThermodynamicEquilibrium(Numeric I0, Numeric r, Numeric drdSELFVMR, Numeric drdT, Numeric QT0, Numeric QT, Numeric dQTdT, Numeric br, Numeric dbr_dT_rat, Numeric stim, Numeric dstim_dT, Numeric dstim_dF0) noexcept
Definition: lineshape.h:360
LineShape::RosenkranzQuadratic::dNdT
Numeric dNdT(Numeric, Numeric f) const noexcept
Definition: lineshape.cc:1808
LineShape::RosenkranzQuadratic::dNdF0
Numeric dNdF0() const noexcept
Definition: lineshape.cc:1816
LineShape::RosenkranzQuadratic::operator()
Numeric operator()(Numeric f) noexcept
Definition: lineshape.cc:1820
LineShape::RosenkranzQuadratic::dNdf
Numeric dNdf(Numeric f) const noexcept
Definition: lineshape.cc:1812
LineShape::RosenkranzQuadratic::RosenkranzQuadratic
RosenkranzQuadratic(Numeric F0, Numeric T) noexcept
Definition: lineshape.cc:1797
LineShape::SimpleFrequencyScaling::operator()
Numeric operator()(Numeric f) noexcept
Definition: lineshape.cc:1844
LineShape::SimpleFrequencyScaling::dNdT
Numeric dNdT(Numeric t_, Numeric f) const ARTS_NOEXCEPT
Definition: lineshape.cc:1825
LineShape::SimpleFrequencyScaling::dNdf
Numeric dNdf(Numeric f) const noexcept
Definition: lineshape.cc:1838
LineShape::SparseLimitRange
Struct to keep the cutoff limited range values and the sparse limits.
Definition: lineshape.cc:2048
LineShape::SpeedDependentVoigt::dFdD2
Complex dFdD2(Numeric dD2dD2) const noexcept
Definition: lineshape.cc:141
LineShape::SpeedDependentVoigt::dFdD0
Complex dFdD0(Numeric dD0dD0) const noexcept
Definition: lineshape.cc:92
LineShape::SpeedDependentVoigt::dFdH
Complex dFdH(Numeric dZ) const noexcept
Definition: lineshape.cc:217
LineShape::SpeedDependentVoigt::update_calcs
void update_calcs() noexcept
Definition: lineshape.cc:428
LineShape::SpeedDependentVoigt::dFdG2
Complex dFdG2(Numeric dG2dG2) const noexcept
Definition: lineshape.cc:186
LineShape::SpeedDependentVoigt::dFdT
Complex dFdT(const Output &d, Numeric T) const noexcept
Definition: lineshape.cc:307
LineShape::SpeedDependentVoigt::dFdf
Complex dFdf() const noexcept
Definition: lineshape.cc:47
LineShape::SpeedDependentVoigt::dFdF0
Complex dFdF0() const noexcept
Definition: lineshape.cc:64
LineShape::SpeedDependentVoigt::init
CalcType init(const Complex c2) const noexcept
Definition: lineshape.cc:416
LineShape::SpeedDependentVoigt::dFdG0
Complex dFdG0(Numeric dG0dG0) const noexcept
Definition: lineshape.cc:123
LineShape::SpeedDependentVoigt::SpeedDependentVoigt
SpeedDependentVoigt(Numeric F0_noshift, const Output &ls, Numeric GD_div_F0, Numeric dZ) noexcept
Definition: lineshape.cc:37
LineShape::SpeedDependentVoigt::dFdVMR
Complex dFdVMR(const Output &d) const noexcept
Definition: lineshape.cc:248
LineShape::SpeedDependentVoigt::operator()
Complex operator()(Numeric f) noexcept
Definition: lineshape.cc:403
LineShape::VanVleckHuber::VanVleckHuber
VanVleckHuber(Numeric F0, Numeric T) noexcept
Definition: lineshape.cc:1768
LineShape::VanVleckHuber::operator()
Numeric operator()(Numeric f) noexcept
Definition: lineshape.cc:1791
LineShape::VanVleckHuber::dNdf
Numeric dNdf(Numeric f) const noexcept
Definition: lineshape.cc:1779
LineShape::VanVleckHuber::dNdF0
Numeric dNdF0() const noexcept
Definition: lineshape.cc:1785
LineShape::VanVleckHuber::dNdT
Numeric dNdT(Numeric T, Numeric f) const noexcept
Definition: lineshape.cc:1772
LineShape::VibrationalTemperaturesNonLocalThermodynamicEquilibriumInitializer::VibrationalTemperaturesNonLocalThermodynamicEquilibriumInitializer
VibrationalTemperaturesNonLocalThermodynamicEquilibriumInitializer(Numeric T, Numeric T0, Numeric F0, Numeric E0, Numeric Tl, Numeric Evl, Numeric Tu, Numeric Evu, Numeric gamma, Numeric gamma_ref, Numeric r_low, Numeric r_upp) noexcept
Definition: lineshape.cc:1901
LineShape::VibrationalTemperaturesNonLocalThermodynamicEquilibriumInitializer::operator()
constexpr VibrationalTemperaturesNonLocalThermodynamicEquilibrium operator()(Numeric I0, Numeric QT0, Numeric QT, Numeric dQTdT, Numeric r, Numeric drdSELFVMR, Numeric drdT) &&noexcept
Definition: lineshape.cc:1921
LineShape::VibrationalTemperaturesNonLocalThermodynamicEquilibrium::VibrationalTemperaturesNonLocalThermodynamicEquilibrium
constexpr VibrationalTemperaturesNonLocalThermodynamicEquilibrium(Numeric I0, Numeric QT0, Numeric QT, Numeric dQTdT, Numeric r, Numeric drdSELFVMR, Numeric drdT, Numeric K1, Numeric dK1dT, Numeric K2, Numeric dK2dT, Numeric dK2dF0, Numeric K3, Numeric dK3dT, Numeric dK3dF0, Numeric dK3dTl, Numeric dK3dTu, Numeric K4, Numeric dK4dT, Numeric dK4dTu, Numeric B, Numeric dBdT, Numeric dBdF0) noexcept
Definition: lineshape.h:475
LineShape::Voigt::operator()
Complex operator()(Numeric f) noexcept
Definition: lineshape.cc:30
LineShape::Derivatives::Values::Values
constexpr Values() noexcept
Definition: lineshape.cc:2179