cva1 Subroutine

subroutine cva1 ( kd, m, q, cv )

  !*****************************************************************************80
  !
  !! CVA1 computes a sequence of characteristic values of Mathieu 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:
  !
  !    25 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 ) KD, the case code.
  !    1, for cem(x,q)  ( m = 0,2,4,... )
  !    2, for cem(x,q)  ( m = 1,3,5,... )
  !    3, for sem(x,q)  ( m = 1,3,5,... )
  !    4, for sem(x,q)  ( m = 2,4,6,... )
  !
  !    Input, integer ( kind = 4 ) M, the maximum order of the Mathieu functions.
  !
  !    Input, real ( kind = 8 ) Q, the parameter of the Mathieu functions.
  !
  !    Output, real ( kind = 8 ) CV(*), characteristic values.
  !    For KD = 1, CV(1), CV(2), CV(3),..., correspond to
  !    the characteristic values of cem for m = 0,2,4,...
  !    For KD = 2, CV(1), CV(2), CV(3),..., correspond to
  !    the characteristic values of cem for m = 1,3,5,...
  !    For KD = 3, CV(1), CV(2), CV(3),..., correspond to
  !    the characteristic values of sem for m = 1,3,5,...
  !    For KD = 4, CV(1), CV(2), CV(3),..., correspond to
  !    the characteristic values of sem for m = 0,2,4,...
  !       
  implicit none

  real ( kind = 8 ) cv(200)
  real ( kind = 8 ) d(500)
  real ( kind = 8 ) e(500)
  real ( kind = 8 ) eps
  real ( kind = 8 ) f(500)
  real ( kind = 8 ) g(200)
  real ( kind = 8 ) h(200)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ic
  integer ( kind = 4 ) icm
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) k1
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) m
  integer ( kind = 4 ) nm
  integer ( kind = 4 ) nm1
  real ( kind = 8 ) q
  real ( kind = 8 ) s
  real ( kind = 8 ) t
  real ( kind = 8 ) t1
  real ( kind = 8 ) x1
  real ( kind = 8 ) xa
  real ( kind = 8 ) xb

  eps = 1.0D-14

  if ( kd == 4 ) then
     icm = m / 2
  else
     icm = int ( m / 2 ) + 1
  end if

  if ( q == 0.0D+00 ) then

     if ( kd == 1 ) then
        do ic = 1, icm
           cv(ic) = 4.0D+00 * ( ic - 1.0D+00 ) ** 2
        end do
     else if ( kd /= 4 ) then
        do ic = 1, icm
           cv(ic) = ( 2.0D+00 * ic - 1.0D+00 ) ** 2
        end do
     else
        do ic = 1, icm
           cv(ic) = 4.0D+00 * ic * ic
        end do
     end if

  else

     nm = int ( 10D+00 + 1.5D+00 * m + 0.5D+00 * q )
     e(1) = 0.0D+00
     f(1) = 0.0D+00

     if ( kd == 1 ) then

        d(1) = 0.0D+00
        do i = 2, nm
           d(i) = 4.0D+00 * ( i - 1.0D+00 ) ** 2
           e(i) = q
           f(i) = q * q
        end do
        e(2) = sqrt ( 2.0D+00 ) * q
        f(2) = 2.0D+00 * q * q

     else if ( kd /= 4 ) then

        d(1) = 1.0D+00 + ( -1.0D+00 ) ** kd * q
        do i = 2, nm
           d(i) = ( 2.0D+00 * i - 1.0D+00 ) ** 2
           e(i) = q
           f(i) = q * q
        end do

     else

        d(1) = 4.0D+00
        do i = 2, nm
           d(i) = 4.0D+00 * i * i
           e(i) = q
           f(i) = q * q
        end do

     end if

     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
        xa = max ( xa, t1 )
        t1 = d(i) - t
        xb = min ( xb, t1 )
     end do

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

     do k = 1, icm

        do k1 = k, icm
           if ( g(k1) < g(k) ) then
              g(k) = g(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 = ( g(k) + h(k) ) /2.0D+00
           cv(k) = x1

           if ( abs ( ( g(k) - h(k) ) / x1 ) < eps ) 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
              g(k) = x1
              if ( icm <= j ) then
                 g(icm) = x1
              else
                 h(j+1) = max ( h(j+1), x1 )
                 g(j) = min ( g(j), x1 )
              end if
           end if

        end do

        cv(k) = x1

     end do

  end if

  return
end subroutine cva1