{*******************************************************************
*           Program to demonstrate INT_BESSEL Function             *
* ---------------------------------------------------------------- *
*         Reference: BASIC Scientific Subroutines, Vol. II         *
*      By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [BIBLI 01]     *
*                                                                  *
*                       Pascal Version by J.-P. Moreau, Paris.     *
*                                (www.jpmoreau.fr)                 *
* ---------------------------------------------------------------- *
* 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 *
*******************************************************************}
PROGRAM INTBESSL;
Uses WinCrtMy, Typedef;

VAR
        i,m : integer;
        x   : ar_reel;
        Y   : Array[0..4] of ar_reel;


{************************************************************
*        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.                  *
************************************************************}
PROCEDURE Int_Bessel;
Label 100,200,fin;
Var j,n : integer;
    c   : ar_reel;
Begin
  {Test for range}
  if x<=0 then goto fin;
  if m<=0 then goto fin;
  Y[0]:=1.0;
  Y[1]:=0.0;
  c:=0.0 ; n:=m;
  {Update results}
100: for j:=4 downto 1 do
    Y[j]:=Y[j-1];
  {Apply recursion relation}
  Y[0] := 2.0*n*Y[1]/x - Y[2];
  n := n - 1;
  if n=0 then goto 200;
  if (n MOD 2)<>0 then goto 100;  {if n is odd}
  c:=c+2.0*Y[0];
  goto 100;
200: c:=c+Y[0];
  {Scale the results}
  for j:=0 to 4 do
    Y[j]:=Y[j]/c;
fin: End;


{main program}
BEGIN
  clrscr;
  m:=20;
  writeln(' X       J0(X)       J1(X)       J2(X)       J3(X)       J4(X)   ');
  writeln('-----------------------------------------------------------------');
  x:=0.1;
  for i:=1 to 30 do
  begin
    Int_Bessel;
    writeln(x:3:1,'   ',Y[0]:10:7,'   ',Y[1]:9:7,'   ',Y[2]:9:7,'   ',Y[3]:9:7,'   ',Y[4]:9:7);
    x := x + 0.1;
    if i=20 then Readkey
  end;
  ReadKey; DoneWinCrt
END.

{End of file intbessl.pas}