aswfa Subroutine

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

  !*****************************************************************************80
  !
  !! ASWFA: 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:
  !
  !    13 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, with N = M, M+1, ...
  !
  !    Input, real ( kind = 8 ) C, the spheroidal parameter.
  !
  !    Input, real ( kind = 8 ) X, the argument of the angular function.
  !    |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 ) a0
  real ( kind = 8 ) c
  real ( kind = 8 ) ck(200)
  real ( kind = 8 ) cv
  real ( kind = 8 ) d0
  real ( kind = 8 ) d1
  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 ) n
  integer ( kind = 4 ) nm
  integer ( kind = 4 ) nm2
  real ( kind = 8 ) r
  real ( kind = 8 ) s1d
  real ( kind = 8 ) s1f
  real ( kind = 8 ) su1
  real ( kind = 8 ) su2
  real ( kind = 8 ) x
  real ( kind = 8 ) x0
  real ( kind = 8 ) x1

  eps = 1.0D-14
  x0 = x
  x = abs ( x )

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

  nm = 10 + int ( ( n - m ) / 2 + c )
  nm2 = nm / 2 - 2 
  call sdmn ( m, n, c, cv, kd, df )
  call sckb ( m, n, c, df, ck )
  x1 = 1.0D+00 - x * x

  if ( m == 0 .and. x1 == 0.0D+00 ) then
     a0 = 1.0D+00
  else
     a0 = x1 ** ( 0.5D+00 * m )
  end if

  su1 = ck(1)
  do k = 1, nm2
     r = ck(k+1) * x1 ** k
     su1 = su1 + r
     if ( 10 <= k .and. abs ( r / su1 ) < eps ) then
        exit
     end if
  end do

  s1f = a0 * x ** ip * su1

  if ( x == 1.0D+00 ) then

     if ( m == 0 ) then
        s1d = ip * ck(1) - 2.0D+00 * ck(2)
     else if ( m == 1 ) then
        s1d = -1.0D+100
     else if ( m == 2 ) then
        s1d = -2.0D+00 * ck(1)
     else if ( 3 <= m ) then
        s1d = 0.0D+00
     end if

  else

     d0 = ip - m / x1 * x ** ( ip + 1.0D+00 )
     d1 = -2.0D+00 * a0 * x ** ( ip + 1.0D+00 )
     su2 = ck(2)
     do k = 2, nm2
        r = k * ck(k+1) * x1 ** ( k - 1.0D+00 )
        su2 = su2 + r
        if ( 10 <= k .and. abs ( r / su2 ) < eps ) then
           exit
        end if
     end do

     s1d = d0 * a0 * su1 + d1 * su2

  end if

  if ( x0 < 0.0D+00 ) then
     if ( ip == 0 ) then
        s1d = -s1d
     else if ( ip == 1 ) then
        s1f = -s1f
     end if
  end if

  x = x0

  return
end subroutine aswfa