13#include <Faddeeva/Faddeeva.hh>
25 2.772588722239781237668928485832706272302000537441021016482720037973574487879;
29 x = (f - mF0) * invGD;
30 F = invGD * inv_sqrt_pi * std::exp(-pow2(x));
35 real_val(z) = invGD * (f - mF0);
36 F = inv_sqrt_pi * invGD * Faddeeva::w(z);
37 dF = 2 * invGD * (Complex(0, inv_pi * invGD) - z * F);
42 Numeric GD_div_F0, Numeric dZ) noexcept
43 : mF0(F0_noshift + dZ + ls.D0 - 1.5 * ls.D2),
45 invc2(1.0 / Complex(ls.G2, ls.D2)), dx(Complex(ls.G0 - 1.5 * ls.G2, mF0)),
46 x(dx * invc2), sqrty(invc2 / (2 * invGD)),
47 calcs(init(Complex(ls.G2, ls.D2))) {
56 return dw1 * pow2(
invGD) * inv_sqrt_pi;
58 return dw1 * pow2(
invGD) * inv_sqrt_pi -
63 return 1i * pow2(
invc2) * (
x - 3) / (pi * pow3(
x));
91 return 1i * pow2(
invc2) * (3 -
x) / (pi * pow3(
x));
100 (4 * pow2(invGD) * (w1 - w2) * sq +
102 (-dw1 * (Complex(0, 2 * mF0 * pow2(invGD)) - 2 * invGD * sq +
104 dw2 * (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
106 (4 * sqrt_pi * invGD * mF0 * sq);
107 case CalcType::Voigt:
108 return -dD0dD0 * invGD *
109 (Complex(0, invGD) * dw1 * (dx - Complex(0, mF0)) + w1) /
111 case CalcType::LowXandHighY:
113 (4 * pow2(invGD) * (w1 - w2) * sq +
114 1i * (4 * dw1 * pow3(invGD) * (dx - Complex(0, mF0)) * sq +
116 (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
118 (4 * sqrt_pi * invGD * mF0 * sq);
119 case CalcType::LowYandLowX:
120 return dD0dD0 * pow2(invc2) * (dw1 * sq - 1i * w1) / (sqrt_pi * sq);
121 case CalcType::LowYandHighX:
122 return -Complex(0, dD0dD0) * pow2(invc2) * (x - 3) / (pi * pow3(x));
130 return Complex(0, dG0dG0) * invGD * invc2 * (dw1 - dw2) /
132 case CalcType::Voigt:
133 return Complex(0, dG0dG0) * dw1 * pow2(invGD) * inv_sqrt_pi;
134 case CalcType::LowXandHighY:
135 return Complex(0, dG0dG0) * invGD * (2 * dw1 * invGD * sq - dw2 * invc2) /
137 case CalcType::LowYandLowX:
138 return -dG0dG0 * pow2(invc2) * (1i * dw1 * sq + w1) / (sqrt_pi * sq);
139 case CalcType::LowYandHighX:
140 return -dG0dG0 * pow2(invc2) * (x - 3) / (pi * pow3(x));
149 (12 * pow2(invGD) * (w1 - w2) * sq +
151 (dw1 * (-Complex(0, 2 * mF0 * pow2(invGD)) *
152 (2 * dx * invc2 + 3) +
153 4 * Complex(0, invGD) * invc2 * mF0 * sq +
154 6 * invGD * sq - Complex(0, 2 * mF0) * pow2(invc2) -
157 (Complex(0, 2 * mF0 * pow2(invGD)) * (2 * dx * invc2 + 3) +
158 4 * Complex(0, invGD) * invc2 * mF0 * sq +
159 6 * invGD * sq + Complex(0, 2 * mF0) * pow2(invc2) +
161 (8 * sqrt_pi * invGD * mF0 * sq);
162 case CalcType::Voigt:
163 return 3 * dD2dD2 * invGD *
164 (Complex(0, invGD) * dw1 * (dx - Complex(0, mF0)) + w1) /
166 case CalcType::LowXandHighY:
168 (12 * pow2(invGD) * (w1 - w2) * sq +
170 (12 * dw1 * pow3(invGD) * (dx - Complex(0, mF0)) * sq +
172 (Complex(0, 2 * mF0 * pow2(invGD)) * (2 * dx * invc2 + 3) +
173 4 * Complex(0, invGD) * invc2 * mF0 * sq +
174 6 * invGD * sq + Complex(0, 2 * mF0) * pow2(invc2) +
176 (8 * sqrt_pi * invGD * mF0 * sq);
177 case CalcType::LowYandLowX:
178 return dD2dD2 * pow2(invc2) *
179 (4 * 1i * sq * (sqrt_pi * w1 * sq - 1) -
180 sqrt_pi * (dw1 * sq - 1i * w1) * (2 * dx * invc2 + 3)) /
182 case CalcType::LowYandHighX:
183 return Complex(0, dD2dD2) * pow2(invc2) *
184 (-x * (2 * x - 3) + (x - 3) * (2 * dx * invc2 + 3)) /
193 return Complex(0, dG2dG2) * invc2 *
194 (dw1 * (-pow2(invGD) * (2 * dx * invc2 + 3) +
195 2 * invGD * invc2 * sq - pow2(invc2)) +
196 dw2 * (pow2(invGD) * (2 * dx * invc2 + 3) + 2 * invGD * invc2 * sq +
198 (4 * sqrt_pi * invGD * sq);
199 case CalcType::Voigt:
200 return -3 * Complex(0, dG2dG2) * dw1 * pow2(invGD) / (2 * sqrt_pi);
201 case CalcType::LowXandHighY:
202 return Complex(0, dG2dG2) *
203 (-6 * dw1 * pow3(invGD) * sq +
205 (pow2(invGD) * (2 * dx * invc2 + 3) + 2 * invGD * invc2 * sq +
207 (4 * sqrt_pi * invGD * sq);
208 case CalcType::LowYandLowX:
209 return dG2dG2 * pow2(invc2) *
210 (4 * sq * (sqrt_pi * w1 * sq - 1) +
211 sqrt_pi * (2 * dx * invc2 + 3) * (1i * dw1 * sq + w1)) /
213 case CalcType::LowYandHighX:
214 return dG2dG2 * pow2(invc2) *
215 (-x * (2 * x - 3) + (x - 3) * (2 * dx * invc2 + 3)) /
225 (4 * pow2(invGD) * (w1 - w2) * sq +
227 (-dw1 * (Complex(0, 2 * mF0 * pow2(invGD)) - 2 * invGD * sq +
229 dw2 * (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
231 (4 * sqrt_pi * invGD * mF0 * sq);
232 case CalcType::Voigt:
234 (Complex(0, invGD) * dw1 * (dx - Complex(0, mF0)) + w1) /
236 case CalcType::LowXandHighY:
238 (4 * pow2(invGD) * (w1 - w2) * sq +
239 1i * (4 * dw1 * pow3(invGD) * (dx - Complex(0, mF0)) * sq +
241 (Complex(0, 2 * mF0 * pow2(invGD)) + 2 * invGD * sq +
243 (4 * sqrt_pi * invGD * mF0 * sq);
244 case CalcType::LowYandLowX:
245 return dZ * pow2(invc2) * (dw1 * sq - 1i * w1) / (sqrt_pi * sq);
246 case CalcType::LowYandHighX:
247 return -Complex(0, dZ) * pow2(invc2) * (x - 3) / (pi * pow3(x));
255 return (-4 * pow2(invGD) * (2 *
d.D0 - 3 *
d.D2) * (w1 - w2) * sq +
257 (dw1 * (-2 * pow2(invGD) * mF0 *
258 (Complex(3 *
d.G2 - 2 *
d.G0, 3 *
d.D2 - 2 *
d.D0) +
259 2 * dx * invc2 * Complex(
d.G2,
d.D2)) +
260 4 * invGD * invc2 * mF0 * sq * Complex(
d.G2,
d.D2) -
261 2 * invGD * (2 *
d.D0 - 3 *
d.D2) * sq -
262 2 * pow2(invc2) * mF0 * Complex(
d.G2,
d.D2) +
263 invc2 * (2 *
d.D0 - 3 *
d.D2)) -
264 dw2 * (-2 * pow2(invGD) * mF0 *
265 (Complex(3 *
d.G2 - 2 *
d.G0, 3 *
d.D2 - 2 *
d.D0) +
266 2 * dx * invc2 * Complex(
d.G2,
d.D2)) -
267 4 * invGD * invc2 * mF0 * sq * Complex(
d.G2,
d.D2) +
268 2 * invGD * (2 *
d.D0 - 3 *
d.D2) * sq -
269 2 * pow2(invc2) * mF0 * Complex(
d.G2,
d.D2) +
270 invc2 * (2 *
d.D0 - 3 *
d.D2)))) /
271 (8 * sqrt_pi * invGD * mF0 * sq);
272 case CalcType::Voigt:
274 (Complex(0, invGD) * dw1 *
275 (dx * (2 *
d.D0 - 3 *
d.D2) -
276 mF0 * Complex(2 *
d.G0 - 3 *
d.G2, 2 *
d.D0 - 3 *
d.D2)) +
277 w1 * (2 *
d.D0 - 3 *
d.D2)) /
279 case CalcType::LowXandHighY:
280 return -(4 * pow2(invGD) * (2 *
d.D0 - 3 *
d.D2) * (w1 - w2) * sq +
281 1i * (4 * dw1 * pow3(invGD) * sq *
282 (dx * (2 *
d.D0 - 3 *
d.D2) -
284 Complex(2 *
d.G0 - 3 *
d.G2, 2 *
d.D0 - 3 *
d.D2)) +
286 (2 * pow2(invGD) * mF0 *
287 (Complex(2 *
d.G0 - 3 *
d.G2, 2 *
d.D0 - 3 *
d.D2) -
288 2 * dx * invc2 * Complex(
d.G2,
d.D2)) -
289 4 * invGD * invc2 * mF0 * sq * Complex(
d.G2,
d.D2) +
290 2 * invGD * (2 *
d.D0 - 3 *
d.D2) * sq -
291 2 * pow2(invc2) * mF0 * Complex(
d.G2,
d.D2) +
292 invc2 * (2 *
d.D0 - 3 *
d.D2)))) /
293 (8 * sqrt_pi * invGD * mF0 * sq);
294 case CalcType::LowYandLowX:
296 (4 * sq * Complex(
d.G2,
d.D2) * (sqrt_pi * w1 * sq - 1) -
297 sqrt_pi * (1i * dw1 * sq + w1) *
298 (Complex(2 *
d.G0 - 3 *
d.G2, 2 *
d.D0 - 3 *
d.D2) -
299 2 * dx * invc2 * Complex(
d.G2,
d.D2))) /
301 case CalcType::LowYandHighX:
302 return -pow2(invc2) *
303 (x * (2 * x - 3) * Complex(
d.G2,
d.D2) +
304 (x - 3) * (Complex(2 *
d.G0 - 3 *
d.G2, 2 *
d.D0 - 3 *
d.D2) -
305 2 * dx * invc2 * Complex(
d.G2,
d.D2))) /
314 return (-pow2(invGD) * (w1 - w2) * sq *
315 (T * (2 *
d.D0 - 3 *
d.D2) * sqrt_ln_2 - pow3(invGD) * mF0) *
319 (T * invGD * invc2 * mF0 * sq * Complex(
d.G2,
d.D2) *
322 (T * (2 *
d.D0 - 3 *
d.D2) * sqrt_ln_2 -
325 (2 * T * pow2(invGD) * mF0 *
326 (Complex(3 *
d.G2 - 2 *
d.G0, 3 *
d.D2 - 2 *
d.D0) +
327 2 * dx * invc2 * Complex(
d.G2,
d.D2)) *
329 2 * T * pow2(invc2) * mF0 * Complex(
d.G2,
d.D2) *
331 invc2 * (T * (-2 *
d.D0 + 3 *
d.D2) * sqrt_ln_2 +
332 pow3(invGD) * mF0)) *
335 (T * invGD * invc2 * mF0 * sq * Complex(
d.G2,
d.D2) *
338 (T * (-2 *
d.D0 + 3 *
d.D2) * sqrt_ln_2 +
341 (-2 * T * pow2(invGD) * mF0 *
342 (Complex(2 *
d.G0 - 3 *
d.G2, 2 *
d.D0 - 3 *
d.D2) -
343 2 * dx * invc2 * Complex(
d.G2,
d.D2)) *
345 2 * T * pow2(invc2) * mF0 * Complex(
d.G2,
d.D2) *
347 invc2 * (T * (-2 *
d.D0 + 3 *
d.D2) * sqrt_ln_2 +
348 pow3(invGD) * mF0)) *
351 (8 * sqrt_pi * T * invGD * mF0 * sq * pow3(sqrt_ln_2));
352 case CalcType::Voigt:
354 (-Complex(0, invGD) * dw1 *
355 (T * mF0 * Complex(2 *
d.G0 - 3 *
d.G2, 2 *
d.D0 - 3 *
d.D2) *
357 dx * (T * (2 *
d.D0 - 3 *
d.D2) * sqrt_ln_2 -
358 pow3(invGD) * mF0)) +
359 w1 * (T * (2 *
d.D0 - 3 *
d.D2) * sqrt_ln_2 - pow3(invGD) * mF0)) /
360 (2 * sqrt_pi * T * mF0 * sqrt_ln_2);
361 case CalcType::LowXandHighY:
362 return (-pow2(invGD) * (w1 - w2) * sq *
363 (T * (2 *
d.D0 - 3 *
d.D2) * sqrt_ln_2 - pow3(invGD) * mF0) *
365 1i * (dw1 * pow3(invGD) * sq *
367 Complex(2 *
d.G0 - 3 *
d.G2, 2 *
d.D0 - 3 *
d.D2) *
369 dx * (T * (2 *
d.D0 - 3 *
d.D2) * sqrt_ln_2 -
370 pow3(invGD) * mF0)) *
373 (T * invGD * invc2 * mF0 * sq * Complex(
d.G2,
d.D2) *
376 (T * (-2 *
d.D0 + 3 *
d.D2) * sqrt_ln_2 +
379 (-2 * T * pow2(invGD) * mF0 *
380 (Complex(2 *
d.G0 - 3 *
d.G2, 2 *
d.D0 - 3 *
d.D2) -
381 2 * dx * invc2 * Complex(
d.G2,
d.D2)) *
383 2 * T * pow2(invc2) * mF0 * Complex(
d.G2,
d.D2) *
385 invc2 * (T * (-2 *
d.D0 + 3 *
d.D2) * sqrt_ln_2 +
386 pow3(invGD) * mF0)) *
389 (8 * sqrt_pi * T * invGD * mF0 * sq * pow3(sqrt_ln_2));
390 case CalcType::LowYandLowX:
392 (4 * sq * Complex(
d.G2,
d.D2) * (sqrt_pi * w1 * sq - 1) -
393 sqrt_pi * (1i * dw1 * sq + w1) *
394 (Complex(2 *
d.G0 - 3 *
d.G2, 2 *
d.D0 - 3 *
d.D2) -
395 2 * dx * invc2 * Complex(
d.G2,
d.D2))) /
397 case CalcType::LowYandHighX:
398 return -pow2(invc2) *
399 (x * (2 * x - 3) * Complex(
d.G2,
d.D2) +
400 (x - 3) * (Complex(2 *
d.G0 - 3 *
d.G2, 2 *
d.D0 - 3 *
d.D2) -
401 2 * dx * invc2 * Complex(
d.G2,
d.D2))) /
408 imag_val(dx) = mF0 - f;
416 return pow2(z.real()) + pow2(z.imag());
422 return CalcType::Voigt;
424 return CalcType::LowXandHighY;
427 return CalcType::LowYandLowX;
429 return CalcType::LowYandHighX;
430 return CalcType::Full;
463 w1 = Faddeeva::w(1i *
sq);
464 F = 2 * inv_pi *
invc2 * (1 - sqrt_pi *
sq *
w1);
465 dw1 = 2i * (inv_sqrt_pi -
sq *
w1);
468 F = inv_pi *
invc2 * (1 /
x - 1.5 / pow2(
x));
474 Numeric GD_div_F0, Numeric dZ) noexcept
475 :
G0(ls.G0),
D0(ls.D0),
G2(ls.G2),
D2(ls.D2),
FVC(ls.FVC),
ETA(ls.ETA),
476 mF0(F0_noshift + dZ + (1 - ls.ETA) * (ls.D0 - 1.5 * ls.D2)),
478 deltax(ls.FVC + (1 - ls.ETA) * (ls.G0 - 3 * ls.G2 / 2), mF0),
479 sqrty(1 / (2 * (1 - ls.ETA) * Complex(ls.G2, ls.D2) * invGD)) {
484 constexpr Complex ddeltax{0, -1};
485 const Complex dx = -ddeltax / ((
ETA - 1) * Complex(
G2,
D2));
486 const Complex dsqrtxy = dx / (2 *
sqrtxy);
490 const Complex dz1 = dsqrtxy;
491 const Complex dz2 = dsqrtxy;
492 const Complex dA = Complex(0, sqrt_pi *
invGD) * (
dw1 * dz1 -
dw2 * dz2);
495 ((pow2(
z1) - 1) * 1i *
dw1 * dz1 - (pow2(
z2) - 1) * 1i *
dw2 * dz2 +
496 2 *
w1 *
z1 * dz1 - 2 *
w2 *
z2 * dz2) /
498 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
500 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
503 const Complex dz1 =
invGD * ddeltax;
504 const Complex dA = Complex(0, sqrt_pi *
invGD) *
dw1 * dz1;
507 (sqrt_pi * ((pow2(
z1) - 1) * 1i *
dw1 * dz1 + 2 *
w1 *
z1 * dz1) - dz1);
508 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
510 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
513 const Complex dz1 =
invGD * ddeltax;
514 const Complex dA = Complex(0, sqrt_pi *
invGD) *
dw1 * dz1;
515 const Complex dB = Complex(0, sqrt_pi *
invGD) *
dw1 * dz1 -
517 9 *
invGD * dz1 / (4 * pow4(
z1));
518 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
520 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
523 const Complex dz1 =
invGD * ddeltax;
524 const Complex dz2 = dsqrtxy;
525 const Complex dA = Complex(0, sqrt_pi *
invGD) * (
dw1 * dz1 -
dw2 * dz2);
526 const Complex dB = Complex(0, sqrt_pi *
invGD) *
dw1 * dz1 -
528 9 *
invGD * dz1 / (4 * pow4(
z1));
529 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
531 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
534 const Complex dz1 = dsqrtxy;
535 const Complex dz2 = dx / (2 *
sqrtx);
536 const Complex dA = 2 * sqrt_pi * (
w2 * dz2 +
z2 * 1i *
dw2 * dz2) /
539 -(2 * sqrt_pi * (
w2 * dz2 +
z2 * 1i *
dw2 * dz2) *
540 (2 * pow2(
sqrty) +
x - 1) +
541 2 * sqrt_pi *
w1 * dz1 + Complex(0, 2 * sqrt_pi) *
z1 *
dw1 * dz1 +
542 2 * (sqrt_pi *
w2 *
z2 - 1) * dx) /
544 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
546 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
549 const Complex dz1 = dsqrtxy;
550 const Complex dA = (
x - 3) * dx / ((
ETA - 1) * Complex(
G2,
D2) * pow3(
x));
552 (-2 * sqrt_pi * (
w1 * dz1 +
z1 * 1i *
dw1 * dz1) * pow3(
x) -
553 (
x - 3) * (2 * pow2(
sqrty) +
x - 1) * dx + (2 *
x - 3) *
x * dx / 2) /
554 ((
ETA - 1) * Complex(
G2,
D2) * pow3(
x));
555 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
557 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
564 const Numeric dGD = (1 / (
invGD *
mF0));
565 const Numeric dinvGD = -dGD * pow2(
invGD);
566 const Complex dsqrty =
567 dinvGD / (2 * (
ETA - 1) * Complex(
G2,
D2) * pow2(
invGD));
568 constexpr Complex ddeltax = 1i;
569 const Complex dx = -ddeltax / ((
ETA - 1) * Complex(
G2,
D2));
570 const Complex dsqrtxy = (
sqrty * dsqrty + dx / 2) /
sqrtxy;
574 const Complex dz1 = dsqrtxy - dsqrty;
575 const Complex dz2 = dsqrtxy + dsqrty;
581 ((-(pow2(
z1) - 1) *
w1 + (pow2(
z2) - 1) *
w2) * dsqrty +
582 ((pow2(
z1) - 1) * 1i *
dw1 * dz1 - (pow2(
z2) - 1) * 1i *
dw2 * dz2 +
583 2 *
w1 *
z1 * dz1 - 2 *
w2 *
z2 * dz2) *
586 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
588 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
591 const Complex dz1 =
deltax * dinvGD +
invGD * ddeltax;
592 const Complex dA = sqrt_pi * (Complex(0,
invGD) *
dw1 * dz1 +
w1 * dinvGD);
594 -(sqrt_pi * ((pow2(
z1) - 1) * 1i *
dw1 * dz1 + 2 *
w1 *
z1 * dz1) -
597 (sqrt_pi * (pow2(
z1) - 1) *
w1 -
z1) * dinvGD;
598 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
600 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
603 const Complex dz1 =
deltax * dinvGD +
invGD * ddeltax;
604 const Complex dA = sqrt_pi * (Complex(0,
invGD) *
dw1 * dz1 +
w1 * dinvGD);
606 ((4 * sqrt_pi *
w1 * pow3(
z1) + 2 * pow2(
z1) - 3) *
z1 * dinvGD +
607 (Complex(0, 4 * sqrt_pi) * pow4(
z1) *
dw1 * dz1 - 2 * pow2(
z1) * dz1 +
611 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
613 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
616 const Complex dz1 =
deltax * dinvGD +
invGD * ddeltax;
617 const Complex dz2 = dsqrtxy + dsqrty;
622 ((4 * sqrt_pi *
w1 * pow3(
z1) + 2 * pow2(
z1) - 3) *
z1 * dinvGD +
623 (Complex(0, 4 * sqrt_pi) * pow4(
z1) *
dw1 * dz1 - 2 * pow2(
z1) * dz1 +
627 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
629 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
632 const Complex dz1 = dsqrtxy;
633 const Complex dz2 = dx / (2 *
sqrtx);
634 const Complex dA = 2 * sqrt_pi * (
w2 * dz2 +
z2 * 1i *
dw2 * dz2) /
637 -(2 * sqrt_pi * (
w2 * dz2 +
z2 * 1i *
dw2 * dz2) *
638 (2 * pow2(
sqrty) +
x - 1) +
639 2 * sqrt_pi *
w1 * dz1 + Complex(0, 2 * sqrt_pi) *
z1 *
dw1 * dz1 +
640 2 * (4 *
sqrty * dsqrty + dx) * (sqrt_pi *
w2 *
z2 - 1)) /
642 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
644 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
647 const Complex dz1 = dsqrtxy;
648 const Complex dA = (
x - 3) * dx / ((
ETA - 1) * Complex(
G2,
D2) * pow3(
x));
650 (-2 * sqrt_pi * (
w1 * dz1 +
z1 * 1i *
dw1 * dz1) * pow3(
x) +
651 (4 *
sqrty * dsqrty + dx) * (2 *
x - 3) *
x / 2 -
652 (
x - 3) * (2 * pow2(
sqrty) +
x - 1) * dx) /
653 ((
ETA - 1) * Complex(
G2,
D2) * pow3(
x));
654 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
656 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
663 const Numeric dmF0 = (1 -
ETA) * dD0;
664 const Numeric dGD = (dmF0 / (invGD * mF0));
665 const Numeric dinvGD = -dGD * pow2(invGD);
666 const Complex dsqrty =
667 dinvGD / (2 * (
ETA - 1) * Complex(
G2,
D2) * pow2(invGD));
668 const Complex ddeltax = Complex(0, 1 -
ETA) * dD0;
669 const Complex dx = -ddeltax / ((
ETA - 1) * Complex(
G2,
D2));
670 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
673 case CalcType::Full: {
674 const Complex dz1 = dsqrtxy - dsqrty;
675 const Complex dz2 = dsqrtxy + dsqrty;
678 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
681 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
682 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
683 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) *
685 (2 * (
ETA - 1) * Complex(
G2,
D2) * pow2(sqrty));
686 const Complex dK =
ETA * Complex(
G2,
D2) * dB + Complex(0,
ETA * dD0) * A +
688 return inv_pi * (-A * dK + K * dA) / pow2(K);
690 case CalcType::Noc2tLowZ: {
691 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
692 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
694 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
697 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
698 const Complex dK =
ETA * Complex(
G2,
D2) * dB + Complex(0,
ETA * dD0) * A +
700 return inv_pi * (-A * dK + K * dA) / pow2(K);
702 case CalcType::Noc2tHighZ: {
703 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
704 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
706 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
707 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
711 const Complex dK =
ETA * Complex(
G2,
D2) * dB + Complex(0,
ETA * dD0) * A +
713 return inv_pi * (-A * dK + K * dA) / pow2(K);
715 case CalcType::LowXandHighY: {
716 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
717 const Complex dz2 = dsqrtxy + dsqrty;
720 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
722 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
723 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
727 const Complex dK =
ETA * Complex(
G2,
D2) * dB + Complex(0,
ETA * dD0) * A +
729 return inv_pi * (-A * dK + K * dA) / pow2(K);
731 case CalcType::LowYandLowX: {
732 const Complex dz1 = dsqrtxy;
733 const Complex dz2 = dx / (2 * sqrtx);
734 const Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
737 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
738 (2 * pow2(sqrty) + x - 1) +
739 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
740 2 * (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) /
742 const Complex dK =
ETA * Complex(
G2,
D2) * dB + Complex(0,
ETA * dD0) * A +
744 return inv_pi * (-A * dK + K * dA) / pow2(K);
746 case CalcType::LowYandHighX: {
747 const Complex dz1 = dsqrtxy;
748 const Complex dA = (x - 3) * dx / ((
ETA - 1) * Complex(
G2,
D2) * pow3(x));
750 (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
751 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x / 2 -
752 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) /
753 ((
ETA - 1) * Complex(
G2,
D2) * pow3(x));
754 const Complex dK =
ETA * Complex(
G2,
D2) * dB + Complex(0,
ETA * dD0) * A +
756 return inv_pi * (-A * dK + K * dA) / pow2(K);
763 const Numeric ddeltax = (1 -
ETA) * dG0;
764 const Complex dx = -ddeltax / ((
ETA - 1) * Complex(
G2,
D2));
765 const Complex dsqrtxy = dx / (2 * sqrtxy);
768 case CalcType::Full: {
769 const Complex dz1 = dsqrtxy;
770 const Complex dz2 = dsqrtxy;
771 const Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
774 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
775 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) /
776 (2 * sqrty * (
ETA - 1) * Complex(
G2,
D2));
777 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
ETA * A * dG0 +
779 return inv_pi * (-A * dK + K * dA) / pow2(K);
781 case CalcType::Noc2tLowZ: {
782 const Complex dz1 = invGD * ddeltax;
783 const Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
786 (sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) - dz1);
787 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
ETA * A * dG0 +
789 return inv_pi * (-A * dK + K * dA) / pow2(K);
791 case CalcType::Noc2tHighZ: {
792 const Complex dz1 = invGD * ddeltax;
793 const Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
794 const Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
795 invGD * dz1 / (2 * pow2(z1)) +
796 9 * invGD * dz1 / (4 * pow4(z1));
797 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
ETA * A * dG0 +
799 return inv_pi * (-A * dK + K * dA) / pow2(K);
801 case CalcType::LowXandHighY: {
802 const Complex dz1 = invGD * ddeltax;
803 const Complex dz2 = dsqrtxy;
804 const Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
805 const Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
806 invGD * dz1 / (2 * pow2(z1)) +
807 9 * invGD * dz1 / (4 * pow4(z1));
808 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
ETA * A * dG0 +
810 return inv_pi * (-A * dK + K * dA) / pow2(K);
812 case CalcType::LowYandLowX: {
813 const Complex dz1 = dsqrtxy;
814 const Complex dz2 = dx / (2 * sqrtx);
815 const Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
818 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
819 (2 * pow2(sqrty) + x - 1) +
820 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
821 2 * (sqrt_pi * w2 * z2 - 1) * dx) /
823 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
ETA * A * dG0 +
825 return inv_pi * (-A * dK + K * dA) / pow2(K);
827 case CalcType::LowYandHighX: {
828 const Complex dz1 = dsqrtxy;
829 const Complex dA = (x - 3) * dx / ((
ETA - 1) * Complex(
G2,
D2) * pow3(x));
831 (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) -
832 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx + (2 * x - 3) * x * dx / 2) /
833 ((
ETA - 1) * Complex(
G2,
D2) * pow3(x));
834 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
ETA * A * dG0 +
836 return inv_pi * (-A * dK + K * dA) / pow2(K);
843 const Numeric dmF0 = -3 * (1 -
ETA) * dD2 / 2;
844 const Numeric dGD = (dmF0 / (invGD * mF0));
845 const Numeric dinvGD = -dGD * pow2(invGD);
846 const Complex dsqrty = (Complex(
G2,
D2) * dinvGD + Complex(0, invGD) * dD2) /
847 (2 * (
ETA - 1) * pow2(Complex(
G2,
D2)) * pow2(invGD));
848 const Complex ddeltax = 1.5 * Complex(0,
ETA - 1) * dD2;
849 const Complex dx = (-Complex(
G2,
D2) * ddeltax + Complex(0, dD2) * deltax) /
850 ((
ETA - 1) * pow2(Complex(
G2,
D2)));
851 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
854 case CalcType::Full: {
855 const Complex dz1 = dsqrtxy - dsqrty;
856 const Complex dz2 = dsqrtxy + dsqrty;
859 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
861 (sqrt_pi * Complex(
G2,
D2) *
862 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
863 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
864 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
868 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
871 (2 * (
ETA - 1) * pow2(Complex(
G2,
D2)) * pow2(sqrty));
872 const Complex dK =
ETA * Complex(
G2,
D2) * dB -
873 Complex(0, 1.5 * dD2 *
ETA) * A +
874 Complex(0, dD2 *
ETA) * B +
876 return inv_pi * (-A * dK + K * dA) / pow2(K);
878 case CalcType::Noc2tLowZ: {
879 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
880 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
882 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
885 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
886 const Complex dK =
ETA * Complex(
G2,
D2) * dB -
887 Complex(0, 1.5 * dD2 *
ETA) * A +
888 Complex(0, dD2 *
ETA) * B +
890 return inv_pi * (-A * dK + K * dA) / pow2(K);
892 case CalcType::Noc2tHighZ: {
893 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
894 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
896 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
897 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
901 const Complex dK =
ETA * Complex(
G2,
D2) * dB -
902 Complex(0, 1.5 * dD2 *
ETA) * A +
903 Complex(0, dD2 *
ETA) * B +
905 return inv_pi * (-A * dK + K * dA) / pow2(K);
907 case CalcType::LowXandHighY: {
908 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
909 const Complex dz2 = dsqrtxy + dsqrty;
912 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
914 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
915 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
919 const Complex dK =
ETA * Complex(
G2,
D2) * dB -
920 Complex(0, 1.5 * dD2 *
ETA) * A +
921 Complex(0, dD2 *
ETA) * B +
923 return inv_pi * (-A * dK + K * dA) / pow2(K);
925 case CalcType::LowYandLowX: {
926 const Complex dz1 = dsqrtxy;
927 const Complex dz2 = dx / (2 * sqrtx);
930 (sqrt_pi * Complex(
G2,
D2) * (w2 * dz2 + z2 * 1i * dw2 * dz2) -
931 1i * (sqrt_pi * w2 * z2 - 1) * dD2) /
932 ((
ETA - 1) * pow2(Complex(
G2,
D2)));
934 (-2 * Complex(
G2,
D2) *
935 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
936 (2 * pow2(sqrty) + x - 1) +
937 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
938 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
940 (2 * sqrt_pi * w1 * z1 +
941 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
943 ((
ETA - 1) * pow2(Complex(
G2,
D2)));
944 const Complex dK =
ETA * Complex(
G2,
D2) * dB -
945 Complex(0, 1.5 * dD2 *
ETA) * A +
946 Complex(0, dD2 *
ETA) * B +
948 return inv_pi * (-A * dK + K * dA) / pow2(K);
950 case CalcType::LowYandHighX: {
951 const Complex dz1 = dsqrtxy;
953 (Complex(
G2,
D2) * (x - 3) * dx + 1i * (2 * x - 3) * x * dD2 / 2) /
954 ((
ETA - 1) * pow2(Complex(
G2,
D2)) * pow3(x));
957 (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
958 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
959 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) -
961 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
962 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
964 (2 * (
ETA - 1) * pow2(Complex(
G2,
D2)) * pow3(x));
965 const Complex dK =
ETA * Complex(
G2,
D2) * dB -
966 Complex(0, 1.5 * dD2 *
ETA) * A +
967 Complex(0, dD2 *
ETA) * B +
969 return inv_pi * (-A * dK + K * dA) / pow2(K);
976 const Complex dsqrty = dG2 / (2 * invGD * (
ETA - 1) * pow2(Complex(
G2,
D2)));
977 const Numeric ddeltax = 3 * (
ETA - 1) * dG2 / 2;
978 const Complex dx = (-Complex(
G2,
D2) * ddeltax + deltax * dG2) /
979 ((
ETA - 1) * pow2(Complex(
G2,
D2)));
980 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
983 case CalcType::Full: {
984 const Complex dz1 = dsqrtxy - dsqrty;
985 const Complex dz2 = dsqrtxy + dsqrty;
986 const Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
988 (sqrt_pi * Complex(
G2,
D2) *
989 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
990 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
991 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
994 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
997 (2 * (
ETA - 1) * pow2(Complex(
G2,
D2)) * pow2(sqrty));
998 const Complex dK =
ETA * Complex(
G2,
D2) * dB - 3 *
ETA * A * dG2 / 2 +
1001 return inv_pi * (-A * dK + K * dA) / pow2(K);
1003 case CalcType::Noc2tLowZ: {
1004 const Complex dz1 = invGD * ddeltax;
1005 const Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
1008 (sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) - dz1);
1009 const Complex dK =
ETA * Complex(
G2,
D2) * dB - 3 *
ETA * A * dG2 / 2 +
1012 return inv_pi * (-A * dK + K * dA) / pow2(K);
1014 case CalcType::Noc2tHighZ: {
1015 const Complex dz1 = invGD * ddeltax;
1016 const Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
1017 const Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
1018 invGD * dz1 / (2 * pow2(z1)) +
1019 9 * invGD * dz1 / (4 * pow4(z1));
1020 const Complex dK =
ETA * Complex(
G2,
D2) * dB - 3 *
ETA * A * dG2 / 2 +
1023 return inv_pi * (-A * dK + K * dA) / pow2(K);
1025 case CalcType::LowXandHighY: {
1026 const Complex dz1 = invGD * ddeltax;
1027 const Complex dz2 = dsqrtxy + dsqrty;
1028 const Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
1029 const Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
1030 invGD * dz1 / (2 * pow2(z1)) +
1031 9 * invGD * dz1 / (4 * pow4(z1));
1032 const Complex dK =
ETA * Complex(
G2,
D2) * dB - 3 *
ETA * A * dG2 / 2 +
1035 return inv_pi * (-A * dK + K * dA) / pow2(K);
1037 case CalcType::LowYandLowX: {
1038 const Complex dz1 = dsqrtxy;
1039 const Complex dz2 = dx / (2 * sqrtx);
1042 (sqrt_pi * Complex(
G2,
D2) * (w2 * dz2 + z2 * 1i * dw2 * dz2) -
1043 (sqrt_pi * w2 * z2 - 1) * dG2) /
1044 ((
ETA - 1) * pow2(Complex(
G2,
D2)));
1046 (-2 * Complex(
G2,
D2) *
1047 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1048 (2 * pow2(sqrty) + x - 1) +
1049 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
1050 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
1051 (2 * sqrt_pi * w1 * z1 +
1052 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
1054 ((
ETA - 1) * pow2(Complex(
G2,
D2)));
1055 const Complex dK =
ETA * Complex(
G2,
D2) * dB - 3 *
ETA * A * dG2 / 2 +
1058 return inv_pi * (-A * dK + K * dA) / pow2(K);
1060 case CalcType::LowYandHighX: {
1061 const Complex dz1 = dsqrtxy;
1063 (Complex(
G2,
D2) * (x - 3) * dx + (2 * x - 3) * x * dG2 / 2) /
1064 ((
ETA - 1) * pow2(Complex(
G2,
D2)) * pow3(x));
1067 (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
1068 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
1069 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) -
1070 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1071 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1073 (2 * (
ETA - 1) * pow2(Complex(
G2,
D2)) * pow3(x));
1074 const Complex dK =
ETA * Complex(
G2,
D2) * dB - 3 *
ETA * A * dG2 / 2 +
1077 return inv_pi * (-A * dK + K * dA) / pow2(K);
1084 const Numeric ddeltax = dFVC;
1085 const Complex dx = -ddeltax / ((
ETA - 1) * Complex(
G2,
D2));
1086 const Complex dsqrtxy = dx / (2 * sqrtxy);
1089 case CalcType::Full: {
1090 const Complex dz1 = dsqrtxy;
1091 const Complex dz2 = dsqrtxy;
1092 const Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
1095 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
1096 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) /
1097 (2 * sqrty * (
ETA - 1) * Complex(
G2,
D2));
1100 (
ETA * Complex(
G0 - 1.5 *
G2,
D0 - 1.5 *
D2) -
FVC) * dA - A * dFVC;
1101 return inv_pi * (-A * dK + K * dA) / pow2(K);
1103 case CalcType::Noc2tLowZ: {
1104 const Complex dz1 = invGD * ddeltax;
1105 const Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
1108 (sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) - dz1);
1111 (
ETA * Complex(
G0 - 1.5 *
G2,
D0 - 1.5 *
D2) -
FVC) * dA - A * dFVC;
1112 return inv_pi * (-A * dK + K * dA) / pow2(K);
1114 case CalcType::Noc2tHighZ: {
1115 const Complex dz1 = invGD * ddeltax;
1116 const Complex dA = Complex(0, sqrt_pi * invGD) * dw1 * dz1;
1117 const Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
1118 invGD * dz1 / (2 * pow2(z1)) +
1119 9 * invGD * dz1 / (4 * pow4(z1));
1122 (
ETA * Complex(
G0 - 1.5 *
G2,
D0 - 1.5 *
D2) -
FVC) * dA - A * dFVC;
1123 return inv_pi * (-A * dK + K * dA) / pow2(K);
1125 case CalcType::LowXandHighY: {
1126 const Complex dz1 = invGD * ddeltax;
1127 const Complex dz2 = dsqrtxy;
1128 const Complex dA = Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
1129 const Complex dB = Complex(0, sqrt_pi * invGD) * dw1 * dz1 -
1130 invGD * dz1 / (2 * pow2(z1)) +
1131 9 * invGD * dz1 / (4 * pow4(z1));
1134 (
ETA * Complex(
G0 - 1.5 *
G2,
D0 - 1.5 *
D2) -
FVC) * dA - A * dFVC;
1135 return inv_pi * (-A * dK + K * dA) / pow2(K);
1137 case CalcType::LowYandLowX: {
1138 const Complex dz1 = dsqrtxy;
1139 const Complex dz2 = dx / (2 * sqrtx);
1140 const Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
1141 ((
ETA - 1) * Complex(
G2,
D2));
1143 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1144 (2 * pow2(sqrty) + x - 1) +
1145 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
1146 2 * (sqrt_pi * w2 * z2 - 1) * dx) /
1147 ((
ETA - 1) * Complex(
G2,
D2));
1150 (
ETA * Complex(
G0 - 1.5 *
G2,
D0 - 1.5 *
D2) -
FVC) * dA - A * dFVC;
1151 return inv_pi * (-A * dK + K * dA) / pow2(K);
1153 case CalcType::LowYandHighX: {
1154 const Complex dz1 = dsqrtxy;
1155 const Complex dA = (x - 3) * dx / ((
ETA - 1) * Complex(
G2,
D2) * pow3(x));
1157 (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) -
1158 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx + (2 * x - 3) * x * dx / 2) /
1159 ((
ETA - 1) * Complex(
G2,
D2) * pow3(x));
1162 (
ETA * Complex(
G0 - 1.5 *
G2,
D0 - 1.5 *
D2) -
FVC) * dA - A * dFVC;
1163 return inv_pi * (-A * dK + K * dA) / pow2(K);
1170 const Numeric dmF0 = -(
D0 - 3 *
D2 / 2) * dETA;
1171 const Numeric dGD = (dmF0 / (invGD * mF0));
1172 const Numeric dinvGD = -dGD * pow2(invGD);
1173 const Complex dsqrty = ((
ETA - 1) * dinvGD + invGD * dETA) /
1174 (2 * Complex(
G2,
D2) * pow2(
ETA - 1) * pow2(invGD));
1175 const Complex ddeltax = -dETA * Complex(
G0 - 1.5 *
G2,
D0 - 1.5 *
D2);
1176 const Complex dx = (-(
ETA - 1) * ddeltax + deltax * dETA) /
1177 (Complex(
G2,
D2) * pow2(
ETA - 1));
1178 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
1181 case CalcType::Full: {
1182 const Complex dz1 = dsqrtxy - dsqrty;
1183 const Complex dz2 = dsqrtxy + dsqrty;
1186 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1189 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
1190 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
1191 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
1192 2 * w2 * z2 * dz2) *
1195 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
1198 (2 * Complex(
G2,
D2) * pow2(
ETA - 1) * pow2(sqrty));
1201 Complex(
G2,
D2) * B * dETA + Complex(
G2,
D2) *
ETA * dB -
1202 Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) * A * dETA;
1203 return inv_pi * (-A * dK + K * dA) / pow2(K);
1205 case CalcType::Noc2tLowZ: {
1206 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1207 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1209 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
1212 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
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);
1219 case CalcType::Noc2tHighZ: {
1220 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1221 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1223 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1224 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1230 Complex(
G2,
D2) * B * dETA + Complex(
G2,
D2) *
ETA * dB -
1231 Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) * A * dETA;
1232 return inv_pi * (-A * dK + K * dA) / pow2(K);
1234 case CalcType::LowXandHighY: {
1235 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1236 const Complex dz2 = dsqrtxy + dsqrty;
1239 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1241 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1242 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1248 Complex(
G2,
D2) * B * dETA + Complex(
G2,
D2) *
ETA * dB -
1249 Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) * A * dETA;
1250 return inv_pi * (-A * dK + K * dA) / pow2(K);
1252 case CalcType::LowYandLowX: {
1253 const Complex dz1 = dsqrtxy;
1254 const Complex dz2 = dx / (2 * sqrtx);
1255 const Complex dA = 2 *
1256 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) * (
ETA - 1) -
1257 (sqrt_pi * w2 * z2 - 1) * dETA) /
1258 (Complex(
G2,
D2) * pow2(
ETA - 1));
1261 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1262 (2 * pow2(sqrty) + x - 1) +
1263 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
1264 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
1265 (2 * sqrt_pi * w1 * z1 +
1266 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
1268 (Complex(
G2,
D2) * pow2(
ETA - 1));
1271 Complex(
G2,
D2) * B * dETA + Complex(
G2,
D2) *
ETA * dB -
1272 Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) * A * dETA;
1273 return inv_pi * (-A * dK + K * dA) / pow2(K);
1275 case CalcType::LowYandHighX: {
1276 const Complex dz1 = dsqrtxy;
1277 const Complex dA = ((
ETA - 1) * (x - 3) * dx + (2 * x - 3) * x * dETA / 2) /
1278 (Complex(
G2,
D2) * pow2(
ETA - 1) * pow3(x));
1280 (-(2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1281 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1284 (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
1285 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
1286 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx)) /
1287 (2 * Complex(
G2,
D2) * pow2(
ETA - 1) * pow3(x));
1290 Complex(
G2,
D2) * B * dETA + Complex(
G2,
D2) *
ETA * dB -
1291 Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) * A * dETA;
1292 return inv_pi * (-A * dK + K * dA) / pow2(K);
1299 const Numeric dmF0 = dZ;
1300 const Numeric dGD = (dmF0 / (invGD * mF0));
1301 const Numeric dinvGD = -dGD * pow2(invGD);
1302 const Complex dsqrty =
1303 dinvGD / (2 * (
ETA - 1) * Complex(
G2,
D2) * pow2(invGD));
1304 const Complex ddeltax = Complex(0, dZ);
1305 const Complex dx = -ddeltax / ((
ETA - 1) * Complex(
G2,
D2));
1306 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
1309 case CalcType::Full: {
1310 const Complex dz1 = dsqrtxy - dsqrty;
1311 const Complex dz2 = dsqrtxy + dsqrty;
1314 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1317 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
1318 ((pow2(z1) - 1) * 1i * dw1 * dz1 - (pow2(z2) - 1) * 1i * dw2 * dz2 +
1319 2 * w1 * z1 * dz1 - 2 * w2 * z2 * dz2) *
1321 (2 * (
ETA - 1) * Complex(
G2,
D2) * pow2(sqrty));
1322 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
1324 return inv_pi * (-A * dK + K * dA) / pow2(K);
1326 case CalcType::Noc2tLowZ: {
1327 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1328 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1330 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
1333 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
1334 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
1336 return inv_pi * (-A * dK + K * dA) / pow2(K);
1338 case CalcType::Noc2tHighZ: {
1339 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1340 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1342 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1343 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1347 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
1349 return inv_pi * (-A * dK + K * dA) / pow2(K);
1351 case CalcType::LowXandHighY: {
1352 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1353 const Complex dz2 = dsqrtxy + dsqrty;
1356 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1358 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1359 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1363 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
1365 return inv_pi * (-A * dK + K * dA) / pow2(K);
1367 case CalcType::LowYandLowX: {
1368 const Complex dz1 = dsqrtxy;
1369 const Complex dz2 = dx / (2 * sqrtx);
1370 const Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
1371 ((
ETA - 1) * Complex(
G2,
D2));
1373 -(2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1374 (2 * pow2(sqrty) + x - 1) +
1375 2 * sqrt_pi * w1 * dz1 + Complex(0, 2 * sqrt_pi) * z1 * dw1 * dz1 +
1376 2 * (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) /
1377 ((
ETA - 1) * Complex(
G2,
D2));
1378 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
1380 return inv_pi * (-A * dK + K * dA) / pow2(K);
1382 case CalcType::LowYandHighX: {
1383 const Complex dz1 = dsqrtxy;
1384 const Complex dA = (x - 3) * dx / ((
ETA - 1) * Complex(
G2,
D2) * pow3(x));
1386 (-2 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
1387 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x / 2 -
1388 (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) /
1389 ((
ETA - 1) * Complex(
G2,
D2) * pow3(x));
1390 const Complex dK =
ETA * Complex(
G2,
D2) * dB +
1392 return inv_pi * (-A * dK + K * dA) / pow2(K);
1399 const Numeric dmF0 =
1400 (1 -
ETA) * (
d.D0 - 3 *
d.D2 / 2) - (
D0 - 3 *
D2 / 2) *
d.ETA;
1401 const Numeric dGD = (dmF0 / (invGD * mF0));
1402 const Numeric dinvGD = -dGD * pow2(invGD);
1403 const Complex dsqrty =
1404 (Complex(
G2,
D2) * (
ETA - 1) * dinvGD + Complex(
G2,
D2) * invGD *
d.ETA +
1405 Complex(
d.G2,
d.D2) * (
ETA - 1) * invGD) /
1406 (2 * pow2(Complex(
G2,
D2)) * pow2(
ETA - 1) * pow2(invGD));
1407 const Complex ddeltax =
1408 -(
ETA - 1) * Complex(
d.G0 - 1.5 *
d.G2,
d.D0 - 1.5 *
d.D2) -
1409 Complex(
G0 - 1.5 *
G2,
D0 - 1.5 *
D2) *
d.ETA +
d.FVC;
1410 const Complex dx = (-Complex(
G2,
D2) * (
ETA - 1) * ddeltax +
1411 Complex(
G2,
D2) * deltax *
d.ETA +
1412 Complex(
d.G2,
d.D2) * (
ETA - 1) * deltax) /
1413 (pow2(Complex(
G2,
D2)) * pow2(
ETA - 1));
1414 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
1417 case CalcType::Full: {
1418 const Complex dz1 = dsqrtxy - dsqrty;
1419 const Complex dz2 = dsqrtxy + dsqrty;
1422 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1424 (sqrt_pi * Complex(
G2,
D2) *
1425 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
1426 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
1427 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
1428 2 * w2 * z2 * dz2) *
1432 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
1435 Complex(
d.G2,
d.D2) * (
ETA - 1) *
1436 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
1439 (2 * pow2(Complex(
G2,
D2)) * pow2(
ETA - 1) * pow2(sqrty));
1441 Complex(
G2,
D2) * B *
d.ETA + Complex(
G2,
D2) *
ETA * dB +
1442 Complex(
d.G2,
d.D2) * B *
ETA +
1444 (-Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) *
d.ETA -
1445 Complex(1.5 *
d.G2 -
d.G0, 1.5 *
d.D2 -
d.D0) *
ETA -
d.FVC) *
1447 return inv_pi * (-A * dK + K * dA) / pow2(K);
1449 case CalcType::Noc2tLowZ: {
1450 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1451 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1453 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
1456 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
1458 Complex(
G2,
D2) * B *
d.ETA + Complex(
G2,
D2) *
ETA * dB +
1459 Complex(
d.G2,
d.D2) * B *
ETA +
1461 (-Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) *
d.ETA -
1462 Complex(1.5 *
d.G2 -
d.G0, 1.5 *
d.D2 -
d.D0) *
ETA -
d.FVC) *
1464 return inv_pi * (-A * dK + K * dA) / pow2(K);
1466 case CalcType::Noc2tHighZ: {
1467 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1468 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1470 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1471 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1476 Complex(
G2,
D2) * B *
d.ETA + Complex(
G2,
D2) *
ETA * dB +
1477 Complex(
d.G2,
d.D2) * B *
ETA +
1479 (-Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) *
d.ETA -
1480 Complex(1.5 *
d.G2 -
d.G0, 1.5 *
d.D2 -
d.D0) *
ETA -
d.FVC) *
1482 return inv_pi * (-A * dK + K * dA) / pow2(K);
1484 case CalcType::LowXandHighY: {
1485 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1486 const Complex dz2 = dsqrtxy + dsqrty;
1489 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1491 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1492 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1497 Complex(
G2,
D2) * B *
d.ETA + Complex(
G2,
D2) *
ETA * dB +
1498 Complex(
d.G2,
d.D2) * B *
ETA +
1500 (-Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) *
d.ETA -
1501 Complex(1.5 *
d.G2 -
d.G0, 1.5 *
d.D2 -
d.D0) *
ETA -
d.FVC) *
1503 return inv_pi * (-A * dK + K * dA) / pow2(K);
1505 case CalcType::LowYandLowX: {
1506 const Complex dz1 = dsqrtxy;
1507 const Complex dz2 = dx / (2 * sqrtx);
1510 (sqrt_pi * Complex(
G2,
D2) * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1512 Complex(
G2,
D2) * (sqrt_pi * w2 * z2 - 1) *
d.ETA -
1513 Complex(
d.G2,
d.D2) * (sqrt_pi * w2 * z2 - 1) * (
ETA - 1)) /
1514 (pow2(Complex(
G2,
D2)) * pow2(
ETA - 1));
1516 (-2 * Complex(
G2,
D2) * (
ETA - 1) *
1517 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1518 (2 * pow2(sqrty) + x - 1) +
1519 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
1520 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
1522 (2 * sqrt_pi * w1 * z1 +
1523 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
1525 Complex(
d.G2,
d.D2) * (
ETA - 1) *
1526 (2 * sqrt_pi * w1 * z1 +
1527 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1)) /
1528 (pow2(Complex(
G2,
D2)) * pow2(
ETA - 1));
1530 Complex(
G2,
D2) * B *
d.ETA + Complex(
G2,
D2) *
ETA * dB +
1531 Complex(
d.G2,
d.D2) * B *
ETA +
1533 (-Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) *
d.ETA -
1534 Complex(1.5 *
d.G2 -
d.G0, 1.5 *
d.D2 -
d.D0) *
ETA -
d.FVC) *
1536 return inv_pi * (-A * dK + K * dA) / pow2(K);
1538 case CalcType::LowYandHighX: {
1539 const Complex dz1 = dsqrtxy;
1540 const Complex dA = (2 * Complex(
G2,
D2) * (
ETA - 1) * (x - 3) * dx +
1541 Complex(
G2,
D2) * (2 * x - 3) * x *
d.ETA +
1542 Complex(
d.G2,
d.D2) * (
ETA - 1) * (2 * x - 3) * x) /
1543 (2 * pow2(Complex(
G2,
D2)) * pow2(
ETA - 1) * pow3(x));
1546 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1547 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1550 (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
1551 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
1552 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) -
1553 Complex(
d.G2,
d.D2) *
1554 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1555 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1557 (2 * pow2(Complex(
G2,
D2)) * pow2(
ETA - 1) * pow3(x));
1559 Complex(
G2,
D2) * B *
d.ETA + Complex(
G2,
D2) *
ETA * dB +
1560 Complex(
d.G2,
d.D2) * B *
ETA +
1562 (-Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) *
d.ETA -
1563 Complex(1.5 *
d.G2 -
d.G0, 1.5 *
d.D2 -
d.D0) *
ETA -
d.FVC) *
1565 return inv_pi * (-A * dK + K * dA) / pow2(K);
1572 const Numeric dmF0 =
1573 (1 -
ETA) * (
d.D0 - 3 *
d.D2 / 2) - (
D0 - 3 *
D2 / 2) *
d.ETA;
1575 (dmF0 / (invGD * mF0)) - invGD * invGD / (2 * T * sqrt_ln_2);
1576 const Numeric dinvGD = -dGD * pow2(invGD);
1577 const Complex dsqrty =
1578 (Complex(
G2,
D2) * (
ETA - 1) * dinvGD + Complex(
G2,
D2) * invGD *
d.ETA +
1579 Complex(
d.G2,
d.D2) * (
ETA - 1) * invGD) /
1580 (2 * pow2(Complex(
G2,
D2)) * pow2(
ETA - 1) * pow2(invGD));
1581 const Complex ddeltax =
1582 -(
ETA - 1) * Complex(
d.G0 - 1.5 *
d.G2,
d.D0 - 1.5 *
d.D2) -
1583 Complex(
G0 - 1.5 *
G2,
D0 - 1.5 *
D2) *
d.ETA +
d.FVC;
1584 const Complex dx = (-Complex(
G2,
D2) * (
ETA - 1) * ddeltax +
1585 Complex(
G2,
D2) * deltax *
d.ETA +
1586 Complex(
d.G2,
d.D2) * (
ETA - 1) * deltax) /
1587 (pow2(Complex(
G2,
D2)) * pow2(
ETA - 1));
1588 const Complex dsqrtxy = (sqrty * dsqrty + dx / 2) / sqrtxy;
1591 case CalcType::Full: {
1592 const Complex dz1 = dsqrtxy - dsqrty;
1593 const Complex dz2 = dsqrtxy + dsqrty;
1596 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1598 (sqrt_pi * Complex(
G2,
D2) *
1599 ((-(pow2(z1) - 1) * w1 + (pow2(z2) - 1) * w2) * dsqrty +
1600 ((pow2(z1) - 1) * 1i * dw1 * dz1 -
1601 (pow2(z2) - 1) * 1i * dw2 * dz2 + 2 * w1 * z1 * dz1 -
1602 2 * w2 * z2 * dz2) *
1606 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
1609 Complex(
d.G2,
d.D2) * (
ETA - 1) *
1610 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
1613 (2 * pow2(Complex(
G2,
D2)) * pow2(
ETA - 1) * pow2(sqrty));
1615 Complex(
G2,
D2) * B *
d.ETA + Complex(
G2,
D2) *
ETA * dB +
1616 Complex(
d.G2,
d.D2) * B *
ETA +
1618 (-Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) *
d.ETA -
1619 Complex(1.5 *
d.G2 -
d.G0, 1.5 *
d.D2 -
d.D0) *
ETA -
d.FVC) *
1621 return inv_pi * (-A * dK + K * dA) / pow2(K);
1623 case CalcType::Noc2tLowZ: {
1624 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1625 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1627 -(sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) -
1630 (sqrt_pi * (pow2(z1) - 1) * w1 - z1) * dinvGD;
1632 Complex(
G2,
D2) * B *
d.ETA + Complex(
G2,
D2) *
ETA * dB +
1633 Complex(
d.G2,
d.D2) * B *
ETA +
1635 (-Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) *
d.ETA -
1636 Complex(1.5 *
d.G2 -
d.G0, 1.5 *
d.D2 -
d.D0) *
ETA -
d.FVC) *
1638 return inv_pi * (-A * dK + K * dA) / pow2(K);
1640 case CalcType::Noc2tHighZ: {
1641 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1642 const Complex dA = sqrt_pi * (Complex(0, invGD) * dw1 * dz1 + w1 * dinvGD);
1644 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1645 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1650 Complex(
G2,
D2) * B *
d.ETA + Complex(
G2,
D2) *
ETA * dB +
1651 Complex(
d.G2,
d.D2) * B *
ETA +
1653 (-Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) *
d.ETA -
1654 Complex(1.5 *
d.G2 -
d.G0, 1.5 *
d.D2 -
d.D0) *
ETA -
d.FVC) *
1656 return inv_pi * (-A * dK + K * dA) / pow2(K);
1658 case CalcType::LowXandHighY: {
1659 const Complex dz1 = deltax * dinvGD + invGD * ddeltax;
1660 const Complex dz2 = dsqrtxy + dsqrty;
1663 ((w1 - w2) * dinvGD + (Complex(0, invGD) * (dw1 * dz1 - dw2 * dz2)));
1665 ((4 * sqrt_pi * w1 * pow3(z1) + 2 * pow2(z1) - 3) * z1 * dinvGD +
1666 (Complex(0, 4 * sqrt_pi) * pow4(z1) * dw1 * dz1 - 2 * pow2(z1) * dz1 +
1671 Complex(
G2,
D2) * B *
d.ETA + Complex(
G2,
D2) *
ETA * dB +
1672 Complex(
d.G2,
d.D2) * B *
ETA +
1674 (-Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) *
d.ETA -
1675 Complex(1.5 *
d.G2 -
d.G0, 1.5 *
d.D2 -
d.D0) *
ETA -
d.FVC) *
1677 return inv_pi * (-A * dK + K * dA) / pow2(K);
1679 case CalcType::LowYandLowX: {
1680 const Complex dz1 = dsqrtxy;
1681 const Complex dz2 = dx / (2 * sqrtx);
1684 (sqrt_pi * Complex(
G2,
D2) * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1686 Complex(
G2,
D2) * (sqrt_pi * w2 * z2 - 1) *
d.ETA -
1687 Complex(
d.G2,
d.D2) * (sqrt_pi * w2 * z2 - 1) * (
ETA - 1)) /
1688 (pow2(Complex(
G2,
D2)) * pow2(
ETA - 1));
1690 (-2 * Complex(
G2,
D2) * (
ETA - 1) *
1691 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1692 (2 * pow2(sqrty) + x - 1) +
1693 sqrt_pi * w1 * dz1 + Complex(0, sqrt_pi) * z1 * dw1 * dz1 +
1694 (4 * sqrty * dsqrty + dx) * (sqrt_pi * w2 * z2 - 1)) +
1696 (2 * sqrt_pi * w1 * z1 +
1697 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1) *
1699 Complex(
d.G2,
d.D2) * (
ETA - 1) *
1700 (2 * sqrt_pi * w1 * z1 +
1701 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1)) /
1702 (pow2(Complex(
G2,
D2)) * pow2(
ETA - 1));
1704 Complex(
G2,
D2) * B *
d.ETA + Complex(
G2,
D2) *
ETA * dB +
1705 Complex(
d.G2,
d.D2) * B *
ETA +
1707 (-Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) *
d.ETA -
1708 Complex(1.5 *
d.G2 -
d.G0, 1.5 *
d.D2 -
d.D0) *
ETA -
d.FVC) *
1710 return inv_pi * (-A * dK + K * dA) / pow2(K);
1712 case CalcType::LowYandHighX: {
1713 const Complex dz1 = dsqrtxy;
1714 const Complex dA = (2 * Complex(
G2,
D2) * (
ETA - 1) * (x - 3) * dx +
1715 Complex(
G2,
D2) * (2 * x - 3) * x *
d.ETA +
1716 Complex(
d.G2,
d.D2) * (
ETA - 1) * (2 * x - 3) * x) /
1717 (2 * pow2(Complex(
G2,
D2)) * pow2(
ETA - 1) * pow3(x));
1720 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1721 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1724 (-4 * sqrt_pi * (w1 * dz1 + z1 * 1i * dw1 * dz1) * pow3(x) +
1725 (4 * sqrty * dsqrty + dx) * (2 * x - 3) * x -
1726 2 * (x - 3) * (2 * pow2(sqrty) + x - 1) * dx) -
1727 Complex(
d.G2,
d.D2) *
1728 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1729 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1731 (2 * pow2(Complex(
G2,
D2)) * pow2(
ETA - 1) * pow3(x));
1733 Complex(
G2,
D2) * B *
d.ETA + Complex(
G2,
D2) *
ETA * dB +
1734 Complex(
d.G2,
d.D2) * B *
ETA +
1736 (-Complex(1.5 *
G2 -
G0, 1.5 *
D2 -
D0) *
d.ETA -
1737 Complex(1.5 *
d.G2 -
d.G0, 1.5 *
d.D2 -
d.D0) *
ETA -
d.FVC) *
1739 return inv_pi * (-A * dK + K * dA) / pow2(K);
1746 imag_val(deltax) = mF0 - f;
1747 x = deltax / ((1 -
ETA) * Complex(
G2,
D2));
1748 sqrtxy = std::sqrt(x + sqrty * sqrty);
1756 return CalcType::Noc2tHighZ;
1758 return CalcType::LowXandHighY;
1761 return CalcType::LowYandLowX;
1763 return CalcType::LowYandHighX;
1764 return CalcType::Full;
1776 w1 = Faddeeva::w(1i *
z1);
1777 w2 = Faddeeva::w(1i *
z2);
1779 B = (-1 + sqrt_pi / (2 *
sqrty) * (1 - pow2(
z1)) *
w1 -
1780 sqrt_pi / (2 *
sqrty) * (1 - pow2(
z2)) *
w2) /
1781 ((1 -
ETA) * Complex(
G2,
D2));
1786 w1 = Faddeeva::w(1i *
z1);
1790 B = sqrt_pi *
invGD * ((1 - pow2(
z1)) *
w1 +
z1 / sqrt_pi);
1793 B =
invGD * (sqrt_pi *
w1 + 1 /
z1 / 2 - 3 / pow3(
z1) / 4);
1799 w1 = Faddeeva::w(1i *
z1);
1800 w2 = Faddeeva::w(1i *
z2);
1802 B =
invGD * (sqrt_pi *
w1 + 1 /
z1 / 2 - 3 / pow3(
z1) / 4);
1808 w1 = Faddeeva::w(1i *
z1);
1809 w2 = Faddeeva::w(1i *
z2);
1810 A = (2 * sqrt_pi / ((1 -
ETA) * Complex(
G2,
D2))) * (inv_sqrt_pi -
z2 *
w2);
1811 B = (1 / ((1 -
ETA) * Complex(
G2,
D2))) *
1814 2 * sqrt_pi *
z1 *
w1);
1818 w1 = Faddeeva::w(1i *
z1);
1819 A = (1 / ((1 -
ETA) * Complex(
G2,
D2))) * (1 /
x - 3 / pow2(
x) / 2);
1820 B = (1 / ((1 -
ETA) * Complex(
G2,
D2))) *
1821 (-1 + (1 -
x - 2 *
sqrty *
sqrty) * (1 /
x - 3 / pow2(
x) / 2) +
1822 2 * sqrt_pi *
z1 *
w1);
1826 dw1 = 2i * (inv_sqrt_pi -
z1 *
w1);
1827 dw2 = 2i * (inv_sqrt_pi -
z2 *
w2);
1835 inv_denom(1.0 / (F0 * tanh_c1f0)) {}
1838 const Numeric part1 =
1839 (1 - pow2(tanh_c1f0)) * c1 * f * tanh_c1f / (T * pow2(tanh_c1f0));
1840 const Numeric part2 = (pow2(tanh_c1f) - 1) * c1 * pow2(f) * inv_denom / T;
1841 return part1 + part2;
1845 const Numeric part1 = tanh_c1f * inv_denom;
1846 const Numeric part2 = c1 * f * inv_denom / (1 - pow2(tanh_c1f));
1847 return part1 + part2;
1853 return part1 + part2;
1857 tanh_c1f = std::tanh(c1 * f);
1858 N = f * tanh_c1f * inv_denom;
1869 dfacdF0(-
fac * (1 / F0 +
fac * F0 *
1874 return f * f * dfacdT;
1878 return 2.0 * f *
fac;
1893 "Temperature is stored internally in SimpleFrequencyScaling\n"
1894 "but for interface reasons this function nevertheless takes "
1895 "temprature as an input\n"
1896 "For some reason, the two temperatures disagree, so regardless, "
1897 "you have encountered\n"
1898 "a path of the code that should never be encountered. The two "
1899 "temperatures are: ",
1900 T,
" and ", t_,
" K")
1921 Numeric I0, Numeric T0, Numeric T, Numeric F0, Numeric E0, Numeric QT,
1922 Numeric QT0, Numeric dQTdT, Numeric r, Numeric drdSELFVMR,
1923 const Numeric drdT) noexcept
1925 I0, r, drdSELFVMR, drdT, QT0, QT, dQTdT,
boltzman_ratio(T, T0, E0),
1935 Numeric T, Numeric r1,
1936 Numeric r2, Numeric c2,
1939 :
k(c3 * (r1 * x - r2) * (A21 / c2)),
dkdF0(-2.0 *
k / F0),
1940 dkdr1(c3 * x * (A21 / c2)),
dkdr2(-c3 * (A21 / c2)),
e(c3 * r2 * A21),
1952 operator()(Numeric r, Numeric drdSELFVMR, Numeric drdT) &&
noexcept {
1959 Numeric F0, Numeric A21, Numeric T, Numeric g1, Numeric g2, Numeric r1,
1960 Numeric r2, Numeric r, Numeric drdSELFVMR, Numeric drdT) noexcept
1963 F0, A21, T, r1, r2, c0 * F0 * F0 * F0, c1 * F0,
1964 g2 / g1)(r, drdSELFVMR, drdT)) {}
1967 Numeric
K1,
dK1dT,
K2,
dK2dT,
dK2dF0,
K3,
dK3dT,
dK3dF0,
dK3dTl,
dK3dTu,
K4,
1971 Numeric T, Numeric T0, Numeric F0, Numeric E0, Numeric Tl, Numeric Evl,
1972 Numeric Tu, Numeric Evu, Numeric gamma, Numeric gamma_ref, Numeric r_low,
1973 Numeric r_upp) noexcept
1995 operator()(Numeric I0, Numeric QT0, Numeric QT, Numeric dQTdT, Numeric r,
1996 Numeric drdSELFVMR, Numeric drdT) &&
noexcept {
1997 return {I0, QT0, QT, dQTdT, r, drdSELFVMR, drdT,
K1,
2005 Numeric I0, Numeric T0, Numeric T, Numeric Tl, Numeric Tu, Numeric F0,
2006 Numeric E0, Numeric Evl, Numeric Evu, Numeric QT, Numeric QT0,
2007 Numeric dQTdT, Numeric r, Numeric drdSELFVMR, Numeric drdT) noexcept
2010 T, T0, F0, E0, Tl, Evl, Tu, Evu,
stimulated_emission(T, F0),
2034 const Numeric *it0 = f_grid.data_handle();
2035 const Numeric *itn = it0 + f_grid.size();
2036 const Numeric *itl = std::lower_bound(it0, itn, fl);
2038 std::distance(itl, std::upper_bound(itl, itn, fu))};
2049 const Numeric flc,
const Numeric fuc,
const Numeric fls,
const Numeric fus,
2050 const Vector &f_grid,
const Vector &sparse_f_grid)
ARTS_NOEXCEPT {
2055 const Index nvs = sparse_f_grid.size();
2058 const Numeric *it0s = sparse_f_grid.data_handle();
2059 const Numeric *itns = it0s + nvs;
2060 const Numeric *itlc =
2061 std::lower_bound(it0s, itns, std::nextafter(flc, fuc));
2062 const Numeric *ituc =
2063 std::upper_bound(itlc, itns, std::nextafter(fuc, flc));
2064 const Numeric *itls = std::upper_bound(itlc, ituc, fls);
2065 const Numeric *itus = std::lower_bound(itls, ituc, fus);
2074 Index
const beg_lr =
2075 std::distance(it0s, itlc);
2076 Index end_lr = std::distance(it0s, itls);
2079 Index beg_ur = std::distance(it0s, itus);
2082 Index
const end_ur =
2083 std::distance(it0s, ituc);
2087 (end_lr <= 0 or end_lr >= nvs) ? flc : sparse_f_grid[end_lr];
2089 (beg_ur <= 0 or beg_ur >= nvs) ? fuc : sparse_f_grid[beg_ur];
2092 const Numeric *it0 = f_grid.data_handle();
2093 const Numeric *itn = it0 + f_grid.size();
2094 const Numeric *itl = std::lower_bound(it0, itn, fl);
2095 const Numeric *itu =
2096 std::upper_bound(itl, itn, std::nextafter(fu, fl));
2098 return {std::distance(it0, itl),
2099 std::distance(itl, itu),
2107 const Numeric flc,
const Numeric fuc,
const Numeric fls,
const Numeric fus,
2108 const Vector &f_grid,
const Vector &sparse_f_grid)
ARTS_NOEXCEPT {
2113 const Index nvs = sparse_f_grid.size();
2114 const Index nv = f_grid.size();
2117 const Numeric *
const it0s = sparse_f_grid.data_handle();
2118 const Numeric *
const itns = it0s + nvs;
2119 const Numeric *itlc =
2120 std::lower_bound(it0s, itns, std::nextafter(flc, fuc));
2121 const Numeric *ituc =
2122 std::upper_bound(itlc, itns, std::nextafter(fuc, flc));
2123 const Numeric *itls =
2124 std::upper_bound(itlc, ituc, std::nextafter(fls, flc));
2125 const Numeric *itus =
2126 std::lower_bound(itls, ituc, std::nextafter(fus, fuc));
2143 while (std::distance(it0s, itls) % 3)
2145 while (std::distance(itlc, itls) % 3 == 1)
2147 while (std::distance(it0s, itus) % 3)
2149 while (std::distance(itus, ituc) % 3 == 1)
2153 const Numeric *
const it0 = f_grid.data_handle();
2154 const Numeric *
const itn = it0 + nv;
2158 if (itls not_eq itns) {
2159 itl = std::lower_bound(it0, itn, *itls);
2161 itl = std::lower_bound(it0, itn, flc);
2164 if (itus not_eq itns and itl not_eq itn) {
2165 itu = std::upper_bound(itl, itn,
2166 std::nextafter(*itus, *itl));
2168 itu = std::lower_bound(itl, itn,
2169 std::nextafter(fuc, flc));
2172 return {std::distance(it0, itl), std::distance(itl, itu),
2173 std::distance(it0s, itlc), std::distance(itlc, itls),
2174 std::distance(it0s, itus), std::distance(itus, ituc)};
2199 return std::distance(
2201 std::find_if(derivs.cbegin(), derivs.cend(),
2202 [](
auto &deriv) { return deriv.deriv == nullptr; }));
2210 const Numeric *
const f;
2217 [[nodiscard]] Index
jac_pos(Index iv, Index ij)
const noexcept {
2222 ComplexMatrix &dN_,
const Vector &f_grid,
const Index start,
2224 const bool do_nlte_) noexcept
2225 :
F(F_.data_handle() + start),
2226 dF(dF_.data_handle() + start * derivs_.size()),
2227 N(N_.data_handle() + start),
2228 dN(dN_.data_handle() + start * derivs_.size()),
2229 f(f_grid.data_handle() + start),
size(nv),
derivs(derivs_),
2234 std::vector<Complex> &dN_,
const Numeric &f_lim,
2236 :
F(&F_),
dF(dF_.data()),
N(&N_),
dN(dN_.data()),
f(&f_lim),
size(1),
2244 for (Index iv = 0; iv <
size; iv++) {
2251 for (Index iv = 0; iv <
size; iv++) {
2263 return std::visit([&](
auto &&LS) {
return LS.dFdT(dXdT, T); }, ls);
2267 return std::visit([](
auto &&LS) {
return LS.dFdf(); },
ls);
2271 return std::visit([](
auto &&LS) {
return LS.dFdF0(); },
ls);
2275 return std::visit([dfdH](
auto &&LS) {
return LS.dFdH(dfdH); }, ls);
2279 return std::visit([&](
auto &&LS) {
return LS.dFdVMR(dXdVMR); }, ls);
2283 return std::visit([
d](
auto &&LS) {
return LS.dFdFVC(
d); }, ls);
2287 return std::visit([
d](
auto &&LS) {
return LS.dFdETA(
d); }, ls);
2291 return std::visit([
d](
auto &&LS) {
return LS.dFdDV(
d); }, ls);
2295 return std::visit([
d](
auto &&LS) {
return LS.dFdD0(
d); }, ls);
2299 return std::visit([
d](
auto &&LS) {
return LS.dFdG0(
d); }, ls);
2303 return std::visit([
d](
auto &&LS) {
return LS.dFdD2(
d); }, ls);
2307 return std::visit([
d](
auto &&LS) {
return LS.dFdG2(
d); }, ls);
2311 return std::visit([](
auto &&LS) {
return LS.F; },
ls);
2315 return std::visit([f](
auto &&LS) {
return LS(f); }, ls);
2319 const Numeric DC,
const Numeric DZ,
2320 bool manually_mirrored) noexcept
2322 if (not manually_mirrored) {
2327 case Type::LP: [[fallthrough]];
2331 case Type::VP: [[fallthrough]];
2333 ls =
Voigt(F0, X, DC, DZ);
2335 case Type::SDVP: [[fallthrough]];
2336 case Type::SplitSDVP:
2339 case Type::HTP: [[fallthrough]];
2340 case Type::SplitHTP:
2350 const Numeric F0,
const Output &X,
const Numeric DC,
2354 case Absorption::MirroringType::Lorentz:
2357 case Absorption::MirroringType::SameAsLineShape:
2360 case Absorption::MirroringType::Manual:
2363 case Absorption::MirroringType::None:
2365 case Absorption::MirroringType::FINAL: {
2371 return std::visit([T, f](
auto &&LSN) {
return LSN.dNdT(T, f); }, ls_norm);
2375 return std::visit([f](
auto &&LS) {
return LS.dNdf(f); }, ls_norm);
2379 return std::visit([](
auto &&LS) {
return LS.dNdF0(); },
ls_norm);
2383 return std::visit([f](
auto &&LSN) {
return LSN(f); }, ls_norm);
2387 const Numeric F0,
const Numeric T) noexcept
2390 case Absorption::NormalizationType::None:
2392 case Absorption::NormalizationType::RQ:
2395 case Absorption::NormalizationType::VVH:
2398 case Absorption::NormalizationType::VVW:
2401 case Absorption::NormalizationType::SFS:
2404 case Absorption::NormalizationType::FINAL: {
2410 return scale * std::visit([](
auto &&
S) {
return S.S; },
ls_str);
2414 return scale * std::visit([](
auto &&LSN) {
return LSN.dSdT(); },
ls_str);
2418 return scale * std::visit([](
auto &&LS) {
return LS.dSdI0(); },
ls_str);
2422 return scale * std::visit([](
auto &&LS) {
return LS.dSdF0(); },
ls_str);
2426 return scale * std::visit([](
auto &&LS) {
return LS.dSdNLTEu(); },
ls_str);
2430 return scale * std::visit([](
auto &&LS) {
return LS.dSdNLTEl(); },
ls_str);
2435 ? std::visit([](
auto &&
S) {
return S.S; },
ls_str)
2437 scale * std::visit([](
auto &&LS) {
return LS.dSdSELFVMR(); },
ls_str);
2443 : std::visit([](
auto &&
S) {
return S.S; },
ls_str);
2447 return scale * std::visit([](
auto &&
S) {
return S.N; },
ls_str);
2451 return scale * std::visit([](
auto &&LSN) {
return LSN.dNdT(); },
ls_str);
2455 return scale * std::visit([](
auto &&LS) {
return LS.dNdI0(); },
ls_str);
2459 return scale * std::visit([](
auto &&LS) {
return LS.dNdF0(); },
ls_str);
2463 return scale * std::visit([](
auto &&LS) {
return LS.dNdNLTEu(); },
ls_str);
2467 return scale * std::visit([](
auto &&LS) {
return LS.dNdNLTEl(); },
ls_str);
2472 ? std::visit([](
auto &&
S) {
return S.N; },
ls_str)
2474 scale * std::visit([](
auto &&LS) {
return LS.dNdSELFVMR(); },
ls_str);
2480 : std::visit([](
auto &&
S) {
return S.N; },
ls_str);
2485 Species::Species other)
noexcept {
2487 self_species = self;
2488 scaling_species = other;
2493 const Numeric T,
const Numeric QT,
const Numeric QT0,
const Numeric dQTdT,
2494 const Numeric r,
const Numeric drdSELFVMR,
const Numeric drdT,
2496 const Index line_index) noexcept
2498 const auto &line = band.lines[line_index];
2499 switch (band.population) {
2500 case Absorption::PopulationType::ByHITRANFullRelmat:
2501 case Absorption::PopulationType::ByHITRANRosenkranzRelmat:
2502 case Absorption::PopulationType::ByMakarovFullRelmat:
2503 case Absorption::PopulationType::ByRovibLinearDipoleLineMixing:
2504 case Absorption::PopulationType::LTE:
2507 QT0, dQTdT, r, drdSELFVMR, drdT);
2509 case Absorption::PopulationType::NLTE: {
2510 const auto [r_low, r_upp] = nlte.get_ratio_params(band, line_index);
2512 line.gupp, r_low, r_upp, r,
2515 case Absorption::PopulationType::VibTemps: {
2516 const auto [E_low, E_upp, T_low, T_upp] = nlte.get_vibtemp_params(band, T);
2517 ls_str = VibrationalTemperaturesNonLocalThermodynamicEquilibrium(
2518 line.I0, band.T0, T, T_low, T_upp, line.F0, line.E0, E_low, E_upp, QT,
2519 QT0, dQTdT, r, drdSELFVMR, drdT);
2521 case Absorption::PopulationType::FINAL: {
2526#define CutInternalDerivativesImpl(X, Y) \
2527 else if (deriv == Jacobian::Line::Shape##X##Y) { \
2528 const Numeric d = value.n; \
2529 const Complex dFm = std::conj(ls_mirr.dFd##X(d) - ls_mirr_cut.dFd##X(d)); \
2530 const Complex dFls = ls.dFd##X(d) - ls_cut.dFd##X(d) + dFm; \
2531 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls; \
2533 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls; \
2537#define CutInternalDerivatives(X) \
2538 CutInternalDerivativesImpl(X, X0) CutInternalDerivativesImpl(X, X1) \
2539 CutInternalDerivativesImpl(X, X2) CutInternalDerivativesImpl(X, X3)
2541#define InternalDerivativesImpl(X, Y) \
2542 else if (deriv == Jacobian::Line::Shape##X##Y) { \
2543 const Numeric d = value.n; \
2544 const Complex dFm = std::conj(ls_mirr.dFd##X(d)); \
2545 const Complex dFls = ls.dFd##X(d) + dFm; \
2546 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls; \
2548 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls; \
2552#define InternalDerivatives(X) \
2553 InternalDerivativesImpl(X, X0) InternalDerivativesImpl(X, X1) \
2554 InternalDerivativesImpl(X, X2) InternalDerivativesImpl(X, X3)
2556#define InternalDerivativesGImpl(X) \
2557 else if (deriv == Jacobian::Line::ShapeG##X) { \
2558 const Numeric dLM = value.n; \
2559 com.dF[com.jac_pos(iv, ij)] += dLM * S * Fls; \
2561 com.dN[com.jac_pos(iv, ij)] += dLM * DS * Fls; \
2565#define InternalDerivativesG \
2566 InternalDerivativesGImpl(X0) InternalDerivativesGImpl(X1) \
2567 InternalDerivativesGImpl(X2) InternalDerivativesGImpl(X3)
2569#define InternalDerivativesYImpl(X) \
2570 else if (deriv == Jacobian::Line::ShapeY##X) { \
2571 const Complex dLM = Complex(0, -value.n); \
2572 com.dF[com.jac_pos(iv, ij)] += dLM * S * Fls; \
2574 com.dN[com.jac_pos(iv, ij)] += dLM * DS * Fls; \
2578#define InternalDerivativesY \
2579 InternalDerivativesYImpl(X0) InternalDerivativesYImpl(X1) \
2580 InternalDerivativesYImpl(X2) InternalDerivativesYImpl(X3)
2582#define InternalDerivativesSetupImpl(X, Y, ...) \
2583 else if (deriv == Jacobian::Line::Shape##X##Y) { \
2585 band.BroadeningSpeciesPosition(deriv.Target().species_id); \
2586 if (pos >= 0 __VA_OPT__(and pos == __VA_ARGS__)) { \
2587 derivs[ij].value.n = (__VA_OPT__(1 ? 1 :) vmrs[pos]) * \
2588 band.lines[i].lineshape[pos].dX( \
2589 T, band.T0, P, Jacobian::Line::Shape##X##Y); \
2591 derivs[ij].value.n = 0; \
2595#define InternalDerivativesSetup(X, ...) \
2596 InternalDerivativesSetupImpl(X, X0 __VA_OPT__(, __VA_ARGS__)) \
2597 InternalDerivativesSetupImpl(X, X1 __VA_OPT__(, __VA_ARGS__)) \
2598 InternalDerivativesSetupImpl(X, X2 __VA_OPT__(, __VA_ARGS__)) \
2599 InternalDerivativesSetupImpl(X, X3 __VA_OPT__(, __VA_ARGS__))
2669 const Numeric &T,
const Numeric &dfdH,
2672 const Index nv = com.size;
2673 const bool do_nlte = com.do_nlte;
2675 const Numeric Si = ls_str.S();
2676 const Numeric DNi = ls_str.N();
2678 for (Index iv = 0; iv < nv; iv++) {
2679 const Numeric f = com.f[iv];
2681 const Numeric Sn = ls_norm(f);
2682 const Numeric S = Sz * Sn * Si;
2683 const Numeric DS = Sz * Sn * DNi;
2684 const Complex Fm = std::conj(ls_mirr(f) - ls_mirr_cut.F());
2685 const Complex Fls = ls(f) - ls_cut.F() + Fm;
2686 com.F[iv] += S * LM * Fls;
2688 com.N[iv] += DS * LM * Fls;
2691 for (
const auto &[lt, value, ij, deriv_ptr] : derivs) {
2694 const auto &deriv = *deriv_ptr;
2696 if (deriv == Jacobian::Atm::Temperature) {
2697 const auto &dXdT = value.o;
2698 const Numeric dSn = ls_norm.dNdT(T, f);
2699 const Numeric dSi = ls_str.dSdT();
2700 const Complex dLM(dXdT.G, -dXdT.Y);
2704 const Complex dFls = ls.dFdT(dXdT, T) - ls_cut.dFdT(dXdT, T) + dFm;
2705 com.dF[com.jac_pos(iv, ij)] +=
2706 S * (LM * dFls + dLM * Fls) + Sz * (dSn * Si + Sn * dSi) * LM * Fls;
2708 const Numeric dDSi = ls_str.dNdT();
2709 com.dN[com.jac_pos(iv, ij)] += DS * (LM * dFls + dLM * Fls) +
2710 Sz * (dSn * Si + Sn * dDSi) * LM * Fls;
2712 }
else if (deriv.is_wind()) {
2713 const Complex dFm = std::conj(ls_mirr.dFdf() - ls_mirr_cut.dFdf());
2714 const Complex dFls = ls.dFdf() - ls_cut.dFdf() + dFm;
2715 const Numeric dS = Sz * ls_norm.dNdf(f) * Si;
2716 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dS * LM * Fls;
2718 const Numeric dDS = Sz * ls_norm.dNdf(f) * DNi;
2719 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + dDS * LM * Fls;
2721 }
else if (deriv.is_mag()) {
2723 std::conj(ls_mirr.dFdH(-dfdH) - ls_mirr_cut.dFdH(-dfdH));
2724 const Complex dFls = ls.dFdH(dfdH) - ls_cut.dFdH(dfdH) + dFm;
2725 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls;
2727 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls;
2729 }
else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR) {
2730 com.dF[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dSdSELFVMR() * LM * Fls;
2732 com.dN[com.jac_pos(iv, ij)] +=
2733 Sz * Sn * ls_str.dNdSELFVMR() * LM * Fls;
2735 }
else if (deriv.Target().needQuantumIdentity()) {
2736 if (deriv == Jacobian::Line::VMR) {
2737 const auto &dXdVMR = value.o;
2741 const Complex dLM(dXdVMR.G, -dXdVMR.Y);
2742 const Complex dFls = ls.dFdVMR(dXdVMR) - ls_cut.dFdVMR(dXdVMR) + dFm;
2743 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dLM * S * Fls;
2745 if (self_species == deriv.Target().species_id) {
2746 com.dF[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dSdSELFVMR() * LM * Fls;
2749 if (ls_str.scaler() == deriv.Target().species_id) {
2750 com.dF[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dSdOTHERVMR_if() * LM * Fls;
2753 if (ls_str.scaler() == Species::Species::Bath) {
2754 com.dF[com.jac_pos(iv, ij)] -= Sz * Sn * ls_str.dSdOTHERVMR_if() * LM * Fls;
2758 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + dLM * DS * Fls;
2760 if (self_species == deriv.QuantumIdentity().Species()) {
2761 com.dN[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dNdSELFVMR() * LM * Fls;
2764 if (ls_str.scaler() == deriv.QuantumIdentity().Species()) {
2765 com.dN[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dNdOTHERVMR_if() * LM * Fls;
2768 if (ls_str.scaler() == Species::Species::Bath) {
2769 com.dN[com.jac_pos(iv, ij)] -= Sz * Sn * ls_str.dNdOTHERVMR_if() * LM * Fls;
2773 if (lt == Quantum::Number::StateMatchType::Full) {
2774 if (deriv == Jacobian::Line::Center) {
2775 const Numeric
d = ls_norm.dNdF0() * Si + ls_str.dSdF0() * Sn;
2777 std::conj(ls_mirr.dFdF0() - ls_mirr_cut.dFdF0());
2778 const Complex dFls = ls.dFdF0() - ls_cut.dFdF0() + dFm;
2779 const Numeric dS = Sz *
d;
2780 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dS * LM * Fls;
2782 const Numeric Dd = ls_norm.dNdF0() * DNi + ls_str.dNdF0() * Sn;
2783 const Numeric DdS = Sz * Dd;
2784 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + DdS * LM * Fls;
2786 }
else if (deriv == Jacobian::Line::Strength) {
2787 const Numeric dS = Sz * Sn * ls_str.dSdI0();
2788 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
2790 const Numeric DdS = Sz * Sn * ls_str.dNdI0();
2791 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
2800 }
else if (lt == Quantum::Number::StateMatchType::Level) {
2801 if (deriv == Jacobian::Line::NLTE) {
2803 const Numeric dS = Sz * Sn * ls_str.dSdNLTEu();
2804 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
2807 const Numeric DdS = Sz * Sn * ls_str.dNdNLTEu();
2808 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
2813 const Numeric dS = Sz * Sn * ls_str.dSdNLTEl();
2814 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
2817 const Numeric DdS = Sz * Sn * ls_str.dNdNLTEl();
2818 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
2893 const Numeric &T,
const Numeric &dfdH,
const Numeric &Sz,
2895 const Index nv = com.size;
2896 const bool do_nlte = com.do_nlte;
2898 const Numeric Si = ls_str.S();
2899 const Numeric DNi = ls_str.N();
2901 for (Index iv = 0; iv < nv; iv++) {
2902 const Numeric f = com.f[iv];
2904 const Numeric Sn = ls_norm(f);
2905 const Numeric S = Sz * Sn * Si;
2906 const Numeric DS = Sz * Sn * DNi;
2907 const Complex Fm = std::conj(ls_mirr(f));
2908 const Complex Fls = ls(f) + Fm;
2909 com.F[iv] += S * LM * Fls;
2911 com.N[iv] += DS * LM * Fls;
2914 for (
const auto &[lt, value, ij, deriv_ptr] : derivs) {
2917 const auto &deriv = *deriv_ptr;
2919 if (deriv == Jacobian::Atm::Temperature) {
2920 const auto &dXdT = value.o;
2921 const Numeric dSn = ls_norm.dNdT(T, f);
2922 const Numeric dSi = ls_str.dSdT();
2923 const Complex dLM(dXdT.G, -dXdT.Y);
2924 const Complex dFm = std::conj(ls_mirr.dFdT(
mirroredOutput(dXdT), T));
2925 const Complex dFls = ls.dFdT(dXdT, T) + dFm;
2926 com.dF[com.jac_pos(iv, ij)] +=
2927 S * (LM * dFls + dLM * Fls) + Sz * (dSn * Si + Sn * dSi) * LM * Fls;
2929 const Numeric dDSi = ls_str.dNdT();
2930 com.dN[com.jac_pos(iv, ij)] += DS * (LM * dFls + dLM * Fls) +
2931 Sz * (dSn * Si + Sn * dDSi) * LM * Fls;
2933 }
else if (deriv.is_wind()) {
2934 const Complex dFm = std::conj(ls_mirr.dFdf());
2935 const Complex dFls = ls.dFdf() + dFm;
2936 const Numeric dS = Sz * ls_norm.dNdf(f) * Si;
2937 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dS * LM * Fls;
2939 const Numeric dDS = Sz * ls_norm.dNdf(f) * DNi;
2940 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + dDS * LM * Fls;
2942 }
else if (deriv.is_mag()) {
2943 const Complex dFm = std::conj(ls_mirr.dFdH(-dfdH));
2944 const Complex dFls = ls.dFdH(dfdH) + dFm;
2945 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls;
2947 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls;
2949 }
else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR) {
2950 com.dF[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dSdSELFVMR() * LM * Fls;
2952 com.dN[com.jac_pos(iv, ij)] +=
2953 Sz * Sn * ls_str.dNdSELFVMR() * LM * Fls;
2955 }
else if (deriv.Target().needQuantumIdentity()) {
2956 if (deriv == Jacobian::Line::VMR) {
2957 const auto &dXdVMR = value.o;
2958 const Complex dFm = std::conj(ls_mirr.dFdVMR(
mirroredOutput(dXdVMR)));
2959 const Complex dLM(dXdVMR.G, -dXdVMR.Y);
2960 const Complex dFls = ls.dFdVMR(dXdVMR) + dFm;
2961 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dLM * S * Fls;
2963 if (self_species == deriv.Target().species_id) {
2964 com.dF[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dSdSELFVMR() * LM * Fls;
2967 if (ls_str.scaler() == deriv.Target().species_id) {
2968 com.dF[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dSdOTHERVMR_if() * LM * Fls;
2971 if (ls_str.scaler() == Species::Species::Bath) {
2972 com.dF[com.jac_pos(iv, ij)] -= Sz * Sn * ls_str.dSdOTHERVMR_if() * LM * Fls;
2976 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + dLM * DS * Fls;
2978 if (self_species == deriv.QuantumIdentity().Species()) {
2979 com.dN[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dNdSELFVMR() * LM * Fls;
2982 if (ls_str.scaler() == deriv.QuantumIdentity().Species()) {
2983 com.dN[com.jac_pos(iv, ij)] += Sz * Sn * ls_str.dNdOTHERVMR_if() * LM * Fls;
2986 if (ls_str.scaler() == Species::Species::Bath) {
2987 com.dN[com.jac_pos(iv, ij)] -= Sz * Sn * ls_str.dNdOTHERVMR_if() * LM * Fls;
2991 if (lt == Quantum::Number::StateMatchType::Full) {
2992 if (deriv == Jacobian::Line::Center) {
2993 const Numeric
d = ls_norm.dNdF0() * Si + ls_str.dSdF0() * Sn;
2994 const Complex dFm = std::conj(ls_mirr.dFdF0());
2995 const Complex dFls = ls.dFdF0() + dFm;
2996 const Numeric dS = Sz *
d;
2997 com.dF[com.jac_pos(iv, ij)] += S * LM * dFls + dS * LM * Fls;
2999 const Numeric Dd = ls_norm.dNdF0() * DNi + ls_str.dNdF0() * Sn;
3000 const Numeric DdS = Sz * Dd;
3001 com.dN[com.jac_pos(iv, ij)] += DS * LM * dFls + DdS * LM * Fls;
3003 }
else if (deriv == Jacobian::Line::Strength) {
3004 const Numeric dS = Sz * Sn * ls_str.dSdI0();
3005 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
3007 const Numeric DdS = Sz * Sn * ls_str.dNdI0();
3008 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
3017 }
else if (lt == Quantum::Number::StateMatchType::Level) {
3018 if (deriv == Jacobian::Line::NLTE) {
3020 const Numeric dS = Sz * Sn * ls_str.dSdNLTEu();
3021 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
3024 const Numeric DdS = Sz * Sn * ls_str.dNdNLTEu();
3025 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
3030 const Numeric dS = Sz * Sn * ls_str.dSdNLTEl();
3031 com.dF[com.jac_pos(iv, ij)] += dS * LM * Fls;
3034 const Numeric DdS = Sz * Sn * ls_str.dNdNLTEl();
3035 com.dN[com.jac_pos(iv, ij)] += DdS * LM * Fls;
3088 const Numeric &H,
const Numeric &sparse_lim,
const Numeric &DC,
3092 const bool do_nlte = ls_str.do_nlte();
3093 const bool do_cutoff = band.cutoff not_eq Absorption::CutoffType::None;
3094 const Numeric fu = band.CutoffFreq(i, X.D0);
3095 const Numeric fl = band.CutoffFreqMinus(i, X.D0);
3096 const Numeric fus = band.lines[i].F0 + sparse_lim;
3097 const Numeric fls = band.lines[i].F0 - sparse_lim;
3100 const auto [dense_start, dense_size, sparse_low_start, sparse_low_size,
3101 sparse_upp_start, sparse_upp_size] =
3104 if ((dense_size + sparse_low_size + sparse_upp_size) == 0)
3108 ComputeValues comval(com.F, com.dF, com.N, com.dN, com.f_grid, dense_start,
3109 dense_size, derivs, do_nlte);
3113 sparse_com.F, sparse_com.dF, sparse_com.N, sparse_com.dN,
3114 sparse_com.f_grid, sparse_low_start, sparse_low_size, derivs, do_nlte);
3116 sparse_com.F, sparse_com.dF, sparse_com.N, sparse_com.dN,
3117 sparse_com.f_grid, sparse_upp_start, sparse_upp_size, derivs, do_nlte);
3119 const Index nz = band.ZeemanCount(i, zeeman_polarization);
3120 for (Index iz = 0; iz < nz; iz++) {
3121 const Numeric dfdH = band.ZeemanSplitting(i, zeeman_polarization, iz);
3122 const Numeric Sz = band.ZeemanStrength(i, zeeman_polarization, iz);
3123 const Complex LM = Complex(1 + X.G, -X.Y);
3124 Calculator ls(band.lineshapetype, band.lines[i].F0, X, DC, dfdH * H,
3125 band.mirroring == Absorption::MirroringType::Manual);
3126 Calculator ls_mirr(band.mirroring, band.lineshapetype, band.lines[i].F0, X,
3137 ls_str, derivs, LM, T, dfdH, Sz, band.Species());
3139 if (sparse_low_size) {
3141 ls_mirr_cut, ls_str, derivs, LM, T, dfdH, Sz,
3145 if (sparse_upp_size) {
3147 ls_mirr_cut, ls_str, derivs, LM, T, dfdH, Sz,
3152 frequency_loop(comval, ls, ls_mirr, ls_norm, ls_str, derivs, LM, T, dfdH,
3153 Sz, band.Species());
3155 if (sparse_low_size) {
3156 frequency_loop(sparse_low_range, ls, ls_mirr, ls_norm, ls_str, derivs,
3157 LM, T, dfdH, Sz, band.Species());
3160 if (sparse_upp_size) {
3161 frequency_loop(sparse_upp_range, ls, ls_mirr, ls_norm, ls_str, derivs,
3162 LM, T, dfdH, Sz, band.Species());
3172 const Numeric &H,
const Numeric &sparse_lim,
const Numeric &DC,
3176 const bool do_nlte = ls_str.do_nlte();
3177 const bool do_cutoff = band.cutoff not_eq Absorption::CutoffType::None;
3178 const Numeric fu = band.CutoffFreq(i, X.D0);
3179 const Numeric fl = band.CutoffFreqMinus(i, X.D0);
3180 const Numeric fus = band.lines[i].F0 + sparse_lim;
3181 const Numeric fls = band.lines[i].F0 - sparse_lim;
3184 const auto [dense_start, dense_size, sparse_low_start, sparse_low_size,
3185 sparse_upp_start, sparse_upp_size] =
3188 if ((dense_size + sparse_low_size + sparse_upp_size) == 0)
3192 ComputeValues comval(com.F, com.dF, com.N, com.dN, com.f_grid, dense_start,
3193 dense_size, derivs, do_nlte);
3197 sparse_com.F, sparse_com.dF, sparse_com.N, sparse_com.dN,
3198 sparse_com.f_grid, sparse_low_start, sparse_low_size, derivs, do_nlte);
3200 sparse_com.F, sparse_com.dF, sparse_com.N, sparse_com.dN,
3201 sparse_com.f_grid, sparse_upp_start, sparse_upp_size, derivs, do_nlte);
3203 const Index nz = band.ZeemanCount(i, zeeman_polarization);
3204 for (Index iz = 0; iz < nz; iz++) {
3205 const Numeric dfdH = band.ZeemanSplitting(i, zeeman_polarization, iz);
3206 const Numeric Sz = band.ZeemanStrength(i, zeeman_polarization, iz);
3207 const Complex LM = Complex(1 + X.G, -X.Y);
3208 Calculator ls(band.lineshapetype, band.lines[i].F0, X, DC, dfdH * H,
3209 band.mirroring == Absorption::MirroringType::Manual);
3210 Calculator ls_mirr(band.mirroring, band.lineshapetype, band.lines[i].F0, X,
3221 ls_str, derivs, LM, T, dfdH, Sz, band.Species());
3223 if (sparse_low_size) {
3225 ls_mirr_cut, ls_str, derivs, LM, T, dfdH, Sz,
3229 if (sparse_upp_size) {
3231 ls_mirr_cut, ls_str, derivs, LM, T, dfdH, Sz,
3236 frequency_loop(comval, ls, ls_mirr, ls_norm, ls_str, derivs, LM, T, dfdH,
3237 Sz, band.Species());
3239 if (sparse_low_size) {
3240 frequency_loop(sparse_low_range, ls, ls_mirr, ls_norm, ls_str, derivs,
3241 LM, T, dfdH, Sz, band.Species());
3244 if (sparse_upp_size) {
3245 frequency_loop(sparse_upp_range, ls, ls_mirr, ls_norm, ls_str, derivs,
3246 LM, T, dfdH, Sz, band.Species());
3293 const Numeric &T,
const Numeric &H,
const Numeric &DC,
3297 const bool do_nlte = ls_str.do_nlte();
3298 const bool do_cutoff = band.cutoff not_eq Absorption::CutoffType::None;
3299 const Numeric fu = band.CutoffFreq(i, X.D0);
3300 const Numeric fl = band.CutoffFreqMinus(i, X.D0);
3303 const auto [cutstart, cutsize] =
limited_range(fl, fu, com.f_grid);
3308 ComputeValues comval(com.F, com.dF, com.N, com.dN, com.f_grid, cutstart,
3309 cutsize, derivs, do_nlte);
3311 const Index nz = band.ZeemanCount(i, zeeman_polarization);
3312 for (Index iz = 0; iz < nz; iz++) {
3313 const Numeric dfdH = band.ZeemanSplitting(i, zeeman_polarization, iz);
3314 const Numeric Sz = band.ZeemanStrength(i, zeeman_polarization, iz);
3315 const Complex LM = Complex(1 + X.G, -X.Y);
3316 Calculator ls(band.lineshapetype, band.lines[i].F0, X, DC, dfdH * H,
3317 band.mirroring == Absorption::MirroringType::Manual);
3318 Calculator ls_mirr(band.mirroring, band.lineshapetype, band.lines[i].F0, X,
3329 ls_str, derivs, LM, T, dfdH, Sz, band.Species());
3332 frequency_loop(comval, ls, ls_mirr, ls_norm, ls_str, derivs, LM, T, dfdH,
3333 Sz, band.Species());
3374 const Numeric &T,
const Numeric &H,
const Numeric &sparse_lim,
3375 const Numeric QT,
const Numeric QT0,
const Numeric dQTdT,
3376 const Numeric r,
const Numeric drdSELFVMR,
const Numeric drdT,
3379 const Index nj = jacobian_quantities.nelem();
3380 const Index nl = band.NumLines();
3386 const Numeric DC = band.DopplerConstant(T);
3388 if (not independent_per_broadener(band.lineshapetype)) {
3389 for (Index i = 0; i < nl; i++) {
3391 for (Index ij = 0; ij < nj; ij++) {
3392 const auto &deriv = jacobian_quantities[ij];
3393 derivs[ij].jac_pos = -1;
3394 derivs[ij].deriv =
nullptr;
3396 if (not deriv.propmattype())
3398 derivs[ij].jac_pos = ij;
3399 derivs[ij].deriv = &deriv;
3401 if (deriv == Jacobian::Atm::Temperature) {
3402 derivs[ij].value.o = band.ShapeParameters_dT(i, T, P,
vmrs);
3403 }
else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR) {
3404 if (not(deriv == self_tag)) {
3405 derivs[ij].jac_pos = -1;
3406 derivs[ij].deriv =
nullptr;
3408 }
else if (deriv.Target().needQuantumIdentity()) {
3409 if (deriv == Jacobian::Line::VMR) {
3410 derivs[ij].value.o =
3411 band.ShapeParameters_dVMR(i, T, P, deriv.QuantumIdentity());
3413 auto < = derivs[ij].target = {deriv.Target().qid,
3414 band.lines[i].localquanta,
3415 band.quantumidentity};
3416 if (lt == Quantum::Number::StateMatchType::Full) {
3417 if constexpr (
false) {
3431 std::remove_if(derivs.begin(), derivs.end(),
3432 [](
Derivatives &dd) { return dd.deriv == nullptr; });
3435 switch (speedup_type) {
3436 case Options::LblSpeedup::None:
3440 band, derivs, band.ShapeParameters(i, T, P,
vmrs), T, H, DC,
3441 i, zeeman_polarization);
3443 case Options::LblSpeedup::QuadraticIndependent:
3446 Normalizer(band.normalization, band.lines[i].F0, T),
3449 band, derivs, band.ShapeParameters(i, T, P,
vmrs), T, H, sparse_lim,
3450 DC, i, zeeman_polarization);
3452 case Options::LblSpeedup::LinearIndependent:
3455 Normalizer(band.normalization, band.lines[i].F0, T),
3458 band, derivs, band.ShapeParameters(i, T, P,
vmrs), T, H, sparse_lim,
3459 DC, i, zeeman_polarization);
3461 case Options::LblSpeedup::FINAL: {
3466 for (Index i = 0; i < nl; i++) {
3467 for (Index ib=0; ib<band.NumBroadeners(); ib++) {
3469 for (Index ij = 0; ij < nj; ij++) {
3470 const auto &deriv = jacobian_quantities[ij];
3471 derivs[ij].jac_pos = -1;
3472 derivs[ij].deriv =
nullptr;
3474 if (not deriv.propmattype())
3476 derivs[ij].jac_pos = ij;
3477 derivs[ij].deriv = &deriv;
3479 if (deriv == Jacobian::Atm::Temperature) {
3480 derivs[ij].value.o = band.ShapeParameters_dT(i, T, P, ib);
3481 }
else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR) {
3482 if (not(deriv == self_tag)) {
3483 derivs[ij].jac_pos = -1;
3484 derivs[ij].deriv =
nullptr;
3486 }
else if (deriv.Target().needQuantumIdentity()) {
3487 if (deriv == Jacobian::Line::VMR) {
3491 auto < = derivs[ij].target = {deriv.Target().qid,
3492 band.lines[i].localquanta,
3493 band.quantumidentity};
3494 if (lt == Quantum::Number::StateMatchType::Full) {
3495 if constexpr (
false) {
3511 std::remove_if(derivs.begin(), derivs.end(),
3512 [](
Derivatives &dd) { return dd.deriv == nullptr; });
3516 drdSELFVMR, drdT, nlte, band, i)
3518 band.broadeningspecies[ib]);
3521 switch (speedup_type) {
3522 case Options::LblSpeedup::None:
3524 ls_str, band, derivs, band.ShapeParameters(i, T, P, ib),
3525 T, H, DC, i, zeeman_polarization);
3527 case Options::LblSpeedup::QuadraticIndependent:
3530 Normalizer(band.normalization, band.lines[i].F0, T), ls_str, band,
3531 derivs, band.ShapeParameters(i, T, P, ib), T, H, sparse_lim, DC,
3532 i, zeeman_polarization);
3534 case Options::LblSpeedup::LinearIndependent:
3537 Normalizer(band.normalization, band.lines[i].F0, T), ls_str, band,
3538 derivs, band.ShapeParameters(i, T, P, ib), T, H, sparse_lim, DC,
3539 i, zeeman_polarization);
3541 case Options::LblSpeedup::FINAL: {
3554 const Numeric &isot_ratio,
const Numeric &P,
const Numeric &T,
3555 const Numeric &H,
const Numeric &sparse_lim,
3557 const Options::LblSpeedup speedup_type,
3559 [[maybe_unused]]
const Index nj = jacobian_quantities.nelem();
3560 const Index nl = band.NumLines();
3561 const Index nv = com.f_grid.nelem();
3564 ARTS_ASSERT(H >= 0,
"Only for positive H. You provided: ", H)
3565 ARTS_ASSERT(P > 0,
"Only for abs positive P. You provided: ", P)
3566 ARTS_ASSERT(T > 0,
"Only for abs positive T. You provided: ", T)
3567 ARTS_ASSERT(band.OK(),
"Band is poorly constructed. You need to use "
3568 "a detailed debugger to find out why.")
3569 ARTS_ASSERT(com.F.size() == nv,
"F is wrong size. Size is (", com.F.size(),
3570 ") but should be: (", nv,
')')
3571 ARTS_ASSERT(not com.do_nlte or com.N.size() == nv,
3572 "N is wrong size. Size is (", com.N.size(),
") but should be (",
3574 ARTS_ASSERT(nj == 0 or (com.dF.nrows() == nv and com.dF.ncols() == nj),
3575 "dF is wrong size. Size is (", com.dF.nrows(),
" x ",
3576 com.dF.ncols(),
") but should be: (", nv,
" x ", nj,
")")
3578 (com.dN.nrows() == nv and com.dN.ncols() == nj),
3579 "dN is wrong size. Size is (", com.dN.nrows(),
" x ",
3580 com.dN.ncols(),
") but should be: (", nv,
" x ", nj,
")")
3581 ARTS_ASSERT((sparse_lim > 0 and sparse_com.f_grid.size() > 1) or
3583 "Sparse limit is either 0, or the sparse frequency grid has to "
3584 "have upper and lower values")
3587 if (nv == 0 or nl == 0 or
3589 band.DoLineMixing(P))) {
3595 if (robust and band.DoLineMixing(P) and band.AnyLinemixing()) {
3596 ComputeData com_safe(com.f_grid, jacobian_quantities, com.do_nlte);
3597 ComputeData sparse_com_safe(sparse_com.f_grid, jacobian_quantities,
3598 sparse_com.do_nlte);
3600 line_loop(com_safe, sparse_com_safe, band, jacobian_quantities, nlte,
vmrs,
3601 self_tag, P, T, H, sparse_lim,
3605 self_vmr * dnumdensdVMR, dnumdensdVMR,
3607 zeeman_polarization, speedup_type);
3613 sparse_com += sparse_com_safe;
3615 line_loop(com, sparse_com, band, jacobian_quantities, nlte,
vmrs, self_tag,
3616 P, T, H, sparse_lim,
3620 self_vmr * dnumdensdVMR, dnumdensdVMR,
3622 zeeman_polarization, speedup_type);
3626#undef InternalDerivatives
3627#undef InternalDerivativesG
3628#undef InternalDerivativesY
3631 const Numeric &sparse_df)
noexcept {
3633 return f_grid.nelem() /
3635 std::abs(f_grid[f_grid.nelem() - 1] - f_grid[0]) / sparse_df);
3642 const Index nv = f_grid.nelem();
3645 std::vector<Numeric> sparse_f_grid;
3646 for (Index iv = 0; iv < nv - n; iv += n) {
3647 sparse_f_grid.emplace_back(f_grid[iv]);
3648 sparse_f_grid.emplace_back(f_grid[iv + n]);
3651 const Numeric f0 = sparse_f_grid.back();
3652 if (f0 not_eq f_grid[nv - 1]) {
3653 sparse_f_grid.emplace_back(f0);
3654 sparse_f_grid.emplace_back(f_grid[nv - 1]);
3657 return sparse_f_grid;
3663 const Vector &f_grid_sparse)
noexcept {
3664 const Index nf_sparse = f_grid_sparse.nelem();
3665 const Index nf_dense = f_grid_dense.nelem();
3670 if (nf_sparse and nf_dense)
3671 return f_grid_dense[0] >= f_grid_sparse[0] and
3672 f_grid_dense[nf_dense - 1] <= f_grid_sparse[nf_sparse - 1];
3678 const Numeric &sparse_df)
noexcept {
3680 const Index nv = f_grid.nelem();
3683 std::vector<Numeric> sparse_f_grid;
3684 for (Index iv = 0; iv < nv - n; iv += n) {
3685 sparse_f_grid.emplace_back(f_grid[iv]);
3686 sparse_f_grid.emplace_back(f_grid[iv] +
3687 0.5 * (f_grid[iv + n] - f_grid[iv]));
3688 sparse_f_grid.emplace_back(f_grid[iv + n]);
3691 const Numeric f0 = sparse_f_grid.back();
3692 if (f0 not_eq f_grid[nv - 1]) {
3693 sparse_f_grid.emplace_back(f0);
3694 sparse_f_grid.emplace_back(f0 + 0.5 * (f_grid[nv - 1] - f0));
3695 sparse_f_grid.emplace_back(f_grid[nv - 1]);
3698 return sparse_f_grid;
3704 const Index nv = f_grid.nelem();
3705 const Index sparse_nv = sparse.f_grid.nelem();
3706 const Index nj = dF.ncols();
3708 ARTS_ASSERT(do_nlte == sparse.do_nlte,
"Must have the same NLTE status")
3709 ARTS_ASSERT(sparse_nv > 1,
"Must have at least two sparse grid-points")
3710 ARTS_ASSERT(nv == 0 or (f_grid[0] == sparse.f_grid[0] and
3711 f_grid[nv - 1] >= sparse.f_grid[sparse_nv - 1]),
3712 "If there are any dense frequency points, then the sparse "
3713 "frequency points must fully cover them")
3714 ARTS_ASSERT(not(sparse_nv % 2),
"Must be multiple of to")
3716 Index sparse_iv = 0;
3717 Numeric f1 = sparse.f_grid[sparse_iv + 1];
3718 Numeric f0 = sparse.f_grid[sparse_iv];
3719 Numeric inv = 1.0 / (f1 - f0);
3720 for (Index iv = 0; iv < nv; iv++) {
3721 if (sparse_iv < (sparse_nv - 2) and f1 == f_grid[iv]) {
3723 f1 = sparse.f_grid[sparse_iv + 1];
3724 f0 = sparse.f_grid[sparse_iv];
3725 inv = 1.0 / (f1 - f0);
3728 const Numeric xm0 = f_grid[iv] - f0;
3729 const Numeric xm1 = f_grid[iv] - f1;
3730 const Numeric l0 = -xm1 * inv;
3731 const Numeric l1 = xm0 * inv;
3733 F[iv] += l0 * sparse.F[sparse_iv] + l1 * sparse.F[sparse_iv + 1];
3734 for (Index ij = 0; ij < nj; ij++) {
3736 l0 * sparse.dF(sparse_iv, ij) + l1 * sparse.dF(sparse_iv + 1, ij);
3739 N[iv] += l0 * sparse.N[sparse_iv] + l1 * sparse.N[sparse_iv + 1];
3740 for (Index ij = 0; ij < nj; ij++) {
3742 l0 * sparse.dN(sparse_iv, ij) + l1 * sparse.dN(sparse_iv + 1, ij);
3750 const Index nv = f_grid.nelem();
3751 const Index sparse_nv = sparse.f_grid.nelem();
3752 const Index nj = dF.ncols();
3754 ARTS_ASSERT(do_nlte == sparse.do_nlte,
"Must have the same NLTE status")
3755 ARTS_ASSERT(sparse_nv > 2,
"Must have at least three sparse grid-points")
3756 ARTS_ASSERT(nv == 0 or (f_grid[0] == sparse.f_grid[0] and
3757 f_grid[nv - 1] >= sparse.f_grid[sparse_nv - 1]),
3758 "If there are any dense frequency points, then the sparse "
3759 "frequency points must fully cover them")
3760 ARTS_ASSERT(not(sparse_nv % 3),
"Must be multiple of three")
3762 Index sparse_iv = 0;
3763 Numeric f2 = sparse.f_grid[sparse_iv + 2];
3764 Numeric f1 = sparse.f_grid[sparse_iv + 1];
3765 Numeric f0 = sparse.f_grid[sparse_iv];
3767 for (Index iv = 0; iv < nv; iv++) {
3768 if (sparse_iv < (sparse_nv - 3) and f2 == f_grid[iv]) {
3770 f2 = sparse.f_grid[sparse_iv + 2];
3771 f1 = sparse.f_grid[sparse_iv + 1];
3772 f0 = sparse.f_grid[sparse_iv];
3776 f_grid[iv] >= f0 and (f_grid[iv] < f2 or (f2 == f_grid[iv] and
3777 sparse_iv == sparse_nv - 3)),
3778 "Out of range frequency grid. Must be caught earlier.\n"
3779 "The sparse range is from: ",
3780 f0,
" to ", f2,
" with ", f1,
3781 " as the half-way grid point.\n"
3782 "The dense frequency is ",
3783 f_grid[iv],
" and the indices are: sparse_iv=", sparse_iv,
"; iv=", iv)
3785 const Numeric xm0 = f_grid[iv] - f0;
3786 const Numeric xm1 = f_grid[iv] - f1;
3787 const Numeric xm2 = f_grid[iv] - f2;
3788 const Numeric l0 = 0.5 * xm1 * xm2 * inv;
3789 const Numeric l1 = -xm0 * xm2 * inv;
3790 const Numeric l2 = 0.5 * xm0 * xm1 * inv;
3792 F[iv] += l0 * sparse.F[sparse_iv] + l1 * sparse.F[sparse_iv + 1] +
3793 l2 * sparse.F[sparse_iv + 2];
3794 for (Index ij = 0; ij < nj; ij++) {
3795 dF(iv, ij) += l0 * sparse.dF(sparse_iv, ij) +
3796 l1 * sparse.dF(sparse_iv + 1, ij) +
3797 l2 * sparse.dF(sparse_iv + 2, ij);
3800 N[iv] += l0 * sparse.N[sparse_iv] + l1 * sparse.N[sparse_iv + 1] +
3801 l2 * sparse.N[sparse_iv + 2];
3802 for (Index ij = 0; ij < nj; ij++) {
3803 dN(iv, ij) += l0 * sparse.dN(sparse_iv, ij) +
3804 l1 * sparse.dN(sparse_iv + 1, ij) +
3805 l2 * sparse.dN(sparse_iv + 2, ij);
This can be used to make arrays out of anything.
Index nelem() const ARTS_NOEXCEPT
Complex dFdH(Numeric dfdH) const noexcept
Complex dFdF0() const noexcept
Complex F() const noexcept
Complex dFdT(const Output &dXdT, Numeric T) const noexcept
Complex dFdD0(Numeric d) const noexcept
Complex dFdVMR(const Output &dXdVMR) const noexcept
Complex dFdG0(Numeric d) const noexcept
Complex dFdG2(Numeric d) const noexcept
Complex dFdFVC(Numeric d) const noexcept
Calculator(const Type type, const Numeric F0, const Output &X, const Numeric DC, const Numeric DZ, bool manually_mirrored) noexcept
Complex dFdD2(Numeric d) const noexcept
Complex operator()(Numeric f) noexcept
Call operator on frequency. Must call this before any of the derivatives.
Complex dFdETA(Numeric d) const noexcept
Complex dFdf() const noexcept
Complex dFdDV(Numeric d) const noexcept
Class encapsulating all supported types of intensity calculations of individual absorption lines.
Numeric N() const noexcept
The line strength source offset.
IntensityCalculator & adaptive_scaling(Numeric x, Species::Species self, Species::Species other) noexcept
Rescale the line strength parameters by x.
Numeric dSdI0() const noexcept
The line strength absorption derivative wrt the reference line strength.
Numeric dSdNLTEl() const noexcept
The line strength absorption derivative wrt either the lower state number density distribution or its...
Numeric dNdF0() const noexcept
The line strength source offset derivative wrt the line center.
Numeric dSdF0() const noexcept
The line strength absorption derivative wrt the line center.
Numeric dSdNLTEu() const noexcept
The line strength absorption derivative wrt either the upper state number density distribution or its...
Numeric dNdOTHERVMR_if() const noexcept
The line source offset derivative wrt the VMR of the rescaling species from adaptive_scaling (if self...
Numeric S() const noexcept
The line strength absorption.
Numeric dSdOTHERVMR_if() const noexcept
The line strength derivative wrt the VMR of the rescaling species from adaptive_scaling (if self not_...
Numeric dNdSELFVMR() const noexcept
The line source offset derivative wrt the VMR of the band's species.
Species::Species self_species
Species::Species scaling_species
Numeric dNdT() const noexcept
The line strength source offset derivative wrt temperature.
Numeric dNdI0() const noexcept
The line strength source offset derivative wrt the reference line strength.
Numeric dNdNLTEu() const noexcept
The line strength source offset derivative wrt either the upper state number density distribution or ...
Numeric dSdSELFVMR() const noexcept
The line strength derivative wrt the VMR of the band's species.
Numeric dNdNLTEl() const noexcept
The line strength source offset derivative wrt either the lower state number density distribution or ...
Numeric dSdT() const noexcept
The line strength absorption derivative wrt temperature.
IntensityCalculator(const Numeric T, const Numeric QT, const Numeric QT0, const Numeric dQTdT, const Numeric r, const Numeric drdSELFVMR, const Numeric drdT, const EnergyLevelMap &nlte, const Absorption::Lines &band, const Index line_index) noexcept
Numeric dNdF0() const noexcept
Numeric operator()(Numeric f) noexcept
Numeric dNdf(Numeric f) const noexcept
Numeric dNdT(Numeric T, Numeric f) const noexcept
Normalizer(const Absorption::NormalizationType type, const Numeric F0, const Numeric T) noexcept
Deals with internal derivatives, Jacobian definition, and OEM calculations.
Helper macros for debugging.
#define ARTS_ASSERT(condition,...)
Numeric stimulated_emission(Numeric T, Numeric F0)
Computes exp(-hf/kT)
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)
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)
Numeric dsingle_partition_function_dT(const Numeric &T, const Species::IsotopeRecord &ir)
Computes the partition function temperature derivative.
Numeric dstimulated_relative_emission_dF0(const Numeric F0, const Numeric T0, const Numeric T) noexcept
Computes.
Numeric boltzman_ratio(const Numeric &T, const Numeric &T0, const Numeric &E0)
Computes exp(E0/c (T - T0) / (T * T0))
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)
Numeric single_partition_function(const Numeric &T, const Species::IsotopeRecord &ir)
Computes the partition function at one temperature.
Numeric dboltzman_ratio_dT(const Numeric &boltzmann_ratio, const Numeric &T, const Numeric &E0)
Computes temperature derivatives exp(E0/k (T - T0) / (T * T0))
Numeric stimulated_relative_emission(const Numeric F0, const Numeric T0, const Numeric T) noexcept
Computes.
Numeric absorption_nlte_ratio(const Numeric &gamma, const Numeric &r_upp, const Numeric &r_low) noexcept
Computes (r_low - r_upp * gamma) / (1 - gamma)
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)
Numeric dstimulated_relative_emission_dT(const Numeric F0, const Numeric T0, const Numeric T) noexcept
Computes.
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))
#define CutInternalDerivatives(X)
#define InternalDerivatives(X)
#define InternalDerivativesSetup(X,...)
#define InternalDerivativesY
#define InternalDerivativesG
Contains the line shape namespace.
Numeric fac(const Index n)
fac
constexpr bool relaxationtype_relmat(PopulationType in) noexcept
constexpr Numeric inv_sqrt_pi
Inverse of the square root of pi.
constexpr Numeric sqrt_ln_2
Square root of natural logarithm of 2.
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric k
Boltzmann constant convenience name [J/K].
constexpr Numeric sqrt_pi
Square root of pi.
constexpr Numeric c
Speed of light convenience name [m/s].
constexpr Numeric inv_pi
Inverse of pi.
constexpr Numeric h
Planck constant convenience name [J s].
Computations of line shape derived parameters.
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.
Index sparse_f_grid_red(const Vector &f_grid, const Numeric &sparse_df) noexcept
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.
Index active_nelem(const ArrayOfDerivatives &derivs) noexcept
Helper function to find the last relevant derivative.
constexpr Numeric abs_squared(Complex z) noexcept
constexpr Output mirroredOutput(Output x) noexcept
Output to be used by mirroring calls.
void compute(ComputeData &com, ComputeData &sparse_com, const AbsorptionLines &band, const ArrayOfRetrievalQuantity &jacobian_quantities, const EnergyLevelMap &nlte, const Vector &vmrs, const ArrayOfSpeciesTag &self_tag, const Numeric &self_vmr, const Numeric &isot_ratio, const Numeric &P, const Numeric &T, const Numeric &H, const Numeric &sparse_lim, const Zeeman::Polarization zeeman_polarization, const Options::LblSpeedup speedup_type, const bool robust) ARTS_NOEXCEPT
Compute the absorption of an absorption band.
void cutoff_loop_sparse_linear(ComputeData &com, ComputeData &sparse_com, Normalizer ls_norm, const IntensityCalculator ls_str, const AbsorptionLines &band, const ArrayOfDerivatives &derivs, const Output X, const Numeric &T, const Numeric &H, const Numeric &sparse_lim, const Numeric &DC, const Index i, const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT
Cutoff considerations of the line shape.
void line_loop(ComputeData &com, ComputeData &sparse_com, const AbsorptionLines &band, const ArrayOfRetrievalQuantity &jacobian_quantities, const EnergyLevelMap &nlte, const Vector &vmrs, const ArrayOfSpeciesTag &self_tag, const Numeric &P, const Numeric &T, const Numeric &H, const Numeric &sparse_lim, const Numeric QT, const Numeric QT0, const Numeric dQTdT, const Numeric r, const Numeric drdSELFVMR, const Numeric drdT, const Zeeman::Polarization zeeman_polarization, const Options::LblSpeedup speedup_type) ARTS_NOEXCEPT
Loop all the lines of the band.
void cutoff_loop_sparse_triple(ComputeData &com, ComputeData &sparse_com, Normalizer ls_norm, const IntensityCalculator ls_str, const AbsorptionLines &band, const ArrayOfDerivatives &derivs, const Output X, const Numeric &T, const Numeric &H, const Numeric &sparse_lim, const Numeric &DC, const Index i, const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT
bool good_linear_sparse_f_grid(const Vector &f_grid_dense, const Vector &f_grid_sparse) noexcept
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
void cutoff_loop(ComputeData &com, Normalizer ls_norm, const IntensityCalculator ls_str, const AbsorptionLines &band, const ArrayOfDerivatives &derivs, const Output X, const Numeric &T, const Numeric &H, const Numeric &DC, const Index i, const Zeeman::Polarization zeeman_polarization) ARTS_NOEXCEPT
Cutoff considerations of the line shape.
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.
Vector triple_sparse_f_grid(const Vector &f_grid, const Numeric &sparse_df) noexcept
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
Vector linear_sparse_f_grid(const Vector &f_grid, const Numeric &sparse_df) ARTS_NOEXCEPT
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.
constexpr auto pow3(auto x) noexcept
power of three
constexpr auto pow4(auto x) noexcept
power of four
constexpr auto pow2(auto x) noexcept
power of two
Polarization
Zeeman polarization selection.
constexpr T abs(T x) noexcept
This file contains declerations of functions of physical character.
constexpr Numeric dnumber_density_dt(Numeric p, Numeric t) noexcept
dnumber_density_dT
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Main computational data for the line shape and strength calculations.
void interp_add_even(const ComputeData &sparse) ARTS_NOEXCEPT
Add a sparse grid to this grid via linear interpolation.
void interp_add_triplequad(const ComputeData &sparse) ARTS_NOEXCEPT
Add a sparse grid to this grid via square interpolation.
void enforce_positive_absorption() noexcept
All four fields are set to zero at i if F[i].real() < 0.
const ArrayOfDerivatives & derivs
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
ComputeValues & operator-=(const ComputeValues &cut) ARTS_NOEXCEPT
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
Index jac_pos(Index iv, Index ij) const noexcept
Struct to keep the cutoff limited range values.
Data struct for keeping derivative keys and values.
const RetrievalQuantity * deriv
union LineShape::Derivatives::Values value
Quantum::Number::StateMatch target
Complex operator()(Numeric f) noexcept
constexpr FullNonLocalThermodynamicEquilibrium operator()(Numeric r, Numeric drdSELFVMR, Numeric drdT) &&noexcept
FullNonLocalThermodynamicEquilibriumInitialization(Numeric F0, Numeric A21, Numeric T, Numeric r1, Numeric r2, Numeric c2, Numeric c3, Numeric x) noexcept
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
Complex dFdG0(Numeric dG0) const noexcept
Complex dFdH(Numeric dZ) const noexcept
Complex dFdT(const Output &d, Numeric T) const noexcept
CalcType init(const Complex c2t) const noexcept
HartmannTran(Numeric F0_noshift, const Output &ls, Numeric GD_div_F0, Numeric dZ) noexcept
Complex dFdf() const noexcept
Complex dFdD2(Numeric dD2) const noexcept
Complex dFdD0(Numeric dD0) const noexcept
void update_calcs() noexcept
Complex dFdG2(Numeric dG2) const noexcept
Complex dFdF0() const noexcept
Complex dFdETA(Numeric dETA) const noexcept
Complex dFdFVC(Numeric dFVC) const noexcept
Complex operator()(Numeric f) noexcept
Complex dFdVMR(const Output &d) const noexcept
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
Numeric dNdT(Numeric, Numeric f) const noexcept
Numeric dNdF0() const noexcept
Numeric operator()(Numeric f) noexcept
Numeric dNdf(Numeric f) const noexcept
RosenkranzQuadratic(Numeric F0, Numeric T) noexcept
Numeric operator()(Numeric f) noexcept
Numeric dNdT(Numeric t_, Numeric f) const ARTS_NOEXCEPT
Numeric dNdf(Numeric f) const noexcept
Struct to keep the cutoff limited range values and the sparse limits.
Complex dFdD2(Numeric dD2dD2) const noexcept
Complex dFdD0(Numeric dD0dD0) const noexcept
Complex dFdH(Numeric dZ) const noexcept
void update_calcs() noexcept
Complex dFdG2(Numeric dG2dG2) const noexcept
Complex dFdT(const Output &d, Numeric T) const noexcept
Complex dFdf() const noexcept
Complex dFdF0() const noexcept
CalcType init(const Complex c2) const noexcept
Complex dFdG0(Numeric dG0dG0) const noexcept
SpeedDependentVoigt(Numeric F0_noshift, const Output &ls, Numeric GD_div_F0, Numeric dZ) noexcept
Complex dFdVMR(const Output &d) const noexcept
Complex operator()(Numeric f) noexcept
VanVleckHuber(Numeric F0, Numeric T) noexcept
Numeric operator()(Numeric f) noexcept
Numeric dNdf(Numeric f) const noexcept
Numeric dNdF0() const noexcept
Numeric dNdT(Numeric T, Numeric f) const noexcept
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
constexpr VibrationalTemperaturesNonLocalThermodynamicEquilibrium operator()(Numeric I0, Numeric QT0, Numeric QT, Numeric dQTdT, Numeric r, Numeric drdSELFVMR, Numeric drdT) &&noexcept
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
Complex operator()(Numeric f) noexcept
StateMatchType operates so that a check less than a level should be 'better', bar None.
constexpr Values() noexcept