c-----------------------------------------------------------------------
c Approximate the Bessel function of the first kind of order 0
c via Taylor's series expansion.
c Link with SFUN.LIB (IMSL's Special Functions)
c          +MATH.LIB (IMSL's QDAGS for integration)
c Instructor: Nam Sun Wang
c-----------------------------------------------------------------------

c Print out a program header -------------------------------------------
      print *, 'Approximate the Bessel function via Taylor''s series ',
     &         'expansion, which fails for x far away from the point ',
     &         'of expansion around x=0.  This program is accurate to ',
     &         'about 5 significant figures for 0<x<10.'
      print *, ' '

c Input data -----------------------------------------------------------
      print *, 'Enter x: '
      read *, x

c Print out the results ------------------------------------------------
      print *, ' '
      print *, 'Taylor''s series expansion gives: ', besJ0(x)
      print *, 'Integral definition gives:       ', besintJ0(x)
      print *, 'Compare to IMSL''s J0(x) =        ', BSJ0(x)

      end



c-----------------------------------------------------------------------
      function besJ0(x)
c-----------------------------------------------------------------------
c Approximate the Bessel's function of the first kind of order zero
c via 10-term Taylor's series expansion.
c   J0(x) = sum(sign*x^i/i!) for i=even integers
c   J0(x) = 1 - (x/2)**2/(1!)**2 + (x/2)**4/(2!)**2 - (x/2)**6/(3!)**2 + ...
c Method used: iterate on the following relationship based on Taylor's series:
c   delta = -delta*(x/2)**2/i**2
c   bessel = bessel + delta
c Programming Note:
c   Because Taylor's series is alternating in sign, there is subtractive
c   cancellation.  Nevertheless, convergence is quick enough for a small x.
c   For large values of x, apply the asymptotic approximation formula.
c     J0(x)=sqrt(2/pi/x)*cos(x-pi/4),  for x>25
c Instructor: Nam Sun Wang
c-----------------------------------------------------------------------

c Calculate the Taylor's series expansion
c     Initialize  (The first term is contained in the initialization step.)
      delta = 1.
      besJ0 = 1.

      do i=1, 99999
        delta = -delta*(x/2.)*(x/2.)/float(i)/float(i)
        besJ0 = besJ0 + delta
c       Converge to 5 digits.
c       Since J0 is O(1), this is almost identical to 5 significant digits.
        if(abs(delta) .lt. 1.E-5) exit
      enddo

      end



c-----------------------------------------------------------------------
      function besintJ0(x)
c-----------------------------------------------------------------------
c Approximate the Bessel's function of the first kind of order zero
c via the following integral definition.
c   J0(x) = (1/pi)*integral (from 0 to pi) of cos(x*sin(q))*dq
c Programming Note:
c   For large values of x, apply the asymptotic approximation formula.
c     J0(x)=sqrt(2/pi/x)*cos(x-pi/4),  for x>25
c Instructor: Nam Sun Wang
c-----------------------------------------------------------------------

      external bintegrand
      common xx
      pi = 3.141593

c Pass x into the integrand through "common"
      xx = x

c Call an integration routine to do the grunt work.
c The integral is returned in xintegral
      call QDAGS(bintegrand, 0, pi, errabs, 1.e-4, xintegral, errest)
      besintJ0 = 1./pi*xintegral

      end



c-----------------------------------------------------------------------
      function bintegrand(q)
c-----------------------------------------------------------------------
c Specify the integrand function
c-----------------------------------------------------------------------
      common xx

      bintegrand = cos(xx*sin(q))

      end