C     Last change:  GL    1 Sep 2004    7:37 pm
      SUBROUTINE QUAD (NZ, ANGLE, W)
C        Generates the abscissas and weights for an even 2*NZ point
C      Gauss-Legendre quadrature.  Only the NZ positive points are 
C  returned.
      INTEGER  NZ
      REAL     ANGLE(*), W(*)

      INTEGER  N, I, J, K, NUM
      PARAMETER (NUM=50)  
      REAL*8   ABSCISSAS(NUM), WEIGHTS(NUM)
      REAL*8   Z, Z1, P1, P2, P3, PP, EPS
      PARAMETER (EPS=3.0D-14)
 
      N=2*NZ
      DO 12 I=1,NZ
        Z=DCOS(3.141592654D0*(I-.25D0)/(N+.5D0))
        K=0
1       CONTINUE
          P1=1.D0
          P2=0.D0
          DO 11 J=1,N
            P3=P2
            P2=P1
            P1=((2.D0*J-1.D0)*Z*P2-(J-1.D0)*P3)/J
11        CONTINUE
          PP=N*(Z*P1-P2)/(Z*Z-1.D0)
          Z1=Z
          Z=Z1-P1/PP
          K=K+1
        IF (ABS(Z-Z1).GT.EPS .AND. K.LT.10) GO TO 1
        ABSCISSAS(NZ+1-I)=Z
        WEIGHTS(NZ+1-I)=2.D0/((1.D0-Z*Z)*PP*PP)
12    CONTINUE
 
      DO 15 I=1,NZ
        ANGLE(I)=SNGL(ABSCISSAS(I))
        W(I)=SNGL(WEIGHTS(I))
 15   CONTINUE

      RETURN
      END