13#include <Faddeeva/Faddeeva.hh>
25 2.772588722239781237668928485832706272302000537441021016482720037973574487879;
29 x = (f - mF0) * invGD;
30 F = invGD * inv_sqrt_pi * std::exp(-pow2(x));
36 F = inv_sqrt_pi * invGD * Faddeeva::w(z);
37 dF = 2 * invGD * (
Complex(0, inv_pi * invGD) - z * F);
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 *
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 *
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:
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:
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) +
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) +
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:
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) -
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) -
301 case CalcType::LowYandHighX:
302 return -pow2(invc2) *
304 (x - 3) * (
Complex(2 *
d.G0 - 3 *
d.G2, 2 *
d.D0 - 3 *
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) +
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) -
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:
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) -
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) -
397 case CalcType::LowYandHighX:
398 return -pow2(invc2) *
400 (x - 3) * (
Complex(2 *
d.G0 - 3 *
d.G2, 2 *
d.D0 - 3 *
d.D2) -
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));
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 = -1i;
495 ((pow2(
z1) - 1) * 1i *
dw1 * dz1 - (pow2(
z2) - 1) * 1i *
dw2 * dz2 +
496 2 *
w1 *
z1 * dz1 - 2 *
w2 *
z2 * dz2) /
500 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
507 (sqrt_pi * ((pow2(
z1) - 1) * 1i *
dw1 * dz1 + 2 *
w1 *
z1 * dz1) - dz1);
510 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
517 9 *
invGD * dz1 / (4 * pow4(
z1));
520 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
528 9 *
invGD * dz1 / (4 * pow4(
z1));
531 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
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) /
546 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
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) /
557 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
568 constexpr Complex ddeltax = 1i;
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) *
588 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
594 -(sqrt_pi * ((pow2(
z1) - 1) * 1i *
dw1 * dz1 + 2 *
w1 *
z1 * dz1) -
597 (sqrt_pi * (pow2(
z1) - 1) *
w1 -
z1) * dinvGD;
600 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
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 +
613 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
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 +
629 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
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)) /
644 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
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) /
656 return inv_pi * (-
A * dK +
K * dA) / pow2(
K);
664 const Numeric dGD = (dmF0 / (invGD * mF0));
665 const Numeric dinvGD = -dGD * pow2(invGD);
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) *
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;
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 +
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 +
729 return inv_pi * (-A * dK + K * dA) / pow2(K);
731 case CalcType::LowYandLowX: {
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)) /
744 return inv_pi * (-A * dK + K * dA) / pow2(K);
746 case CalcType::LowYandHighX: {
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) /
756 return inv_pi * (-A * dK + K * dA) / pow2(K);
765 const Complex dsqrtxy = dx / (2 * sqrtxy);
768 case CalcType::Full: {
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) /
779 return inv_pi * (-A * dK + K * dA) / pow2(K);
781 case CalcType::Noc2tLowZ: {
782 const Complex dz1 = invGD * ddeltax;
786 (sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) - dz1);
789 return inv_pi * (-A * dK + K * dA) / pow2(K);
791 case CalcType::Noc2tHighZ: {
792 const Complex dz1 = invGD * ddeltax;
795 invGD * dz1 / (2 * pow2(z1)) +
796 9 * invGD * dz1 / (4 * pow4(z1));
799 return inv_pi * (-A * dK + K * dA) / pow2(K);
801 case CalcType::LowXandHighY: {
802 const Complex dz1 = invGD * ddeltax;
804 const Complex dA =
Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
806 invGD * dz1 / (2 * pow2(z1)) +
807 9 * invGD * dz1 / (4 * pow4(z1));
810 return inv_pi * (-A * dK + K * dA) / pow2(K);
812 case CalcType::LowYandLowX: {
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) /
825 return inv_pi * (-A * dK + K * dA) / pow2(K);
827 case CalcType::LowYandHighX: {
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) /
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);
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)));
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 +
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;
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 +
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 +
923 return inv_pi * (-A * dK + K * dA) / pow2(K);
925 case CalcType::LowYandLowX: {
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) /
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) *
948 return inv_pi * (-A * dK + K * dA) / pow2(K);
950 case CalcType::LowYandHighX: {
953 (
Complex(
G2,
D2) * (x - 3) * dx + 1i * (2 * x - 3) * x * dD2 / 2) /
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)) *
969 return inv_pi * (-A * dK + K * dA) / pow2(K);
977 const Numeric ddeltax = 3 * (
ETA - 1) * dG2 / 2;
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);
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 +
1001 return inv_pi * (-A * dK + K * dA) / pow2(K);
1003 case CalcType::Noc2tLowZ: {
1004 const Complex dz1 = invGD * ddeltax;
1008 (sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) - dz1);
1012 return inv_pi * (-A * dK + K * dA) / pow2(K);
1014 case CalcType::Noc2tHighZ: {
1015 const Complex dz1 = invGD * ddeltax;
1018 invGD * dz1 / (2 * pow2(z1)) +
1019 9 * invGD * dz1 / (4 * pow4(z1));
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);
1030 invGD * dz1 / (2 * pow2(z1)) +
1031 9 * invGD * dz1 / (4 * pow4(z1));
1035 return inv_pi * (-A * dK + K * dA) / pow2(K);
1037 case CalcType::LowYandLowX: {
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) /
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) *
1058 return inv_pi * (-A * dK + K * dA) / pow2(K);
1060 case CalcType::LowYandHighX: {
1063 (
Complex(
G2,
D2) * (x - 3) * dx + (2 * x - 3) * x * dG2 / 2) /
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)) *
1077 return inv_pi * (-A * dK + K * dA) / pow2(K);
1086 const Complex dsqrtxy = dx / (2 * sqrtxy);
1089 case CalcType::Full: {
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) /
1101 return inv_pi * (-A * dK + K * dA) / pow2(K);
1103 case CalcType::Noc2tLowZ: {
1104 const Complex dz1 = invGD * ddeltax;
1108 (sqrt_pi * ((pow2(z1) - 1) * 1i * dw1 * dz1 + 2 * w1 * z1 * dz1) - dz1);
1112 return inv_pi * (-A * dK + K * dA) / pow2(K);
1114 case CalcType::Noc2tHighZ: {
1115 const Complex dz1 = invGD * ddeltax;
1118 invGD * dz1 / (2 * pow2(z1)) +
1119 9 * invGD * dz1 / (4 * pow4(z1));
1123 return inv_pi * (-A * dK + K * dA) / pow2(K);
1125 case CalcType::LowXandHighY: {
1126 const Complex dz1 = invGD * ddeltax;
1128 const Complex dA =
Complex(0, sqrt_pi * invGD) * (dw1 * dz1 - dw2 * dz2);
1130 invGD * dz1 / (2 * pow2(z1)) +
1131 9 * invGD * dz1 / (4 * pow4(z1));
1135 return inv_pi * (-A * dK + K * dA) / pow2(K);
1137 case CalcType::LowYandLowX: {
1139 const Complex dz2 = dx / (2 * sqrtx);
1140 const Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
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) /
1151 return inv_pi * (-A * dK + K * dA) / pow2(K);
1153 case CalcType::LowYandHighX: {
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) /
1163 return inv_pi * (-A * dK + K * dA) / pow2(K);
1171 const Numeric dGD = (dmF0 / (invGD * mF0));
1172 const Numeric dinvGD = -dGD * pow2(invGD);
1173 const Complex dsqrty = ((
ETA - 1) * dinvGD + invGD * dETA) /
1176 const Complex dx = (-(
ETA - 1) * ddeltax + deltax * dETA) /
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 +
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;
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 +
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 +
1250 return inv_pi * (-A * dK + K * dA) / pow2(K);
1252 case CalcType::LowYandLowX: {
1254 const Complex dz2 = dx / (2 * sqrtx);
1256 (sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) * (
ETA - 1) -
1257 (sqrt_pi * w2 * z2 - 1) * dETA) /
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) *
1273 return inv_pi * (-A * dK + K * dA) / pow2(K);
1275 case CalcType::LowYandHighX: {
1277 const Complex dA = ((
ETA - 1) * (x - 3) * dx + (2 * x - 3) * x * dETA / 2) /
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)) /
1292 return inv_pi * (-A * dK + K * dA) / pow2(K);
1300 const Numeric dGD = (dmF0 / (invGD * mF0));
1301 const Numeric dinvGD = -dGD * pow2(invGD);
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) *
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;
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 +
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 +
1365 return inv_pi * (-A * dK + K * dA) / pow2(K);
1367 case CalcType::LowYandLowX: {
1369 const Complex dz2 = dx / (2 * sqrtx);
1370 const Complex dA = 2 * sqrt_pi * (w2 * dz2 + z2 * 1i * dw2 * dz2) /
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)) /
1380 return inv_pi * (-A * dK + K * dA) / pow2(K);
1382 case CalcType::LowYandHighX: {
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) /
1392 return inv_pi * (-A * dK + K * dA) / pow2(K);
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);
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)));
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 +
1436 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
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;
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 +
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 +
1503 return inv_pi * (-A * dK + K * dA) / pow2(K);
1505 case CalcType::LowYandLowX: {
1507 const Complex dz2 = dx / (2 * sqrtx);
1510 (sqrt_pi *
Complex(
G2,
D2) * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1513 Complex(
d.G2,
d.D2) * (sqrt_pi * w2 * z2 - 1) * (
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) *
1526 (2 * sqrt_pi * w1 * z1 +
1527 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1)) /
1536 return inv_pi * (-A * dK + K * dA) / pow2(K);
1538 case CalcType::LowYandHighX: {
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) -
1554 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1555 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1565 return inv_pi * (-A * dK + K * dA) / pow2(K);
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);
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)));
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 +
1610 (sqrt_pi * (pow2(z1) - 1) * w1 - sqrt_pi * (pow2(z2) - 1) * w2 +
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;
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 +
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 +
1677 return inv_pi * (-A * dK + K * dA) / pow2(K);
1679 case CalcType::LowYandLowX: {
1681 const Complex dz2 = dx / (2 * sqrtx);
1684 (sqrt_pi *
Complex(
G2,
D2) * (w2 * dz2 + z2 * 1i * dw2 * dz2) *
1687 Complex(
d.G2,
d.D2) * (sqrt_pi * w2 * z2 - 1) * (
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) *
1700 (2 * sqrt_pi * w1 * z1 +
1701 2 * (sqrt_pi * w2 * z2 - 1) * (2 * pow2(sqrty) + x - 1) - 1)) /
1710 return inv_pi * (-A * dK + K * dA) / pow2(K);
1712 case CalcType::LowYandHighX: {
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) -
1728 (2 * (-2 * sqrt_pi * w1 * z1 + 1) * pow2(x) +
1729 (2 * x - 3) * (2 * pow2(sqrty) + x - 1)) *
1739 return inv_pi * (-A * dK + K * dA) / pow2(K);
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) /
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);
1814 2 * sqrt_pi *
z1 *
w1);
1818 w1 = Faddeeva::w(1i *
z1);
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)) {}
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")
1925 I0, r, drdSELFVMR, drdT, QT0, QT, dQTdT,
boltzman_ratio(T, T0, E0),
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),
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,
1997 return {I0, QT0, QT, dQTdT, r, drdSELFVMR, drdT,
K1,
2010 T, T0, F0, E0, Tl, Evl, Tu, Evu,
stimulated_emission(T, F0),
2034 const Numeric *it0 = f_grid.get_c_array();
2035 const Numeric *itn = it0 + f_grid.size();
2036 const Numeric *itl = std::lower_bound(it0, itn, fl);
2038 std::distance(itl, std::upper_bound(itl, itn, fu))};
2055 const Index nvs = sparse_f_grid.size();
2058 const Numeric *it0s = sparse_f_grid.get_c_array();
2059 const Numeric *itns = it0s + nvs;
2061 std::lower_bound(it0s, itns, std::nextafter(flc, fuc));
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.get_c_array();
2093 const Numeric *itn = it0 + f_grid.size();
2094 const Numeric *itl = std::lower_bound(it0, itn, fl);
2096 std::upper_bound(itl, itn, std::nextafter(fu, fl));
2098 return {std::distance(it0, itl),
2099 std::distance(itl, itu),
2113 const Index nvs = sparse_f_grid.size();
2114 const Index nv = f_grid.size();
2117 const Numeric *
const it0s = sparse_f_grid.get_c_array();
2118 const Numeric *
const itns = it0s + nvs;
2120 std::lower_bound(it0s, itns, std::nextafter(flc, fuc));
2122 std::upper_bound(itlc, itns, std::nextafter(fuc, flc));
2124 std::upper_bound(itlc, ituc, std::nextafter(fls, flc));
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.get_c_array();
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; }));
2224 const bool do_nlte_) noexcept
2225 :
F(F_.get_c_array() + start),
2226 dF(dF_.get_c_array() + start * derivs_.size()),
2227 N(N_.get_c_array() + start),
2228 dN(dN_.get_c_array() + start * derivs_.size()),
2229 f(f_grid.get_c_array() + start),
size(nv),
derivs(derivs_),
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),
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);
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:
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);
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;
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__))
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++) {
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;
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;
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;
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++) {
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);
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;
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;
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;
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;
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);
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());
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);
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());
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);
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());
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: {
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];
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;
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;
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,...)
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
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
NUMERIC Numeric
The type to use for all floating point numbers.
INDEX Index
The type to use for all integer numbers and indices.
Numeric & real_val(Complex &c) noexcept
Return a reference to the real value of c.
Numeric & imag_val(Complex &c) noexcept
Return a reference to the imaginary value of c.
std::complex< Numeric > Complex
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