qstar Subroutine

subroutine qstar ( m, n, c, ck, ck1, qs, qt )

  !*****************************************************************************80
  !
  !! QSTAR computes Q*mn(-ic) for oblate radial functions with a 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:
  !
  !    18 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 ) CK(*), ?
  !
  !    Input, real ( kind = 8 ) CK1, ?
  !
  !    Output, real ( kind = 8 ) QS, ?
  !
  !    Output, real ( kind = 8 ) QT, ?
  !
  implicit none

  real ( kind = 8 ) ap(200)
  real ( kind = 8 ) c
  real ( kind = 8 ) ck(200)
  real ( kind = 8 ) ck1
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  real ( kind = 8 ) qs
  real ( kind = 8 ) qs0
  real ( kind = 8 ) qt
  real ( kind = 8 ) r
  real ( kind = 8 ) s
  real ( kind = 8 ) sk

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

  r = 1.0D+00 / ck(1) ** 2
  ap(1) = r
  do i = 1, m
     s = 0.0D+00
     do l = 1, i
        sk = 0.0D+00
        do k = 0, l
           sk = sk + ck(k+1) * ck(l-k+1)
        end do
        s = s + sk * ap(i-l+1)
     end do
     ap(i+1) = -r * s
  end do

  qs0 = ap(m+1)     
  do l = 1, m
     r = 1.0D+00
     do k = 1, l
        r = r * ( 2.0D+00 * k + ip ) &
             * ( 2.0D+00 * k - 1.0D+00 + ip ) / ( 2.0D+00 * k ) ** 2
     end do
     qs0 = qs0 + ap(m-l+1) * r
  end do

  qs = ( -1.0D+00 ) ** ip * ck1 * ( ck1 * qs0 ) / c
  qt = - 2.0D+00 / ck1 * qs

  return
end subroutine qstar