mtu12 Subroutine

subroutine mtu12 ( kf, kc, m, q, x, f1r, d1r, f2r, d2r )

  !*****************************************************************************80
  !
  !! MTU12 computes modified Mathieu functions of the first and second kind.
  !
  !  Discussion:
  !
  !    This procedure computes modified Mathieu functions of the first and
  !    second kinds, Mcm(1)(2)(x,q) and Msm(1)(2)(x,q),
  !    and their 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:
  !
  !    31 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 Mcm(x,q);
  !    2 for computing Msm(x,q).
  !
  !    Input, integer ( kind = 4 ) KC, the function code.
  !    1, for computing the first kind
  !    2, for computing the second kind or Msm(2)(x,q) and Msm(2)'(x,q)
  !    3, for computing both the first and second kinds.
  !
  !    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.
  !
  !    Output, real ( kind = 8 ) F1R, D1R, F2R, D2R, the values of 
  !    Mcm(1)(x,q) or Msm(1)(x,q), Derivative of Mcm(1)(x,q) or Msm(1)(x,q),
  !    Mcm(2)(x,q) or Msm(2)(x,q), Derivative of Mcm(2)(x,q) or Msm(2)(x,q).
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) bj1(0:251)
  real ( kind = 8 ) bj2(0:251)
  real ( kind = 8 ) by1(0:251)
  real ( kind = 8 ) by2(0:251)
  real ( kind = 8 ) c1
  real ( kind = 8 ) c2
  real ( kind = 8 ) d1r
  real ( kind = 8 ) d2r
  real ( kind = 8 ) dj1(0:251)
  real ( kind = 8 ) dj2(0:251)
  real ( kind = 8 ) dy1(0:251)
  real ( kind = 8 ) dy2(0:251)
  real ( kind = 8 ) eps
  real ( kind = 8 ) f1r
  real ( kind = 8 ) f2r
  real ( kind = 8 ) fg(251)
  integer ( kind = 4 ) ic
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kc
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) kf
  integer ( kind = 4 ) km
  integer ( kind = 4 ) m
  integer ( kind = 4 ) nm
  real ( kind = 8 ) q
  real ( kind = 8 ) qm
  real ( kind = 8 ) u1
  real ( kind = 8 ) u2
  real ( kind = 8 ) w1
  real ( kind = 8 ) w2
  real ( kind = 8 ) x

  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 )

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

  c1 = exp ( - x )
  c2 = exp ( x )
  u1 = sqrt ( q ) * c1
  u2 = sqrt ( q ) * c2

  call jynb ( km, u1, nm, bj1, dj1, by1, dy1 )
  call jynb ( km, u2, nm, bj2, dj2, by2, dy2 )

  if ( kc == 1 ) then

     f1r = 0.0D+00

     do k = 1, km

        if ( kd == 1 ) then
           f1r = f1r + ( - 1.0D+00 ) ** ( ic + k ) * fg(k) * bj1(k-1) * bj2(k-1)
        else if ( kd == 2 .or. kd == 3 ) then
           f1r = f1r + ( -1.0D+00 ) ** ( ic + k ) * fg(k) * ( bj1(k-1) * bj2(k) &
                + ( - 1.0D+00 ) ** kd * bj1(k) * bj2(k-1) )
        else
           f1r = f1r + ( -1.0D+00 ) ** ( ic + k ) * fg(k) &
                * ( bj1(k-1) * bj2(k+1) - bj1(k+1) * bj2(k-1) )
        end if

        if ( 5 <= k .and. abs ( f1r - w1 ) < abs ( f1r ) * eps ) then
           exit
        end if

        w1 = f1r

     end do

     f1r = f1r / fg(1)
     d1r = 0.0D+00
     do k = 1, km
        if ( kd == 1 ) then
           d1r = d1r + ( - 1.0D+00 ) ** ( ic + k ) * fg(k) &
                * ( c2 * bj1(k-1) * dj2(k-1) - c1 * dj1(k-1) * bj2(k-1) )
        else if ( kd == 2 .or. kd == 3 ) then
           d1r = d1r + ( -1.0D+00 ) ** ( ic + k ) * fg(k) &
                * ( c2 * ( bj1(k-1) * dj2(k) &
                + ( -1.0D+00 ) ** kd * bj1(k) * dj2(k-1) ) &
                - c1 * ( dj1(k-1) * bj2(k) &
                + ( -1.0D+00 ) ** kd * dj1(k) * bj2(k-1) ) )
        else
           d1r = d1r + ( -1.0D+00 ) ** ( ic + k ) * fg(k) &
                * ( c2 * ( bj1(k-1) * dj2(k+1) - bj1(k+1) * dj2(k-1) ) &
                - c1 * ( dj1(k-1) * bj2(k+1) - dj1(k+1) * bj2(k-1) ) )
        end if
        if ( 5 <= k .and. abs ( d1r - w2 ) < abs ( d1r ) * eps ) then
           exit
        end if
        w2 = d1r
     end do

     d1r = d1r * sqrt ( q ) / fg(1)

  else

     f2r = 0.0D+00

     do k = 1, km
        if ( kd == 1 ) then
           f2r = f2r + ( -1.0D+00 ) ** ( ic + k ) * fg(k) &
                * bj1(k-1) * by2(k-1)
        else if ( kd == 2 .or. kd == 3 ) then
           f2r = f2r + ( -1.0D+00 ) ** ( ic + k ) * fg(k) * ( bj1(k-1) * by2(k) &
                + ( -1.0D+00 ) ** kd * bj1(k) * by2(k-1) )
        else
           f2r = f2r + ( -1.0D+00 ) ** ( ic + k ) * fg(k) &
                * ( bj1(k-1) * by2(k+1) - bj1(k+1) * by2(k-1) )
        end if
        if ( 5 <= k .and. abs ( f2r - w1 ) < abs ( f2r ) * eps ) then
           exit
        end if
        w1 = f2r
     end do

     f2r = f2r / fg(1)
     d2r = 0.0D+00

     do k = 1, km
        if ( kd == 1 ) then
           d2r = d2r + ( -1.0D+00 ) ** ( ic + k ) * fg(k) &
                * ( c2 * bj1(k-1) * dy2(k-1) - c1 * dj1(k-1) * by2(k-1) )
        else if ( kd == 2 .or. kd == 3 ) then
           d2r = d2r + ( -1.0D+00 ) ** ( ic + k ) * fg(k) &
                * ( c2 * ( bj1(k-1) * dy2(k) &
                + ( -1.0D+00 ) ** kd * bj1(k) * dy2(k-1) ) &
                - c1 * ( dj1(k-1) * by2(k) + ( -1.0D+00 ) ** kd &
                * dj1(k) * by2(k-1) ) )
        else
           d2r = d2r + ( -1.0D+00 ) ** ( ic + k ) * fg(k) &
                * ( c2 * ( bj1(k-1) * dy2(k+1) - bj1(k+1) * dy2(k-1) ) &
                - c1 * ( dj1(k-1) * by2(k+1) - dj1(k+1) * by2(k-1) ) )
        end if

        if ( 5 <= k .and. abs ( d2r - w2 ) < abs ( d2r ) * eps ) then
           exit
        end if

        w2 = d2r

     end do

     d2r = d2r * sqrt ( q ) / fg(1)

  end if

  return
end subroutine mtu12