!************************************************************** !* CALCULATE THE Kth ZERO OF THE DERIVATIVE OF BESSEL * !* FUNCTION OF ORDER N, J(N,X) * !* ---------------------------------------------------------- * !* SAMPLE RUN: * !* (Calculate the 10th zero of the derivative of J(1,X) ). * !* * !* X= 30.6019229941596 * !* * !* BESSJP(X)= -1.2229801E-16 * !* * !* ---------------------------------------------------------- * !* Reference: From Numath Library By Tuan Dang Trong in * !* Fortran 77 [BIBLI 18]. * !* * !* F90 Release By J-P Moreau, Paris * !* (www.jpmoreau.fr) * !************************************************************** PROGRAM TEST_ZEROJP REAL*8 RES, ZEROJP N=1 K=10 RES = ZEROJP(N,K) print *,' ' print *,' X=', RES print *,' ' print *,' BESSJP(X)=', BESSJP(N,RES) print *,' ' END REAL*8 FUNCTION ZEROJP(N,K) !-------------------------------------------------------------------- ! CALCULATE THE Kth ZERO OF THE DERIVATIVE OF BESSEL FUNCTION ! OF ORDER N, J(N,X) !-------------------------------------------------------------------- ! CALLING MODE: ! RES = ZEROJP(N,K) ! ! INPUTS: ! N ORDER OF BESSEL FUNCTION J (INTEGER >= 0) I*4 ! K RANK OF ZERO (INTEGER > 0) I*4 ! OUTPUT: ! ZEROJP R*8 ! REFERENCE: ! ABRAMOWITZ M. & STEGUN IRENE A. ! HANDBOOK OF MATHEMATICAL FUNCTIONS !--------------------------------------------------------------------- REAL*8 BESSJP,B0,B1,B2,B3,B5,B7,T0,T1,T3,T5,T7,PI,FN,FK, & C1,C2,C3,C4,F1,F2,F3,P,DP,P0,P1,Q0,Q1,TOL LOGICAL IMPROV DATA TOL/1.D-7/,NITMX/15/ DATA C1,C2,C3,C4 /0.8086165D0,0.072490D0,.05097D0,.0094D0/ DATA IMPROV/.TRUE./ PI = 4.d0*ATAN(1.d0) FN = DFLOAT(N) FK = DFLOAT(K) IF (K.GT.1) GO TO 10 ! SI N = 0 ET K = 1 IF (N.EQ.0) THEN ZEROJP= 0.D0 RETURN ! TCHEBYCHEV'S SERIES FOR K <= N ELSE F1 = FN**(1.D0/3.D0) F2 = F1*F1*FN ZEROJP = FN+C1*F1+(C2/F1)-(C3/FN)+(C4/F2) GO TO 20 ENDIF ! MAC MAHON'S SERIES FOR K >> N 10 B0 = (FK+.5D0*FN-.75D0)*PI B1 = 8.D0*B0 B2 = B1*B1 B3 = 3.D0*B1*B2 B5 = 5.D0*B3*B2 B7 = 7.D0*B5*B2 T0 = 4.D0*FN*FN T1 = T0+3.D0 T3 = 4.D0*((7.D0*T0+82.D0)*T0-9.D0) T5 = 32.D0*(((83.D0*T0+2075.D0)*T0-3039.D0)*T0+3537.D0) T7 = 64.D0*((((6949.D0*T0+296492.D0)*T0-1248002.D0)*T0 & +7414380.D0)*T0-5853627.D0) ZEROJP = B0-(T1/B1)-(T3/B3)-(T5/B5)-(T7/B7) 20 IF (IMPROV) THEN ! IMPROVE SOLUTION BY SECANT METHOD WHEN K > N ! AND IMPROV = .TRUE. P0 = 0.9D0*ZEROJP P1 = ZEROJP IER = 0 NEV = 2 Q0 = BESSJP(N,P0) Q1 = BESSJP(N,P1) DO 30 IT = 1,NITMX P = P1-Q1*(P1-P0)/(Q1-Q0) DP = P-P1 IF (IT.EQ.1) GO TO 25 IF (ABS(DP).LT.TOL) GO TO 40 25 NEV = NEV+1 P0 = P1 Q0 = Q1 P1 = P Q1 = BESSJP(N,P1) 30 CONTINUE IER = 1 WRITE(*,'(1X,A)') '** ZEROJP ** NITMX EXCEEDED' RETURN 40 ZEROJP = P ENDIF RETURN END FUNCTION BESSJP (N,X) ! ---------------------------------------------------------------------- ! NAME : BESSJP ! DATE : 06/01/1982 ! IV : 1 ! IE : 1 ! AUTHOR: DANG TRONG TUAN ! ...................................................................... ! ! FIRST DERIVATIVE OF FIRST KIND BESSEL FUNCTION OF ORDER N, FOR REAL X ! ! MODULE BESSJP . ! ...................................................................... ! ! THIS SUBROUTINE CALCULATES THE FIRST DERIVATIVE OF FIRST KIND BESSEL ! FUNCTION OF ORDER N, FOR REAL X. ! . ! ...................................................................... ! ! I VARIABLE DIMENSION/TYPE DESCRIPTION (INPUTS) ! N I*4 ORDER OF FUNCTION . ! X R*8 ABSCISSA OF FUNCTION BESSJP(N,X) . ! ! O VARIABLE,DIMENSION/TYPE DESCRIPTION (OUTPUT) ! ! BESSJP R*8 FUNCTION EVALUATION AT X . !....................................................................... ! CALLED SUBROUTINE ! ! BESSJ FIRST KIND BESSEL FUNCTION ! ! ---------------------------------------------------------------------- DOUBLE PRECISION X,BESSJP,BESSJ IF (N.EQ.0) THEN BESSJP=-BESSJ(1,X) ELSE IF(X.EQ.0.D0) THEN X=1.D-30 ELSE BESSJP=BESSJ(N-1,X)-( FLOAT(N)/X)*BESSJ(N,X) ENDIF RETURN END REAL*8 FUNCTION BESSJ (N,X) !----------------------------------------------------------------------- ! THIS SUBROUTINE CALCULATES THE 1ST KIND BESSEL FUNCTION OF ORDER ! N, INTEGER FOR ANY REAL X. THE CLASSICAL RECURRENCE FORMULA IS USED ! HERE, WHEN X IS > N, FOR X < N, THE MILLER'S ALGORITHM ALLOWS ! AVOIDING 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,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 REAL*8 FUNCTION BESSJ0 (X) REAL *8 X,AX,FR,FS,Z,FP,FQ,XX !---------------------------------------------------------------------- ! THIS SUBROUTINE CALCULATES THE 1ST KIND BESSEL FUNCTION OF ORDER ! ZERO FOR ANY REAL X. THE CHEBYSHEV'S POLYNOMIAL SERIES IS USED ! FOR 0