SUBROUTINE MFCVAL ( N, R, QQ, CV, J )

c*********************************************************************72
C
CC MFCVAL computes characteristic values of Mathieu's differential equation.
C
C  Modified:
C
C    06 October 2006
C
C  Author:
C
C    Donald Clemm
C
C  Reference:
C
C    Donald Clemm,
C    Algorithm 352: Characteristic Values and Associated 
C    Solutions of Matthieu's Differential Equation,
C    Communications of the ACM,
C    Volume 12, Number 7, pages 399-407, June 1969.
C
C  Parameters:
C
C    Input, integer N, the number of characteristic values desired.
C
C    Input, integer R, given as N-1 or N according as the characteristic
C    values are to be associated with the even or odd solutions, respectively.
C
C    Input, double precision QQ, the nonnegative parameter.
C
C    Output, double precision CV(6,N), the computed array of characteristic
C    values and bounds.
C
C    Output, integer J, the number of characteristic values successfully
C    computed.  If J is not equal to N then an error occurred while trying 
C    to compute value J+1.
C
      implicit none

      integer n

      double precision a
      double precision a0
      double precision a1
      double precision a2
      double precision cv(6,n)
      double precision dl
      double precision dr
      double precision dtm
      integer j
      integer k
      integer kk
      integer l
      integer m
      integer m0
      integer m1
      integer m2s
      integer mf
      double precision q
      double precision qq
      integer r
      double precision t
      double precision tm
      double precision tol
      double precision tola
      integer type

      equivalence ( dl, dr, t )

      common / mf1 / q, tol, type, m1, m0, m2s, mf

      common / mf2 / a0, a2, a1

      tol = 1.0d-13

      if ( n - r ) 10, 10, 20

10    l = 1
      go to 30

20    l = 2

30    q = qq

      do 500 k = 1, n

        j = k

        if ( q .eq. 0.0d+00 ) then
          go to 490
        end if

        if ( q .lt. 0.0d+00 ) then
          write ( 6, 961 )
          q = -q
        end if

        kk = min0 ( k, 4 )
        type = 2 * mod ( l, 2 ) + mod ( k - l + 1, 2 )
c
c  first approximation.
c
        go to ( 100, 200, 300, 400 ), kk

100     if ( q - 1.0d+00 ) 110, 140, 140

110     go to ( 120, 130 ), l

120     a = 1.0d+00 - q - 0.125d+00 * q * q
        go to 420

130     a = q * q
        a = a * ( -0.5d+00 + 0.0546875d+00 * a )
        go to 420

140     if ( q - 2.0d+00 ) 150, 180, 180

150     go to ( 160, 170 ), l

160     a = 1.033d+00 - 1.0746d+00 * q - 0.0688d+00 * q * q
        go to 420

170     a = 0.23d+00 - 0.495d+00 * q - 0.191d+00 * q * q
        go to 420

180     a = -0.25d+00 - 2.0d+00 * q + 2.0d+00 * dsqrt ( q )
        go to 420

200     dl = l
        if ( q * dl - 6.0d+00 ) 210, 350, 350

210     go to ( 220, 230 ), l

220     a = 4.01521d+00 - q * ( 0.046d+00 + 0.0667857d+00 * q )
        go to 420

230     a = 1.0d+00 + 1.05007d+00 * q - 0.180143d+00 * q * q
        go to 420

300     if ( q - 8.0d+00 ) 310, 350, 350

310     go to ( 320, 330 ), l

320     a = 8.93867d+00 + 0.178156d+00 * q - 0.0252132d+00 * q * q
        go to 420

330     a = 3.70017d+00 + 0.953485d+00 * q - 0.0475065d+00 * q * q
        go to 420

350     dr = k - 1
        a = cv(1,k-1) - dr + 4.0d+00 * dsqrt ( q )
        go to 420

400     a = cv(1,k-1) - cv(1,k-2)
        a = 3.0d+00 * a + cv(1,k-3)

420     if ( q .ge. 1.0d+00 ) go to 440
        if ( k .ne. 1 ) go to 430

        tola = dmax1 ( dmin1 ( tol, dabs ( a ) ), 1.0d-14 )
        go to 450

430     tola = tol * dabs ( a )
        go to 450

440     tola = tol * dmax1 ( q, dabs ( a ) )

        tola = dmax1 ( 
     &    dmin1 ( tola, dabs ( a ), 0.4d+00 * dsqrt ( q ) ), 
     &    1.0d-14 )
c
c  crude upper and lower bounds.
c
450     continue

        call bounds ( k, a, tola, cv, n, m )
        if ( m .ne. 0 ) if ( m - 1 ) 470, 910, 900
c
c  iterate.
c
        call mfitr8 ( tola, cv(1,k), cv(2,k), m )
        if ( m .gt. 0 ) go to 920
c
c  final bounds and functions.
c
470     t = cv(1,k) - tola
        call tmofa ( t, tm, dtm, m )

        if ( m .gt. 0 ) then
          write ( 6, 941 ) k
          cv(3,k) = 0.0d+00
          cv(4,k) = 0.0d+00
        else
          cv(3,k) = t
          cv(4,k) = -tm / dtm
        end if

        t = cv(1,k) + tola
        call tmofa ( t, tm, dtm, m )
        if ( m .gt. 0 ) go to 950
        cv(5,k) = t
        cv(6,k) = -tm / dtm
        go to 500
