!************************************************************************ !* * !* Program to calculate the first kind modified Bessel function of * !* integer order N, for any REAL X, using the function BESSI(N,X). * !* * !* -------------------------------------------------------------------- * !* * !* SAMPLE RUN: * !* * !* (Calculate Bessel function for N=2, X=0.75). * !* * !* Bessel function of order 2 for X = 0.7500: * !* * !* Y = 0.73666878E-01 * !* * !* -------------------------------------------------------------------- * !* Reference: From Numath Library By Tuan Dang Trong in Fortran 77. * !* * !* F90 Release 1.2 By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !* * !* Version 1.1: corected value of P4 in BESSIO (P4=1.2067492 and not * !* 1.2067429) Aug. 2011. * !* Version 2: all variables are declared. * !************************************************************************ PROGRAM TBESSI IMPLICIT NONE REAL*8 BESSI, X, Y INTEGER N N=2 X=0.75D0 Y = BESSI(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/) END ! ---------------------------------------------------------------------- FUNCTION BESSI(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 N, M, J REAL *8 X,BESSI,BESSI0,BESSI1,TOX,BIM,BI,BIP IF (N.EQ.0) THEN BESSI = BESSI0(X) RETURN ENDIF IF (N.EQ.1) THEN BESSI = BESSI1(X) RETURN ENDIF IF(X.EQ.0.D0) THEN BESSI=0.D0 RETURN ENDIF TOX = 2.D0/X BIP = 0.D0 BI = 1.D0 BESSI = 0.D0 M = 2*((N+INT(SQRT(FLOAT(IACC*N))))) DO 12 J = M,1,-1 BIM = BIP+DFLOAT(J)*TOX*BI BIP = BI BI = BIM IF (ABS(BI).GT.BIGNO) THEN BI = BI*BIGNI BIP = BIP*BIGNI BESSI = BESSI*BIGNI ENDIF IF (J.EQ.N) BESSI = BIP 12 CONTINUE BESSI = BESSI*BESSI0(X)/BI RETURN END ! ---------------------------------------------------------------------- ! Auxiliary Bessel functions for N=0, N=1 FUNCTION BESSI0(X) IMPLICIT NONE REAL *8 X,BESSI0,Y,P1,P2,P3,P4,P5,P6,P7, & Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX DATA P1,P2,P3,P4,P5,P6,P7/1.D0,3.5156229D0,3.0899424D0,1.2067492D0, & 0.2659732D0,0.360768D-1,0.45813D-2/ DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,0.1328592D-1, & 0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1, & 0.2635537D-1,-0.1647633D-1,0.392377D-2/ IF(ABS(X).LT.3.75D0) THEN Y=(X/3.75D0)**2 BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))) ELSE AX=ABS(X) Y=3.75D0/AX BX=EXP(AX)/SQRT(AX) AX=Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))) BESSI0=AX*BX ENDIF RETURN END ! ---------------------------------------------------------------------- FUNCTION BESSI1(X) IMPLICIT NONE REAL *8 X,BESSI1,Y,P1,P2,P3,P4,P5,P6,P7, & Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0, & 0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/ DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,-0.3988024D-1, & -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1, & -0.2895312D-1,0.1787654D-1,-0.420059D-2/ IF(ABS(X).LT.3.75D0) THEN Y=(X/3.75D0)**2 BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))) ELSE AX=ABS(X) Y=3.75D0/AX BX=EXP(AX)/SQRT(AX) AX=Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))) BESSI1=AX*BX ENDIF RETURN END ! ---------------------------------------------------------------------- ! end of file tbessi.f90