!********************************************************************
!*           Program to demonstrate INT_BESSEL Function             *
!* ---------------------------------------------------------------- *
!*         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)       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 *
!********************************************************************
PROGRAM INTBESSL

implicit none
integer i, m
real*8 x
real*8 Y(0:4)

  m=20
  print *,' X       J0(X)       J1(X)       J2(X)       J3(X)       J4(X)  '
  print *,'-----------------------------------------------------------------'
  x = 0.1d0
  do i = 1, 30
    call Int_Bessel(m,x,Y)
    write(*,50) x,Y(0),Y(1),Y(2),Y(3),Y(4)
    x = x + 0.1d0
    if (i.eq.20) pause '<Enter> to continue...'
  end do
50 format(F3.1,'   ',F10.7,'   ',F9.7,'   ',F9.7,'   ',F9.7,'   ',F9.7)
  stop
end


!*************************************************************
!*        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.                  *
!*************************************************************
Subroutine Int_Bessel(m,x,Y)  
! Labels: 100, 200
  implicit none
  integer j,m,n
  real*8 c,x,xn
  real*8 Y(0:4)
  ! 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
100 do j =4, 1, -1
    Y(j)=Y(j-1)
  end do
  ! Apply recursion relation
  xn = n
  Y(0) = 2.d0*xn*Y(1)/x - Y(2)
  n = n - 1
  if (n.eq.0) goto 200
  if (int(xn/2).ne.n/2) goto 100
  c = c + 2*Y(0)
  goto 100
200 c = c + Y(0)
  ! Scale the results
  do j = 0, 4
    Y(j)=Y(j)/c
  end do
  return
end

! End of file intbessl.f90