mtu0 Subroutine

subroutine mtu0 ( kf, m, q, x, csf, csd )

  !*****************************************************************************80
  !
  !! MTU0 computes Mathieu functions CEM(x,q) and SEM(x,q) and derivatives.
  !
  !  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:
  !
  !    20 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 ) KF, the function code.
  !    1 for computing cem(x,q) and cem'(x,q)
  !    2 for computing sem(x,q) and sem'(x,q).
  !
  !    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 ) X, the argument of the Mathieu functions,
  !    in degrees.
  !
  !    Output, real ( kind = 8 ) CSF, CSD, the values of cem(x,q) and cem'(x,q),
  !    or of sem(x,q) and sem'(x,q).
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) csd
  real ( kind = 8 ) csf
  real ( kind = 8 ) eps
  real ( kind = 8 ) fg(251)
  integer ( kind = 4 ) ic
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) kf
  integer ( kind = 4 ) km
  integer ( kind = 4 ) m
  real ( kind = 8 ) q
  real ( kind = 8 ) qm
  real ( kind = 8 ) rd
  real ( kind = 8 ) x
  real ( kind = 8 ) xr

  eps = 1.0D-14

  if ( kf == 1 ) then

     if ( m == 2 * int ( m / 2 ) ) then
        kd = 1
     else
        kd = 2
     end if

  else

     if ( m /= 2 * int ( m / 2 ) ) then
        kd = 3
     else
        kd = 4
     end if

  end if

  call cva2 ( kd, m, q, a )

  if ( q <= 1.0D+00 ) then
     qm = 7.5D+00 + 56.1D+00 * sqrt ( q ) - 134.7D+00 * q &
          + 90.7D+00 * sqrt ( q ) * q
  else
     qm = 17.0D+00 + 3.1D+00 * sqrt ( q ) - 0.126D+00 * q &
          + 0.0037D+00 * sqrt ( q ) * q
  end if

  km = int ( qm + 0.5D+00 * m )
  call fcoef ( kd, m, q, a, fg )
  ic = int ( m / 2 ) + 1
  rd = 1.74532925199433D-02
  xr = x * rd

  csf = 0.0D+00

  do k = 1, km

     if ( kd == 1 ) then
        csf = csf + fg(k) * cos ( ( 2.0D+00 * k - 2.0D+00 ) * xr )
     else if ( kd == 2 ) then
        csf = csf + fg(k) * cos ( ( 2.0D+00 * k - 1.0D+00 ) * xr )
     else if ( kd == 3 ) then
        csf = csf + fg(k) * sin ( ( 2.0D+00 * k - 1.0D+00 ) * xr )
     else if ( kd == 4 ) then
        csf = csf + fg(k) * sin ( 2.0D+00 * k * xr )
     end if

     if ( ic <= k .and. abs ( fg(k) ) < abs ( csf ) * eps ) then
        exit
     end if

  end do

  csd = 0.0D+00

  do k = 1, km

     if ( kd == 1 ) then
        csd = csd - ( 2 * k - 2 ) * fg(k) * sin ( ( 2 * k - 2 ) * xr )
     else if ( kd == 2 ) then
        csd = csd - ( 2 * k - 1 ) * fg(k) * sin ( ( 2 * k - 1 ) * xr )
     else if ( kd == 3 ) then
        csd = csd + ( 2 * k - 1 ) * fg(k) * cos ( ( 2 * k - 1 ) * xr )
     else if ( kd == 4 ) then
        csd = csd + 2.0D+00 * k * fg(k) * cos ( 2 * k * xr )
     end if

     if ( ic <= k .and. abs ( fg(k) ) < abs ( csd ) * eps ) then
        exit
     end if

  end do

  return
end subroutine mtu0