{ -------------------------------------------------------------------------- ! 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 (4 items per line) ! 5.1356223E+0000 8.4172442E+0000 1.1619841E+0001 1.4795952E+0001 ! 1.7959819E+0001 2.1116997E+0001 2.4270112E+0001 2.7420574E+0001 ! 3.0569204E+0001 3.3716520E+0001 ! ! Table of error codes (4 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]. ! ! TPW Release 1.0 By J-P Moreau, Paris ! (www.jpmoreau.fr) ! -------------------------------------------------------------------------} PROGRAM TROOTJ; Uses Wincrt; Type Tab10 = Array[1..10] of Double; ITab10 = Array[1..10] of Integer; Var JZ: Tab10; IE: ITab10; I,J,N,NR: Integer; Procedure SECANT(N,NITMX:Integer;TOL:Double;Var ZEROJ:Double;Var IER:Integer); Forward; Function BESSJ(N:Integer;X:Double): Double; Forward; Function BESSJ0(X:Double): Double; Forward; Function BESSJ1(X:Double): Double; Forward; 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.0 then Sign:=ABS(a) else Sign:=-ABS(a) End; { ---------------------------------------------------------------------} PROCEDURE ROOTJ(N,NK:Integer; Var JZERO:Tab10; Var IER:ITab10); { ---------------------------------------------------------------------- ! 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 ! ---------------------------------------------------------------------} Var ZEROJ,B0,B1,B2,B3,B5,B7,T0,T1,T3,T5,T7,FN,FK: Double; C1,C2,C3,C4,C5,F1,F2,F3,TOL,ERRJ: Double; IERROR,K,NITMX: Integer; Begin TOL:=1E-8; NITMX:=10; C1:=1.8557571; C2:=1.033150; C3:=0.00397; C4:=0.0908; C5:=0.043; FN := 1.0*N; { FIRST ZERO } IF N=0 THEN begin ZEROJ := C1+C2-C3-C4+C5; { WRITE(*,'(1X,A,I5,E15.6)') 'K=1,N=0,ZEROJ',K,ZEROJ } SECANT(N,NITMX,TOL,ZEROJ,IERROR); IER[1]:=IERROR; JZERO[1]:=ZEROJ end ELSE begin F1 := Power(FN,(1.0/3.0)); F2 := F1*F1*FN; F3 := F1*FN*FN; ZEROJ := FN+C1*F1+(C2/F1)-(C3/FN)-(C4/F2)+(C5/F3); SECANT(N,NITMX,TOL,ZEROJ,IERROR); IER[1]:=IERROR; JZERO[1]:=ZEROJ END; T0 := 4.0*FN*FN; T1 := T0-1.0; T3 := 4.0*T1*(7.0*T0-31.0); T5 := 32.0*T1*((83.0*T0-982.0)*T0+3779.0); T7 := 64.0*T1*(((6949.0*T0-153855.0)*T0+1585743.0)*T0-6277237.0); { OTHER ZEROES } For K := 2 to NK do begin JZERO[K] := 0.0; FK := 1.0*K; { MAC MAHON'S SERIES FOR K>>N } B0 := (FK+0.5*FN-0.25)*PI; B1 := 8.0*B0; B2 := B1*B1; B3 := 3.0*B1*B2; B5 := 5.0*B3*B2; B7 := 7.0*B5*B2; ZEROJ := B0-(T1/B1)-(T3/B3)-(T5/B5)-(T7/B7); ERRJ:=ABS(BESSJ(N,ZEROJ)); { IMPROVE SOLUTION USING PROCEDURE SECANT } IF ERRJ > TOL THEN SECANT(N,NITMX,TOL,ZEROJ,IERROR); JZERO[K]:=ZEROJ; IER[K]:=IERROR; end End; { ------------------------------------------------------------------------------ } PROCEDURE SECANT(N,NITMX:Integer;TOL:Double;Var ZEROJ:Double;Var IER:Integer); Label 5,10,15,20; Var P0,P1,Q0,Q1,DP,P: Double; C: array[1..2] of Double; IT,NEV,NTRY: Integer; Begin C[1]:=0.95; C[2]:=0.999; 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 } For IT := 1 to NITMX do begin IF Q1=Q0 THEN GOTO 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=1 THEN GOTO 10; IF ABS(DP) < TOL THEN GOTO 20; 10: NEV := NEV+1; P0 := P1; Q0 := Q1; P1 := P; Q1 := BESSJ(N,P1) end; 15: NTRY:=NTRY+1; IF NTRY <= 2 THEN GOTO 5; IER := NTRY; 20: ZEROJ := P End; { ---------------------------------------------------------------------- } FUNCTION BESSJ(N:Integer;X:Double): Double; { 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. } Label EndProc; Const IACC = 40; BIGNO = 1E10; BIGNI = 1E-10; Var TOX,BJM,BJ,BJP,SUM,TMP: Double; J,JSUM,M: Integer; Begin IF N=0 THEN begin BESSJ := BESSJ0(X); goto EndProc end; IF N=1 THEN begin BESSJ := BESSJ1(X); goto EndProc end; IF X=0.0 THEN begin BESSJ := 0.0; goto EndProc end; TOX := 2.0/X; IF X > 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; BESSJ := BJ end end ELSE begin M := 2*((N+Round(SQRT(IACC*N))) div 2); TMP := 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; TMP := TMP * BIGNI; SUM := SUM * BIGNI end; IF JSUM <> 0 THEN SUM := SUM+BJ; JSUM := 1-JSUM; IF J = N THEN TMP := BJP; SUM := 2.0*SUM-BJ end; BESSJ := TMP/SUM end; EndProc: End; { ---------------------------------------------------------------------- } FUNCTION BESSJ0(X:Double): Double; { 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