lqmns Subroutine

subroutine lqmns ( m, n, x, qm, qd )

  !*****************************************************************************80
  !
  !! LQMNS computes associated Legendre functions Qmn(x) and derivatives Qmn'(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:
  !
  !    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 order.
  !
  !    Input, integer ( kind = 4 ) N, the degree.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) QM(0:N), QD(0:N), the values of Qmn(x) 
  !    and Qmn'(x).
  !
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) k
  integer ( kind = 4 ) km
  integer ( kind = 4 ) l
  integer ( kind = 4 ) ls
  integer ( kind = 4 ) m
  real ( kind = 8 ) q0
  real ( kind = 8 ) q00
  real ( kind = 8 ) q01
  real ( kind = 8 ) q0l
  real ( kind = 8 ) q10
  real ( kind = 8 ) q11
  real ( kind = 8 ) q1l
  real ( kind = 8 ) qd(0:n)
  real ( kind = 8 ) qf0
  real ( kind = 8 ) qf1
  real ( kind = 8 ) qf2
  real ( kind = 8 ) qg0
  real ( kind = 8 ) qg1
  real ( kind = 8 ) qh0
  real ( kind = 8 ) qh1
  real ( kind = 8 ) qh2
  real ( kind = 8 ) qm(0:n)
  real ( kind = 8 ) qm0
  real ( kind = 8 ) qm1
  real ( kind = 8 ) qmk
  real ( kind = 8 ) x
  real ( kind = 8 ) xq

  do k = 0, n
     qm(k) = 0.0D+00
     qd(k) = 0.0D+00
  end do

  if ( abs ( x ) == 1.0D+00 ) then
     do k = 0, n
        qm(k) = 1.0D+300
        qd(k) = 1.0D+300
     end do
     return
  end if

  if ( 1.0D+00 < abs ( x ) ) then
     ls = -1
  else
     ls = +1
  end if

  xq = sqrt ( ls * ( 1.0D+00 - x * x ) )
  q0 = 0.5D+00 * log ( abs ( ( x + 1.0D+00 ) / ( x - 1.0D+00 ) ) )
  q00 = q0
  q10 = -1.0D+00 / xq
  q01 = x * q0 - 1.0D+00
  q11 = - ls * xq * ( q0 + x / ( 1.0D+00 - x * x ) )
  qf0 = q00
  qf1 = q10
  do k = 2, m
     qm0 = -2.0D+00 * ( k - 1.0D+00 ) / xq * x * qf1 &
          - ls * ( k - 1.0D+00 ) * ( 2.0D+00 - k ) * qf0
     qf0 = qf1
     qf1 = qm0
  end do

  if ( m == 0 ) then
     qm0 = q00
  else if ( m == 1 ) then
     qm0 = q10
  end if

  qm(0) = qm0

  if ( abs ( x ) < 1.0001D+00 ) then

     if ( m == 0 .and. 0 < n ) then

        qf0 = q00
        qf1 = q01
        do k = 2, n
           qf2 = ( ( 2.0D+00 * k - 1.0D+00 ) * x * qf1 &
                - ( k - 1.0D+00 ) * qf0 ) / k
           qm(k) = qf2
           qf0 = qf1
           qf1 = qf2
        end do

     end if
     qg0 = q01
     qg1 = q11
     do k = 2, m
        qm1 = - 2.0D+00 * ( k - 1.0D+00 ) / xq * x * qg1 &
             - ls * k * ( 3.0D+00 - k ) * qg0
        qg0 = qg1
        qg1 = qm1
     end do

     if ( m == 0 ) then
        qm1 = q01
     else if ( m == 1 ) then
        qm1 = q11
     end if
     qm(1) = qm1

     if ( m == 1 .and. 1 < n ) then

        qh0 = q10
        qh1 = q11
        do k = 2, n
           qh2 = ( ( 2.0D+00 * k - 1.0D+00 ) * x * qh1 - k * qh0 ) &
                / ( k - 1.0D+00 )
           qm(k) = qh2
           qh0 = qh1
           qh1 = qh2
        end do

     else if ( 2 <= m ) then

        qg0 = q00
        qg1 = q01
        qh0 = q10
        qh1 = q11

        do l = 2, n
           q0l = ( ( 2.0D+00 * l - 1.0D+00 ) * x * qg1 &
                - ( l - 1.0D+00 ) * qg0 ) / l
           q1l = ( ( 2.0D+00 * l - 1.0D+00 ) * x * qh1 - l * qh0 ) &
                / ( l - 1.0D+00 )
           qf0 = q0l
           qf1 = q1l
           do k = 2, m
              qmk = - 2.0D+00 * ( k - 1.0D+00 ) / xq * x * qf1 &
                   - ls * ( k + l - 1.0D+00 ) * ( l + 2.0D+00 - k ) * qf0
              qf0 = qf1
              qf1 = qmk
           end do
           qm(l) = qmk
           qg0 = qg1
           qg1 = q0l
           qh0 = qh1
           qh1 = q1l
        end do

     end if

  else

     if ( 1.1D+00 < abs ( x ) ) then
        km = 40 + m + n
     else
        km = ( 40 + m + n ) * int ( - 1.0D+00 - 1.8D+00 * log ( x - 1.0D+00 ) )
     end if

     qf2 = 0.0D+00
     qf1 = 1.0D+00
     do k = km, 0, -1
        qf0 = ( ( 2.0D+00 * k + 3.0D+00 ) * x * qf1 &
             - ( k + 2.0D+00 - m ) * qf2 ) / ( k + m + 1.0D+00 )
        if ( k <= n ) then
           qm(k) = qf0
        end if
        qf2 = qf1
        qf1 = qf0
     end do

     do k = 0, n
        qm(k) = qm(k) * qm0 / qf0
     end do

  end if

  if ( abs ( x ) < 1.0D+00 ) then
     do k = 0, n
        qm(k) = ( -1 ) ** m * qm(k)
     end do
  end if

  qd(0) = ( ( 1.0D+00 - m ) * qm(1) - x * qm(0) )  / ( x * x - 1.0D+00 )
  do k = 1, n
     qd(k) = ( k * x * qm(k) - ( k + m ) * qm(k-1) ) / ( x * x - 1.0D+00 )
  end do

  return
end subroutine lqmns