kmn Subroutine

subroutine kmn ( m, n, c, cv, kd, df, dn, ck1, ck2 )

  !*****************************************************************************80
  !
  !! KMN: expansion coefficients of prolate or oblate spheroidal 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:
  !
  !    02 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 ) M, the mode parameter;  M = 0, 1, 2, ...
  !
  !    Input, integer ( kind = 4 ) N, mode parameter, N = M, M + 1, M + 2, ...
  !
  !    Input, real ( kind = 8 ) C, spheroidal parameter.
  !
  !    Input, real ( kind = 8 ) CV, the characteristic value.
  !
  !    Input, integer ( kind = 4 ) KD, the function code.
  !    1, the prolate function.
  !    -1, the oblate function.
  !
  !    Input, real ( kind = 8 ) DF(*), the expansion coefficients.
  !
  implicit none

  real ( kind = 8 ) c
  real ( kind = 8 ) ck1
  real ( kind = 8 ) ck2
  real ( kind = 8 ) cs
  real ( kind = 8 ) cv
  real ( kind = 8 ) df(200)
  real ( kind = 8 ) dn(200)
  real ( kind = 8 ) dnp
  real ( kind = 8 ) g0
  real ( kind = 8 ) gk0
  real ( kind = 8 ) gk1
  real ( kind = 8 ) gk2
  real ( kind = 8 ) gk3
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nm
  integer ( kind = 4 ) nm1
  integer ( kind = 4 ) nn
  real ( kind = 8 ) r
  real ( kind = 8 ) r1
  real ( kind = 8 ) r2
  real ( kind = 8 ) r3
  real ( kind = 8 ) r4
  real ( kind = 8 ) r5
  real ( kind = 8 ) rk(200)
  real ( kind = 8 ) sa0
  real ( kind = 8 ) sb0
  real ( kind = 8 ) su0
  real ( kind = 8 ) sw
  real ( kind = 8 ) t
  real ( kind = 8 ) tp(200)
  real ( kind = 8 ) u(200)
  real ( kind = 8 ) v(200)
  real ( kind = 8 ) w(200)

  nm = 25 + int ( 0.5D+00 * ( n - m ) + c )
  nn = nm + m
  cs = c * c * kd

  if ( n - m == 2 * int ( ( n - m ) / 2 ) ) then
     ip = 0
  else
     ip = 1
  end if

  do i = 1, nn + 3

     if ( ip == 0 ) then
        k = - 2 * ( i - 1 )
     else
        k = - ( 2 * i - 3 )
     end if

     gk0 = 2.0D+00 * m + k
     gk1 = ( m + k ) * ( m + k + 1.0D+00 )
     gk2 = 2.0D+00 * ( m + k ) - 1.0D+00
     gk3 = 2.0D+00 * ( m + k ) + 3.0D+00
     u(i) = gk0 * ( gk0 - 1.0D+00 ) * cs / ( gk2 * ( gk2 + 2.0D+00 ) )
     v(i) = gk1 - cv + ( 2.0D+00 * ( gk1 - m * m ) - 1.0D+00 ) * cs &
          / ( gk2 * gk3 )
     w(i) = ( k + 1.0D+00 ) * ( k + 2.0D+00 ) * cs / ( ( gk2 + 2.0D+00 ) * gk3 )

  end do

  do k = 1, m
     t = v(m+1)
     do l = 0, m - k - 1
        t = v(m-l) - w(m-l+1) * u(m-l) / t
     end do
     rk(k) = -u(k) / t
  end do

  r = 1.0D+00
  do k = 1, m
     r = r * rk(k)
     dn(k) = df(1) * r
  end do

  tp(nn) = v(nn+1)
  do k = nn - 1, m + 1,-1
     tp(k) = v(k+1) - w(k+2) * u(k+1) / tp(k+1)
     if ( m + 1 < k ) then
        rk(k) = -u(k) / tp(k)
     end if
  end do

  if ( m == 0 ) then
     dnp = df(1)
  else
     dnp = dn(m)
  end if

  dn(m+1) = ( - 1.0D+00 ) ** ip * dnp * cs &
       / ( ( 2.0D+00 * m - 1.0D+00 ) &
       * ( 2.0D+00 * m + 1.0D+00 - 4.0D+00 * ip ) * tp(m+1) )
  do k = m + 2, nn
     dn(k) = rk(k) * dn(k-1)
  end do

  r1 = 1.0D+00
  do j = 1, ( n + m + ip ) / 2
     r1 = r1 * ( j + 0.5D+00 * ( n + m + ip ) )
  end do
  nm1 = ( n - m ) / 2
  r = 1.0D+00
  do j = 1, 2 * m + ip
     r = r * j
  end do
  su0 = r * df(1)

  do k = 2, nm
     r = r * ( m + k - 1.0D+00 ) * ( m + k + ip - 1.5D+00 ) &
          / ( k - 1.0D+00 ) / ( k + ip - 1.5D+00 )
     su0 = su0 + r * df(k)
     if ( nm1 < k .and. &
          abs ( ( su0 - sw ) / su0 ) < 1.0D-14 ) then
        exit
     end if
     sw = su0
  end do

  if ( kd /= 1 ) then

     r2 = 1.0D+00
     do j = 1,m
        r2 = 2.0D+00 * c * r2 * j
     end do
     r3 = 1.0D+00
     do j = 1, ( n - m - ip ) / 2
        r3 = r3 * j
     end do
     sa0 = ( 2.0D+00 * ( m + ip ) + 1.0D+00 ) * r1 &
          / ( 2.0D+00 ** n * c ** ip * r2 * r3 * df(1) )
     ck1 = sa0 * su0

     if ( kd == -1 ) then
        return
     end if

  end if

  r4 = 1.0D+00
  do j = 1, ( n - m - ip ) / 2
     r4 = 4.0D+00 * r4 * j
  end do
  r5 = 1.0D+00
  do j = 1, m
     r5 = r5 * ( j + m ) / c
  end do

  if ( m == 0 ) then
     g0 = df(1)
  else
     g0 = dn(m)
  end if

  sb0 = ( ip + 1.0D+00 ) * c ** ( ip + 1 ) &
       / ( 2.0D+00 * ip * ( m - 2.0D+00 ) + 1.0D+00 ) &
       / ( 2.0D+00 * m - 1.0D+00 )

  ck2 = ( -1 ) ** ip * sb0 * r4 * r5 * g0 / r1 * su0

  return
end subroutine kmn