On Fri, 19 Jun 2009, Klein, Christoph wrote:
 thank you *so* much!!! this function is much more accurate than
 my version. 
I can't claim credit for the revised version; it comes from the
netlib "specfun" repository.
 Additionaly, I propose to dump my bessel_I function, too. With
 the attached patch bessel_I for negative orders is calculated
 using the new besselK function. bessel_I for negative orders
 can't be calculated with the cephes function, because cephes
 assumes symmetry around 0 for the order (see code), which is
 false. 
Thanks.  I've applied your patch.  I've also added (minimal)
docmentation for the bessel function.  One little refinement: you
don't need to quote the first argument.  For example,
x = bessel(J, 2, 5.6)
will do OK.
Allin.