'*********************************************************** '* This program tests the subroutine BESSK to calculate the* '* modified Bessel function of the third kind of order N * '* for any real positive argument X. * '* ------------------------------------------------------- * '* SAMPLE RUN: * '* * '* X = 1.2097 * '* * '* N = 0 * '* Y = 0.31432449 * '* * '* N = 1 * '* Y = 0.42805075 * '* * '* N = 2 * '* Y = 1.02202186 * '* * '* N = 3 * '* Y = 3.80747328 * '* * '* N = 4 * '* Y = 19.90673679 * '* * '* ------------------------------------------------------- * '* Reference: From Numath Library By Tuan Dang Trong * '* in Fortran 77 [BIBLI 18]. * '* * '* Basic Version By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '*********************************************************** 'PROGRAM TBESSK; DEFDBL A-H, O-Z DEFINT I-N BIG = 1D+30 ' X : Double ' i, NO: Integer NO = 5 'Number of integer orders X = 1.2097# 'value of argument F$ = "##.########" CLS PRINT PRINT " X = "; X PRINT FOR i = 0 TO NO - 1 N = i GOSUB 2000 'call BESSK(i,X) PRINT " N = "; i PRINT " Y = "; USING F$; BESSK PRINT NEXT i END 'of main program ' Bessel Function of the 1st kind of order zero 500 'FUNCTION BESSI0(X) ' Y,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX: Double P1 = 1#: P2 = 3.5156229#: P3 = 3.0899424#: P4 = 1.2067429# P5 = .2659732#: P6 = .0360768#: P7 = .0045813# Q1 = .39894228#: Q2 = .01328592#: Q3 = .00225319#: Q4 = -.00157565# Q5 = .00916281#: Q6 = -.02057706#: Q7 = .02635537# Q8 = -.01647633#: Q9 = .00392377# IF ABS(X) < 3.75 THEN Y = (X / 3.75) ^ 2 TMP = P2 + Y * (P3 + Y * (P4 + Y * (P5 + Y * (P6 + Y * P7)))) BESSI0 = P1 + Y * TMP ELSE AX = ABS(X) Y = 3.75 / AX BX = EXP(AX) / SQR(AX) TMP0 = Q7 + Y * (Q8 + Y * Q9) TMP = Q2 + Y * (Q3 + Y * (Q4 + Y * (Q5 + Y * (Q6 + Y * TMP0)))) AX = Q1 + Y * TMP BESSI0 = AX * BX END IF RETURN ' Bessel Function of the 1st kind of order one 600 'FUNCTION BESSI1(X) ' Y,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX: Double P1 = .5: P2 = .87890594#: P3 = .51498869#: P4 = .15084934# P5 = .02658733#: P6 = .00301532#: P7 = .00032411# Q1 = .39894228#: Q2 = -.03988024#: Q3 = -.00362018# Q4 = .00163801#: Q5 = -.01031555#: Q6 = .02282967# Q7 = -.02895312#: Q8 = .01787654#: Q9 = -.00420059# IF ABS(X) < 3.75 THEN Y = (X / 3.75) ^ 2 TMP0 = P6 + Y * P7 TMP = P1 + Y * (P2 + Y * (P3 + Y * (P4 + Y * (P5 + Y * TMP0)))) BESSI1 = X * TMP ELSE AX = ABS(X) Y = 3.75 / AX BX = EXP(AX) / SQR(AX) TMP0 = Q7 + Y * (Q8 + Y * Q9) TMP = Q2 + Y * (Q3 + Y * (Q4 + Y * (Q5 + Y * (Q6 + Y * TMP0)))) AX = Q1 + Y * TMP BESSI1 = AX * BX END IF RETURN 1000 'FUNCTION BESSK0(X) ' ---------------------------------------------------------------------- ' CALCUL DE LA FONCTION BESSEL MODIFIEE DU 3EME ESPECE D'ORDRE 0 ' POUR TOUT X REEL NON NUL. ' ' CALCULATES THE THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ' ORDER ZERO FOR ANY POSITIVE REAL ARGUMENT, X. ' ---------------------------------------------------------------------} ' Y,AX,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7: Double IF X <= 2# THEN GOSUB 500'Calculate BESSIO(X) P1 = -.57721566#: P2 = .4227842#: P3 = .23069756#: P4 = .0348859# P5 = .00262698#: P6 = .0001075#: P7 = .0000074# Q1 = 1.25331414#: Q2 = -.07832358#: Q3 = .02189568#: Q4 = -.01062446# Q5 = .00587872#: Q6 = -.0025154#: Q7 = .00053208# IF X = 0# THEN BESSK0 = BIG 'arbitrary big value RETURN END IF IF X <= 2# THEN Y = X * X / 4# AX = -LOG(X / 2#) * BESSI0 TMP0 = P6 + Y * P7 TMP = P1 + Y * (P2 + Y * (P3 + Y * (P4 + Y * (P5 + Y * TMP0)))) BESSK0 = AX + TMP ELSE Y = 2# / X AX = EXP(-X) / SQR(X) TMP0 = Q6 + Y * Q7 TMP = Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + Y * (Q5 + Y * TMP0)))) BESSK0 = AX + TMP END IF RETURN 1500 'FUNCTION BESSK1(X) ' ------------------------------------------------------------------------- ' CALCUL DE LA FONCTION BESSEL MODIFIEE DE 3EME ESPECE D'ORDRE 1 ' POUR TOUT X REEL POSITF NON NUL. ' ' CALCULATES THE THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ' ORDER ONE FOR ANY POSITIVE REAL ARGUMENT, X. ' ------------------------------------------------------------------------} ' Y,AX,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7: Double IF X <= 2# THEN GOSUB 600'Calculate BESSI1(X) P1 = 1#: P2 = .15443144#: P3 = -.6727857900000001#: P4 = -.18156897# P5 = -.01919402#: P6 = -.00110404#: P7 = -.00004686# Q1 = 1.25331414#: Q2 = .23498619#: Q3 = -.0365562#: Q4 = .01504268# Q5 = -.00780353#: Q6 = .00325614#: Q7 = -.00068245# IF X = 0# THEN BESSK1 = BIG RETURN END IF IF X <= 2# THEN Y = X * X / 4# AX = LOG(X / 2#) * BESSI1 TMP0 = P6 + Y * P7 TMP = P1 + Y * (P2 + Y * (P3 + Y * (P4 + Y * (P5 + Y * TMP0)))) BESSK1 = AX + (1# / X) * TMP ELSE Y = 2# / X AX = EXP(-X) / SQR(X) TMP0 = Q6 + Y * Q7 TMP = Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + Y * (Q5 + Y * TMP0)))) BESSK1 = AX * TMP END IF RETURN 2000 'FUNCTION BESSK(N, X) ' ------------------------------------------------------------------------ ' CE SOUS-PROGRAMME CALCULE LA FONCTION BESSEL MODIFIFIEE 3E ESPECE ' D'ORDRE N ENTIER POUR TOUT X REEL POSITIF > 0. ON UTILISE ICI LA ' FORMULE DE RECURRENCE CLASSIQUE EN PARTANT DE BESSK0 ET BESSK1. ' ' THIS ROUTINE CALCULATES THE MODIFIED BESSEL FUNCTION OF THE THIRD ' KIND OF INTEGER ORDER, N FOR ANY POSITIVE REAL ARGUMENT, X. THE ' CLASSICAL RECURSION FORMULA IS USED, STARTING FROM BESSK0 AND BESSK1. ' ------------------------------------------------------------------------ ' REFERENCE: ' C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS, ' MATHEMATICAL TABLES, VOL.5, 1962. ' -----------------------------------------------------------------------} ' TOX,BK,BKM,BKP: Double ' J: Integer IF N = 0 THEN GOSUB 1000: BESSK = BESSK0 RETURN END IF IF N = 1 THEN GOSUB 1500: BESSK = BESSK1 RETURN END IF IF X = 0# THEN BESSK = BIG 'arbitrary big value RETURN END IF TOX = 2# / X GOSUB 1500: BK = BESSK1 GOSUB 1000: BKM = BESSK0 FOR J = 1 TO N - 1 BKP = BKM + J * TOX * BK BKM = BK BK = BKP NEXT J BESSK = BK RETURN 'End of file Tbessk.bas