fcoef Subroutine

subroutine fcoef(kd, m, q, a, fc)

************80

! FCOEF: expansion coefficients for Mathieu and modified 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:

01 August 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 function.

Input, real ( kind = 8 ) Q, the parameter of the Mathieu functions.

Input, real ( kind = 8 ) A, the characteristic value of the Mathieu
functions for given m and q.

Output, real ( kind = 8 ) FC(*), the expansion coefficients of Mathieu
functions ( k =  1,2,...,KM ).  FC(1),FC(2),FC(3),... correspond to
A0,A2,A4,... for KD = 1 case,
A1,A3,A5,... for KD = 2 case,
B1,B3,B5,... for KD = 3 case,
B2,B4,B6,... for KD = 4 case.

Arguments

Type IntentOptional Attributes Name
integer(kind=4) :: kd
integer(kind=4) :: m
real(kind=8) :: q
real(kind=8) :: a
real(kind=8) :: fc(251)

Source Code

subroutine fcoef ( kd, m, q, a, fc )

  !*****************************************************************************80
  !
  !! FCOEF: expansion coefficients for Mathieu and modified 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:
  !
  !    01 August 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 function.
  !
  !    Input, real ( kind = 8 ) Q, the parameter of the Mathieu functions.
  !
  !    Input, real ( kind = 8 ) A, the characteristic value of the Mathieu
  !    functions for given m and q.
  !
  !    Output, real ( kind = 8 ) FC(*), the expansion coefficients of Mathieu
  !    functions ( k =  1,2,...,KM ).  FC(1),FC(2),FC(3),... correspond to
  !    A0,A2,A4,... for KD = 1 case, 
  !    A1,A3,A5,... for KD = 2 case,
  !    B1,B3,B5,... for KD = 3 case,
  !    B2,B4,B6,... for KD = 4 case.
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) f
  real ( kind = 8 ) f1
  real ( kind = 8 ) f2
  real ( kind = 8 ) f3
  real ( kind = 8 ) fc(251)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kb
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) km
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
  real ( kind = 8 ) q
  real ( kind = 8 ) qm
  real ( kind = 8 ) s
  real ( kind = 8 ) s0
  real ( kind = 8 ) sp
  real ( kind = 8 ) ss
  real ( kind = 8 ) u
  real ( kind = 8 ) v

  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 )

  if ( q == 0.0D+00 ) then

     do k = 1, km
        fc(k) = 0.0D+00
     end do

     if ( kd == 1 ) then
        fc((m+2)/2) = 1.0D+00
        if (m == 0 ) then
           fc(1) = 1.0D+00 / sqrt ( 2.0D+00 )
        end if
     else if ( kd == 4 ) then
        fc(m/2) = 1.0D+00
     else
        fc((m+1)/2) = 1.0D+00
     end if

     return

  end if

  kb = 0
  s = 0.0D+00
  f = 1.0D-100
  u = 0.0D+00
  fc(km) = 0.0D+00

  if ( kd == 1 ) then

     l = 0

     do k = km, 3, -1

        v = u
        u = f
        f = ( a - 4.0D+00 * k * k ) * u / q - v

        if ( abs ( f ) < abs ( fc(k+1) ) ) then

           kb = k
           fc(1) = 1.0D-100
           sp = 0.0D+00
           f3 = fc(k+1)
           fc(2) = a / q * fc(1)
           fc(3) = ( a - 4.0D+00 ) * fc(2) / q - 2.0D+00 * fc(1)
           u = fc(2)
           f1 = fc(3)

           do i = 3, kb
              v = u
              u = f1
              f1 = ( a - 4.0D+00 * ( i - 1.0D+00 ) ** 2 ) * u / q - v
              fc(i+1) = f1
              if ( i == kb ) then
                 f2 = f1
              else
                 sp = sp + f1 * f1
              end if
           end do

           sp = sp + 2.0D+00 * fc(1) ** 2 + fc(2) ** 2 + fc(3) ** 2
           ss = s + sp * ( f3 / f2 ) ** 2
           s0 = sqrt ( 1.0D+00 / ss )
           do j = 1, km
              if ( j <= kb + 1 ) then
                 fc(j) = s0 * fc(j) * f3 / f2
              else
                 fc(j) = s0 * fc(j)
              end if
           end do
           l = 1
           exit
        else
           fc(k) = f
           s = s + f * f
        end if

     end do

     if ( l == 0 ) then
        fc(2) = q * fc(3) / ( a - 4.0D+00 - 2.0D+00 * q * q / a )
        fc(1) = q / a * fc(2)
        s = s + 2.0D+00 * fc(1) ** 2 + fc(2) ** 2
        s0 = sqrt ( 1.0D+00 / s )
        do k = 1, km
           fc(k) = s0 * fc(k)
        end do
     end if

  else if ( kd == 2 .or. kd == 3 ) then

     l = 0

     do k = km, 3, -1

        v = u
        u = f
        f = ( a - ( 2.0D+00 * k - 1 ) ** 2 ) * u / q - v

        if ( abs ( fc(k) ) <= abs ( f ) ) then
           fc(k-1) = f
           s = s + f * f
        else
           kb = k
           f3 = fc(k)
           l = 1
           exit
        end if

     end do

     if ( l == 0 ) then

        fc(1) = q / ( a - 1.0D+00 - ( - 1 ) ** kd * q ) * fc(2)
        s = s + fc(1) * fc(1)
        s0 = sqrt ( 1.0D+00 / s )
        do k = 1, km
           fc(k) = s0 * fc(k)
        end do

     else

        fc(1) = 1.0D-100
        fc(2) = ( a - 1.0D+00 - ( - 1 ) ** kd * q ) / q * fc(1)
        sp = 0.0D+00
        u = fc(1)
        f1 = fc(2)
        do i = 2, kb - 1
           v = u
           u = f1
           f1 = ( a - ( 2.0D+00 * i - 1.0D+00 ) ** 2 ) * u / q - v
           if ( i /= kb - 1 ) then
              fc(i+1) = f1
              sp = sp + f1 * f1
           else
              f2 = f1
           end if
        end do

        sp = sp + fc(1) ** 2 + fc(2) ** 2
        ss = s + sp * ( f3 / f2 ) ** 2
        s0 = 1.0D+00 / sqrt ( ss )
        do j = 1, km
           if ( j < kb ) then
              fc(j) = s0 * fc(j) * f3 / f2
           else
              fc(j) = s0 * fc(j)
           end if
        end do

     end if

  else if ( kd == 4 ) then

     l = 0

     do k = km, 3, -1
        v = u
        u = f
        f = ( a - 4.0D+00 * k * k ) * u / q - v
        if ( abs ( fc(k) ) <= abs ( f ) ) then
           fc(k-1) = f
           s = s + f * f
        else
           kb = k
           f3 = fc(k)
           l = 1
           exit
        end if
     end do

     if ( l == 0 ) then

        fc(1) = q / ( a - 4.0D+00 ) * fc(2)
        s = s + fc(1) * fc(1)
        s0 = sqrt ( 1.0D+00 / s )
        do k = 1, km
           fc(k) = s0 * fc(k)
        end do

     else

        fc(1) = 1.0D-100
        fc(2) = ( a - 4.0D+00 ) / q * fc(1)
        sp = 0.0D+00
        u = fc(1)
        f1 = fc(2)

        do i = 2, kb - 1
           v = u
           u = f1
           f1 = ( a - 4.0D+00 * i * i ) * u / q - v
           if ( i /= kb - 1 ) then
              fc(i+1) = f1
              sp = sp + f1 * f1
           else
              f2 = f1
           end if
        end do

        sp = sp + fc(1) ** 2 + fc(2) ** 2
        ss = s + sp * ( f3 / f2 ) ** 2
        s0 = 1.0D+00 / sqrt ( ss )

        do j = 1, km
           if ( j < kb ) then
              fc(j) = s0 * fc(j) * f3 / f2
           else
              fc(j) = s0 * fc(j)
           end if
        end do

     end if

  end if

  if ( fc(1) < 0.0D+00 ) then
     do j = 1, km
        fc(j) = -fc(j)
     end do
  end if

  return
end subroutine fcoef