*DECK XPQNU
      SUBROUTINE XPQNU (NU1, NU2, MU, THETA, ID, PQA, IPQA, IERROR)
C***BEGIN PROLOGUE  XPQNU
C***SUBSIDIARY
C***PURPOSE  To compute the values of Legendre functions for XLEGF.
C            This subroutine calculates initial values of P or Q using
C            power series, then performs forward nu-wise recurrence to
C            obtain P(-MU,NU,X), Q(0,NU,X), or Q(1,NU,X). The nu-wise
C            recurrence is stable for P for all mu and for Q for mu=0,1.
C***LIBRARY   SLATEC
C***CATEGORY  C3A2, C9
C***TYPE      SINGLE PRECISION (XPQNU-S, DXPQNU-D)
C***KEYWORDS  LEGENDRE FUNCTIONS
C***AUTHOR  Smith, John M., (NBS and George Mason University)
C***ROUTINES CALLED  XADD, XADJ, XPSI
C***COMMON BLOCKS    XBLK1
C***REVISION HISTORY  (YYMMDD)
C   820728  DATE WRITTEN
C   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
C   901019  Revisions to prologue.  (DWL and WRB)
C   901106  Changed all specific intrinsics to generic.  (WRB)
C           Corrected order of sections in prologue and added TYPE
C           section.  (WRB)
C   920127  Revised PURPOSE section of prologue.  (DWL)
C***END PROLOGUE  XPQNU
      REAL A,NU,NU1,NU2,PQ,PQA,XPSI,R,THETA,W,X,X1,X2,XS,
     1 Y,Z
      REAL DI,DMU,PQ1,PQ2,FACTMU,FLOK
      DIMENSION PQA(*),IPQA(*)
      COMMON /XBLK1/ NBITSF
      SAVE /XBLK1/
C
C        J0, IPSIK, AND IPSIX ARE INITIALIZED IN THIS SUBROUTINE.
C        J0 IS THE NUMBER OF TERMS USED IN SERIES EXPANSION
C        IN SUBROUTINE XPQNU.
C        IPSIK, IPSIX ARE VALUES OF K AND X RESPECTIVELY
C        USED IN THE CALCULATION OF THE XPSI FUNCTION.
C
C***FIRST EXECUTABLE STATEMENT  XPQNU
      IERROR=0
      J0=NBITSF
      IPSIK=1+(NBITSF/10)
      IPSIX=5*IPSIK
      IPQ=0
C        FIND NU IN INTERVAL [-.5,.5) IF ID=2  ( CALCULATION OF Q )
      NU=MOD(NU1,1.)
      IF(NU.GE..5) NU=NU-1.
C        FIND NU IN INTERVAL (-1.5,-.5] IF ID=1,3, OR 4  ( CALC. OF P )
      IF(ID.NE.2.AND.NU.GT.-.5) NU=NU-1.
C        CALCULATE MU FACTORIAL
      K=MU
      DMU=MU
      IF(MU.LE.0) GO TO 60
      FACTMU=1.
      IF=0
      DO 50 I=1,K
      FACTMU=FACTMU*I
   50 CALL XADJ(FACTMU,IF,IERROR)
      IF (IERROR.NE.0) RETURN
   60 IF(K.EQ.0) FACTMU=1.
      IF(K.EQ.0) IF=0
