lqmn Subroutine

subroutine lqmn ( mm, m, n, x, qm, qd )

  !*****************************************************************************80
  !
  !! LQMN computes associated Legendre functions Qmn(x) and derivatives.
  !
  !  Discussion:
  !
  !    This routine computes the associated Legendre functions of the
  !    second kind, Qmn(x) and 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:
  !
  !    13 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 ) MM, determines the leading dimension 
  !    of QM and QD.
  !
  !    Input, integer ( kind = 4 ) M, the order of Qmn(x).
  !
  !    Input, integer ( kind = 4 ) N, the degree of Qmn(x).
  !
  !    Output, real ( kind = 8 ) QM(0:MM,0:N), QD(0:MM,0:N), contains the values
  !    of Qmn(x) and Qmn'(x).
  !
  implicit none

  integer ( kind = 4 ) mm
  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) km
  integer ( kind = 4 ) ls
  integer ( kind = 4 ) m
  real ( kind = 8 ) q0
  real ( kind = 8 ) q1
  real ( kind = 8 ) q10
  real ( kind = 8 ) qd(0:mm,0:n)
  real ( kind = 8 ) qf
  real ( kind = 8 ) qf0
  real ( kind = 8 ) qf1
  real ( kind = 8 ) qf2
  real ( kind = 8 ) qm(0:mm,0:n)
  real ( kind = 8 ) x
  real ( kind = 8 ) xq
  real ( kind = 8 ) xs

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

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

  xs = ls * ( 1.0D+00 - x * x )
  xq = sqrt ( xs )
  q0 = 0.5D+00 * log ( abs ( ( x + 1.0D+00 ) / ( x - 1.0D+00 ) ) )

  if ( abs ( x ) < 1.0001D+00 ) then
     qm(0,0) = q0
     qm(0,1) = x * q0 - 1.0D+00
     qm(1,0) = -1.0D+00 / xq
     qm(1,1) = -xq * ( q0 + x / ( 1.0D+00 - x * x ) )
     do i = 0, 1
        do j = 2, n
           qm(i,j) = ( ( 2.0D+00 * j - 1.0D+00 ) * x * qm(i,j-1) &
                - ( j + i - 1.0D+00 ) * qm(i,j-2))/ ( j - i )
        end do
     end do

     do j = 0, n
        do i = 2, m
           qm(i,j) = -2.0D+00 * ( i - 1.0D+00 ) * x / xq * qm(i-1,j) &
                - ls * ( j + i - 1.0D+00 ) * ( j - i + 2.0D+00 ) * qm(i-2,j)
        end do
     end do

  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 * k + 3.0D+00 ) * x * qf1 &
             - ( k + 2.0D+00 ) * qf2 ) / ( k + 1.0D+00 )
        if ( k <= n ) then
           qm(0,k) = qf0
        end if
        qf2 = qf1
        qf1 = qf0
     end do

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

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

     q10 = -1.0D+00 / xq
     do k = 0, n
        qm(1,k) = q10 * qm(1,k) / qf0
     end do

     do j = 0, n
        q0 = qm(0,j)
        q1 = qm(1,j)
        do i = 0, m - 2
           qf = -2.0D+00 * ( i + 1 ) * x / xq * q1 &
                + ( j - i ) * ( j + i + 1.0D+00 ) * q0
           qm(i+2,j) = qf
           q0 = q1
           q1 = qf
        end do
     end do

  end if

  qd(0,0) = ls / xs
  do j = 1, n
     qd(0,j) = ls * j * ( qm(0,j-1) - x * qm(0,j) ) / xs
  end do

  do j = 0, n
     do i = 1, m
        qd(i,j) = ls * i * x / xs * qm(i,j) &
             + ( i + j ) * ( j - i + 1.0D+00 ) / xq * qm(i-1,j)
     end do
  end do

  return
end subroutine lqmn