{*****************************************************************
*  Program to demonstrate the Bessel function asymptotic series  *
* -------------------------------------------------------------- *
*       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)        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  *
* -------------------------------------------------------------- *
*                                                                *
*****************************************************************}
PROGRAM BESSEL1;
Uses WinCrtMy, Typedef, Utilit;

Var
        e, e3, J0, J1, x  : ar_reel;
        m1, m2, n1, n2 : ar_reel;
        i,n : integer;


{*********************************************************
* 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.                                   *
*********************************************************}
PROCEDURE Calcul_Bessel;
Label 100, 200;
Var   a,a1,a2,b,c,e1,e2,e4,x1 : ar_reel;
      m : integer;
Begin
  {Calculate P and Q polynomials; P0(x):=m1, P1(x):=m2,
   Q0(x):=n1, Q1(x):=n2  }
  a := 1; a1 := 1; a2 := 1; b := 1; c := 1;
  e1 := 1e6;
  m := -1;
  x1 := 1.0 / (8.0 * x);
  x1 := x1 * x1;
  m1 := 1.0; m2 := 1.0;
  n1 := -1.0 / (8.0 * x);
  n2 := -3.0 * n1;
  n := 0;
100:  m := m + 2; a := a * m * m;
  m := m + 2; a := a * m * m;
  c := c * x1;
  a1 := a1 * a2; a2 := a2 + 1.0; a1 := a1 * a2;
  a2 := a2 + 1.0; e2 := a * c / a1;
  e4 := 1.0 + (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 200;
  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 200;
  GOTO 100;
200: a := PI;
  e := e2;
  b := SQRT(2.0 / (a * x));
  J0 := b * (m1 * COS(x - a / 4) - n1 * SIN(x - a / 4.0));
  J1 := b * (m2 * COS(x - 3 * a / 4.0) - n2 * SIN(x - 3 * a / 4.0))
End;


{main program}
BEGIN
  CLRSCR;
  e3:=1e-10;
  Writeln;
  Writeln('   X            J0(X)            J1(X)        N            E      ');
  Writeln('------------------------------------------------------------------');
  FOR i := 1 TO 15 DO
  begin
    x := i;
    Calcul_Bessel;
    Write('  ',x:2:0,'        '); Aff_reel(J0); Write('       '); Aff_reel(J1);
    Write('       ',n); Write('       '); Aff_reel(e); Writeln
  end;
  Writeln('------------------------------------------------------------------');
  Readkey; DoneWinCrt
END.

{End of file bessel1.pas}