rmn2l Subroutine

subroutine rmn2l ( m, n, c, x, df, kd, r2f, r2d, id )

  !*****************************************************************************80
  !
  !! RMN2L: prolate and oblate spheroidal functions, second kind, large CX.
  !
  !  Discussion:
  !
  !    This procedure computes prolate and oblate spheroidal radial functions 
  !    of the second kind for given m, n, c and a large cx.
  !
  !  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:
  !
  !    30 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 ) X, the argument.
  !
  !    Input, real ( kind = 8 ) DF(*), the expansion coefficients.
  !
  !    Input, integer ( kind = 4 ) KD, the function code.
  !    1, the prolate function.
  !    -1, the oblate function.
  !
  !    Output, real ( kind = 8 ) R2F, R2D, the function and derivative values.
  !
  implicit none

  real ( kind = 8 ) a0
  real ( kind = 8 ) b0
  real ( kind = 8 ) c
  real ( kind = 8 ) cx
  real ( kind = 8 ) df(200)
  real ( kind = 8 ) dy(0:251)
  real ( kind = 8 ) eps
  real ( kind = 8 ) eps1
  real ( kind = 8 ) eps2
  integer ( kind = 4 ) id
  integer ( kind = 4 ) id1
  integer ( kind = 4 ) id2
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) l
  integer ( kind = 4 ) lg
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nm
  integer ( kind = 4 ) nm1
  integer ( kind = 4 ) nm2
  integer ( kind = 4 ) np
  real ( kind = 8 ) r
  real ( kind = 8 ) r0
  real ( kind = 8 ) r2d
  real ( kind = 8 ) r2f
  real ( kind = 8 ) reg
  real ( kind = 8 ) sw
  real ( kind = 8 ) suc
  real ( kind = 8 ) sud
  real ( kind = 8 ) sy(0:251)
  real ( kind = 8 ) x

  eps = 1.0D-14

  nm1 = int ( ( n - m ) / 2 )

  if ( n - m == 2 * nm1 ) then
     ip = 0
  else
     ip = 1
  end if
  nm = 25 + nm1 + int ( c )

  if ( 80 < m + nm ) then
     reg = 1.0D-200
  else
     reg = 1.0D+00
  end if
  nm2 = 2 * nm + m
  cx = c * x
  call sphy ( nm2, cx, nm2, sy, dy )
  r0 = reg
  do j = 1, 2 * m + ip
     r0 = r0 * j
  end do
  r = r0    
  suc = r * df(1)
  do k = 2, nm
     r = r * ( m + k - 1.0D+00 ) * ( m + k + ip - 1.5D+00 ) &
          / ( k - 1.0D+00 ) / ( k + ip - 1.5D+00 )
     suc = suc + r * df(k)
     if ( nm1 < k .and. abs ( suc - sw ) < abs ( suc ) * eps ) then
        exit
     end if
     sw = suc
  end do

  a0 = ( 1.0D+00 - kd / ( x * x ) ) ** ( 0.5D+00 * m ) / suc
  r2f = 0.0D+00
  do k = 1, nm
     l = 2 * k + m - n - 2 + ip
     if ( l == 4 * int ( l / 4 ) ) then
        lg = 1
     else
        lg = -1
     end if

     if ( k == 1 ) then
        r = r0
     else
        r = r * ( m + k - 1.0D+00 ) * ( m + k + ip - 1.5D+00 ) &
             / ( k - 1.0D+00 ) / ( k + ip - 1.5D+00 )
     end if

     np = m + 2 * k - 2 + ip
     r2f = r2f + lg * r * ( df(k) * sy(np) )
     eps1 = abs ( r2f - sw )
     if ( nm1 < k .and. eps1 < abs ( r2f ) * eps ) then
        exit
     end if
     sw = r2f
  end do

  id1 = int ( log10 ( eps1 / abs ( r2f ) + eps ) )
  r2f = r2f * a0

  if ( nm2 <= np ) then
     id = 10
     return
  end if

  b0 = kd * m / x ** 3.0D+00 / ( 1.0D+00 - kd / ( x * x ) ) * r2f                
  sud = 0.0D+00
  do k = 1, nm
     l = 2 * k + m - n - 2 + ip
     if ( l == 4 * int ( l / 4 ) ) then
        lg = 1
     else
        lg = -1
     end if
     if (k == 1) then
        r = r0
     else
        r = r * ( m + k - 1.0D+00 ) * ( m + k + ip - 1.5D+00 ) &
             / ( k - 1.0D+00 ) / ( k + ip - 1.5D+00 )
     end if
     np = m + 2 * k - 2 + ip
     sud = sud + lg * r * ( df(k) * dy(np) )
     eps2 = abs ( sud - sw )
     if ( nm1 < k .and. eps2 < abs ( sud ) * eps ) then
        exit
     end if
     sw = sud
  end do

  r2d = b0 + a0 * c * sud
  id2 = int ( log10 ( eps2 / abs ( sud ) + eps ) )
  id = max ( id1, id2 )

  return
end subroutine rmn2l