cvf Subroutine

subroutine cvf ( kd, m, q, a, mj, f )

  !*****************************************************************************80
  !
  !! CVF computes F for the characteristic equation 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:
  !
  !    16 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.
  !
  !    Input, real ( kind = 8 ) A, the characteristic value.
  !
  !    Input, integer ( kind = 4 ) MJ, ?
  !
  !    Output, real ( kind = 8 ) F, the value of the function for the
  !    characteristic equation.
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  real ( kind = 8 ) f
  integer ( kind = 4 ) ic
  integer ( kind = 4 ) j
  integer ( kind = 4 ) j0
  integer ( kind = 4 ) jf
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) l
  integer ( kind = 4 ) l0
  integer ( kind = 4 ) m
  integer ( kind = 4 ) mj
  real ( kind = 8 ) q
  real ( kind = 8 ) t0
  real ( kind = 8 ) t1
  real ( kind = 8 ) t2

  b = a
  ic = int ( m / 2 )
  l = 0
  l0 = 0
  j0 = 2
  jf = ic

  if ( kd == 1 ) then
     l0 = 2
     j0 = 3
  else if ( kd == 2 .or. kd == 3 ) then
     l = 1
  else if ( kd == 4 ) then
     jf = ic - 1
  end if

  t1 = 0.0D+00
  do j = mj, ic + 1, -1
     t1 = - q * q / ( ( 2.0D+00 * j + l ) ** 2 - b + t1 )
  end do

  if ( m <= 2 ) then

     t2 = 0.0D+00

     if ( kd == 1 ) then
        if ( m == 0 ) then
           t1 = t1 + t1
        else if ( m == 2 ) then
           t1 = - 2.0D+00 * q * q / ( 4.0D+00 - b + t1 ) - 4.0D+00
        end if
     else if ( kd == 2 ) then
        if ( m == 1 ) then
           t1 = t1 + q
        end if
     else if ( kd == 3 ) then
        if ( m == 1 ) then
           t1 = t1 - q
        end if
     end if

  else

     if ( kd == 1 ) then
        t0 = 4.0D+00 - b + 2.0D+00 * q * q / b
     else if ( kd == 2 ) then
        t0 = 1.0D+00 - b + q
     else if ( kd == 3 ) then
        t0 = 1.0D+00 - b - q
     else if ( kd == 4 ) then
        t0 = 4.0D+00 - b
     end if

     t2 = - q * q / t0
     do j = j0, jf
        t2 = - q * q / ( ( 2.0D+00 * j - l - l0 ) ** 2 - b + t2 )
     end do

  end if

  f = ( 2.0D+00 * ic + l ) ** 2 + t1 + t2 - b

  return
end subroutine cvf