!************************************************************************ !* * !* Program to calculate the first kind Bessel function of integer * !* order N, for any REAL X, using the function BESSJ(N,X). * !* * !* -------------------------------------------------------------------- * !* * !* SAMPLE RUN: * !* * !* (Calculate Bessel function for N=2, X=0.75). * !* * !* Bessel function of order 2 for X = 0.7500: * !* * !* Y = 0.67073997E-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 TBESSJ IMPLICIT NONE REAL*8 BESSJ, X, Y INTEGER N N=2 X=0.75D0 Y = BESSJ(N,X) write(*,10) N, X write(*,20) Y stop 10 format (/' Bessel Function of order ',I2,' for X=',F8.4,':') 20 format(/' Y = ',E15.8/) STOP END FUNCTION BESSJ (N,X) ! This subroutine calculates the first kind modified Bessel function ! of integer order N, for any REAL X. We use here the classical ! recursion formula, when X > N. For X < N, the Miller's algorithm ! is used to avoid overflows. ! REFERENCE: ! C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS, ! MATHEMATICAL TABLES, VOL.5, 1962. IMPLICIT NONE INTEGER, PARAMETER :: IACC = 40 REAL*8, PARAMETER :: BIGNO = 1.D10, BIGNI = 1.D-10 INTEGER M, N, J, JSUM REAL *8 X,BESSJ,BESSJ0,BESSJ1,TOX,BJM,BJ,BJP,SUM IF (N.EQ.0) THEN BESSJ = BESSJ0(X) RETURN ENDIF IF (N.EQ.1) THEN BESSJ = BESSJ1(X) RETURN ENDIF IF (X.EQ.0.) THEN BESSJ = 0. RETURN ENDIF TOX = 2./X IF (X.GT.FLOAT(N)) THEN BJM = BESSJ0(X) BJ = BESSJ1(X) DO 11 J = 1,N-1 BJP = J*TOX*BJ-BJM BJM = BJ BJ = BJP 11 CONTINUE BESSJ = BJ ELSE M = 2*((N+INT(SQRT(FLOAT(IACC*N))))/2) BESSJ = 0. JSUM = 0 SUM = 0. BJP = 0. BJ = 1. DO 12 J = M,1,-1 BJM = J*TOX*BJ-BJP BJP = BJ BJ = BJM IF (ABS(BJ).GT.BIGNO) THEN BJ = BJ*BIGNI BJP = BJP*BIGNI BESSJ = BESSJ*BIGNI SUM = SUM*BIGNI ENDIF IF (JSUM.NE.0) SUM = SUM+BJ JSUM = 1-JSUM IF (J.EQ.N) BESSJ = BJP 12 CONTINUE SUM = 2.*SUM-BJ BESSJ = BESSJ/SUM ENDIF RETURN END FUNCTION BESSJ0 (X) IMPLICIT NONE REAL *8 X,BESSJ0,AX,FR,FS,Z,FP,FQ,XX ! This subroutine calculates the First Kind Bessel Function of ! order 0, for any real number X. The polynomial approximation by ! series of Chebyshev polynomials is used for 0