c
c  q equals zero.
c
490     cv(1,k) = ( k - l + 1 )**2
        cv(2,k) = 0.0d+00
        cv(3,k) = cv(1,k)
        cv(4,k) = 0.0d+00
        cv(5,k) = cv(1,k)
        cv(6,k) = 0.0d+00

500   continue

550   return
c
c  print error messages.
c
900   write ( 6, 901 ) k
901   format ( 20h0crude bounds cannot,
     &  22h be located, no output,
     &  7h for k=, i2 )
      go to 930

910   write ( 6, 911 ) k
911   format ( ' error in subprogram tmofa, via subprogram',
     &  ' bounds, no output, for k = ', i2 )
      go to 930

920   write ( 6, 921 ) k
921   format ( ' error in subprogram tmofa, via subprogram',
     &  ' mfitr8, no output for k = ', i2 )

930   j = j - 1
      go to 550
941   format ( ' error in subprogram tmofa, no lower bound',
     &  ' for k = ', i2 )

950   write ( 6, 951 ) k
951   format ( ' error in subprogram tmofa, no upper bound',
     &  ' for k = ', i2 )
      cv(5,k) = 0.0d+00
      cv(6,k) = 0.0d+00
      go to 500
961   format ( 20h0q given negatively,,
     &  20h used absolute value )
      end
      subroutine math ( xx, qq, r, cv, sol, fnc, norm, f, k, m )
 
c*********************************************************************72

      implicit none

      double precision a
      double precision ab(200)
      double precision cv
      external dc
      double precision dc
      external ddc
      double precision ddc
      external dds
      double precision dds
      double precision dlast
      double precision dmax
      external dpc
      double precision dpc
      external dps
      double precision dps
      external ds
      double precision ds
      real dum1(578)
      real dum2(6)
      double precision f(3)
      integer fnc
      double precision g
      integer i
      double precision j(250)
      integer k(2)
      integer klast
      integer kmax
      integer l
      integer ll
      integer m
      integer mf
      integer ml
      integer mm
      integer m0
      integer m1
      integer m2s
      integer n
      integer norm
      integer p
      external pc
      double precision pc
      external ps
      double precision ps
      double precision q
      double precision qq
      integer r
      integer s
      integer sol
      double precision t
      double precision tol
      integer type
      double precision u1
      double precision u2
      double precision x
      double precision xx
      double precision y(250)

      common j, y, u1, u2, n, p, s, l, x, t, i, ll, g, dmax, dlast,
     &  kmax, klast, dum1, a, dum2, mm, ml, ab

      common / mf1 / q, tol, type, m1, m0, m2s, mf

      m = 0
      if ( sol .lt. 1 .or. 
     &  sol .gt. 3 .or. 
     &  fnc .lt. 1 .or. 
     &  fnc .gt. 4 ) go to 400

      a = cv
      q = qq
      tol = 1.0d-13
      type = 2 * mod ( fnc, 2 ) + mod ( r, 2 )
      call coef ( m )

      if ( m ) 410, 10, 420

10    n = r / 2
      p = mod ( r, 2 )
      s = mm / 2
      l = ml / 2
      x = xx
      t = 1.0d+00

      if ( sol .eq. 3 ) go to ( 150, 160, 170, 180 ), fnc

      u1 = dsqrt ( q ) * dexp ( -x )
      u2 = q / u1
      ll = l + s + p
c
c  compute bessel functions.
c
      call bessel ( 1, u1, j, ll )
      call bessel ( sol, u2, y, ll )
c
c  evaluate selected function.
c
      go to ( 50, 60, 70, 80 ), fnc

50    call sum_series ( ds )
      go to 300
60    call sum_series ( dc )
      go to 300
70    call sum_series ( dds )
      go to 300
80    call sum_series ( ddc )
      go to 300
150   call sum_series ( ps )
      go to 200
160   call sum_series ( pc )
      go to 200
170   call sum_series ( dps )
      go to 200
180   call sum_series ( dpc )

200   if ( norm - 2 ) 300, 210, 250
c
c  ince normalization.
c
210   t = ab(1)**2
      if ( type .eq. 0 ) t = t + t
      do 220 i = 1, l
        t = t + ab(i+1)**2
220   continue
      t = dsqrt ( t )
      i = m0 / 2
      if ( ab(i) .lt. 0.0d+00 ) t = -t
      go to 300
c
c  stratton normalization.
c
250   if ( type .gt. 1 ) go to 270
      t = ab(1)
      do 260 i = 1, l
        t = t + ab(i+1)
260   continue
      go to 300
270   t = dble ( float ( p ) ) * ab(1)
      do 280 i = 1, l
        t = t + ab(i+1) * dble ( float ( 2 * i + p ) )
280   continue

300   f(1) = g / t
      f(2) = dmax / t
      f(3) = dlast / t
      k(1) = kmax
      k(2) = klast
350   return
c
c  print error messages.
c
400   write ( 6, 401 )
401   format ( 18h0sol or fnc out of,
     &  17h range, no output )
      go to 450
410   write ( 6, 411 )
411   format ( 15h0more than 200 ,
     &  22hcoefficients required,,
     &  20h qq and r too large,,
     &  10h no output )
      go to 450
420   write ( 6, 421 )
421   format ( 20h0error in subprogram,
     &  22h tmofa, via subprogram,
     &  13h coef, verify,
     &  21h arguments, no output )

450   m = 1
      f(1) = 0.0d+00
      f(2) = 0.0d+00
      f(3) = 0.0d+00
      k(1) = 0
      k(2) = 0
      go to 350
      end
      subroutine bounds ( k, approx, tola, cv, n, mm )

