lqnb Subroutine

subroutine lqnb ( n, x, qn, qd )

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

  integer ( kind = 4 ) n

  real ( kind = 8 ) eps
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  integer ( kind = 4 ) nl
  real ( kind = 8 ) q0
  real ( kind = 8 ) q1
  real ( kind = 8 ) qc1
  real ( kind = 8 ) qc2
  real ( kind = 8 ) qd(0:n)
  real ( kind = 8 ) qf
  real ( kind = 8 ) qf0
  real ( kind = 8 ) qf1
  real ( kind = 8 ) qf2
  real ( kind = 8 ) qn(0:n)
  real ( kind = 8 ) qr
  real ( kind = 8 ) x
  real ( kind = 8 ) x2

  eps = 1.0D-14

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

  if ( x <= 1.021D+00 ) then

     x2 = abs ( ( 1.0D+00 + x ) / ( 1.0D+00 - x ) )
     q0 = 0.5D+00 * log ( x2 )
     q1 = x * q0 - 1.0D+00
     qn(0) = q0
     qn(1) = q1
     qd(0) = 1.0D+00 / ( 1.0D+00 - x * x )
     qd(1) = qn(0) + x * qd(0)
     do k = 2, n
        qf = ( ( 2.0D+00 * k - 1.0D+00 ) * x * q1 &
             - ( k - 1.0D+00 ) * q0 ) / k
        qn(k) = qf
        qd(k) = ( qn(k-1) - x * qf ) * k / ( 1.0D+00 - x * x )
        q0 = q1
        q1 = qf
     end do

  else

     qc2 = 1.0D+00 / x
     do j = 1, n
        qc2 = qc2 * j / ( ( 2.0D+00 * j + 1.0D+00 ) * x )
        if ( j == n - 1 ) then
           qc1 = qc2
        end if
     end do

     do l = 0, 1

        nl = n + l
        qf = 1.0D+00
        qr = 1.0D+00
        do k = 1, 500
           qr = qr * ( 0.5D+00 * nl + k - 1.0D+00 ) &
                * ( 0.5D+00 * ( nl - 1 ) + k ) &
                / ( ( nl + k - 0.5D+00 ) * k * x * x )
           qf = qf + qr
           if ( abs ( qr / qf ) < eps ) then
              exit
           end if
        end do

        if ( l == 0 ) then
           qn(n-1) = qf * qc1
        else
           qn(n) = qf * qc2
        end if

     end do

     qf2 = qn(n)
     qf1 = qn(n-1)
     do k = n, 2, -1
        qf0 = ( ( 2.0D+00 * k - 1.0D+00 ) * x * qf1 - k * qf2 ) / ( k - 1.0D+00 )
        qn(k-2) = qf0
        qf2 = qf1
        qf1 = qf0
     end do

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

  end if

  return
end subroutine lqnb