| Re: bessel [message #17826 is a reply to message #17819] |
Fri, 05 November 1999 00:00  |
Michael Kueppers
Messages: 4
Registered: June 1999
|
Junior Member |
|
|
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
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
