'*********************************************************
'* Program to demonstrate Bessel Coefficients Subroutine *
'* ----------------------------------------------------- *
'*    Reference: BASIC Scientific Subroutines, Vol. II   *
'*    by F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].  *
'*                                                       *
'* ----------------------------------------------------- *
'* SAMPLE RUN:                                           *
'*                                                       *
'* BESSEL COEFFICIENTS                                   *
'*                                                       *
'* What is the order of the Bessel function ?  0         *
'* What degree is desired ? 41                           *
'*                                                       *
'* The coefficients are:                                 *
'*                                                       *
'*  A( 0) =  1                        A( 1) =  0         *
'*  A( 2) = -.25                      A( 3) =  0         *
'*  A( 4) =  .015625                  A( 5) =  0         *
'*  A( 6) = -4.340277777777778D-04    A( 7) =  0         *
'*  A( 8) =  6.781684027777777D-06    A( 9) =  0         *
'*  A(10) = -6.781684027777778D-08    A(11) =  0         *
'*  A(12) =  4.709502797067901D-10    A(13) =  0         *
'*  A(14) = -2.40280754952444D-12     A(15) =  0         *
'*  A(16) =  9.385966990329842D-15    A(17) =  0         *
'*  A(18) = -2.896903392077112D-17    A(19) =  0         *
'*  A(20) =  7.242258480192779D-20    A(21) =  0         *
'*  A(22) = -1.496334396734045D-22    A(23) =  0         *
'*  A(24) =  2.597802772107717D-25    A(25) =  0         *
'*  A(26) = -3.842903509035085D-28    A(27) =  0         *
'*  A(28) =  4.901662639075364D-31    A(29) =  0         *
'*  A(30) = -5.446291821194849D-34    A(31) =  0         *
'*  A(32) =  5.318644356635594D-37    A(33) =  0         *
'*  A(34) = -4.60090342269515D-40     A(35) =  0         *
'*  A(36) =  3.550079801462307D-43    A(37) =  0         *
'*  A(38) = -2.458504017633177D-46    A(39) =  0         *
'*  A(40) =  1.536565011020736D-49    A(41) =  0         *
'*                                                       *
'* Argument ?  1                                         *
'* Y=  0.765198                                          *
'*                                                       *
'*********************************************************
DEFINT I-N
DEFDBL A-H, O-Z
CLS
PRINT "BESSEL COEFFICIENTS"
PRINT
INPUT "What is the order of the Bessel function ? ", n
INPUT "What degree is desired ? ", m
PRINT
DIM A(m + n), B(m + n)
GOSUB 1000
PRINT "The coefficients are:"
PRINT
F$ = " A(##) = "
FOR i = 0 TO m STEP 2
  PRINT USING F$; i; : PRINT A(i);
  PRINT TAB(35); : PRINT USING F$; i + 1; : PRINT A(i + 1)
NEXT
PRINT
INPUT "Argument ? ", x
y = A(0) + A(1) * x
IF m > 1 THEN y = y + A(2) * x * x
IF m > 2 THEN y = y + A(3) * x * x * x
IF m > 3 THEN y = y + A(4) * x ^ 4
FOR i = 4 TO m - 1
  IF m > i THEN y = y + A(i + 1) * x ^ (i + 1)
NEXT
PRINT
PRINT USING " Y=##.#######"; y
PRINT
END

'*************************************************************
'* Bessel function series coefficient evaluation subroutine  *
'* m+1 is the number of coefficients desired, n is the order *
'* of the Bessel function. The coefficients are returned in  *
'* A(i).                                                     *
'*************************************************************
1000 a1 = 1: b1 = 1
  FOR i = 1 TO n
    B(i - 1) = 0: b1 = b1 * i: a1 = a1 / 2
  NEXT
  b1 = a1 / b1: a1 = 1
  FOR i = 0 TO m STEP 2
    A(i) = a1 * b1: A(i + 1) = 0
    a1 = -a1 / ((i + 2) * (n + n + i + 2))
  NEXT
  a1 = a1 / 2
  FOR i = 0 TO m
    B(i + n) = A(i)
  NEXT
  FOR i = 0 TO n + m
    A(i) = B(i)
  NEXT
RETURN

'End file Bessel.bas