1 #include <complex>
2
3 #include "_faddeeva.h"
4
5 using namespace std;
6
7 EXTERN_C_START
8
faddeeva_w(npy_cdouble zp)9 npy_cdouble faddeeva_w(npy_cdouble zp)
10 {
11 complex<double> z(zp.real, zp.imag);
12 std::complex<double> w = Faddeeva::w(z);
13 return npy_cpack(real(w), imag(w));
14 }
15
faddeeva_erf(npy_cdouble zp)16 npy_cdouble faddeeva_erf(npy_cdouble zp)
17 {
18 complex<double> z(zp.real, zp.imag);
19 complex<double> w = Faddeeva::erf(z);
20 return npy_cpack(real(w), imag(w));
21 }
22
faddeeva_erfc(npy_cdouble zp)23 npy_cdouble faddeeva_erfc(npy_cdouble zp)
24 {
25 complex<double> z(zp.real, zp.imag);
26 complex<double> w = Faddeeva::erfc(z);
27 return npy_cpack(real(w), imag(w));
28 }
29
faddeeva_erfcx(double x)30 double faddeeva_erfcx(double x)
31 {
32 return Faddeeva::erfcx(x);
33 }
34
faddeeva_erfcx_complex(npy_cdouble zp)35 npy_cdouble faddeeva_erfcx_complex(npy_cdouble zp)
36 {
37 complex<double> z(zp.real, zp.imag);
38 complex<double> w = Faddeeva::erfcx(z);
39 return npy_cpack(real(w), imag(w));
40 }
41
faddeeva_erfi(double x)42 double faddeeva_erfi(double x)
43 {
44 return Faddeeva::erfi(x);
45 }
46
faddeeva_erfi_complex(npy_cdouble zp)47 npy_cdouble faddeeva_erfi_complex(npy_cdouble zp)
48 {
49 complex<double> z(zp.real, zp.imag);
50 complex<double> w = Faddeeva::erfi(z);
51 return npy_cpack(real(w), imag(w));
52 }
53
faddeeva_dawsn(double x)54 double faddeeva_dawsn(double x)
55 {
56 return Faddeeva::Dawson(x);
57 }
58
faddeeva_dawsn_complex(npy_cdouble zp)59 npy_cdouble faddeeva_dawsn_complex(npy_cdouble zp)
60 {
61 complex<double> z(zp.real, zp.imag);
62 complex<double> w = Faddeeva::Dawson(z);
63 return npy_cpack(real(w), imag(w));
64 }
65
66 /*
67 * A wrapper for a normal CDF for complex argument
68 */
69
faddeeva_ndtr(npy_cdouble zp)70 npy_cdouble faddeeva_ndtr(npy_cdouble zp)
71 {
72 complex<double> z(zp.real, zp.imag);
73 z *= NPY_SQRT1_2;
74 complex<double> w = 0.5 * Faddeeva::erfc(-z);
75 return npy_cpack(real(w), imag(w));
76 }
77
78 /*
79 * Log of the normal CDF for complex arguments.
80 *
81 * This is equivalent to log(ndtr(z)), but is more robust to overflow at $z\to\infty$.
82 * This implementation uses the Faddeva computation, $\erfc(z) = \exp(-z^2) w(iz)$,
83 * taking special care to select the principal branch of the log function
84 * log( exp(-z^2) w(i z) )
85 */
faddeeva_log_ndtr(npy_cdouble zp)86 npy_cdouble faddeeva_log_ndtr(npy_cdouble zp)
87 {
88
89 complex<double> z(zp.real, zp.imag);
90 if (zp.real > 6) {
91 // Underflow. Close to the real axis, expand the log in log(1 - ndtr(-z)).
92 complex<double> w = -0.5 * Faddeeva::erfc(z*NPY_SQRT1_2);
93 if (abs(w) < 1e-8) {
94 return npy_cpack(real(w), imag(w));
95 }
96 }
97
98 z *= -NPY_SQRT1_2;
99 double x = real(z), y = imag(z);
100
101 /* Compute the principal branch of $log(exp(-z^2))$, using the fact that
102 * $log(e^t) = log|e^t| + i Arg(e^t)$, and that if $t = r + is$, then
103 * $e^t = e^r (\cos(s) + i \sin(s))$.
104 */
105 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
106 double mIm_z2 = -2*x*y; // Im(-z^2)
107
108 double im = fmod(mIm_z2, 2.0*NPY_PI);
109 if (im > NPY_PI) {im -= 2.0*NPY_PI;}
110
111 complex<double> val1 = complex<double>(mRe_z2, im);
112
113 complex<double> val2 = log(Faddeeva::w(complex<double>(-y, x)));
114 complex<double> result = val1 + val2 - NPY_LOGE2;
115
116 /* Again, select the principal branch: log(z) = log|z| + i arg(z), thus
117 * the imaginary part of the result should belong to [-pi, pi].
118 */
119 im = imag(result);
120 if (im >= NPY_PI){ im -= 2*NPY_PI; }
121 if (im < -NPY_PI){ im += 2*NPY_PI; }
122
123 return npy_cpack(real(result), im);
124 }
125
faddeeva_voigt_profile(double x,double sigma,double gamma)126 double faddeeva_voigt_profile(double x, double sigma, double gamma)
127 {
128 const double INV_SQRT_2 = 0.707106781186547524401;
129 const double SQRT_2PI = 2.5066282746310002416123552393401042;
130
131 if(sigma == 0){
132 if (gamma == 0){
133 if (std::isnan(x))
134 return x;
135 if (x == 0)
136 return NPY_INFINITY;
137 return 0;
138 }
139 return gamma / NPY_PI / (x*x + gamma*gamma);
140 }
141 if (gamma == 0){
142 return 1 / SQRT_2PI / sigma * exp(-(x/sigma)*(x/sigma) / 2);
143 }
144
145 double zreal = x / sigma * INV_SQRT_2;
146 double zimag = gamma / sigma * INV_SQRT_2;
147 std::complex<double> z(zreal, zimag);
148 std::complex<double> w = Faddeeva::w(z);
149 return real(w) / sigma / SQRT_2PI;
150 }
151
152 EXTERN_C_END
153