sckb Subroutine

subroutine sckb ( m, n, c, df, ck )

  !*****************************************************************************80
  !
  !! SCKB: 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:
  !
  !    15 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 ) DF(*), the expansion coefficients DK.
  !
  !    Output, real ( kind = 8 ) CK(*), the expansion coefficients CK.
  !
  implicit none

  real ( kind = 8 ) c
  real ( kind = 8 ) ck(200)
  real ( kind = 8 ) d1
  real ( kind = 8 ) d2
  real ( kind = 8 ) d3
  real ( kind = 8 ) df(200)
  real ( kind = 8 ) fac
  integer ( kind = 4 ) i
  integer ( kind = 4 ) i1
  integer ( kind = 4 ) i2
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nm
  real ( kind = 8 ) r
  real ( kind = 8 ) r1
  real ( kind = 8 ) reg
  real ( kind = 8 ) sum
  real ( kind = 8 ) sw

  c = max ( c, 1.0D-10 )

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

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

  if ( 80 < m + nm ) then
     reg = 1.0D-200
  else
     reg = 1.0D+00
  end if

  fac = - 0.5D+00 ** m

  do k = 0, nm - 1

     fac = - fac
     i1 = 2 * k + ip + 1
     r = reg
     do i = i1, i1 + 2 * m - 1
        r = r * i
     end do

     i2 = k + m + ip
     do i = i2, i2 + k - 1
        r = r * ( i + 0.5D+00 )
     end do

     sum = r * df(k+1)
     do i = k + 1, nm
        d1 = 2.0D+00 * i + ip
        d2 = 2.0D+00 * m + d1
        d3 = i + m + ip - 0.5D+00
        r = r * d2 * ( d2 - 1.0D+00 ) * i * ( d3 + k ) &
             / ( d1 * ( d1 - 1.0D+00 ) * ( i - k ) * d3 )
        sum = sum + r * df(i+1)
        if ( abs ( sw - sum ) < abs ( sum ) * 1.0D-14 ) then
           exit
        end if
        sw = sum
     end do

     r1 = reg
     do i = 2, m + k
        r1 = r1 * i
     end do

     ck(k+1) = fac * sum / r1

  end do

  return
end subroutine sckb