c*********************************************************************72

      implicit none

      integer n

      double precision a
      double precision a0
      double precision a1
      double precision approx
      double precision cv(6,n)
      double precision d0
      double precision d1
      double precision dtm
      integer k
      integer ka
      integer m
      integer m0
      integer m1
      integer m2s
      integer mf
      integer mm
      double precision q
      double precision tm
      double precision tol
      double precision tola
      integer type

      common / mf1 / q, tol, type, m1, m0, m2s, mf

      common / mf2 / a0, a, a1

      ka = 0
      if ( k .eq. 1 ) go to 20

      if ( approx - cv(1,k-1) ) 10, 10, 20

10    a0 = cv(1,k-1) + 1.0d+00
      go to 30

20    a0 = approx

30    call tmofa ( a0, tm, dtm, m )
      if ( m .gt. 0 ) go to 250
      d0 = -tm / dtm
      if ( d0 ) 100, 300, 50
c
c  a0 is lower bound, search for upper bound.
c
50    a1 = a0 + d0 + 0.1d+00
      call tmofa ( a1, tm, dtm, m )

      if ( m .gt. 0 ) go to 250
      d1 = -tm / dtm

      if ( d1 ) 200, 350, 60

60    a0 = a1
      d0 = d1
      ka = ka + 1

      if ( ka - 4 ) 50, 400, 400
c
c  a1 is upper bound, search for lower bound.
c
100   a1 = a0
      d1 = d0
      a0 = dmax1 ( a1 + d1 - 0.1d+00, -2.0d+00 * q )
      if ( k .eq. 1 ) go to 110
      if ( a0 - cv(1,k-1) ) 150, 150, 110

110   call tmofa ( a0, tm, dtm, m )
      if ( m .gt. 0 ) go to 250
      d0 = -tm / dtm
      if ( d0 ) 120, 300, 200

120   ka = ka + 1
      if ( ka - 4 ) 100, 400, 400

150   ka = ka + 1
      if ( ka - 4 ) 160, 400, 400

160   a0 = a1 + dmax1 ( tola, dabs ( d1 ) )
      go to 30

200   a = 0.5d+00 * ( a0 + d0 + a1 + d1 )
      if ( a .le. a0 .or. a .ge. a1 ) a = 0.5d+00 * ( a0 + a1 )

250   mm = m
      return

300   cv(1,k) = a0

310   cv(2,k) = 0.0d+00
      m = -1
      go to 250

350   cv(1,k) = a1
      go to 310

400   m = 2
      go to 250
      end
      subroutine mfitr8 ( tola, cv, dcv, mm )

c*********************************************************************72

      implicit none

      double precision a
      double precision a0
      double precision a1
      double precision a2
      double precision cv
      double precision d
      double precision dcv
      double precision dtm
      logical last
      integer m
      integer mm
      integer n
      double precision tm
      double precision tola

      common / mf2 / a0, a, a1

      n = 0
      last = .false.

50    n = n + 1

      call tmofa ( a, tm, dtm, m )
      if ( m .gt. 0 ) go to 400
      d = -tm / dtm
c
c  is tolerance met?
c
      if ( n .eq. 40 .or. 
     &  a - a0 .le. tola .or.
     &  a1 - a .le. tola .or.
     &  dabs ( d ) .lt. tola ) last = .true.
      if ( d ) 110, 100, 120

100   cv = a
      dcv = 0.0d+00
      go to 320
c
c  replace upper bound by a.
c
110   a1 = a
      go to 200
c
c  replace lower bound by a.
c
120   a0 = a

200   a2 = a + d

      if ( last ) go to 300
      if ( a2 .gt. a0 .and. a2 .lt. a1 ) go to 250

      a = 0.5d+00 * ( a0 + a1 )
      go to 50

250   a = a2
      go to 50

300   if ( a2 .le. a0 .or. a2 .ge. a1 ) go to 350

      call tmofa ( a2, tm, dtm, m )

      if ( m .gt. 0 ) go to 400

      d = - tm / dtm
      cv = a2

310   dcv = d

320   mm = m
      return

350   cv = a
      go to 310

400   cv = 0.0d+00
      dcv = 0.0d+00
      go to 320
      end
      subroutine tmofa ( alfa, tm, dtm, nd )

c*********************************************************************72

      implicit none

      double precision a(3)
      double precision aa
      double precision alfa
      double precision b(3)
      double precision dg(200,2)
      double precision dtm
      double precision dtype
      double precision f
      double precision fl
      double precision g(200,2)
      double precision h(200)
      double precision hp
      integer k
      integer kk
      integer kt
      integer l
      integer mf
      integer m0
      integer m1
      integer m2s
      integer nd
      integer type
      double precision q
      double precision q1
      double precision q2
      double precision qinv
      double precision t
      double precision tm
      double precision tol
      double precision tt
      double precision v

      common g, dg, aa, a, b, dtype, qinv, q1, q2, t, tt, k, l, kk, kt

      common / mf1 / q, tol, type, m1, m0, m2s, mf

      equivalence ( h(1), g(1,1) )
      equivalence ( q1, hp )
      equivalence ( q2, f )

      data fl / 1.0d+30 /
c
c  statement function.
c
      v ( k ) = ( aa - dble ( float ( k ) )**2 ) / q

      nd = 0
      kt = 0
      aa = alfa
      dtype = type
      qinv = 1.0d+00 / q
      do 10 l = 1, 2
        do 5 k = 1, 200
          g(k,l) = 0.0d+00
          dg(k,l) = 0.0d+00
