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