cva2 Subroutine

subroutine cva2 ( kd, m, q, a )

  !*****************************************************************************80
  !
  !! CVA2 computes a specific characteristic value 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:
  !
  !    22 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 order of the Mathieu functions.
  !
  !    Input, real ( kind = 8 ) Q, the parameter of the Mathieu functions.
  !
  !    Output, real ( kind = 8 ) A, the characteristic value.
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) a1
  real ( kind = 8 ) a2
  real ( kind = 8 ) delta
  integer ( kind = 4 ) i
  integer ( kind = 4 ) iflag
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) m
  integer ( kind = 4 ) ndiv
  integer ( kind = 4 ) nn
  real ( kind = 8 ) q
  real ( kind = 8 ) q1
  real ( kind = 8 ) q2
  real ( kind = 8 ) qq

  if ( m <= 12 .or. q <= 3.0D+00 * m .or. m * m < q ) then

     call cv0 ( kd, m, q, a )

     if ( q /= 0.0D+00 ) then
        call refine ( kd, m, q, a, 1 )
     end if

  else

     ndiv = 10
     delta = ( m - 3.0D+00 ) * m / real ( ndiv, kind = 8 )

     if ( ( q - 3.0D+00 * m ) <= ( m * m - q ) ) then

        do

           nn = int ( ( q - 3.0D+00 * m ) / delta ) + 1
           delta = ( q - 3.0D+00 * m ) / nn
           q1 = 2.0D+00 * m
           call cvqm ( m, q1, a1 )
           q2 = 3.0D+00 * m
           call cvqm ( m, q2, a2 )
           qq = 3.0D+00 * m

           do i = 1, nn

              qq = qq + delta
              a = ( a1 * q2 - a2 * q1 + ( a2 - a1 ) * qq ) / ( q2 - q1 )

              if ( i == nn ) then
                 iflag = -1
              else
                 iflag = 1
              end if

              call refine ( kd, m, qq, a, iflag )
              q1 = q2
              q2 = qq
              a1 = a2
              a2 = a

           end do

           if ( iflag /= -10 ) then
              exit
           end if

           ndiv = ndiv * 2
           delta = ( m - 3.0D+00 ) * m / real ( ndiv, kind = 8 )

        end do

     else

        do

           nn = int ( ( m * m - q ) / delta ) + 1
           delta = ( m * m - q ) / nn
           q1 = m * ( m - 1.0D+00 )
           call cvql ( kd, m, q1, a1 )
           q2 = m * m
           call cvql ( kd, m, q2, a2 )
           qq = m * m

           do i = 1, nn

              qq = qq - delta
              a = ( a1 * q2 - a2 * q1 + ( a2 - a1 ) * qq ) / ( q2 - q1 )

              if ( i == nn ) then
                 iflag = -1
              else
                 iflag = 1
              end if

              call refine ( kd, m, qq, a, iflag )
              q1 = q2
              q2 = qq
              a1 = a2
              a2 = a

           end do

           if ( iflag /= -10 ) then
              exit
           end if

           ndiv = ndiv * 2
           delta = ( m - 3.0D+00 ) * m / real ( ndiv, kind = 8 )

        end do

     end if

  end if

  return
end subroutine cva2