'********************************************************************
'*     Program to demonstrate integer Bessel Function Subroutine    *
'* ---------------------------------------------------------------- *
'*         Reference: BASIC Scientific Subroutines, Vol. II         *
'*         By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].        *
'*                                                                  *
'* ---------------------------------------------------------------- *
'* SAMPLE RUN:                                                      *
'*                                                                  * 
'* X       J0(X)       J1(X)       J2(X)       J3(X)       J4(X)    *
'*----------------------------------------------------------------- *
'* 0.1    0.9975016   0.0499375   0.0012490   0.0000208   0.0000003 *
'* 0.2    0.9900250   0.0995008   0.0049834   0.0001663   0.0000042 *
'* 0.3    0.9776262   0.1483188   0.0111659   0.0005593   0.0000210 *
'* 0.4    0.9603982   0.1960266   0.0197347   0.0013201   0.0000661 *
'* 0.5    0.9384698   0.2422685   0.0306040   0.0025637   0.0001607 *
'* 0.6    0.9120049   0.2867010   0.0436651   0.0043997   0.0003315 *
'* 0.7    0.8812009   0.3289957   0.0587869   0.0069297   0.0006101 *
'* 0.8    0.8462874   0.3688420   0.0758178   0.0102468   0.0010330 *
'* 0.9    0.8075238   0.4059495   0.0945863   0.0144340   0.0016406 *
'* 1.0    0.7651977   0.4400506   0.1149035   0.0195634   0.0024766 *
'* 1.1    0.7196220   0.4709024   0.1365642   0.0256945   0.0035878 *
'* 1.2    0.6711327   0.4982891   0.1593490   0.0328743   0.0050227 *
'* 1.3    0.6200860   0.5220232   0.1830267   0.0411358   0.0068310 *
'* 1.4    0.5668551   0.5419477   0.2073559   0.0504977   0.0090629 *
'* 1.5    0.5118277   0.5579365   0.2320877   0.0609640   0.0117681 *
'* 1.6    0.4554022   0.5698959   0.2569678   0.0725234   0.0149952 *
'* 1.7    0.3979849   0.5777652   0.2817389   0.0851499   0.0187902 *
'* 1.8    0.3399864   0.5815170   0.3061435   0.0988020   0.0231965 *
'* 1.9    0.2818186   0.5811571   0.3299257   0.1134234   0.0282535 *
'* 2.0    0.2238908   0.5767248   0.3528340   0.1289432   0.0339957 *
'* 2.1    0.1666070   0.5682921   0.3746236   0.1452767   0.0404526 *
'* 2.2    0.1103623   0.5559630   0.3950587   0.1623255   0.0476471 *
'* 2.3    0.0555398   0.5398725   0.4139146   0.1799789   0.0555957 *
'* 2.4    0.0025077   0.5201853   0.4309800   0.1981148   0.0643070 *
'* 2.5   -0.0483838   0.4970941   0.4460591   0.2166004   0.0737819 *
'* 2.6   -0.0968050   0.4708183   0.4589729   0.2352938   0.0840129 *
'* 2.7   -0.1424494   0.4416014   0.4695615   0.2540453   0.0949836 *
'* 2.8   -0.1850360   0.4097092   0.4776855   0.2726986   0.1066687 *
'* 2.9   -0.2243115   0.3754275   0.4832271   0.2910926   0.1190335 *
'* 3.0   -0.2600520   0.3390590   0.4860913   0.3090627   0.1320342 *
'********************************************************************
defint i-n
defdbl a-h,o-z

cls
print
m=20
print " X       J0(X)       J1(X)       J2(X)       J3(X)       J4(X)   "
print "-----------------------------------------------------------------"
F$ =  "#.#   ##.#######   #.#######   #.#######   #.#######   #.####### "
x=0.1
for i=1 to 30
  gosub 1000
  print using F$; x; Y(0); Y(1); Y(2); Y(3); Y(4)
  x = x + 0.1
  if i=15 then input r$
next i
print
end

'*************************************************************
'*        Integer order Bessel function subroutine           *
'* Calculates Bessel functions of order 0 to 4, for x > 0.   *
'* Argument is x, number of steps is m, returned results are *
'* Y(i).                                                     *
'* --------------------------------------------------------- *
'* Reference: Miller's method, see Henrici.                  *
'*************************************************************
1000 'Test for range
  if x<=0 then return
  if m<=0 then return
  Y(0)=1#
  Y(1)=0#
  c=0# : n=m
  'Update results
1100 for j=4 to 1 step -1
    Y(j)=Y(j-1)
  next j
  'Apply recursion relation
  Y(0) = 2#*n*Y(1)/x - Y(2)
  n = n - 1
  if n=0 then goto 1200
  if INT(n/2)<>n/2 then goto 1100
  c=c+2#*Y(0)
  goto 1100
1200 c=c+Y(0)
  'Scale the results
  for j=0 to 4
    Y(j)=Y(j)/c
  next j
return

'End of file INTBESSL.BAS