{************************************************************* * 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= 3.06019229735146E+0001 * * * * BESSJP(X)= -1.75243069973916E-0016 * * * * ---------------------------------------------------------- * * Reference: From Numath Library By Tuan Dang Trong in * * Fortran 77 [BIBLI 18]. * * * * TPW Release By J-P Moreau, Paris * * (www.jpmoreau.fr) * *************************************************************} PROGRAM TEST_ZEROJP; Uses WinCrt; VAR RES: Double; K,N: Integer; FUNCTION Power(y,x:Double):Double; Begin IF x<0 THEN EXIT; Power:=Exp(x*Ln(y)) End; Function Sign(a,b:Double): Double; Begin if b<0 then Sign:=-ABS(a) else Sign:=ABS(a) End; FUNCTION BESSJ (N:Integer;X:Double):Double; Forward; FUNCTION BESSJ0(X:Double):Double; Forward; FUNCTION BESSJ1(X:Double):Double; Forward; FUNCTION BESSJP(N:Integer;X:Double):Double; Forward; FUNCTION ZEROJP(N,K:Integer): Double; {-------------------------------------------------------------------- ! 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 !--------------------------------------------------------------------} Label 10,20,25, 40,50; Var B0,B1,B2,B3,B5,B7,T0,T1,T3,T5,T7,FN,FK: Double; C1,C2,C3,C4,F1,F2,P,DP,P0,P1,Q0,Q1,TOL: Double; ZJP: Double; IER,IT,NEV,NITMX: Integer; IMPROV: Boolean; Begin TOL:=1.e-7; NITMX:=15; C1:=0.8086165; C2:=0.072490; C3:=0.05097; C4:=0.0094; IMPROV:=TRUE; FN := 1.0*N; FK := 1.0*K; IF K > 1 THEN GOTO 10; { if N = 0 ET K = 1 } IF N = 0 THEN begin ZEROJP:= 0.0; goto 50; end { TCHEBYCHEV'S SERIES FOR K <= N } ELSE begin F1 := Power(FN,1.0/3.0); F2 := F1*F1*FN; ZJP := FN+C1*F1+(C2/F1)-(C3/FN)+(C4/F2); GOTO 20 end; { MAC MAHON'S SERIES FOR K >> N } 10: B0 := (FK+0.5*FN-0.75)*PI; B1 := 8.0*B0; B2 := B1*B1; B3 := 3.0*B1*B2; B5 := 5.0*B3*B2; B7 := 7.0*B5*B2; T0 := 4.0*FN*FN; T1 := T0 + 3.0; T3 := 4.0*((7.0*T0+82.0)*T0-9.0); T5 := 32.0*(((83.0*T0+2075.0)*T0-3039.0)*T0+3537.0); T7 := 64.0*((((6949.0*T0+296492.0)*T0-1248002.0)*T0 +7414380.0)*T0-5853627.0); ZJP := B0-(T1/B1)-(T3/B3)-(T5/B5)-(T7/B7); 20: IF IMPROV THEN begin { IMPROVE SOLUTION BY SECANT METHOD WHEN K > N AND IMPROV = TRUE } P0 := 0.9*ZJP; P1 := ZJP; IER := 0; NEV := 2; Q0 := BESSJP(N,P0); Q1 := BESSJP(N,P1); For IT := 1 to NITMX do begin P := P1-Q1*(P1-P0)/(Q1-Q0); DP := P-P1; IF IT = 1 THEN GOTO 25; IF ABS(DP) < TOL THEN GOTO 40; 25: NEV := NEV+1; P0 := P1; Q0 := Q1; P1 := P ; Q1 := BESSJP(N,P1); end; IER := 1; WRITELN(' ** ZEROJP ** NITMX EXCEEDED.'); ZEROJP:=ZJP; goto 50; 40: ZEROJP := P end; 50: End; FUNCTION BESSJP(N:Integer; X:Double): Double; { ---------------------------------------------------------------------- ! 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 ! ! ----------------------------------------------------------------------} Begin IF N = 0 THEN BESSJP:=-BESSJ(1,X) ELSE IF X = 0.0 THEN X:=1.e-30 ELSE BESSJP:=BESSJ(N-1,X)-(1.0*N/X)*BESSJ(N,X); End; FUNCTION BESSJ(N:Integer; X:Double): Double; {----------------------------------------------------------------------- ! 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. !----------------------------------------------------------------------} Label 10; Const IACC = 40; BIGNO = 1e10; BIGNI = 1e-10; Var TOX,BJM,BJ,BJP,BSJ,SUM: Double; J,JSUM,M:Integer; Begin IF N = 0 THEN begin BESSJ := BESSJ0(X); goto 10 end; IF N = 1 THEN begin BESSJ := BESSJ1(X); goto 10 end; IF X = 0.0 THEN begin BESSJ := 0.0; goto 10 end; TOX := 2.0/X; IF X > 1.0*N THEN begin BJM := BESSJ0(X); BJ := BESSJ1(X); For J := 1 to N-1 do begin BJP := J*TOX*BJ-BJM; BJM := BJ; BJ := BJP end; BESSJ := BJ end ELSE begin M := Round(2*((N+Int(SQRT(IACC*N)))/2)); BSJ := 0.0; JSUM := 0; SUM := 0.0; BJP := 0.0; BJ := 1.0; For J := M Downto 1 do begin BJM := J*TOX*BJ-BJP; BJP := BJ; BJ := BJM; IF ABS(BJ) > BIGNO THEN begin BJ := BJ*BIGNI; BJP := BJP*BIGNI; BSJ := BSJ*BIGNI; SUM := SUM*BIGNI end; IF JSUM <> 0 then SUM := SUM+BJ; JSUM := 1-JSUM; IF J = N then BSJ := BJP end; SUM := 2.0*SUM-BJ; BESSJ := BSJ/SUM end; 10: End; FUNCTION BESSJ0(X:Double):Double; {---------------------------------------------------------------------- ! THIS SUBROUTINE CALCULATES THE 1ST KIND BESSEL FUNCTION OF ORDER ! ZERO FOR ANY REAL X. THE CHEBYSHEV'S POLYNOMIAL SERIES IS USED ! FOR 0