rmn2so Subroutine

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

  !*****************************************************************************80
  !
  !! RMN2SO: oblate radial functions of the second kind with small argument.
  !
  !  Discussion:
  !
  !    This procedure computes oblate radial functions of the second kind
  !    with a small argument, Rmn(-ic,ix) and Rmn'(-ic,ix).
  !
  !  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:
  !
  !    27 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 Rmn(-ic,ix) 
  !    and Rmn'(-ic,ix).
  !
  implicit none

  real ( kind = 8 ) bk(200)
  real ( kind = 8 ) c
  real ( kind = 8 ) ck(200)
  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 ) gd
  real ( kind = 8 ) gf
  real ( kind = 8 ) h0
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) j
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nm
  real ( kind = 8 ) pi
  real ( kind = 8 ) qs
  real ( kind = 8 ) qt
  real ( kind = 8 ) r1d
  real ( kind = 8 ) r1f
  real ( kind = 8 ) r2d
  real ( kind = 8 ) r2f
  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
  pi = 3.141592653589793D+00
  nm = 25 + int ( ( n - m ) / 2 + c )
  if ( n - m == 2 * int ( ( n - m ) / 2 ) ) then
     ip = 0
  else
     ip = 1
  end if

  call sckb ( m, n, c, df, ck )
  call kmn ( m, n, c, cv, kd, df, dn, ck1, ck2 )
  call qstar ( m, n, c, ck, ck1, qs, qt )
  call cbk ( m, n, c, cv, qt, ck, bk )

  if ( x == 0.0D+00 ) then

     sum = 0.0D+00
     do j = 1, nm
        sum = sum + ck(j)
        if ( abs ( sum - sw ) < abs ( sum ) * eps ) then
           exit
        end if
        sw = sum
     end do

     if ( ip == 0 ) then
        r1f = sum / ck1
        r2f = - 0.5D+00 * pi * qs * r1f
        r2d = qs * r1f + bk(1)
     else if ( ip == 1 ) then
        r1d = sum / ck1
        r2f = bk(1)
        r2d = -0.5D+00 * pi * qs * r1d
     end if

     return

  else

     call gmn ( m, n, c, x, bk, gf, gd )
     call rmn1 ( m, n, c, x, df, kd, r1f, r1d )
     h0 = atan ( x ) - 0.5D+00 * pi
     r2f = qs * r1f * h0 + gf
     r2d = qs * ( r1d * h0 + r1f / ( 1.0D+00 + x * x ) ) + gd

  end if

  return
end subroutine rmn2so