program bes1
c
c -- test the bessel function of the first kind
c -- function gamma is required
c -- figure 11.14
c
      integer out
      real x, ans, order
c
      out = 6
10    write(out, 101)
      read(*, 102) order
      if (order .lt. -25.0) goto 99
      write(out, 103)
      read(*, 102) x
      ans = besj(x, order)
      write(out, 104) ans
      goto 10
99    stop
101   format(/' order? ' )
102   format(e10.0)
103   format('  arg? ' )
104   format(' j bessel =', f10.4)
      end
      function besj(x, ord)
c
c -- bessel function of the first kind
c -- gamma function is required
c
      integer i
      real x, ord, tol, pi, term, term2, sum, x2
      data tol/1.0e-5/, pi/3.141593/
c
      x2 = x * x
      besj = 0.0
      if ((x .eq. 0.0) .and. (ord .eq. 1.0)) goto 99
      if (x .gt. 15.0) goto 80
      if (ord .eq. 0.0) sum = 1.0
      if (ord .ne. 0.0) 
     *  sum = exp(ord * alog(x/2)) / gamma(ord + 1)
      term2 = sum
      i = 0
60    i = i + 1
       term = term2
       term2 = -term * x2 * 0.25/(i * (ord + i))
       sum = sum + term2
      if (abs(term2) .gt. abs(sum * tol)) goto 60
      besj = sum
      goto 99
c
c -- asymptotic expansion for large x
c
80    besj = sqrt(2/(pi*x))*cos(x - pi/4 - ord*pi/2)
99    return
      end