'******************************************************************
'*  Program to demonstrate the Bessel function asymptotic series  *
'* -------------------------------------------------------------- *
'*       Reference: BASIC Scientific Subroutines, Vol. II         *
'*       By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].        *
'*                                                                *
'* -------------------------------------------------------------- *
'* SAMPLE RUN:                                                    *
'*                                                                *
'*  X            J0(X)            J1(X)        N            E     *
'* -------------------------------------------------------------- *
'*  1           0.733562         0.402234      1       0.1121521  *
'*  2           0.221488         0.578634      1       0.0070095  *
'*  3          -0.259956         0.340699      2       0.0007853  *
'*  4          -0.396826        -0.065886      2       0.0001398  *
'*  5          -0.177486        -0.327696      3       0.0000155  *
'*  6           0.150635        -0.276760      3       0.0000036  *
'*  7           0.300051        -0.004696      4       0.0000004  *
'*  8           0.171638         0.234650      5       0.0000000  *
'*  9          -0.090332         0.245324      5       0.0000000  *
'* 10          -0.245930         0.043476      6       0.0000000  *
'* 11          -0.171187        -0.176788      5       0.0000000  *
'* 12           0.047689        -0.223450      5       0.0000000  *
'* 13           0.206925        -0.070319      4       0.0000000  *
'* 14           0.171072         0.133376      4       0.0000000  *
'* 15          -0.014225         0.205105      4       0.0000000  *
'* -------------------------------------------------------------- *
'*                                                                *
'******************************************************************
DEFINT I-N
DEFDBL A-H, O-Z
DIM J0 AS DOUBLE, J1 AS DOUBLE
DIM m1 AS DOUBLE, m2 AS DOUBLE, n1 AS DOUBLE, n2 AS DOUBLE
CLS
PRINT
INPUT "What is the desired error bound ? ", e3
CLS
PRINT
PRINT "  X             J0(X)            J1(X)        N            E      "
PRINT "------------------------------------------------------------------"
FOR x = 1 TO 15
  GOSUB 1000
  PRINT USING " ##           #.########"; x; J0;
  PRINT USING "       #.########     #"; J1; n;
  PRINT USING "         #.#######"; e
NEXT x
PRINT
END

'**********************************************************
'* Bessel function asymptotic series subroutine. This     *
'* program calculates the zeroth and first order Bessel   *
'* functions using an asymptotic series expansion.        *
'* The required input are X and a convergence factor e.   *
'* Returned are the two Bessel functions, J0(X) and J1(X) *
'* ------------------------------------------------------ *
'* Reference: Algorithms for RPN calculators, by Ball,    *
'* L.A. Wiley and sons.                                   *
'**********************************************************
' Calculate P and Q polynomials: P0(x)=m1, P1(x)=m2,
' Q0(x)=n1, Q1(x)=n2
1000 a = 1: a1 = 1: a2 = 1: b = 1: c = 1
  e1 = 1000000#
  m = -1
  x1 = 1# / (8# * x)
  x1 = x1 * x1
  m1 = 1#: m2 = 1#
  n1 = -1# / (8# * x)
  n2 = -3# * n1
  n = 0:
1100  m = m + 2: a = a * m * m
  m = m + 2: a = a * m * m
  c = c * x1
  a1 = a1 * a2: a2 = a2 + 1#: a1 = a1 * a2
  a2 = a2 + 1#: e2 = a * c / a1
  e4 = 1# + (m + 2) / m + (m + 2) * (m + 2) / (a2 * 8# * x) + (m + 2) * (m + 4) / (a2 * 8# * x)
  e4 = e4 * e2
  'Test for divergence
  IF ABS(e4) > e1 THEN GOTO 1200
  e1 = ABS(e2)
  m1 = m1 - e2
  m2 = m2 + e2 * (m + 2) / m
  n1 = n1 + e2 * (m + 2) * (m + 2) / (a2 * 8# * x)
  n2 = n2 - e2 * (m + 2) * (m + 4) / (a2 * 8# * x)
  n = n + 1
  'Test for convergence criterion
  IF e1 < e3 THEN GOTO 1200
  GOTO 1100
1200 a = 3.1415926535#
  e = e2
  b = SQR(2# / (a * x))
  J0 = b * (m1 * COS(x - a / 4#) - n1 * SIN(x - a / 4#))
  J1 = b * (m2 * COS(x - 3# * a / 4#) - n2 * SIN(x - 3# * a / 4#))
RETURN

'End of file Bessel1.bas