Re: bessel

[Michael Kueppers <michael.kueppers(at)phim.unibe.ch>]

enea wrote:

> I have to calculate the modified Besell functions K(y).
> I 'm not able to do it in idl.
> Someone can help me?
>
> Excuse me for my bad english
>
> Claudia

     The  IDL-functions below are the  Bessel-functions
      K_0(y)  and K_1(y) taken from "Numerical Recipes in C"
      (Press et al. 1992, Cambridge Univ. Press) and
       translated to the
       Interactive Data Language.  Should your question refer
       to the other idl (I am sufficiently ignorant not to know if  this
       is a possibility), please apologize for bothering.
        You can construct higher order bessel functions by

        -2n / x  * K_n(x)  =  K_(n-1) (x)   - K_(n+1) (x)


        Best wishes,
         Michael



FUNCTION beselk0,X,z
; Compute the modified Bessel-function of second kind and zeroth order
;                                               M.K., 1.5.97
; Completely changed: now taken from "Numerical recipes in C", Press et al.,
; 1992, and translated from C to IDL.
;                                M.K., 2.5.97
; Corection in case that X is a vector. Also, dummy variable z added to allow
function for INT_2d
;                                M.K., 14.7.99


ans = X
FOR I =0,(N_ELEMENTS(X)-1) DO BEGIN
If X(I) LE 2. THEN BEGIN
     Y = X(I)*X(I)/4.
     ans(I) = (-alog(X(I)/2.)*BESELI(X(I),0)) + (-0.57721566 + Y*(0.42278420 +$
           Y*(0.23069756 + Y*(0.03488590 + Y*(0.262698e-2 +$
          Y*(0.10750e-3 + Y*0.74e-5))))))
ENDIF ELSE BEGIN
     Y = 2.0/X(I)
     ans(I) = (EXP(-X(I))/SQRT(X(I)))*(1.25331414 + Y*(-0.07832358 +$
           Y*(0.02189568 + Y*(-0.01062446 + Y*(0.587872e-2 +$
           Y*(-0.251540e-2 + Y*0.53208E-3))))))
ENDELSE
ENDFOR

RETURN, ans

END



 FUNCTION beselk1,X
; Compute the modified Bessel-function of second kind and first order
; Taken from "Numerical recipes in C", Press et al.,
; 1992, and translated from C to IDL.
;                                M.K., 14.7.99

ans = X
FOR I =0,(N_ELEMENTS(X)-1) DO BEGIN
If X(I) LE 2. THEN BEGIN
     Y = X(I)*X(I)/4.
     ans(I) = (alog(X(I)/2.)*BESELI(X(I),1))  + $
(1./X(I))*(1.+Y*(0.15443144+Y*(-0.67278579+Y*(-0.18156897+Y*(-0.01919402+Y*$
        (-0.00110404+Y*(4.686e-5)))))))
ENDIF ELSE BEGIN
     Y = 2.0/X(I)
     ans(I) = (EXP(-X(I))/SQRT(X(I)))*(1.25331414 +
Y*(0.23498619+Y*(-0.03655620+Y*(0.01504268+Y*(-0.00780353+y*(0.00325614+Y*$
        (-0.00068245)))))))
ENDELSE
ENDFOR

RETURN, ans

END

--
Michael Kueppers          email: kueppers@phim.unibe.ch
Physikalisches Institut     Tel: [41] - (31) - 631 4419
Universitaet Bern           FAX: [41] - (31) - 631 4405
CH-3012 Bern, Switzerland