cbk Subroutine

subroutine cbk ( m, n, c, cv, qt, ck, bk )

  !*****************************************************************************80
  !
  !! CBK computes coefficients for oblate radial functions with small argument.
  !
  !  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:
  !
  !    20 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;  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, real ( kind = 8 ) QT, ?
  !
  !    Input, real ( kind = 8 ) CK(*), ?
  !
  !    Output, real ( kind = 8 ) BK(*), the coefficients.
  !
  implicit none

  real ( kind = 8 ) bk(200)
  real ( kind = 8 ) c
  real ( kind = 8 ) ck(200)
  real ( kind = 8 ) cv
  real ( kind = 8 ) eps
  integer ( kind = 4 ) i
  integer ( kind = 4 ) i1
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) n2
  integer ( kind = 4 ) nm
  real ( kind = 8 ) qt
  real ( kind = 8 ) r1
  real ( kind = 8 ) s1
  real ( kind = 8 ) sw
  real ( kind = 8 ) t
  real ( kind = 8 ) u(200)
  real ( kind = 8 ) v(200)
  real ( kind = 8 ) w(200)

  eps = 1.0D-14
  if ( n - m == 2 * int ( ( n - m ) / 2 ) ) then
     ip = 0
  else
     ip = 1
  end if
  nm = 25 + int ( 0.5D+00 * ( n - m ) + c )
  u(1) = 0.0D+00
  n2 = nm - 2
  do j = 2, n2
     u(j) = c * c
  end do

  do j = 1, n2
     v(j) = ( 2.0D+00 * j - 1.0D+00 - ip ) &
          * ( 2.0D+00 * ( j - m ) - ip ) + m * ( m - 1.0D+00 ) - cv
  end do

  do j = 1, nm - 1
     w(j) = ( 2.0D+00 * j - ip ) * ( 2.0D+00 * j + 1.0D+00 - ip )
  end do

  if ( ip == 0 ) then

     do k = 0, n2 - 1

        s1 = 0.0D+00
        i1 = k - m + 1

        do i = i1, nm
           if ( 0 <= i ) then
              r1 = 1.0D+00
              do j = 1, k
                 r1 = r1 * ( i + m - j ) / j
              end do
              s1 = s1 + ck(i+1) * ( 2.0D+00 * i + m ) * r1
              if ( abs ( s1 - sw ) < abs ( s1 ) * eps ) then
                 exit
              end if
              sw = s1
           end if
        end do

        bk(k+1) = qt * s1

     end do

  else if ( ip == 1 ) then

     do k = 0, n2 - 1

        s1 = 0.0D+00
        i1 = k - m + 1

        do i = i1, nm

           if ( 0 <= i ) then

              r1 = 1.0D+00
              do j = 1, k
                 r1 = r1 * ( i + m - j ) / j
              end do

              if ( 0 < i ) then
                 s1 = s1 + ck(i) * ( 2.0D+00 * i + m - 1 ) * r1
              end if
              s1 = s1 - ck(i+1) * ( 2.0D+00 * i + m ) * r1
              if ( abs ( s1 - sw ) < abs ( s1 ) * eps ) then
                 exit
              end if
              sw = s1

           end if

        end do

        bk(k+1) = qt * s1

     end do

  end if

  w(1) = w(1) / v(1)
  bk(1) = bk(1) / v(1)
  do k = 2, n2
     t = v(k) - w(k-1) * u(k)
     w(k) = w(k) / t
     bk(k) = ( bk(k) - bk(k-1) * u(k) ) / t
  end do

  do k = n2 - 1, 1, -1
     bk(k) = bk(k) - w(k) * bk(k+1)
  end do

  return
end subroutine cbk