5       continue
10    continue

      if ( mod ( type, 2 ) ) 20, 30, 20

20    m0 = 3
      go to 40

30    m0 = type + 2
40    k = int ( 0.5d+00 
     &  * dsqrt ( dmax1 ( 3.0d+00 * q + aa, 0.0d+00 ) ) )
      m2s = min0 ( 2 * k + m0 + 4, 398 + mod ( m0, 2 ) )
c
c  evaluation of the tail of a continued fraction.
c
      a(1) = 1.0d+00
      a(2) = v(m2s+2)
      b(1) = v(m2s)
      b(2) = a(2) * b(1) - 1.0d+00
      q1 = a(2) / b(2)

      do 50 k = 1, 200
        mf = m2s + 2 + 2 * k
        t = v(mf)
        a(3) = t * a(2) - a(1)
        b(3) = t * b(2) - b(1)
        q2 = a(3) / b(3)

        if ( dabs ( q1 - q2 ) .lt. tol ) go to 70

        q1 = q2
        a(1) = a(2)
        a(2) = a(3)
        b(1) = b(2)
        b(2) = b(3)
50    continue

      kt = 1

70    t = 1.0d+00 / t
      tt = - t * t * qinv
      l = mf - m2s

      do 80 k = 2, l, 2
        t = 1.0d+00 / ( v(mf-k) - t )
        tt = t * t * ( tt - qinv )
80    continue

      kk = m2s / 2 + 1

      if ( kt .eq. 1 ) q2 = t

      g(kk,2) = 0.5d+00 * ( q2 + t )
      dg(kk,2) = tt
c
c  stage 1.
c
      g(2,1) = 1.0d+00

      do 140 k = m0, m2s, 2

        kk = k / 2 + 1
        if ( k .lt. 5 ) if ( k - 3 ) 100, 110, 120
        g(kk,1) = v(k-2) - 1.0d+00 / g(kk-1,1)
        dg(kk,1) = qinv + dg(kk-1,1) / g(kk-1,1)**2
        go to 130

100     g(2,1) = v ( 0 )
        dg(2,1) = qinv
        go to 130

110     g(2,1) = v ( 1 ) + dtype - 2.0d+00
        dg(2,1) = qinv
        go to 130


120     g(3,1) = v ( 2 ) + ( dtype - 2.0d+00 ) / g(2,1)
        dg(3,1) = qinv + ( 2.0d+00 - dtype ) * dg(2,1) / g(2,1)**2

        if ( type .eq. 2 ) g(2,1) = 0.0d+00

130     if ( dabs ( g(kk,1) ) .lt. 1.0d+00 ) go to 200

140   continue
c
c  backtrack.
c
      tm = g(kk,2) - g(kk,1)
      dtm = dg(kk,2) - dg(kk,1)
      m1 = m2s
      kt = m2s - m0

      do 180 l = 2, kt, 2

        k = m2s - l
        kk = k / 2 + 1
        g(kk,2) = 1.0d+00 / ( v(k) - g(kk+1,2) )
        dg(kk,2) = -g(kk,2)**2 * ( qinv - dg(kk+1,2) )
        if ( k - 2 ) 150, 150, 160

150     g(2,2) = 2.0d+00 * g(2,2)
        dg(2,2) = 2.0d+00 * dg(2,2)

160     t = g(kk,2) - g(kk,1)
        if ( dabs ( t ) - dabs ( tm ) ) 170, 180, 180

170     tm = t
        dtm = dg(kk,2) - dg(kk,1)
        m1 = k

180   continue

      go to 320
c
c  stage 2.
c
200   continue

      m1 = k
      k = m2s
      kk = k / 2 + 1

210   continue

      if ( k .eq. m1 ) if ( k - 2 ) 300, 300, 310
      k = k - 2
      kk = kk - 1
      t = v(k) - g(kk+1,2)
      if ( dabs ( t ) - 1.0d+00 ) 250, 220, 220

220   continue

      g(kk,2) = 1.0d+00 / t
      dg(kk,2) = ( dg(kk+1,2) - qinv ) / t**2
      go to 210
c
c  stage 3.
c
250   continue

      if ( k .eq. m1 ) if ( t ) 220, 290, 220
      hp = dg(kk+1,2) - qinv

      g(kk,2) = fl
      h(kk) = t
      k = k - 2
      kk = kk - 1
      f = v(k) * t - 1.0d+00
      if ( k .eq. m1 ) if ( f ) 280, 290, 280
      if ( dabs ( f ) - dabs ( t ) ) 270, 280, 280

270   continue

      hp = hp / t**2 - qinv
      t = f / t
      go to 270

280   g(kk,2) = t / f
      dg(kk,2) = ( hp - qinv * t * t ) / f**2
      go to 210

290   nd = 1
      go to 320
c
c  chaining m equals 3.
c
300   g(2,2) = 2.0d+00 * g(2,2)
      dg(2,2) = 2.0d+00 * dg(2,2)

310   tm = g(kk,2) - g(kk,1)
      dtm = dg(kk,2) - dg(kk,1)

320   return
      end
      subroutine coef ( m )

