' ---------------------------------------------------------------------
' 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