/*******************************************************************
*           Program to demonstrate INT_BESSEL Function             *
* ---------------------------------------------------------------- *
*         Reference: BASIC Scientific Subroutines, Vol. II         *
*      By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [BIBLI 01]     *
*                                                                  *
*                          C++ 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 *
*******************************************************************/
#include <stdio.h>
#include <math.h>

int i, m;
double x;
double Y[5];


/************************************************************
*        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.                  *
************************************************************/
void Int_Bessel()  {
  // Labels: e100,e200
  int j,n;
  double c,xn;
  // Test for range
  if (x<=0) return;
  if (m<=0) return;
  Y[0]=1.0;
  Y[1]=0.0;
  c=0.0; 
  n=m;
  // Update results
e100: for (j=4; j>0; j--) {
	Y[j]=Y[j-1];
  }
  // Apply recursion relation
  xn = (double) n;
  Y[0] = 2.0*xn*Y[1]/x - Y[2];
  n--;
  if (n==0) goto e200;
  if (n % 2) goto e100; // if n is odd
  c += 2*Y[0];
  goto e100;
e200: c += Y[0];
  // Scale the results
  for (j=0; j<5; j++) {
	Y[j]=Y[j]/c;
  }    
}


void main()  {
  m=20;
  printf(" X       J0(X)       J1(X)       J2(X)       J3(X)       J4(X)   \n");
  printf("-----------------------------------------------------------------\n");
  x = (double) 0.1;
  for (i=1; i<31; i++) {
    Int_Bessel();
    printf("%3.1f   %10.7f   %9.7f   %9.7f   %9.7f   %9.7f\n",x,Y[0],Y[1],Y[2],Y[3],Y[4]);
    x += (double) 0.1;
    if (i==20) getchar();
  }
}


// End of file INTBESSL.CPP