{*********************************************************************** * * * Program to calculate the first kind Bessel function of integer * * order N, for any REAL X, using the function BESSJ(N,X). * * * * -------------------------------------------------------------------- * * * * SAMPLE RUN: * * * * (Calculate Bessel function for N=2, X=0.75). * * * * Bessel function of order 2 for X = 0.7500: * * * * Y = 6.70739973016599E-0002 * * * * -------------------------------------------------------------------- * * Reference: From Numath Library By Tuan Dang Trong in Fortran 77. * * * * TPW Release 1.0 By J-P Moreau, Paris. * * (www.jpmoreau.fr) * ***********************************************************************} PROGRAM TBESSJ; Uses WinCrt; Var X, Y: Double; N: Integer; FUNCTION BESSJ0 (X:Double): Double; { This subroutine calculates the First Kind Bessel Function of order 0, for any real number X. The polynomial approximation by series of Chebyshev polynomials is used for 0 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 10; Const IACC = 40; BIGNO = 1E10; BIGNI = 1E-10; Var B,TOX,BJM,BJ,BJP,SUM,TMP: 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 := 2*((N+Round(SQRT(1.0*(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 end; SUM := 2.*SUM-BJ; BESSJ := TMP/SUM end; 10: End; {main program} BEGIN N:=2; X:=0.75; Y := BESSJ(N,X); writeln; writeln(' Bessel Function of order ', N,' for X=', X:8:4,':'); writeln; writeln(' Y = ', Y); readkey; Donewincrt END. {End of file Tbessj.pas}