c*********************************************************************72

      implicit none

      double precision a
      double precision ab(200)
      real dum1(800)
      double precision fl
      double precision g(200,2)
      double precision h(200)
      integer k
      integer ka
      integer kb
      integer kk
      integer m
      integer mf
      integer ml
      integer mm
      integer m0
      integer m1
      integer m2s
      double precision q
      double precision t
      double precision tol
      integer type
      double precision v
      double precision v2

      save fl
      save v2

      common g, dum1, a, t, k, ka, kb, kk, mm, ml, ab

      common / mf1 / q, tol, type, m1, m0, m2s, mf

      equivalence ( h(1), g(1,1) )

      data fl / 1.0d+30 /
      data v2 / 1.0d-15 /
c
c  statement function:
c
      v ( k ) = ( a - dble ( float ( k ) )**2 ) / q

      call tmofa ( a, t, t, m )

      if ( m .ne. 0 ) go to 300

      do 60 k = 1, 200
        ab(k) = 0.0d+00
60    continue

      ka = m1 - m0 + 2

      do 90 k = 2, ka, 2
        kk = ( m1 - k ) / 2 + 1
        if ( k - 2 ) 70, 70, 80
70      ab(kk) = 1.0d+00
        go to 90
80      ab(kk) = ab(kk+1) / g(kk+1,1)
90    continue

      ka = 0

      do 130 k = m1, m2s, 2

        kk = k / 2 + 1
        ml = k
        if ( g(kk,2) .eq. fl ) go to 100
        ab(kk) = ab(kk-1) * g(kk,2)
        go to 110

100     t = ab(kk-2)
        if ( k .eq. 4 .and. m1 .eq. 2 ) t = t + t
        ab(kk) = t / ( v(k-2) * h(kk) - 1.0d+00 )

110     if ( dabs ( ab(kk) ) .ge. 1.0d-17 ) ka = 0
        if ( ka .eq. 5 ) go to 260
        ka = ka + 1

130   continue

      t = dlog ( dabs ( ab(kk) ) / v2 ) 
     &  / dlog ( 1.0d+00 / dabs ( g(kk,2) ) )
      ka = 2 * idint ( t )
      ml = ka + 2 + m2s
      if ( ml .gt. 399 ) go to 400
      kb = ka + 2 + mf
      t = 1.0d+00 / v(kb)
      kk = mf - m2s

      do 150 k = 2, kk, 2
        t = 1.0d+00 / ( v(kb-k) - t )
150   continue

      kk = ml / 2 + 1
      g(kk,2) = t

      do 200 k = 2, ka, 2
        kk = ( ml - k ) / 2 + 1
        g(kk,2) = 1.0d+00 / ( v(ml-k) - g(kk+1,2) )
200   continue

      ka = m2s + 2
      do 250 k = ka, ml, 2
        kk = k / 2 + 1
        ab(kk) = ab(k-1) * g(kk,2)
250   continue
c
c  neutral normalization.
c
260   t = ab(1)
      mm = mod ( type, 2 )
      ka = mm + 2
      do 280 k = ka, ml, 2
        kk = k / 2 + 1
        if ( dabs ( t ) - dabs ( ab(kk) ) ) 270, 280, 280
270     t = ab(kk)
        mm = k
280   continue

      do 290 k = 1, kk
        ab(k) = ab(k) / t
290   continue

300   return
400   m = -1
      go to 300
      end
      subroutine sum_series ( dum )

c*********************************************************************72
c
c  to avoid conflict with the fortran90 intrinsic sum function,
c  this routine, originally named "sum", has been renamed "sum_series".
c
      implicit none

      double precision dlast
      double precision dmax
      external dum
      double precision dum
      real dum1(1006)
      real dum2(6)
      double precision f
      integer k
      integer klast
      integer kmax
      integer l
      integer s
      double precision t

      common dum1, s, l, dum2, f, dmax, dlast, kmax, klast, t

      k = 0
      f = dum ( 0 )
      dmax = f
      t = dabs ( f )
      kmax = 0
      do 30 klast = 1, l
        dlast = dum ( klast )
        f = f + dlast
        if ( t - dabs ( dlast ) ) 10, 10, 20
10      dmax = dlast
        t = dabs ( dmax )
        kmax = klast
20      if ( klast .le. s ) go to 30
        if ( dabs ( dlast ) / t .gt. 1.0d-13 ) k = 0
        k = k + 1
        if ( k .eq. 3 ) go to 40
30    continue

      klast = l

40    return
      end
      subroutine bessel ( sol, u, jy, n )

c*********************************************************************72
c
      implicit none

      integer n

      double precision jy(n)
      integer k
      integer nn
      integer sol
      double precision u

      nn = min0 ( n, 249 )

      if ( u .eq. 0.0d+00 .and. sol .eq. 2 ) go to 80
      if ( u .ge. 8.0d+00 ) go to 30
      go to ( 10, 20 ), sol

10    continue

      call j0j1 ( u, jy )
      go to 40

20    continue

      call y0y1 ( u, jy )
      go to 40

30    continue

      call luke ( u, sol, jy )

40    continue

      if ( n .lt. 2 ) go to 100
      go to ( 50, 60 ), sol

50    continue

      call jns ( jy, u, nn )
      go to 100
c
c  recurrence formula.
c
60    continue

      do 70 k = 2, nn
        jy(k+1) = 2.0d+00 
     &    * dble ( float ( k - 1 ) ) * jy(k) / u - jy(k-1)
70    continue

      go to 100

80    continue

      nn = nn + 1

      do 90 k = 1, nn
        jy(k) = -1.0d+37
90    continue

100   continue

      return
      end
      subroutine j0j1 ( x, j )

