1import ctypes
2from libc.math cimport sqrt, fabs
3from libc.stdlib cimport free
4from numpy import nan
5
6cdef extern from "Python.h":
7    object PyCapsule_New(void *pointer, char *name, void *destructor)
8
9from scipy._lib._ccallback import LowLevelCallable
10from ._ellip_harm cimport ellip_harmonic, ellip_harm_eval, lame_coefficients
11
12ctypedef struct _ellip_data_t:
13    double *eval
14    double h2, k2
15    int n, p
16
17cdef double _F_integrand(double t, void *user_data) nogil:
18    cdef _ellip_data_t *data = <_ellip_data_t *>user_data
19    cdef double h2, k2, t2, i, a, result
20    cdef int n, p
21    cdef double * eval
22    t2 = t*t
23    h2 = data[0].h2
24    k2 = data[0].k2
25    n = data[0].n
26    p = data[0].p
27    eval = data[0].eval
28    i = ellip_harm_eval( h2, k2, n, p, 1/t, eval, 1, 1)
29    result = 1/(i*i*sqrt(1 - t2*k2)*sqrt(1 - t2*h2))
30    return result
31
32cdef double _F_integrand1(double t, void *user_data) nogil:
33    cdef _ellip_data_t *data = <_ellip_data_t *>user_data
34    cdef double h2, k2, t2, i, a, h, result
35    cdef int n, p
36    cdef double * eval
37    t2 = t*t
38    h2 = data[0].h2
39    k2 =data[0].k2
40    n = data[0].n
41    p = data[0].p
42    eval = data[0].eval
43
44    h = sqrt(h2)
45    k = sqrt(k2)
46    i = ellip_harm_eval( h2, k2, n, p, t, eval, 1, 1)
47    result = i*i/sqrt((t + h)*(t + k))
48    return result
49
50cdef double _F_integrand2(double t, void *user_data) nogil:
51    cdef _ellip_data_t *data = <_ellip_data_t *>user_data
52    cdef double h2, k2, t2, i, a, h, result
53    cdef int n, p
54    cdef double * eval
55    t2 = t*t
56    h2 = data[0].h2
57    k2 =data[0].k2
58    n = data[0].n
59    p = data[0].p
60    eval = data[0].eval
61
62    h = sqrt(h2)
63    k = sqrt(k2)
64    i = ellip_harm_eval( h2, k2, n, p, t, eval, 1, 1)
65    result = t2*i*i/sqrt((t + h)*(t + k))
66    return result
67
68cdef double _F_integrand3(double t, void *user_data) nogil:
69    cdef _ellip_data_t *data = <_ellip_data_t *>user_data
70    cdef double h2, k2, t2, i, a, h, result
71    cdef int n, p
72    cdef double * eval
73    t2 = t*t
74    h2 = data[0].h2
75    k2 =data[0].k2
76    n = data[0].n
77    p = data[0].p
78    eval = data[0].eval
79
80    h = sqrt(h2)
81    k = sqrt(k2)
82    i = ellip_harm_eval( h2, k2, n, p, t, eval, 1, 1)
83    result = i*i/sqrt((t + h)*(k2 - t2))
84    return result
85
86cdef double _F_integrand4(double t, void *user_data) nogil:
87    cdef _ellip_data_t *data = <_ellip_data_t *>user_data
88    cdef double h2, k2, t2, i, a, h, result
89    cdef int n, p
90    cdef double *eval
91    t2 = t*t
92    h2 = data[0].h2
93    k2 =data[0].k2
94    n = data[0].n
95    p = data[0].p
96    eval = data[0].eval
97
98    h = sqrt(h2)
99    k = sqrt(k2)
100    i = ellip_harm_eval( h2, k2, n, p, t, eval, 1, 1)
101    result = i*i*t2/sqrt((t + h)*(k2 - t2))
102    return result
103
104
105def _ellipsoid(double h2, double k2, int n, int p, double s):
106    import scipy.special._ellip_harm_2 as mod
107    from scipy.integrate import quad
108
109    cdef _ellip_data_t data
110
111    cdef double * eval
112    cdef void *bufferp
113    eval = lame_coefficients(h2, k2, n, p, &bufferp, 1, 1)
114    if not eval:
115        return nan
116
117    data.h2 = h2
118    data.k2 = k2
119    data.n = n
120    data.p = p
121    data.eval = eval
122
123    cdef double res, err
124
125    try:
126        capsule = PyCapsule_New(<void*>&data, NULL, NULL)
127        res, err = quad(LowLevelCallable.from_cython(mod, "_F_integrand", capsule), 0, 1/s,
128                                                     epsabs=1e-300, epsrel=1e-15)
129    finally:
130        free(bufferp)
131    if err > 1e-10*fabs(res) + 1e-290:
132        return nan
133    res = res*(2*n + 1)*ellip_harmonic( h2, k2, n, p, s, 1, 1)
134    return res
135
136
137def _ellipsoid_norm(double h2, double k2, int n, int p):
138    import scipy.special._ellip_harm_2 as mod
139    from scipy.integrate import quad
140
141    cdef _ellip_data_t data
142
143    cdef double *eigv
144    cdef void *bufferp
145    eval = lame_coefficients(h2, k2, n, p, &bufferp, 1, 1)
146    if not eval:
147        return nan
148
149    data.h2 = h2
150    data.k2 = k2
151    data.n = n
152    data.p = p
153    data.eval = eval
154
155    cdef double res, res1, res2, res3, err, err1, err2, err3
156
157    h = sqrt(h2)
158    k = sqrt(k2)
159    try:
160        capsule = PyCapsule_New(<void*>&data, NULL, NULL)
161
162        wvar = (-0.5, -0.5)
163
164        res, err = quad(LowLevelCallable.from_cython(mod, "_F_integrand1", capsule), h, k,
165                        epsabs=1e-300, epsrel=1e-15, weight="alg", wvar=wvar)
166
167        res1, err1 = quad(LowLevelCallable.from_cython(mod, "_F_integrand2", capsule), h, k,
168                          epsabs=1e-300, epsrel=1e-15, weight="alg", wvar=wvar)
169
170        wvar = (0, -0.5)
171
172        res2, err2 = quad(LowLevelCallable.from_cython(mod, "_F_integrand3", capsule), 0, h,
173                          epsabs=1e-300, epsrel=1e-15, weight="alg", wvar=wvar)
174
175        res3, err3 = quad(LowLevelCallable.from_cython(mod, "_F_integrand4", capsule), 0, h,
176                          epsabs=1e-300, epsrel=1e-15, weight="alg", wvar=wvar)
177
178    finally:
179        free(bufferp)
180
181    error = 8*(res2*err1 + err2*res1 + res*err3 + res3*err)
182    result = 8*(res1*res2 - res*res3)
183
184    if error > 10e-8*fabs(result):
185        return nan
186    return result
187
188
189# Needed for the _sf_error calls in _ellip_harm.pxd
190
191cimport numpy as np
192
193np.import_array()
194np.import_ufunc()
195
196cdef extern from "numpy/ufuncobject.h":
197    int PyUFunc_getfperr() nogil
198
199cdef public int wrap_PyUFunc_getfperr() nogil:
200    """
201    Call PyUFunc_getfperr in a context where PyUFunc_API array is initialized;
202    this avoids messing with the UNIQUE_SYMBOL #defines
203    """
204    return PyUFunc_getfperr()
205