!******************************************************************
!*  Program to demonstrate the Bessel function asymptotic 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:                                                    *
!*                                                                *
!*  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

implicit none
real*8  e, e3, J0, J1, x
integer i, n

  e3=1.d-8
  print *, ' '
  print *,'   X            J0(X)            J1(X)        N            E      '
  print *,'------------------------------------------------------------------'
  do i = 1, 15
    x = dfloat(i)
    call Calcul_Bessel(n,e,e3,x,J0,J1)
    write(*,100) i,J0,J1,n,e
  end do
  print *,'------------------------------------------------------------------'
  stop
100 format('  ',i2,'          ',f9.6,'        ',f9.6,'       ',i1,'        ',f9.7)

end


!**********************************************************
!* 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.                                   *
!**********************************************************
Subroutine Calcul_Bessel(n,e,e3,x,J0,J1)
! Labels: 100, 200
  implicit none
  real*8  a,a1,a2,b,c,e,e1,e2,e3,e4,J0,J1,pi,x,x1
  real*8  m1, m2, n1, n2
  integer m, n
  pi = 4.d0*atan(1.d0)
  !Calculate P and Q polynomials; P0(x)=m1, P1(x)=m2,
  !Q0(x)=n1, Q1(x)=n2  
  a = 1.d0; a1 = 1.d0; a2 = 1.d0; b = 1.d0; c = 1.d0
  e1 = 1.d6
  m = -1
  x1 = 1.d0 / (8.d0 * x)
  x1 = x1 * x1
  m1 = 1.d0; m2 = 1.d0
  n1 = -1.d0 / (8.d0 * x)
  n2 = -3.d0 * 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.d0; a1 = a1 * a2
  a2 = a2 + 1.d0; e2 = a * c / a1
  e4 = 1.d0 + (m + 2) / m + (m + 2) * (m + 2) / (a2 * 8.d0 * x) + (m + 2) * (m + 4) / (a2 * 8.d0 * x);
  e4 = e4 * e2
  !Test for divergence
  if (dabs(e4) > e1) goto 200
  e1 = dabs(e2)
  m1 = m1 - e2
  m2 = m2 + e2 * (m + 2) / m
  n1 = n1 + e2 * (m + 2) * (m + 2) / (a2 * 8.d0 * x)
  n2 = n2 - e2 * (m + 2) * (m + 4) / (a2 * 8.d0 * x)
  n = n + 1
  !Test for convergence criterion
  if (e1 < e3) goto 200
  goto 100
200  a = pi
  e = e2
  b = dsqrt(2.d0 / (a * x))
  J0 = b * (m1 * dcos(x - a / 4.d0) - n1 * dsin(x - a / 4.d0))
  J1 = b * (m2 * dcos(x - 3.d0 * a / 4.d0) - n2 * dsin(x - 3.d0 * a / 4.d0))
  return
end


! End of file bessel1.f90