/***************************************************
*      Program to demonstrate Bessel Series        *
*              Summation Function                  * 
* ------------------------------------------------ *
* Reference: BASIC Scientific Subroutines, Vol. II *
* by F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981.    *
*                                                  *
*           C++ Version by J.-P. Moreau, Paris.    *
* ------------------------------------------------ *
* SAMPLE RUN:                                      *
*                                                  *
*  BESSEL SERIES SUMMATION                         *
*                                                  *
*  What is the order of the Bessel function ? 2    *
*  Argument ? 1                                    *
*  Convergence criterion ? 1e-6                    *
*                                                  *
*  J(1.000000) of order 2 =  0.114903              *
*                                                  *
*  Number of terms used:  4                        *
*                                                  *
***************************************************/
#include <stdio.h>
#include <math.h>

int    m,n;
double e,x,y;


/************************************
* Bessel function series subroutine *
* The order is n, the argument x.   *
* The returned value is in y.       *
* The number of terms used is in m. *
* e is the convergence criterion    *
* Example: e = 1e-6.                *
************************************/
void Bessel_Series()  {
//Labels: e100, e200, e210
  double a,b0,b1,b2;
  int i;
  a = 1;
  if (n <= 1) goto e100;
  // Calculate N! 
  for (i = 1; i < n+1; i++)  a *= i;
e100: a = 1 / a;
  if (n == 0) goto e200;
  // Calculate multiplying term
  for (i = 1; i < n+1; i++)  a *= (x / 2.0);
e200: b0 = 1.0;
  b2 = 1.0; m = 0;
  //Assemble series sum
e210: m++;
  b1 = -(x * x * b0) / (m * (m + n) * 4);
  b2 = b2 + b1;
  b0 = b1;
  // Test for convergence
  if (fabs(b1) > e)  goto e210;
  // form final answer
  y = a * b2;
}


void main()  {

  printf("\n BESSEL SERIES SUMMATION\n\n");
  printf(" What is the order of the Bessel function ? "); scanf("%d",&n);
  printf(" Argument ? "); scanf("%lf",&x);
  printf(" Convergence criterion ? "); scanf("%lf",&e);
  
  Bessel_Series();

  printf("\n J(%f) of order %d = %f\n\n",x,n,y); 
  printf(" Number of terms used: %d\n\n",m);

}

// End file besslser.cpp