!************************************************************************ !* * !* Program to calculate the second kind Bessel function of integer * !* order N, for any REAL X, using the function BESSY(N,X). * !* * !* -------------------------------------------------------------------- * !* * !* SAMPLE RUN: * !* * !* (Calculate Bessel function for N=2, X=0.75). * !* * !* Second kind Bessel Function of order 2 for X = 0.7500: * !* * !* Y = -0.26297459E+01 * !* * !* -------------------------------------------------------------------- * !* Reference: From Numath Library By Tuan Dang Trong in Fortran 77. * !* * !* F90 Release 1.0 By J-P Moreau, Paris. * !* all variables declared * !* (www.jpmoreau.fr) * !************************************************************************ PROGRAM TBESSY IMPLICIT NONE REAL*8 BESSY, X, Y INTEGER N N=2 X=0.75D0 Y = BESSY(N,X) write(*,10) N, X write(*,20) Y stop 10 format (/' Second Kind Bessel Function of order ',I2,' for X=',F8.4,':') 20 format(/' Y = ',E15.8/) STOP END FUNCTION BESSY (N,X) ! ------------------------------------------------------------------ ! This subroutine calculates the second kind Bessel Function of ! integer order N, for any real X. We use here the classical ! recursive formula. ! ------------------------------------------------------------------ IMPLICIT NONE INTEGER N, J REAL *8 X,BESSY,BESSY0,BESSY1,TOX,BY,BYM,BYP IF (N.EQ.0) THEN BESSY = BESSY0(X) RETURN ENDIF IF (N.EQ.1) THEN BESSY = BESSY1(X) RETURN ENDIF IF (X.EQ.0.) THEN BESSY = -1.E30 RETURN ENDIF TOX = 2./X BY = BESSY1(X) BYM = BESSY0(X) DO 11 J = 1,N-1 BYP = J*TOX*BY-BYM BYM = BY BY = BYP 11 CONTINUE BESSY = BY RETURN END ! --------------------------------------------------------------------------- FUNCTION BESSYP (N,X) IMPLICIT NONE INTEGER N REAL *8 X,BESSYP,BESSY IF (N.EQ.0) THEN BESSYP=-BESSY(1,X) ELSE IF(X.EQ.0.D0) THEN X=1.D-30 ELSE BESSYP=BESSY(N-1,X)-(DFLOAT(N)/X)*BESSY(N,X) ENDIF RETURN END ! --------------------------------------------------------------------------- FUNCTION BESSY0 (X) IMPLICIT NONE REAL *8 X,BESSY0,BESSJ0,FS,FR,Z,FP,FQ,XX ! --------------------------------------------------------------------- ! This subroutine calculates the Second Kind Bessel Function of ! order 0, for any real number X. The polynomial approximation by ! series of Chebyshev polynomials is used for 0