' --------------------------------------------------------------------- ' Program to calculate the first kind modified Bessel function ' of integer order N, for any REAL X, using the function BESSI(N,X). ' --------------------------------------------------------------------- ' SAMPLE RUN: ' ' (Calculate Bessel function for N=2, X=0.75). ' ' Bessel function of order 2 for X = .75 : ' ' Y = 7.366687804794499D-02 ' ' --------------------------------------------------------------------- ' Reference: From Numath Library By Tuan Dang Trong in Fortran 77. ' [BIBLI 18]. ' ' Basic Release 1.0 By J-P Moreau, Paris ' (www.jpmoreau.fr) ' --------------------------------------------------------------------- 'PROGRAM TBESSI; DEFDBL A-H, O-Z DEFINT I-N N = 2 X = .75 GOSUB 500 'Y = BESSI(N,X) CLS PRINT PRINT " Bessel Function of order "; N; " for X="; X; ":" PRINT PRINT " Y = "; Y END ' ---------------------------------------------------------------------- 500 'FUNCTION BESSI(N:Integer; X:Double): Double; ' ' This subroutine calculates the first kind modified Bessel function ' of integer order N, for any REAL X. We use here the classical ' recursion 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 GOSUB 600 RETURN END IF IF N = 1 THEN GOSUB 700 RETURN END IF IF X = 0# THEN Y = 0# RETURN END IF TOX = 2# / X BIP = 0# BI = 1# BSI = 0# M = 2 * ((N + INT(SQR(IACC * N)))) FOR J = M TO 1 STEP -1 BIM = BIP + J * TOX * BI BIP = BI BI = BIM IF ABS(BI) > BIGNO THEN BI = BI * BIGNI BIP = BIP * BIGNI BSI = BSI * BIGNI END IF IF J = N THEN BSI = BIP NEXT J GOSUB 600 Y = BSI * Y / BI RETURN ' ---------------------------------------------------------------------- ' Auxiliary Bessel functions for N=0, N=1. 600 'FUNCTION BESSI0(X:Double): Double; P1 = 1#: P2 = 3.5156229#: P3 = 3.0899424#: P4 = 1.2067429# P5 = .2659732#: P6 = .0360768: P7 = .0045813 Q1 = .39894228#: Q2 = .01328592#: Q3 = 2.25319E-03 Q4 = -1.57565E-03: Q5 = 9.16281E-03: Q6 = -.02057706# Q7 = .02635537#: Q8 = -.01647633#: Q9 = 3.92377E-03 IF ABS(X) < 3.75 THEN Y = (X / 3.75) ^ 2 Y1 = P6 + Y * P7 Y = P1 + Y * (P2 + Y * (P3 + Y * (P4 + Y * (P5 + Y * Y1)))) ELSE AX = ABS(X) Y = 3.75 / AX BX = EXP(AX) / SQR(AX) Y1 = Y * (Q7 + Y * (Q8 + Y * Q9)) AX = Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + Y * (Q5 + Y1)))) Y = AX * BX END IF RETURN ' ---------------------------------------------------------------------- 700 'FUNCTION BESSI1(X:Double): Double; P1 = .5: P2 = .87890594#: P3 = .51498869#: P4 = .15084934# P5 = .02658733#: P6 = 3.01532E-03: P7 = 3.2411E-04 Q1 = .39894228#: Q2 = -.03988024#: Q3 = -3.62018E-03 Q4 = 1.63801E-03: Q5 = -.01031555#: Q6 = .02282967# Q7 = -.02895312#: Q8 = .01787654#: Q9 = -4.20059E-03 IF ABS(X) < 3.75 THEN Y = (X / 3.75) ^ 2 Y1 = P6 + Y * P7 Y = X * (P1 + Y * (P2 + Y * (P3 + Y * (P4 + Y * (P5 + Y * Y1))))) ELSE AX = ABS(X) Y = 3.75 / AX BX = EXP(AX) / SQR(AX) Y1 = Y * (Q7 + Y * (Q8 + Y * Q9)) AX = Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + Y * (Q5 + Y1)))) Y = AX * BX END IF RETURN ' --------------------------------------------------------------------- ' end of file tbessi.bas