!****************************************************
!*  Program to demonstrate Bessel Function Series   *
!* ------------------------------------------------ *
!* Reference: BASIC Scientific Subroutines, Vol. II *
!* by F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981.    *
!*                                                  *
!*           F90 Version by J.-P. Moreau, Paris.    *
!*                  (www.jpmoreau.fr)               *
!* ------------------------------------------------ *
!* SAMPLE RUN:                                      *
!*                                                  *
!*  BESSEL SERIES SUMMATION                         *
!*                                                  *
!*  What is the order of the Bessel function ? 2    *
!*  Argument ? 1                                    *
!*  Convergence criterion ? 1e-6                    *
!*                                                  *
!*  J( 1 ) of order 2 =  0.11490349                 *
!*                                                  *
!*  Number of terms used:  4                        *
!*                                                  *
!*  <(e)nd to exit> _                               *                             
!****************************************************
PROGRAM Besslser

implicit none
integer m,n
real*8  e,x,y
character*3 answ

10 print *,' '
  print *,'BESSEL SERIES SUMMATION'
  print *, ' '
  write(*,'(" What is the order of the Bessel function ? ")',advance='no')
  read *, n
  write(*,'(" Argument ? ")',advance='no')
  read *, x
  write(*,'(" Convergence criterion ? ")',advance='no')
  read *, e
  
  call Bessel_Series(m,n,e,x,y)

  write(*,50) x,n,y
  write(*,60) m

  write (*,'(" <(e)nd> to exit ")',advance='no')
  read *, answ

  
  if (answ.eq.'e') goto 20
  if (answ.eq.'end') goto 20
  goto 10

20 stop

50 format(/' J(',F6.2,') of order ',I2,' = ',F12.8) 
60 format(/' Number of terms used: ',I2/)

end

!*************************************
!* 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.                *
!*************************************
Subroutine Bessel_Series(m,n,e,x,y)  
!Labels: 100, 200, 210
  implicit none
  real*8 a,b0,b1,b2,e,x,y
  integer i,m,n
  a = 1.d0
  if (n <= 1) goto 100
  ! Calculate N! 
  do i = 1, n
    a = a * dfloat(i)
  end do
100  a = 1 / a
  if (n.eq.0) goto 200
  ! Calculate multiplying term
  do i = 1, n
    a = a * x / 2.d0
  end do
200 b0 = 1.d0
  b2 = 1.d0
  m = 0
  !Assemble series sum
210 m = m + 1
  b1 = -(x * x * b0) / (m * (m + n) * 4.d0)
  b2 = b2 + b1
  b0 = b1
  ! Test for convergence
  if (dabs(b1) > e)  goto 210
  ! form final answer
  y = a * b2
  return
end

! End file besslser.f90