[SciPy-dev] some interesting routines in Python SOMNEC program

[Travis Oliphant]

> > neat stuff, and begs to be rewritten in C.  (no not f2c, that becomes
> > unintelligible :)
>
> It seems like you should just leverage the scipy.special library, but...  I just
> checked and there are a few things missing.  The first thing is bessel and
> hankel functions for complex arguments.  The Fortran functions are there in amos
> (zbesh,cbesh,etc.), but they aren't exposed from the cephes interface.  Is there
> a reason for this?

Actually, they are exposed.  Are they not working on your platform?

>>> special.jv(4.0,3+1j)
(0.098910157742492094+0.14280019502387303j)
>>> special.hankel1(4.0,3+1j)
(-0.32517852499177413-0.44754176013594482j)


>
> 2nd, there isn't a function to compute the derivative of the bessel function.
> Rob's somnec.BESSEL(z) returns both the bessel function and its derivative as
> results.  Do you need a special algorithm to compute the derivative, or is it
> possible to compute it analytically by some relationship with the bessel
> function (seems like I remember something like this)?
>
> Also, would the functions in integrate work for calculating the integrals?  This
> is a "real world" problem that SciPy should be able to handle in a few lines,
> correct?  It is just a matter of making sure we have the correct algorithms
> underneath.

Absolutely, you could use integrate.quad to compute the integrals nicely.

The derivatives of any bessel function satisfy the recurrence relation:

d C_v(z)
------- =  ( C_(v-1) (z) - C_(v+1) (z) ) / 2
  dz

          = C_(v-1) (z) - v / z * C_v(z)

          = -C_(v+1) (z) + v / z * C_v(z)


where C is any of the bessel functions or modified bessel functions or any
linear combination of them.

(From Abramowitz and Stegun pg. 361)

We could easily implement these derivatives if it is desired.

-Travis