rmn1 Subroutine

subroutine rmn1 ( m, n, c, x, df, kd, r1f, r1d )

  !*****************************************************************************80
  !
  !! RMN1 computes prolate and oblate spheroidal functions of the first kind.
  !
  !  Discussion:
  !
  !    This procedure computes prolate and oblate spheroidal radial
  !    functions of the first kind for given m, n, c and x.
  !
  !  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:
  !
  !    29 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 ) R1F, R1D, the function and derivative.
  !
  implicit none

  real ( kind = 8 ) a0
  real ( kind = 8 ) b0
  real ( kind = 8 ) c
  real ( kind = 8 ) ck(200)
  real ( kind = 8 ) cx
  real ( kind = 8 ) df(200)
  real ( kind = 8 ) dj(0:251)
  real ( kind = 8 ) eps
  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 ) r1
  real ( kind = 8 ) r1d
  real ( kind = 8 ) r1f
  real ( kind = 8 ) r2
  real ( kind = 8 ) r3
  real ( kind = 8 ) reg
  real ( kind = 8 ) sa0
  real ( kind = 8 ) sj(0:251)
  real ( kind = 8 ) suc
  real ( kind = 8 ) sud
  real ( kind = 8 ) sum
  real ( kind = 8 ) sw
  real ( kind = 8 ) sw1
  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 )
  reg = 1.0D+00
  if ( 80 < m + nm ) then
     reg = 1.0D-200
  end if
  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

  if ( x == 0.0D+00 ) then

     call sckb ( m, n, c, df, ck )
     sum = 0.0D+00
     do j = 1, nm
        sum = sum + ck(j)
        if ( abs ( sum - sw1 ) < abs ( sum ) * eps ) then
           exit
        end if
        sw1 = sum
     end do

     r1 = 1.0D+00
     do j = 1, ( n + m + ip ) / 2
        r1 = r1 * ( j + 0.5D+00 * ( n + m + ip ) )
     end do

     r2 = 1.0D+00
     do j = 1, m
        r2 = 2.0D+00 * c * r2 * j
     end do

     r3 = 1.0D+00
     do j = 1, ( n - m - ip ) / 2
        r3 = r3 * j
     end do

     sa0 = ( 2.0D+00 * ( m + ip ) + 1.0D+00 ) * r1 &
          / ( 2.0D+00 ** n * c ** ip * r2 * r3 )

     if ( ip == 0 ) then
        r1f = sum / ( sa0 * suc ) * df(1) * reg
        r1d = 0.0D+00
     else if ( ip == 1 ) then
        r1f = 0.0D+00
        r1d = sum / ( sa0 * suc ) * df(1) * reg
     end if

     return

  end if

  cx = c * x
  nm2 = 2 * nm + m
  call sphj ( nm2, cx, nm2, sj, dj )
  a0 = ( 1.0D+00 - kd / ( x * x ) ) ** ( 0.5D+00 * m ) / suc  
  r1f = 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
     r1f = r1f + lg * r * df(k) * sj(np)
     if ( nm1 < k .and. abs ( r1f - sw ) < abs ( r1f ) * eps ) then
        exit
     end if
     sw = r1f
  end do

  r1f = r1f * a0
  b0 = kd * m / x ** 3.0D+00 / ( 1.0D+00 - kd / ( x * x ) ) * r1f    
  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) * dj(np)
     if ( nm1 < k .and. abs ( sud - sw ) < abs ( sud ) * eps ) then
        exit
     end if
     sw = sud
  end do

  r1d = b0 + a0 * c * sud

  return
end subroutine rmn1