' -------------------------------------------------------------------------- ' Program to calculate the value of the first kind Bessel function of integer ' order N using the subroutine BESSJ. ' -------------------------------------------------------------------------- ' SAMPLE RUN: ' ' (Calculate the 1st kind Bessel function of order 2 for X=8.41724462111). ' ' ORDER = 2 ' ' X = 8.41724462111 ' ' BESSJ OF ORDER 2 = 1.25443232...D-07 ' ' -------------------------------------------------------------------------- ' Reference: From Numath Library By Tuan Dang Trong in Fortran 77. ' ' Basic Release 1.0 By J-P Moreau, Paris ' (www.jpmoreau.fr) ' -------------------------------------------------------------------------- 'PROGRAM TBESSJ DEFDBL A-H, O-Z DEFINT I-N CLS PRINT PRINT " BESSEL FUNCTIONS OF FIRST KIND" PRINT INPUT " ORDER = ", N PRINT INPUT " X = ", X GOSUB 500 'CALL BESSJ PRINT PRINT " BESSJ OF ORDER "; N; " = "; BESSJ PRINT END ' ---------------------------------------------------------------------- 500 '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. IACC = 40 BIGNO = 1E+10 BIGNI = 1E-10 IF N = 0 THEN 'BESSJ = BESSJ0(X) GOSUB 600: BESSJ = BESSJ0 RETURN END IF IF N = 1 THEN 'BESSJ = BESSJ1(X) GOSUB 700: BESSJ = BESSJ1 RETURN END IF IF X = 0# THEN BESSJ = 0# RETURN END IF TOX = 2# / X IF X > 1# * N THEN 'BJM = BESSJ0(X) GOSUB 600: BJM = BESSJ0 'BJ = BESSJ1(X) GOSUB 700: BJ = BESSJ1 FOR J = 1 TO N - 1 BJP = J * TOX * BJ - BJM BJM = BJ BJ = BJP NEXT J BESSJ = BJ ELSE M = 2 * ((N + INT(SQR(1# * IACC * N))) / 2) BESSJ = 0# JSUM = 0 SUM = 0# BJP = 0# BJ = 1# FOR J = M TO 1 STEP -1 BJM = J * TOX * BJ - BJP BJP = BJ BJ = BJM IF ABS(BJ) > BIGNO THEN BJ = BJ * BIGNI BJP = BJP * BIGNI BESSJ = BESSJ * BIGNI SUM = SUM * BIGNI END IF IF JSUM <> 0 THEN SUM = SUM + BJ JSUM = 1 - JSUM IF J = N THEN BESSJ = BJP NEXT J SUM = 2# * SUM - BJ BESSJ = BESSJ / SUM END IF RETURN ' ---------------------------------------------------------------------- 600 '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= 0# THEN SIGN = ABS(S6) ELSE SIGN = -ABS(S6) END IF BESSJ1 = SQR(P6 / AX) * (COS(XX) * FP - Z * SIN(XX) * FQ) * SIGN END IF RETURN ' end of file tbessj.bas