c*********************************************************************72

      implicit none

      real dum(1014)
      double precision j(2)
      double precision t(5)
      double precision x

      common dum, t

      t(1) = x / 2.0d+00
      j(1) = 1.0d+00
      j(2) = t(1)
      t(2) = -t(1)**2
      t(3) = 1.0d+00
      t(4) = 1.0d+00

10    t(4) = t(4) * t(2) / t(3)**2
      j(1) = j(1) + t(4)
      t(5) = t(4) * t(1) / ( t(3) + 1.0d+00 )
      j(2) = j(2) + t(5)

      if ( dmax1 ( dabs ( t(4) ), dabs ( t(5) ) ) .lt. 1.0d-15 ) return

      t(3) = t(3) + 1.0d+00
      go to 10
      end
      subroutine y0y1 ( x, y )

      implicit none

      real dum(1014)
      double precision t(10)
      double precision x
      double precision y(2)

      common dum, t

      t(1) = x / 2.0d+00
      t(2) = -t(1)**2
      y(1) = 1.0d+00
      y(2) = t(1)
      t(7) = 0.0d+00
      t(10) = -t(1)
      t(3) = 0.0d+00
      t(4) = 0.0d+00
      t(5) = 1.0d+00

10    continue

      t(3) = t(3) + 1.0d+00
      t(4) = t(4) + 1.0d+00 / t(3)
      t(5) = t(5) * t(2) / t(3)**2
      y(1) = y(1) + t(5)
      t(6) = -t(5) * t(4)
      t(7) = t(7) + t(6)
      t(8) = t(5) * t(1) / ( t(3) + 1.0d+00 )
      y(2) = y(2) + t(8)
      t(9) = -t(8) * ( 2.0d+00 * t(4) + 1.0d+00 / ( t(3) + 1.0d+00 ) )
      t(10) = t(10) + t(9)

      if ( dmax1 ( dabs ( t(6) ), dabs ( t(9) ) ) .ge. 1.0d-15 ) 
     &  go to 10

      t(2) = 0.57721566490153286d+00 + dlog ( t(1) )
      y(1) = 0.63661977236758134d+00 * ( y(1) * t(2) + t(7) )
      y(2) = 0.63661977236758134d+00 * ( y(2) * t(2) - 1.0d+00 / x ) 
     &  + t(10) / 3.1415926535897932d+00

      return
      end
      subroutine luke ( u, kind, jy )

c*********************************************************************72

      implicit none

      double precision a(19)
      double precision b(19)
      double precision c(19)
      double precision cs
      double precision d(19)
      real dum(1014)
      double precision g(3)
      double precision jy(2)
      integer k
      integer kind
      double precision r(2)
      double precision s(2)
      double precision sn
      double precision t
      double precision u
      double precision x

      save a
      save b
      save c
      save d

      common dum, r, s, g, x, t, sn, cs

      data a /
     &   0.99959506476867287416d+00,
     &  -0.53807956139606913d-03,
     &  -0.13179677123361570d-03,
     &   0.151422497048644d-05,
     &   0.15846861792063d-06,
     &  -0.856069553946d-08,
     &  -0.29572343355d-09,
     &   0.6573556254d-10,
     &  -0.223749703d-11,
     &  -0.44821140d-12,
     &   0.6954827d-13,
     &  -0.151340d-14,
     &  -0.92422d-15,
     &   0.15558d-15,
     &  -0.476d-17,
     &  -0.274d-17,
     &   0.61d-18,
     &  -0.4d-19,
     &  -0.1d-19 /
      data b /
     &  -0.776935569420532136d-02,
     &  -0.774803230965447670d-02,
     &   0.2536541165430796d-04,
     &   0.394273598399711d-05,
     &  -0.10723498299129d-06,
     &  -0.721389799328d-08,
     &   0.73764602893d-09,
     &   0.150687811d-11,
     &  -0.574589537d-11,
     &   0.45996574d-12,
     &   0.2270323d-13,
     &  -0.887890d-14,
     &   0.74497d-15,
     &   0.5847d-16,
     &  -0.2410d-16,
     &   0.265d-17,
     &   0.13d-18,
     &  -0.10d-18,
     &   0.2d-19 /

      data c /
     &   1.00067753586591346234d+00,
     &   0.90100725195908183d-03,
     &   0.22172434918599454d-03,
     &  -0.196575946319104d-05,
     &  -0.20889531143270d-06,
     &   0.1028144350894d-07,
     &   0.37597054789d-10,
     &  -0.7638891358d-10,
     &   0.238734670d-11,
     &   0.51825489d-12,
     &  -0.7693969d-13,
     &   0.144008d-14,
     &   0.103294d-14,
     &  -0.16821d-15,
     &   0.459d-17,
     &   0.302d-17,
     &  -0.65d-18,
     &   0.4d-19,
     &   0.1d-19 /
      data d /
     &   0.2337682998628580328d-01,
     &   0.2334680122354557533d-01,
     &  -0.3576010590901382d-04,
     &  -0.560863149492627d-05,
     &   0.13273894084340d-06,
     &   0.916975845066d-08,
     &  -0.86838880371d-09,
     &  -0.378073005d-11,
     &   0.663145586d-11,
     &  -0.50584390d-12,
     &  -0.2720782d-13,
     &   0.985381d-14,
     &  -0.79398d-15,
     &  -0.6757d-16,
     &   0.2625d-16,
     &  -0.280d-17,
     &  -0.15d-18,
     &   0.10d-18,
     &  -0.2d-19 /

      x = 8.0d+00 / u
      g(1) = 1.0d+00
      g(2) = 2.0d+00 * x - 1.0d+00
      r(1) = a(1) + a(2) * g(2)
      s(1) = b(1) + b(2) * g(2)
      r(2) = c(1) + c(2) * g(2)
      s(2) = d(1) + d(2) * g(2)
      do 10 k = 3, 19
        g(3) = ( 4.0d+00 * x - 2.0d+00 ) * g(2) - g(1)
        r(1) = r(1) + a(k) * g(3)
        s(1) = s(1) + b(k) * g(3)
        r(2) = r(2) + c(k) * g(3)
        s(2) = s(2) + d(k) * g(3)
        g(1) = g(2)
        g(2) = g(3)
