gmn Subroutine

subroutine gmn ( m, n, c, x, bk, gf, gd )

  !*****************************************************************************80
  !
  !! GMN computes quantities for oblate radial functions with small argument.
  !
  !  Discussion:
  !
  !    This procedure computes Gmn(-ic,ix) and its derivative for oblate
  !    radial functions with a 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:
  !
  !    15 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 ) BK(*), coefficients.
  !
  !    Output, real ( kind = 8 ) GF, GD, the value of Gmn(-C,X) and Gmn'(-C,X).
  !
  implicit none

  real ( kind = 8 ) bk(200)
  real ( kind = 8 ) c
  real ( kind = 8 ) eps
  real ( kind = 8 ) gd
  real ( kind = 8 ) gd0
  real ( kind = 8 ) gd1
  real ( kind = 8 ) gf
  real ( kind = 8 ) gf0
  real ( kind = 8 ) gw
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nm
  real ( kind = 8 ) x
  real ( kind = 8 ) xm

  eps = 1.0D-14

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

  nm = 25 + int ( 0.5D+00 * ( n - m ) + c )
  xm = ( 1.0D+00 + x * x ) ** ( -0.5D+00 * m )
  gf0 = 0.0D+00
  do k = 1, nm
     gf0 = gf0 + bk(k) * x ** ( 2.0D+00 * k - 2.0D+00 )
     if ( abs ( ( gf0 - gw ) / gf0 ) < eps .and. 10 <= k ) then
        exit
     end if
     gw = gf0
  end do

  gf = xm * gf0 * x ** ( 1 - ip )

  gd1 = - m * x / ( 1.0D+00 + x * x ) * gf
  gd0 = 0.0D+00

  do k = 1, nm

     if ( ip == 0 ) then
        gd0 = gd0 + ( 2.0D+00 * k - 1.0D+00 ) * bk(k) &
             * x ** ( 2.0D+00 * k - 2.0D+00 )
     else
        gd0 = gd0 + 2.0D+00 * k * bk(k+1) * x ** ( 2.0D+00 * k - 1.0D+00 )
     end if

     if ( abs ( ( gd0 - gw ) / gd0 ) < eps .and. 10 <= k ) then
        exit
     end if

     gw = gd0

  end do

  gd = gd1 + xm * gd0

  return
end subroutine gmn