program bes2
c
c -- test the bessel function of the second kind
c -- figure 11.17
c
      integer out
      real x, ans, order
c
      out = 6
10    write(out, 101)
      read(*, 102) order
      if (order .lt. 0.0) goto 99
20    write(out, 103)
      read(*, 102) x
      if (x .le. 0.0) goto 20
      ans = besy(x, order)
      write(out, 104) ans
      goto 10
99    stop
101   format(/' order? ' )
102   format(e10.0)
103   format('  arg? ' )
104   format(' y bessel =', f10.4)
      end
      function besy(x, ord)
c
c -- bessel function of the second kind
c
      integer j
      real x, ord, pi, ts, sum, x2, pi2, xx
      real t, y0, y1, ya, yb, yc, small, euler, fj1
      data pi/3.141593/, pi2/0.6366198/
      data small/1.0e-8/, euler/0.5772157/
c
      if (x .ge. 12.0) goto 80
      xx = 0.5 * x
      x2 = xx * xx
      t = alog(xx) + euler
c
c -- calculation of y0
c
      sum = 0.0
      term = t
      y0 = t
      j = 0
60    j = j + 1
       if(j .ne. 1) sum = sum + 1.0/(j - 1)
       ts = t - sum
       term = -x2 * term / (j*j) * (1.0 - 1.0/(j*ts))
       y0 = y0 + term
      if (abs(term) .ge. small) goto 60
      y0 = pi2 * y0
      if (ord .eq. 0.0) goto 90
c
c -- calculation of y1 if needed
c
      term = xx * (t - 0.5)
      sum = 0.0
      y1 = term
      j = 1
70     j = j + 1
       fj1 = j - 1
       sum = sum + 1/(fj1)
       ts = t - sum
       t2 = (-x2*term) / (j * (fj1))
       term = t2 *((ts - 0.5/j) / (ts + 0.5/(fj1)))
       y1 = y1 + term
      if (abs(term) .ge. small) goto 70
      y1 = pi2 * (y1 - 1/x)
      if (ord .eq. 1.0) goto 95
c
c -- recursive calculation from y0 and y1
c
      ts = 2.0 / x
      ya = y0
      yb = y1
      jd = aint(ord + 0.01)
      do 75 j = 2, jd
       yc = ts * (j-1) * yb -ya
       ya = yb
       yb = yc
75    continue
      besy = yc
      return
c
c -- asymptotic expansion for large x
c
80    besy = sqrt(2/(pi*x))*sin(x - pi/4 - ord*pi/2)
      return
90    besy = y0
      return
95    besy = y1
      return
      end