cvql Subroutine

subroutine cvql ( kd, m, q, a0 )

  !*****************************************************************************80
  !
  !! CVQL computes the characteristic value of Mathieu functions for q <= 3*m.
  !
  !  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:
  !
  !    10 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 value.
  !
  !    Output, real ( kind = 8 ) A0, the initial characteristic value.
  !
  implicit none

  real ( kind = 8 ) a0
  real ( kind = 8 ) c1
  real ( kind = 8 ) cv1
  real ( kind = 8 ) cv2
  real ( kind = 8 ) d1
  real ( kind = 8 ) d2
  real ( kind = 8 ) d3
  real ( kind = 8 ) d4
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) m
  real ( kind = 8 ) p1
  real ( kind = 8 ) p2
  real ( kind = 8 ) q
  real ( kind = 8 ) w
  real ( kind = 8 ) w2
  real ( kind = 8 ) w3
  real ( kind = 8 ) w4
  real ( kind = 8 ) w6

  if ( kd == 1 .or. kd == 2 ) then
     w = 2.0D+00 * m + 1.0D+00
  else
     w = 2.0D+00 * m - 1.0D+00
  end if

  w2 = w * w
  w3 = w * w2
  w4 = w2 * w2
  w6 = w2 * w4
  d1 = 5.0D+00 + 34.0D+00 / w2 + 9.0D+00 / w4
  d2 = ( 33.0D+00 + 410.0D+00 / w2 + 405.0D+00 / w4 ) / w
  d3 = ( 63.0D+00 + 1260.0D+00 / w2 + 2943.0D+00 / w4 + 486.0D+00 / w6 ) / w2
  d4 = ( 527.0D+00 + 15617.0D+00 / w2 + 69001.0D+00 / w4 &
       + 41607.0D+00 / w6 ) / w3
  c1 = 128.0D+00
  p2 = q / w4
  p1 = sqrt ( p2 )
  cv1 = - 2.0D+00 * q + 2.0D+00 * w * sqrt ( q ) &
       - ( w2 + 1.0D+00 ) / 8.0D+00
  cv2 = ( w + 3.0D+00 / w ) + d1 / ( 32.0D+00 * p1 ) + d2 &
       / ( 8.0D+00 * c1 * p2 )
  cv2 = cv2 + d3 / ( 64.0D+00 * c1 * p1 * p2 ) + d4 &
       / ( 16.0D+00 * c1 * c1 * p2 * p2 )
  a0 = cv1 - cv2 / ( c1 * p1 )

  return
end subroutine cvql