rmn2sp Subroutine

subroutine rmn2sp ( m, n, c, x, cv, df, kd, r2f, r2d )

  !*****************************************************************************80
  !
  !! RMN2SP: prolate, oblate spheroidal radial functions, kind 2, 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:
  !
  !    28 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 ) CV, the characteristic value.
  !
  !    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 values of the function and 
  !    its derivative.
  !
  implicit none

  real ( kind = 8 ) c
  real ( kind = 8 ) ck1
  real ( kind = 8 ) ck2
  real ( kind = 8 ) cv
  real ( kind = 8 ) df(200)
  real ( kind = 8 ) dn(200)
  real ( kind = 8 ) eps
  real ( kind = 8 ) ga
  real ( kind = 8 ) gb
  real ( kind = 8 ) gc
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) j
  integer ( kind = 4 ) j1
  integer ( kind = 4 ) j2
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) ki
  integer ( kind = 4 ) l1
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nm
  integer ( kind = 4 ) nm1
  integer ( kind = 4 ) nm2
  integer ( kind = 4 ) nm3
  real ( kind = 8 ) pd(0:251)
  real ( kind = 8 ) pm(0:251)
  real ( kind = 8 ) qd(0:251)
  real ( kind = 8 ) qm(0:251)
  real ( kind = 8 ) r1
  real ( kind = 8 ) r2
  real ( kind = 8 ) r2d
  real ( kind = 8 ) r2f
  real ( kind = 8 ) r3
  real ( kind = 8 ) r4
  real ( kind = 8 ) sd
  real ( kind = 8 ) sd0
  real ( kind = 8 ) sd1
  real ( kind = 8 ) sd2
  real ( kind = 8 ) sdm
  real ( kind = 8 ) sf
  real ( kind = 8 ) spd1
  real ( kind = 8 ) spd2
  real ( kind = 8 ) spl
  real ( kind = 8 ) su0
  real ( kind = 8 ) su1
  real ( kind = 8 ) su2
  real ( kind = 8 ) sum
  real ( kind = 8 ) sw
  real ( kind = 8 ) x

  if ( abs ( df(1) ) < 1.0D-280 ) then
     r2f = 1.0D+300
     r2d = 1.0D+300
     return
  end if

  eps = 1.0D-14

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

  if ( n - m .eq. 2 * nm1 ) then
     ip = 0
  else
     ip = 1
  end if

  nm = 25 + nm1 + int ( c )
  nm2 = 2 * nm + m
  call kmn ( m, n, c, cv, kd, df, dn, ck1, ck2 )
  call lpmns ( m, nm2, x, pm, pd )
  call lqmns ( m, nm2, x, qm, qd )

  su0 = 0.0D+00
  do k = 1, nm
     j = 2 * k - 2 + m + ip
     su0 = su0 + df(k) * qm(j)
     if ( nm1 < k .and. abs ( su0 - sw ) < abs ( su0 ) * eps ) then
        exit                                                               
     end if
     sw = su0
  end do

  sd0 = 0.0D+00

  do k = 1, nm
     j = 2 * k - 2 + m + ip
     sd0 = sd0 + df(k) * qd(j)
     if ( nm1 < k .and. abs ( sd0 - sw ) < abs ( sd0 ) * eps ) then
        exit
     end if
     sw = sd0
  end do

  su1 = 0.0D+00
  sd1 = 0.0D+00
  do k = 1, m
     j = m - 2 * k + ip
     if ( j < 0 ) then
        j = - j - 1
     end if
     su1 = su1 + dn(k) * qm(j)
     sd1 = sd1 + dn(k) * qd(j)
  end do

  ga = ( ( x - 1.0D+00 ) / ( x + 1.0D+00 ) ) ** ( 0.5D+00 * m )

  do k = 1, m

     j = m - 2 * k + ip

     if ( 0 <= j ) then
        cycle
     end if

     if ( j < 0 ) then
        j = - j - 1
     end if
     r1 = 1.0D+00
     do j1 = 1, j
        r1 = ( m + j1 ) * r1
     end do
     r2 = 1.0D+00
     do j2 = 1, m - j - 2
        r2 = j2 * r2
     end do
     r3 = 1.0D+00
     sf = 1.0D+00
     do l1 = 1, j
        r3 = 0.5D+00 * r3 * ( - j + l1 - 1.0D+00 ) * ( j + l1 ) &
             / ( ( m + l1 ) * l1 ) * ( 1.0D+00 - x )
        sf = sf + r3
     end do

     if ( m - j <= 1 ) then
        gb = 1.0D+00
     else
        gb = ( m - j - 1.0D+00 ) * r2
     end if

     spl = r1 * ga * gb * sf
     su1 = su1 + ( -1 ) ** ( j + m ) * dn(k) * spl
     spd1 = m / ( x * x - 1.0D+00 ) * spl
     gc = 0.5D+00 * j * ( j + 1.0 ) / ( m + 1.0D+00 )
     sd = 1.0D+00
     r4 = 1.0D+00
     do l1 = 1, j - 1
        r4 = 0.5D+00 * r4 * ( - j + l1 ) * ( j + l1 + 1.0D+00 ) &
             / ( ( m + l1 + 1.0D+00 ) * l1 ) * ( 1.0D+00 - x )
        sd = sd + r4
     end do

     spd2 = r1 * ga * gb * gc * sd
     sd1 = sd1 + ( - 1 ) ** ( j + m ) * dn(k) * ( spd1 + spd2 )

  end do

  su2 = 0.0D+00
  ki = ( 2 * m + 1 + ip ) / 2
  nm3 = nm + ki
  do k = ki, nm3
     j = 2 * k - 1 - m - ip
     su2 = su2 + dn(k) * pm(j)
     if ( m < j .and. &
          abs ( su2 - sw ) < abs ( su2 ) * eps ) then
        exit
     end if
     sw = su2
  end do

  sd2 = 0.0D+00

  do k = ki, nm3
     j = 2 * k - 1 - m - ip
     sd2 = sd2 + dn(k) * pd(j)
     if ( m < j .and. &
          abs ( sd2 - sw ) < abs ( sd2 ) * eps ) then
        exit
     end if
     sw = sd2
  end do

  sum = su0 + su1 + su2
  sdm = sd0 + sd1 + sd2
  r2f = sum / ck2
  r2d = sdm / ck2

  return
end subroutine rmn2sp