C
C        X=COS(THETA)
C        Y=SIN(THETA/2)**2=(1-X)/2=.5-.5*X
C        R=TAN(THETA/2)=SQRT((1-X)/(1+X)
C
      X=COS(THETA)
      Y=SIN(THETA/2.)**2
      R=TAN(THETA/2.)
C
C        USE ASCENDING SERIES TO CALCULATE TWO VALUES OF P OR Q
C        FOR USE AS STARTING VALUES IN RECURRENCE RELATION.
C
      PQ2=0.0
      DO 100 J=1,2
      IPQ1=0
      IF(ID.EQ.2) GO TO 80
C
C        SERIES FOR P ( ID = 1, 3, OR 4 )
C        P(-MU,NU,X)=1./FACTORIAL(MU)*SQRT(((1.-X)/(1.+X))**MU)
C                *SUM(FROM 0 TO J0-1)A(J)*(.5-.5*X)**J
C
      IPQ=0
      PQ=1.
      A=1.
      IA=0
      DO 65 I=2,J0
      DI=I
      A=A*Y*(DI-2.-NU)*(DI-1.+NU)/((DI-1.+DMU)*(DI-1.))
      CALL XADJ(A,IA,IERROR)
      IF (IERROR.NE.0) RETURN
      IF(A.EQ.0.) GO TO 66
      CALL XADD(PQ,IPQ,A,IA,PQ,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
   65 CONTINUE
   66 CONTINUE
      IF(MU.LE.0) GO TO 90
      X2=R
      X1=PQ
      K=MU
      DO 77 I=1,K
      X1=X1*X2
   77 CALL XADJ(X1,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
      PQ=X1/FACTMU
      IPQ=IPQ-IF
      CALL XADJ(PQ,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
      GO TO 90
C
C        Z=-LN(R)=.5*LN((1+X)/(1-X))
C
   80 Z=-LOG(R)
      W=XPSI(NU+1.,IPSIK,IPSIX)
      XS=1./SIN(THETA)
C
C        SERIES SUMMATION FOR Q ( ID = 2 )
C        Q(0,NU,X)=SUM(FROM 0 TO J0-1)((.5*LN((1+X)/(1-X))
C    +XPSI(J+1,IPSIK,IPSIX)-XPSI(NU+1,IPSIK,IPSIX)))*A(J)*(.5-.5*X)**J
C
C        Q(1,NU,X)=-SQRT(1./(1.-X**2))+SQRT((1-X)/(1+X))
C             *SUM(FROM 0 T0 J0-1)(-NU*(NU+1)/2*LN((1+X)/(1-X))
C                 +(J-NU)*(J+NU+1)/(2*(J+1))+NU*(NU+1)*
C     (XPSI(NU+1,IPSIK,IPSIX)-XPSI(J+1,IPSIK,IPSIX))*A(J)*(.5-.5*X)**J
C
C        NOTE, IN THIS LOOP K=J+1
C
      PQ=0.
      IPQ=0
      IA=0
      A=1.
      DO 85 K=1,J0
      FLOK=K
      IF(K.EQ.1) GO TO 81
      A=A*Y*(FLOK-2.-NU)*(FLOK-1.+NU)/((FLOK-1.+DMU)*(FLOK-1.))
      CALL XADJ(A,IA,IERROR)
      IF (IERROR.NE.0) RETURN
   81 CONTINUE
      IF(MU.GE.1) GO TO 83
      X1=(XPSI(FLOK,IPSIK,IPSIX)-W+Z)*A
      IX1=IA
      CALL XADD(PQ,IPQ,X1,IX1,PQ,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
      GO TO 85
   83 X1=(NU*(NU+1.)*(Z-W+XPSI(FLOK,IPSIK,IPSIX))+(NU-FLOK+1.)
     1  *(NU+FLOK)/(2.*FLOK))*A
      IX1=IA
      CALL XADD(PQ,IPQ,X1,IX1,PQ,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
   85 CONTINUE
      IF(MU.GE.1) PQ=-R*PQ
      IXS=0
      IF(MU.GE.1) CALL XADD(PQ,IPQ,-XS,IXS,PQ,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
      IF(J.EQ.2) MU=-MU
      IF(J.EQ.2) DMU=-DMU
   90 IF(J.EQ.1) PQ2=PQ
      IF(J.EQ.1) IPQ2=IPQ
      NU=NU+1.
  100 CONTINUE
      K=0
      IF(NU-1.5.LT.NU1) GO TO 120
      K=K+1
      PQA(K)=PQ2
      IPQA(K)=IPQ2
      IF(NU.GT.NU2+.5) RETURN
  120 PQ1=PQ
      IPQ1=IPQ
      IF(NU.LT.NU1+.5) GO TO 130
      K=K+1
      PQA(K)=PQ
      IPQA(K)=IPQ
      IF(NU.GT.NU2+.5) RETURN
C
C        FORWARD NU-WISE RECURRENCE FOR F(MU,NU,X) FOR FIXED MU
C        USING
C        (NU+MU+1)*F(MU,NU,X)=(2.*NU+1)*F(MU,NU,X)-(NU-MU)*F(MU,NU-1,X)
C        WHERE F(MU,NU,X) MAY BE P(-MU,NU,X) OR IF MU IS REPLACED
C        BY -MU THEN F(MU,NU,X) MAY BE Q(MU,NU,X).
C        NOTE, IN THIS LOOP, NU=NU+1
C
  130 X1=(2.*NU-1.)/(NU+DMU)*X*PQ1
      X2=(NU-1.-DMU)/(NU+DMU)*PQ2
      CALL XADD(X1,IPQ1,-X2,IPQ2,PQ,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
      CALL XADJ(PQ,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
      NU=NU+1.
      PQ2=PQ1
      IPQ2=IPQ1
      GO TO 120
C
      END