segv Subroutine

subroutine segv ( m, n, c, kd, cv, eg )

  !*****************************************************************************80
  !
  !! SEGV computes the characteristic values of spheroidal wave functions.
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    28 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, integer ( kind = 4 ) M, the mode parameter.
  !
  !    Input, integer ( kind = 4 ) N, the mode parameter.
  !
  !    Input, real ( kind = 8 ) C, the spheroidal parameter.
  !
  !    Input, integer ( kind = 4 ) KD, the function code.
  !    1, the prolate function.
  !    -1, the oblate function.
  !
  !    Output, real ( kind = 8 ) CV, the characteristic value.
  !
  !    Output, real ( kind = 8 ) EG(*), the characteristic value for 
  !    mode parameters m and n.  ( L = n - m + 1 )
  !
  implicit none

  real ( kind = 8 ) a(300)
  real ( kind = 8 ) b(100)
  real ( kind = 8 ) c
  real ( kind = 8 ) cs
  real ( kind = 8 ) cv
  real ( kind = 8 ) cv0(100)
  real ( kind = 8 ) d(300)
  real ( kind = 8 ) d2k
  real ( kind = 8 ) dk0
  real ( kind = 8 ) dk1
  real ( kind = 8 ) dk2
  real ( kind = 8 ) e(300)
  real ( kind = 8 ) eg(200)
  real ( kind = 8 ) f(300)
  real ( kind = 8 ) g(300)
  real ( kind = 8 ) h(100)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) icm
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) k1
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nm
  integer ( kind = 4 ) nm1
  real ( kind = 8 ) s
  real ( kind = 8 ) t
  real ( kind = 8 ) t1
  real ( kind = 8 ) x1
  real ( kind = 8 ) xa
  real ( kind = 8 ) xb

  if ( c < 1.0D-10 ) then
     do i = 1, n
        eg(i) = ( i + m ) * ( i + m - 1.0D+00 )
     end do
     cv = eg(n-m+1)
     return
  end if

  icm = ( n - m + 2 ) / 2
  nm = 10 + int ( 0.5D+00 * ( n - m ) + c )
  cs = c * c * kd

  do l = 0, 1

     do i = 1, nm
        if ( l == 0 ) then
           k = 2 * ( i - 1 )
        else
           k = 2 * i - 1
        end if
        dk0 = m + k
        dk1 = m + k + 1
        dk2 = 2 * ( m + k )
        d2k = 2 * m + k
        a(i) = ( d2k + 2.0D+00 ) * ( d2k + 1.0D+00 ) &
             / ( ( dk2 + 3.0D+00 ) * ( dk2 + 5.0D+00 ) ) * cs
        d(i) = dk0 * dk1 + ( 2.0D+00 * dk0 * dk1 &
             - 2.0 * m * m - 1.0D+00 ) &
             / ( ( dk2 - 1.0D+00 ) * ( dk2 + 3.0D+00 ) ) * cs
        g(i) = k * ( k - 1.0D+00 ) / ( ( dk2 - 3.0D+00 ) &
             * ( dk2 - 1.0D+00 ) ) * cs
     end do

     do k = 2, nm
        e(k) = sqrt ( a(k-1) * g(k) )
        f(k) = e(k) * e(k)
     end do

     f(1) = 0.0D+00
     e(1) = 0.0D+00
     xa = d(nm) + abs ( e(nm) )
     xb = d(nm) - abs ( e(nm) )
     nm1 = nm - 1
     do i = 1, nm1
        t = abs ( e(i) ) + abs ( e(i+1) )
        t1 = d(i) + t
        if ( xa < t1 ) then
           xa = t1
        end if
        t1 = d(i) - t
        if ( t1 < xb ) then
           xb = t1
        end if
     end do

     do i = 1, icm
        b(i) = xa
        h(i) = xb
     end do

     do k = 1, icm

        do k1 = k, icm
           if ( b(k1) < b(k) ) then
              b(k) = b(k1)
              exit
           end if
        end do

        if ( k /= 1 .and. h(k) < h(k-1) ) then
           h(k) = h(k-1)
        end if

        do

           x1 = ( b(k) + h(k) ) /2.0D+00
           cv0(k) = x1

           if ( abs ( ( b(k) - h(k) ) / x1 ) < 1.0D-14 ) then
              exit
           end if

           j = 0
           s = 1.0D+00

           do i = 1, nm

              if ( s == 0.0D+00 ) then
                 s = s + 1.0D-30
              end if
              t = f(i) / s
              s = d(i) - t - x1
              if ( s < 0.0D+00 ) then
                 j = j + 1
              end if
           end do

           if ( j < k ) then

              h(k) = x1

           else

              b(k) = x1
              if ( icm <= j ) then
                 b(icm) = x1
              else
                 if ( h(j+1) < x1 ) then
                    h(j+1) = x1
                 end if
                 if ( x1 < b(j) ) then
                    b(j) = x1
                 end if
              end if

           end if

        end do

        cv0(k) = x1

        if ( l == 0 ) then
           eg(2*k-1) = cv0(k)
        else
           eg(2*k) = cv0(k)
        end if

     end do

  end do

  cv = eg(n-m+1)

  return
end subroutine segv