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