lpmv Subroutine

subroutine lpmv ( v, m, x, pmv )

  !*****************************************************************************80
  !
  !! LPMV computes associated Legendre functions Pmv(X) with arbitrary degree.
  !
  !  Discussion:
  !
  !    Compute the associated Legendre function Pmv(x) with an integer order 
  !    and an arbitrary nonnegative degree v.
  !
  !  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:
  !
  !    19 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, real ( kind = 8 ) V, the degree of Pmv(x).
  !
  !    Input, integer ( kind = 4 ) M, the order of Pmv(x).
  !
  !    Input, real ( kind = 8 ) X, the argument of Pm(x).
  !
  !    Output, real ( kind = 8 ) PMV, the value of Pm(x).
  !
  implicit none

  real ( kind = 8 ) c0
  real ( kind = 8 ) el
  real ( kind = 8 ) eps
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) nv
  real ( kind = 8 ) pa
  real ( kind = 8 ) pi
  real ( kind = 8 ) pmv
  real ( kind = 8 ) pss
  real ( kind = 8 ) psv
  real ( kind = 8 ) pv0
  real ( kind = 8 ) qr
  real ( kind = 8 ) r
  real ( kind = 8 ) r0
  real ( kind = 8 ) r1
  real ( kind = 8 ) r2
  real ( kind = 8 ) rg
  real ( kind = 8 ) s
  real ( kind = 8 ) s0
  real ( kind = 8 ) s1
  real ( kind = 8 ) s2
  real ( kind = 8 ) v
  real ( kind = 8 ) v0
  real ( kind = 8 ) vs
  real ( kind = 8 ) x
  real ( kind = 8 ) xq

  pi = 3.141592653589793D+00
  el = 0.5772156649015329D+00
  eps = 1.0D-14
  nv = int ( v )
  v0 = v - nv

  if ( x == -1.0D+00 .and. v /= nv ) then
     if ( m == 0 ) then
        pmv = -1.0D+300
     else
        pmv = 1.0D+300
     end if
     return
  end if

  c0 = 1.0D+00

  if ( m /= 0 ) then

     rg = v * ( v + m )
     do j = 1, m - 1 
        rg = rg * ( v * v - j * j )
     end do
     xq = sqrt ( 1.0D+00 - x * x )
     r0 = 1.0D+00
     do j = 1, m
        r0 = 0.5D+00 * r0 * xq / j
     end do
     c0 = r0 * rg

  end if

  if ( v0 == 0.0D+00 ) then

     pmv = 1.0D+00
     r = 1.0D+00
     do k = 1, nv - m
        r = 0.5D+00 * r * ( - nv + m + k - 1.0D+00 ) &
             * ( nv + m + k ) / ( k * ( k + m ) ) * ( 1.0D+00 + x )
        pmv = pmv + r
     end do
     pmv = ( -1.0D+00 ) ** nv * c0 * pmv

  else

     if ( -0.35D+00 <= x ) then

        pmv = 1.0D+00
        r = 1.0D+00
        do k = 1, 100
           r = 0.5D+00 * r * ( - v + m + k - 1.0D+00 ) &
                * ( v + m + k ) / ( k * ( m + k ) ) * ( 1.0D+00 - x )
           pmv = pmv + r
           if ( 12 < k .and. abs ( r / pmv ) < eps ) then
              exit
           end if
        end do

        pmv = ( -1.0D+00 ) ** m * c0 * pmv

     else

        vs = sin ( v * pi ) / pi
        pv0 = 0.0D+00

        if ( m /= 0 ) then

           qr = sqrt ( ( 1.0D+00 - x ) / ( 1.0D+00 + x ) )
           r2 = 1.0D+00
           do j = 1, m
              r2 = r2 * qr * j
           end do
           s0 = 1.0D+00
           r1 = 1.0D+00
           do k = 1, m - 1 
              r1 = 0.5D+00 * r1 * ( - v + k - 1 ) * ( v + k ) &
                   / ( k * ( k - m ) ) * ( 1.0D+00 + x )
              s0 = s0 + r1
           end do
           pv0 = - vs * r2 / m * s0

        end if

        call psi ( v, psv )
        pa = 2.0D+00 * ( psv + el ) + pi / tan ( pi * v ) &
             + 1.0D+00 / v

        s1 = 0.0D+00
        do j = 1, m
           s1 = s1 + ( j * j + v * v ) / ( j * ( j * j - v * v ) )
        end do

        pmv = pa + s1 - 1.0D+00 / ( m - v ) &
             + log ( 0.5D+00 * ( 1.0D+00 + x ) )
        r = 1.0D+00
        do k = 1, 100
           r = 0.5D+00 * r * ( - v + m + k - 1.0D+00 ) * ( v + m + k ) &
                / ( k * ( k + m ) ) * ( 1.0D+00 + x )
           s = 0.0D+00
           do j = 1, m
              s = s + ( ( k + j ) ** 2 + v * v ) &
                   / ( ( k + j ) * ( ( k + j ) ** 2 - v * v ) )
           end do
           s2 = 0.0D+00
           do j = 1, k
              s2 = s2 + 1.0D+00 / ( j * ( j * j - v * v ) )
           end do
           pss = pa + s + 2.0D+00 * v * v * s2 &
                - 1.0D+00 / ( m + k - v ) &
                + log ( 0.5D+00 * ( 1.0D+00 + x ) )
           r2 = pss * r
           pmv = pmv + r2
           if ( abs ( r2 / pmv ) < eps ) then
              exit
           end if
        end do

        pmv = pv0 + pmv * vs * c0

     end if

  end if

  return
end subroutine lpmv