! -------------------------------------------------------------------------- ! Program to calculate the first zeroes (root abscissas) of the first ! kind Bessel function of integer order N using the subroutine ROOTJ. ! -------------------------------------------------------------------------- ! SAMPLE RUN: ! ! (Calculate the first 10 zeroes of 1st kind Bessel function of order 2). ! ! Zeroes of Bessel Function of order: 2 ! ! Number of calculated zeroes: 10 ! ! Table of root abcissas (5 items per line) ! 0.51356223D+01 0.84172442D+01 0.11619841D+02 0.14795952D+02 0.17959819D+02 ! 0.21116997D+02 0.24270112D+02 0.27420574D+02 0.30569204D+02 0.33716520D+02 ! ! Table of error codes (5 items per line) ! 0 0 0 0 0 ! 0 0 0 0 0 ! ! -------------------------------------------------------------------------- ! Reference: From Numath Library By Tuan Dang Trong in Fortran 77 ! [BIBLI 18]. ! ! F90 Release 1.0 By J-P Moreau, Paris ! (www.jpmoreau.fr) !--------------------------------------------------------------------------- PROGRAM TROOTJ REAL*8 JZ(10) INTEGER IE(10) N=2 NR=10 CALL ROOTJ(N,NR,JZ,IE) print *,' ' write(*,10) N write(*,20) NR print *,' ' print *,'Table of root abcissas (5 items per line)' write(*,50) (JZ(i),i=1,NR) print *,' ' print *,'Table of error codes (5 items per line)' write(*,60) (IE(i),i=1,NR) print *,' ' stop 10 format(/' Zeroes of Bessel Function of order: ', I2) 20 format(/' Number of calculated zeroes: ', I2) 50 format(5D15.8) 60 format(5I4) END SUBROUTINE ROOTJ(N,NK,JZERO,IER) ! ---------------------------------------------------------------------- ! CALCULATE THE FIRST NK ZEROES OF BESSEL FUNCTION J(N,X) ! ! INPUTS: ! N ORDER OF FUNCTION J (INTEGER >= 0) I*4 ! NK NUMBER OF FIRST ZEROES (INTEGER > 0) I*4 ! OUTPUTS: ! JZERO(NK) TABLE OF FIRST ZEROES (ABCISSAS) R*8 ! IER(NK) TABLE OF ERROR CODES (MUST BE ZEROES) I*4 ! ! REFERENCE : ! ABRAMOWITZ M. & STEGUN IRENE A. ! HANDBOOK OF MATHEMATICAL FUNCTIONS ! ---------------------------------------------------------------------- REAL*8 BESSJ,ZEROJ,B0,B1,B2,B3,B5,B7,T0,T1,T3,T5,T7,PI,FN,FK, & C1,C2,C3,C4,C5,F1,F2,F3,TOL,ERRJ,JZERO(NK) INTEGER IER(NK) DATA PI/3.14159265358979D0/,TOL/1.D-8/,NITMX/10/ DATA C1,C2,C3,C4,C5 /1.8557571D0,1.033150D0,.00397D0,.0908D0,.043D0/ FN = FLOAT(N) ! FIRST ZERO IF (N.EQ.0) THEN ZEROJ = C1+C2-C3-C4+C5 ! WRITE(*,'(1X,A,I5,E15.6)') 'K=1,N=0,ZEROJ',K,ZEROJ CALL SECANT(N,NITMX,TOL,ZEROJ,IERROR) IER(1)=IERROR JZERO(1)=ZEROJ ELSE F1 = FN**(1.D0/3.D0) F2 = F1*F1*FN F3 = F1*FN*FN ZEROJ = FN+C1*F1+(C2/F1)-(C3/FN)-(C4/F2)+(C5/F3) ! WRITE(*,'(1X,A,I5,E15.6)') 'K=1,NE.0,ZEROJ',K,ZEROJ CALL SECANT(N,NITMX,TOL,ZEROJ,IERROR) IER(1)=IERROR JZERO(1)=ZEROJ ENDIF T0 = 4.D0*FN*FN T1 = T0-1.D0 T3 = 4.D0*T1*(7.D0*T0-31.D0) T5 = 32.D0*T1*((83.D0*T0-982.D0)*T0+3779.D0) T7 = 64.D0*T1*(((6949.D0*T0-153855.D0)*T0+1585743.D0)*T0 & -6277237.D0) ! OTHER ZEROES DO K = 2,NK JZERO(K) = 0.D0 FK = FLOAT(K) ! MAC MAHON'S SERIES FOR K>>N B0 = (FK+.5D0*FN-.25D0)*PI B1 = 8.D0*B0 B2 = B1*B1 B3 = 3.D0*B1*B2 B5 = 5.D0*B3*B2 B7 = 7.D0*B5*B2 ZEROJ = B0-(T1/B1)-(T3/B3)-(T5/B5)-(T7/B7) ERRJ=ABS(BESSJ(N,ZEROJ)) ! WRITE(*,'(1X,A,2I5,2E15.6)') 'N,K,ZEROJ,ERRJ',N,K,ZEROJ,ERRJ ! IMPROVE SOLUTION USING SUBROUTINE SECANT IF (ERRJ.GT.TOL) CALL SECANT(N,NITMX,TOL,ZEROJ,IERROR) JZERO(K)=ZEROJ IER(K)=IERROR ENDDO RETURN END ! ---------------------------------------------------------------------- SUBROUTINE SECANT(N,NITMX,TOL,ZEROJ,IER) REAL *8 TOL,ZEROJ,P0,P1,Q0,Q1,DP,P,C(2) DATA C/.95D0,.999D0/ NTRY=1 IER=0 5 P0 = C(NTRY)*ZEROJ P1 = ZEROJ NEV = 2 Q0 = BESSJ(N,P0) Q1 = BESSJ(N,P1) ! WRITE(*,'(1X,A,I5,4E15.6)') 'NTRY,P0,Q0,P1,Q1',NTRY,P0,Q0,P1,Q1 DO IT = 1,NITMX IF(Q1.EQ.Q0) GO TO 15 P = P1-Q1*(P1-P0)/(Q1-Q0) DP = P-P1 ! WRITE(*,'(1X,A,I5,4E15.6)') 'IT,P,DP',IT,P,DP IF (IT.EQ.1) GO TO 10 IF (ABS(DP).LT.TOL) GO TO 20 10 NEV = NEV+1 P0 = P1 Q0 = Q1 P1 = P Q1 = BESSJ(N,P1) ENDDO 15 NTRY=NTRY+1 IF(NTRY.LE.2) GO TO 5 IER= NTRY 20 ZEROJ = P END ! ---------------------------------------------------------------------- FUNCTION BESSJ (N,X) ! THIS FUNCTION RETURNS THE VALUE OF THE FIRST KIND BESSEL FUNCTION ! OF ORDER N, INTEGER FOR ANY REAL X. WE USE HERE THE CLASSICAL ! RECURRENT 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. PARAMETER (IACC = 40,BIGNO = 1.D10, BIGNI = 1.D-10) 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) ! THIS FUNCTION RETURNS THE VALUE OF THE FIRST KIND BESSEL FUNCTION ! OF ORDER 0 FOR ANY REAL X. WE USE HERE THE POLYNOMIAL APPROXIMATION ! BY SERIES OF CHEBYSHEV POLYNOMIALS FOR 0