1# -*-cython-*- 2# 3# Implementation of spherical Bessel functions and modified spherical Bessel 4# functions of the first and second kinds, as well as their derivatives. 5# 6# Author: Tadeusz Pudlik 7# 8# Distributed under the same license as SciPy. 9# 10# I attempt to correctly handle the edge cases (0 and infinity), but this is 11# tricky: the values of the functions often depend on the direction in which 12# the limit is taken. At zero, I follow the convention of numpy (1.9.2), 13# which treats zero differently depending on its type: 14# 15# >>> np.cos(0)/0 16# inf 17# >>> np.cos(0+0j)/(0+0j) 18# inf + nan*j 19# 20# So, real zero is assumed to be "positive zero", while complex zero has an 21# unspecified sign and produces nans. Similarly, complex infinity is taken to 22# represent the "point at infinity", an ambiguity which for some functions 23# makes `nan` the correct return value. 24 25import cython 26from libc.math cimport cos, sin, sqrt, M_PI_2 27 28from numpy cimport npy_cdouble 29from ._complexstuff cimport * 30 31from . cimport sf_error 32 33cdef extern from "amos_wrappers.h": 34 npy_cdouble cbesi_wrap( double v, npy_cdouble z) nogil 35 npy_cdouble cbesj_wrap(double v, npy_cdouble z) nogil 36 double cbesj_wrap_real(double v, double x) nogil 37 npy_cdouble cbesk_wrap(double v, npy_cdouble z) nogil 38 double cbesk_wrap_real(double v, double x) nogil 39 npy_cdouble cbesy_wrap(double v, npy_cdouble z) nogil 40 double cbesy_wrap_real(double v, double x) nogil 41 42from ._cephes cimport iv 43 44# Fused type wrappers 45 46cdef inline number_t cbesj(double v, number_t z) nogil: 47 cdef npy_cdouble r 48 if number_t is double: 49 return cbesj_wrap_real(v, z) 50 else: 51 r = cbesj_wrap(v, npy_cdouble_from_double_complex(z)) 52 return double_complex_from_npy_cdouble(r) 53 54cdef inline number_t cbesy(double v, number_t z) nogil: 55 cdef npy_cdouble r 56 if number_t is double: 57 return cbesy_wrap_real(v, z) 58 else: 59 r = cbesy_wrap(v, npy_cdouble_from_double_complex(z)) 60 return double_complex_from_npy_cdouble(r) 61 62cdef inline number_t cbesk(double v, number_t z) nogil: 63 cdef npy_cdouble r 64 if number_t is double: 65 return cbesk_wrap_real(v, z) 66 else: 67 r = cbesk_wrap(v, npy_cdouble_from_double_complex(z)) 68 return double_complex_from_npy_cdouble(r) 69 70 71# Spherical Bessel functions 72 73@cython.cdivision(True) 74cdef inline double spherical_jn_real(long n, double x) nogil: 75 cdef double s0, s1, sn 76 cdef int idx 77 78 if npy_isnan(x): 79 return x 80 if n < 0: 81 sf_error.error("spherical_jn", sf_error.DOMAIN, NULL) 82 return nan 83 if x == inf or x == -inf: 84 return 0 85 if x == 0: 86 if n == 0: 87 return 1 88 else: 89 return 0 90 91 if n > 0 and n >= x: 92 return sqrt(M_PI_2/x)*cbesj(n + 0.5, x) 93 94 s0 = sin(x)/x 95 if n == 0: 96 return s0 97 s1 = (s0 - cos(x))/x 98 if n == 1: 99 return s1 100 101 for idx in range(n - 1): 102 sn = (2*idx + 3)*s1/x - s0 103 s0 = s1 104 s1 = sn 105 if npy_isinf(sn): 106 # Overflow occurred already: terminate recurrence. 107 return sn 108 109 return sn 110 111 112@cython.cdivision(True) 113cdef inline double complex spherical_jn_complex(long n, double complex z) nogil: 114 cdef double complex out 115 if zisnan(z): 116 return z 117 if n < 0: 118 sf_error.error("spherical_jn", sf_error.DOMAIN, NULL) 119 return nan 120 if z.real == inf or z.real == -inf: 121 # https://dlmf.nist.gov/10.52.E3 122 if z.imag == 0: 123 return 0 124 else: 125 return (1+1j)*inf 126 if z.real == 0 and z.imag == 0: 127 if n == 0: 128 return 1 129 else: 130 return 0 131 132 out = zsqrt(M_PI_2/z)*cbesj(n + 0.5, z) 133 134 if z.imag == 0: 135 # Small imaginary part is spurious 136 return out.real 137 else: 138 return out 139 140 141@cython.cdivision(True) 142cdef inline double spherical_yn_real(long n, double x) nogil: 143 cdef double s0, s1, sn 144 cdef int idx 145 146 if npy_isnan(x): 147 return x 148 if n < 0: 149 sf_error.error("spherical_yn", sf_error.DOMAIN, NULL) 150 return nan 151 if x < 0: 152 return (-1)**(n+1)*spherical_yn_real(n, -x) 153 if x == inf or x == -inf: 154 return 0 155 if x == 0: 156 return -inf 157 158 s0 = -cos(x)/x 159 if n == 0: 160 return s0 161 s1 = (s0 - sin(x))/x 162 if n == 1: 163 return s1 164 165 for idx in range(n - 1): 166 sn = (2*idx + 3)*s1/x - s0 167 s0 = s1 168 s1 = sn 169 if npy_isinf(sn): 170 # Overflow occurred already: terminate recurrence. 171 return sn 172 173 return sn 174 175 176@cython.cdivision(True) 177cdef inline double complex spherical_yn_complex(long n, double complex z) nogil: 178 179 if zisnan(z): 180 return z 181 if n < 0: 182 sf_error.error("spherical_yn", sf_error.DOMAIN, NULL) 183 return nan 184 if z.real == 0 and z.imag == 0: 185 # https://dlmf.nist.gov/10.52.E2 186 return nan 187 if z.real == inf or z.real == -inf: 188 # https://dlmf.nist.gov/10.52.E3 189 if z.imag == 0: 190 return 0 191 else: 192 return (1+1j)*inf 193 194 return zsqrt(M_PI_2/z)*cbesy(n + 0.5, z) 195 196 197@cython.cdivision(True) 198cdef inline double spherical_in_real(long n, double z) nogil: 199 200 if npy_isnan(z): 201 return z 202 if n < 0: 203 sf_error.error("spherical_in", sf_error.DOMAIN, NULL) 204 return nan 205 if z == 0: 206 # https://dlmf.nist.gov/10.52.E1 207 if n == 0: 208 return 1 209 else: 210 return 0 211 if npy_isinf(z): 212 # https://dlmf.nist.gov/10.49.E8 213 if z == -inf: 214 return (-1)**n*inf 215 else: 216 return inf 217 218 return sqrt(M_PI_2/z)*iv(n + 0.5, z) 219 220 221@cython.cdivision(True) 222cdef inline double complex spherical_in_complex(long n, double complex z) nogil: 223 cdef npy_cdouble s 224 225 if zisnan(z): 226 return z 227 if n < 0: 228 sf_error.error("spherical_in", sf_error.DOMAIN, NULL) 229 return nan 230 if zabs(z) == 0: 231 # https://dlmf.nist.gov/10.52.E1 232 if n == 0: 233 return 1 234 else: 235 return 0 236 if zisinf(z): 237 # https://dlmf.nist.gov/10.52.E5 238 if z.imag == 0: 239 if z.real == -inf: 240 return (-1)**n*inf 241 else: 242 return inf 243 else: 244 return nan 245 246 s = cbesi_wrap(n + 0.5, npy_cdouble_from_double_complex(z)) 247 return zsqrt(M_PI_2/z)*double_complex_from_npy_cdouble(s) 248 249 250@cython.cdivision(True) 251cdef inline double spherical_kn_real(long n, double z) nogil: 252 253 if npy_isnan(z): 254 return z 255 if n < 0: 256 sf_error.error("spherical_kn", sf_error.DOMAIN, NULL) 257 return nan 258 if z == 0: 259 return inf 260 if npy_isinf(z): 261 # https://dlmf.nist.gov/10.52.E6 262 if z == inf: 263 return 0 264 else: 265 return -inf 266 267 return sqrt(M_PI_2/z)*cbesk(n + 0.5, z) 268 269 270@cython.cdivision(True) 271cdef inline double complex spherical_kn_complex(long n, double complex z) nogil: 272 273 if zisnan(z): 274 return z 275 if n < 0: 276 sf_error.error("spherical_kn", sf_error.DOMAIN, NULL) 277 return nan 278 if zabs(z) == 0: 279 return nan 280 if zisinf(z): 281 # https://dlmf.nist.gov/10.52.E6 282 if z.imag == 0: 283 if z.real == inf: 284 return 0 285 else: 286 return -inf 287 else: 288 return nan 289 290 return zsqrt(M_PI_2/z)*cbesk(n + 0.5, z) 291 292 293# Derivatives 294 295@cython.cdivision(True) 296cdef inline double spherical_jn_d_real(long n, double x) nogil: 297 if n == 0: 298 return -spherical_jn_real(1, x) 299 else: 300 if x == 0: 301 # DLMF 10.51.2 doesn't work, so use 10.51.1 to get the 302 # exact value 303 if n == 1: 304 return 1.0/3 305 else: 306 return 0 307 # DLMF 10.51.2 308 return (spherical_jn_real(n - 1, x) - 309 (n + 1)*spherical_jn_real(n, x)/x) 310 311 312@cython.cdivision(True) 313cdef inline double complex spherical_jn_d_complex(long n, double complex x) nogil: 314 if n == 0: 315 return -spherical_jn_complex(1, x) 316 else: 317 return (spherical_jn_complex(n - 1, x) - 318 (n + 1)*spherical_jn_complex(n, x)/x) 319 320 321@cython.cdivision(True) 322cdef inline double spherical_yn_d_real(long n, double x) nogil: 323 if n == 0: 324 return -spherical_yn_real(1, x) 325 else: 326 return (spherical_yn_real(n - 1, x) - 327 (n + 1)*spherical_yn_real(n, x)/x) 328 329 330@cython.cdivision(True) 331cdef inline double complex spherical_yn_d_complex(long n, double complex x) nogil: 332 if n == 0: 333 return -spherical_yn_complex(1, x) 334 else: 335 return (spherical_yn_complex(n - 1, x) - 336 (n + 1)*spherical_yn_complex(n, x)/x) 337 338 339@cython.cdivision(True) 340cdef inline double spherical_in_d_real(long n, double x) nogil: 341 if n == 0: 342 return spherical_in_real(1, x) 343 else: 344 if x == 0: 345 return 0 346 return (spherical_in_real(n - 1, x) - 347 (n + 1)*spherical_in_real(n, x)/x) 348 349 350@cython.cdivision(True) 351cdef inline double complex spherical_in_d_complex(long n, double complex x) nogil: 352 if n == 0: 353 return spherical_in_complex(1, x) 354 else: 355 if x == 0: 356 return 0 357 return (spherical_in_complex(n - 1, x) - 358 (n + 1)*spherical_in_complex(n, x)/x) 359 360 361@cython.cdivision(True) 362cdef inline double spherical_kn_d_real(long n, double x) nogil: 363 if n == 0: 364 return -spherical_kn_real(1, x) 365 else: 366 return (-spherical_kn_real(n - 1, x) - 367 (n + 1)*spherical_kn_real(n, x)/x) 368 369 370@cython.cdivision(True) 371cdef inline double complex spherical_kn_d_complex(long n, double complex x) nogil: 372 if n == 0: 373 return -spherical_kn_complex(1, x) 374 else: 375 return (-spherical_kn_complex(n - 1, x) - 376 (n + 1)*spherical_kn_complex(n, x)/x) 377