aswfb Subroutine

subroutine aswfb ( m, n, c, x, kd, cv, s1f, s1d )

  !*****************************************************************************80
  !
  !! ASWFB: prolate and oblate spheroidal angular functions of the first kind.
  !
  !  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, the spheroidal parameter.
  !
  !    Input, real ( kind = 8 ) X, the argument, with |X| < 1.0.
  !
  !    Input, integer ( kind = 4 ) KD, the function code.
  !    1, the prolate function.
  !    -1, the oblate function.
  !
  !    Input, real ( kind = 8 ) CV, the characteristic value.
  !
  !    Output, real ( kind = 8 ) S1F, S1D, the angular function of the first
  !    kind and its derivative.
  !
  implicit none

  real ( kind = 8 ) c
  real ( kind = 8 ) cv
  real ( kind = 8 ) df(200)
  real ( kind = 8 ) eps
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) m
  integer ( kind = 4 ) mk
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nm
  integer ( kind = 4 ) nm2
  real ( kind = 8 ) pd(0:251)
  real ( kind = 8 ) pm(0:251)
  real ( kind = 8 ) s1d
  real ( kind = 8 ) s1f
  real ( kind = 8 ) su1
  real ( kind = 8 ) sw
  real ( kind = 8 ) x

  eps = 1.0D-14

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

  nm = 25 + int ( ( n - m ) / 2 + c )
  nm2 = 2 * nm + m
  call sdmn ( m, n, c, cv, kd, df )
  call lpmns ( m, nm2, x, pm, pd )
  su1 = 0.0D+00
  do k = 1, nm
     mk = m + 2 * ( k - 1 ) + ip
     su1 = su1 + df(k) * pm(mk)
     if ( abs ( sw - su1 ) < abs ( su1 ) * eps ) then
        exit
     end if
     sw = su1
  end do

  s1f = ( -1.0D+00 ) ** m * su1

  su1 = 0.0D+00
  do k = 1, nm
     mk = m + 2 * ( k - 1 ) + ip
     su1 = su1 + df(k) * pd(mk)
     if ( abs ( sw - su1 ) < abs ( su1 ) * eps ) then
        exit
     end if
     sw = su1
  end do

  s1d = ( -1.0D+00 ) ** m * su1

  return
end subroutine aswfb