lqnb Subroutine

   subroutine lqnb(n,x,Qn,Qd)

      real(wp),intent(in) :: x !! Argument of Qn(x)
      integer,intent(in) :: n !! Degree of Qn(x)  ( n = 0,1,2,…)
      real(wp),intent(out) :: Qn(0:n) !! Qn(x)
      real(wp),intent(out) :: Qd(0:n) !! Qn'(x)

      real(wp) :: q0 , q1 , qc1 , qc2 , qf , qf0 , qf1 , qf2 , qr , x2
      integer :: j , k , l , nl

      real(wp),parameter :: eps = 1.0e-14_wp

      if ( abs(x)==1.0_wp ) then
         do k = 0 , n
            Qn(k) = 1.0e+300_wp
            Qd(k) = 1.0e+300_wp
         enddo
         return
      endif
      if ( x<=1.021_wp ) then
         x2 = abs((1.0_wp+x)/(1.0_wp-x))
         q0 = 0.5_wp*log(x2)
         q1 = x*q0 - 1.0_wp
         Qn(0) = q0
         Qn(1) = q1
         Qd(0) = 1.0_wp/(1.0_wp-x*x)
         Qd(1) = Qn(0) + x*Qd(0)
         do k = 2 , n
            qf = ((2.0_wp*k-1.0_wp)*x*q1-(k-1.0_wp)*q0)/k
            Qn(k) = qf
            Qd(k) = (Qn(k-1)-x*qf)*k/(1.0_wp-x*x)
            q0 = q1
            q1 = qf
         enddo
      else
         qc1 = 0.0_wp
         qc2 = 1.0_wp/x
         do j = 1 , n
            qc2 = qc2*j/((2.0_wp*j+1.0_wp)*x)
            if ( j==n-1 ) qc1 = qc2
         enddo
         do l = 0 , 1
            nl = n + l
            qf = 1.0_wp
            qr = 1.0_wp
            do k = 1 , 500
               qr = qr*(0.5_wp*nl+k-1.0_wp)*(0.5_wp*(nl-1)+k) &
                    /((nl+k-0.5_wp)*k*x*x)
               qf = qf + qr
               if ( abs(qr/qf)<eps ) exit
            enddo
            if ( l==0 ) then
               Qn(n-1) = qf*qc1
            else
               Qn(n) = qf*qc2
            endif
         enddo
         qf2 = Qn(n)
         qf1 = Qn(n-1)
         do k = n , 2 , -1
            qf0 = ((2*k-1.0_wp)*x*qf1-k*qf2)/(k-1.0_wp)
            Qn(k-2) = qf0
            qf2 = qf1
            qf1 = qf0
         enddo
         Qd(0) = 1.0_wp/(1.0_wp-x*x)
         do k = 1 , n
            Qd(k) = k*(Qn(k-1)-x*Qn(k))/(1.0_wp-x*x)
         enddo
      endif

   end subroutine lqnb