10    continue

      t = 0.7978845608028654d+00 / dsqrt ( u )
      sn = dsin ( u - 0.7853981633974483d+00 )
      cs = dcos ( u - 0.7853981633974483d+00 )

      go to ( 20, 30 ), kind

20    jy(1) = t * ( r(1) * cs - s(1) * sn )
      jy(2) = t * ( r(2) * sn + s(2) * cs )
      go to 40

30    jy(1) = t * ( s(1) * cs + r(1) * sn )
      jy(2) = t * ( s(2) * sn - r(2) * cs )

40    return
      end
      subroutine jns ( jj, u, m )

c*********************************************************************72

      implicit none

      double precision a
      double precision b
      double precision d(2)
      double precision dm
      real dum(1014)
      double precision g(249)
      double precision jj(250)
      integer k
      integer ka
      integer kk
      integer m
      integer m2
      double precision p(3)
      double precision q(3)
      double precision u

      equivalence ( a, g(1) )
      equivalence ( d(1), g(2) )
      equivalence ( p(1), g(4) )
      equivalence ( q(1), g(7) )
      equivalence ( dm, g(10) )
      equivalence ( b, g(11) )

      common dum, g, m2, k, kk, ka

      m2 = m

      dm = 2.0d+00 * m
      p(1) = 0.0d+00
      q(1) = 1.0d+00
      p(2) = 1.0d+00
      q(2) = dm / u
      d(1) = 1.0d+00 / q(2)
      a = 2.0d+00

10    b = ( dm + a ) / u
      p(3) = b * p(2) - p(1)
      q(3) = b * q(2) - q(1)
      d(2) = p(3) / q(3)

      if ( dabs ( d(1) - d(2) ) .lt. 1.0d-15 ) go to 20

      p(1) = p(2)
      p(2) = p(3)
      q(1) = q(2)
      q(2) = q(3)
      d(1) = d(2)
      a = a + 2.0d+00
      go to 10

20    g(m) = d(2)
      ka = m - 2

      do 30 k = 1, ka
        kk = m - k
        a = 2 * kk
        g(kk) = u / ( a - u * g(kk+1) )
        if ( g(kk) .eq. 0.0d+00 ) g(kk) = 1.0d-35
30    continue

      do 40 k = 2, m
        jj(k+1) = g(k) * jj(k)
40    continue

      return
      end
      double precision function ds ( kk )

c*********************************************************************72
c
c  evaluates one term of the radial solution associated with b(q).
c
      implicit none

      double precision ab(200)
      real dum1(1004)
      real dum2(17)
      real dum3(583)
      double precision fj
      double precision fy
      integer k
      integer kk
      integer n
      integer n1
      integer n2
      integer p
      integer s

      common dum1, n, p, s, dum2,
     &  k, n1, n2, dum3, ab

      k = kk
      n1 = k - s
      n2 = k + s + p
      ds = ab(k+1) * ( fj ( n1 ) * fy ( n2 ) - fj ( n2 ) * fy ( n1 ) )
      if ( mod ( k + n, 2 ) .ne. 0 ) ds = -ds

      return
      end
      double precision function dc ( kk )

c*********************************************************************72
c
c  evaluates one term of the radial solution associated with a(q).
c
      implicit none

      double precision ab(200)
      real dum1(1004)
      real dum2(17)
      real dum3(583)
      double precision fj
      double precision fy
      integer k
      integer kk
      integer n
      integer n1
      integer n2
      integer p
      integer s     

      common dum1, n, p, s, dum2,
     &  k, n1, n2, dum3, ab

      k = kk
      n1 = k - s
      n2 = k + s + p     
      dc = ab(k+1) * ( fj ( n1 ) * fy ( n2 ) - fj ( n2 ) * fy ( n1 ) )
      if ( mod ( k + n, 2 ) .ne. 0 ) dc = -dc
      if ( s + p .eq. 0 ) dc = 0.5d+00 * dc

      return
      end
      double precision function dds ( kk )

c*********************************************************************72
c
c  evaluates one term of the derivative of the radial solution
c  associated with b(q).
c
      implicit none

      double precision ab(200)
      real dum1(1000)
      real dum2(17)
      real dum3(583)
      double precision dj
      double precision dy
      double precision fj
      double precision fy
      integer k
      integer kk
      integer n
      integer n1
      integer n2
      integer p
      integer s
      double precision u1
      double precision u2

      common dum1, u1, u2, n, p, s, dum2,
     &  k, n1, n2, dum3, ab

      k = kk
      n1 = k - s
      n2 = k + s + p
      dds = ab(k+1) * ( 
     &    u2 * ( fj ( n1 ) * dy ( n2 ) 
     &         - fj ( n2 ) * dy ( n1 ) )
     &  - u1 * ( fy ( n2 ) * dj ( n1 ) 
     &         - fy ( n1 ) * dj ( n2 ) ) )
 
      if ( mod ( k + n, 2 ) .ne. 0 ) dds = -dds

      return
      end
      double precision function ddc ( kk )

