scka Subroutine

subroutine scka ( m, n, c, cv, kd, ck )

  !*****************************************************************************80
  !
  !! SCKA: expansion coefficients for prolate and 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:
  !
  !    22 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 ) M, the mode parameter.
  !
  !    Input, integer ( kind = 4 ) N, the mode parameter.
  !
  !    Input, real ( kind = 8 ) C, the 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.
  !
  !    Output, real ( kind = 8 ) CK(*), the expansion coefficients.
  !    CK(1), CK(2),... correspond to c0, c2,..., and so on.
  !       
  implicit none

  real ( kind = 8 ) c
  real ( kind = 8 ) ck(200)
  real ( kind = 8 ) cs
  real ( kind = 8 ) cv
  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  real ( kind = 8 ) f2
  real ( kind = 8 ) fl
  real ( kind = 8 ) fs
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) k1
  integer ( kind = 4 ) kb
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nm
  real ( kind = 8 ) r1
  real ( kind = 8 ) r2
  real ( kind = 8 ) s0
  real ( kind = 8 ) su1
  real ( kind = 8 ) su2

  if ( c <= 1.0D-10 ) then
     c = 1.0D-10
  end if

  nm = 25 + int ( ( n - m ) / 2 + c )
  cs = c * c * kd

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

  fs = 1.0D+00
  f1 = 0.0D+00
  f0 = 1.0D-100
  kb = 0
  ck(nm+1) = 0.0D+00

  do k = nm, 1, -1

     f = ((( 2.0D+00 * k + m + ip ) &
          * ( 2.0D+00 * k + m + 1.0D+00 + ip ) - cv + cs ) * f0 &
          - 4.0D+00 * ( k + 1.0D+00 ) * ( k + m + 1.0D+00 ) * f1 ) / cs

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

        ck(k) = f
        f1 = f0
        f0 = f

        if ( 1.0D+100 < abs ( f ) ) then
           do k1 = nm, k, -1
              ck(k1) = ck(k1) * 1.0D-100
           end do
           f1 = f1 * 1.0D-100
           f0 = f0 * 1.0D-100
        end if

     else

        kb = k
        fl = ck(k+1)
        f1 = 1.0D+00
        f2 = 0.25D+00 * ( ( m + ip ) * ( m + ip + 1.0D+00 ) &
             - cv + cs ) / ( m + 1.0D+00 ) * f1
        ck(1) = f1

        if ( kb == 1 ) then
           fs = f2
        else if (kb == 2 ) then
           ck(2) = f2
           fs = 0.125D+00 * ( ( ( m + ip + 2.0D+00 ) &
                * ( m + ip + 3.0D+00 ) - cv + cs ) * f2 &
                - cs * f1 ) / ( m + 2.0D+00 )
        else
           ck(2) = f2
           do j = 3, kb + 1
              f = 0.25D+00 * ( ( ( 2.0D+00 * j + m + ip - 4.0D+00 ) &
                   * ( 2.0D+00 * j + m + ip - 3.0D+00 ) - cv + cs ) * f2 &
                   - cs * f1 ) / ( ( j - 1.0D+00 ) * ( j + m - 1.0D+00 ) )
              if ( j <= kb ) then
                 ck(j) = f
              end if
              f1 = f2
              f2 = f
           end do
           fs = f
        end if

        exit

     end if

  end do

  su1 = 0.0D+00
  do k = 1, kb
     su1 = su1 + ck(k)
  end do

  su2 = 0.0D+00
  do k = kb + 1, nm
     su2 = su2 + ck(k)
  end do

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

  r2 = 1.0D+00
  do j = 1, ( n - m - ip ) / 2
     r2 = - r2 * j
  end do

  if ( kb == 0 ) then
     s0 = r1 / ( 2.0D+00 ** n * r2 * su2 )
  else
     s0 = r1 / ( 2.0D+00 ** n * r2 * ( fl / fs * su1 + su2 ) )
  end if

  do k = 1, kb
     ck(k) = fl / fs * s0 * ck(k)
  end do

  do k = kb + 1, nm
     ck(k) = s0 * ck(k)
  end do

  return
end subroutine scka