c*********************************************************************72
c
c  evaluates one term of the derivative of the radial solution 
c  associated with a(q).
c
      implicit none

      double precision ab(200)
      real dum1(1000)
      real dum2(17)
      real dum3(583)
      double precision dj
      double precision dy
      double precision fj
      double precision fy
      integer k
      integer kk
      integer n
      integer n1
      integer n2
      integer p
      integer s
      double precision u1
      double precision u2

      common dum1, u1, u2, n, p, s, dum2,
     &  k, n1, n2, dum3, ab

      k = kk
      n1 = k - s
      n2 = k + s + p     
      ddc = ab(k+1) * ( 
     &    u2 * ( fj ( n1 ) * dy ( n2 ) 
     &         + fj ( n2 ) * dy ( n1 ) )
     &  - u1 * ( fy ( n2 ) * dj ( n1 ) 
     &         + fy ( n1 ) * dj ( n2 ) ) )

      if ( mod ( k + n, 2 ) .ne. 0 ) ddc = -ddc
      if ( s + p .eq. 0 ) ddc = 0.5d+00 * ddc

      return
      end
      double precision function ps ( k )

c*********************************************************************72
c
c  evaluates one term of the odd periodic solution.
c
      implicit none

      double precision ab(200)
      real dum1(1005)
      real dum2(2)
      real dum3(600)
      integer k
      integer p
      double precision x

      common dum1, p, dum2, x, dum3, ab

      ps = ab(k+1) * dsin ( dble ( float ( 2 * k + p ) ) * x )

      return
      end
      double precision function pc ( k )

c*********************************************************************72
c
c  evaluates one term of the even periodic solution.
c
      implicit none

      double precision ab(200)
      real dum1(1005)
      real dum2(2)
      real dum3(600)
      integer k
      integer p
      double precision x

      common dum1, p, dum2, x, dum3, ab

      pc = ab(k+1) * dcos ( dble ( float ( 2 * k + p ) ) * x )

      return
      end
      double precision function dps ( k )

c*********************************************************************72
c
c  evaluates one term of the derivative of the odd periodic solution.
c
      implicit none

      double precision ab(200)
      real dum1(1005)
      real dum2(2)
      real dum3(14)
      real dum4(584)
      integer k
      integer p
      double precision t
      double precision x

      common dum1, p, dum2, x, dum3, t, dum4, ab

      t = 2 * k + p
      dps = ab(k+1) * t * dcos ( t * x )

      return
      end
      double precision function dpc ( k )

c*********************************************************************72
c
c  evaluates one term of the derivative of the even periodic solution.
c
      implicit none

      double precision ab(200)
      real dum1(1005)
      real dum2(2)
      real dum3(14)
      real dum4(584)
      integer k
      integer p
      double precision t
      double precision x

      common dum1, p, dum2, x, dum3, t, dum4, ab

      t = 2 * k + p
      dpc = -ab(k+1) * t * dsin ( t * x )

      return
      end
      double precision function fj ( n )

c*********************************************************************72
c
cc fj produces bessel functions of the first kind.
c
      implicit none

      real dum(527)
      double precision j(250)
      integer k
      integer n

      common j, dum, k

      k = iabs ( n )

      if ( k .ge. 250 ) go to 20
      fj = j(k+1)
      if ( mod ( n, 2 ) .lt. 0 ) fj = -fj
10    return

20    fj = 0.0d+00
      write ( 6, 99 ) n
99    format ( 2h0j, i3, 7h needed )
      go to 10
      end
      double precision function fy ( n )

c*********************************************************************72
c
cc fy produces bessel functions of the second kind.
c
      implicit none

      real dum1(500)
      real dum2(27)
      integer k
      integer n
      double precision y(250)

      common dum1, y, dum2, k

      k = iabs ( n )

      if ( k .ge. 250 ) go to 20
      fy = y(k+1)
      if ( mod ( n, 2 ) .lt. 0 ) fy = -fy
10    return

20    fy = 0.0d+00
      write ( 6, 99 ) n
99    format ( 2h0y, i3, 7h needed )
      go to 10
      end
      double precision function dj ( n )

c*********************************************************************72
c
cc dj computes derivatives of bessel functions of the first kind.
c
      implicit none

      real dum1(1000)
      real dum2(26)
      double precision fj
      double precision fn
      integer n
      double precision u1

      common dum1, u1, dum2, fn

      fn = n

      if ( n - 249 ) 10, 20, 40

10    dj = fn * fj ( n ) / u1 - fj ( n+1 )
      go to 30

20    dj = fj ( n-1 ) - fn * fj ( n ) / u1

30    return

40    dj = 0.0d+00
      write ( 6, 99 ) n
99    format ( 3h0j@, i3, 7h needed )
      go to 30
      end
      double precision function dy ( n )

c*********************************************************************72
c
c  derivatives of bessel functions of the second kind.
c
      implicit none

      real dum1(1002)
      real dum2(24)
      double precision fn
      double precision fy
      integer n
      double precision u2

      common dum1, u2, dum2, fn

      if ( n .ge. 250 ) go to 20

      fn = n
      dy = fy ( n-1 ) * fn * fy ( n ) / u2
10    return

20    dy = 0.0d+00
      write ( 6, 99 ) n
99    format ( 3h0y@, i3, 7h needed )
      go to 10
      end