PROGRAM MAIN
C This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/gauss.for
C
      DOUBLE PRECISION X1,X2,X(100),W(100),EPS
      INTEGER N,I
C
C     MIT FREUNDLICHEM GRUSS AN HERMANN VON ULRICH
C                               23. SEPTEMBER 1987
      WRITE(*,1)
1     FORMAT(' PLEASE SPECIFY THE LOWER AND THEN UPPER BOUND FOR THE ',/
     *       ' GAUSS-LEGENDRE QUADRATURE. ')
      READ(*,*) X1,X2
      WRITE(*,2)
2     FORMAT(' PLEASE SPECIFY NOW THE NUMBER OF WEIGHTS, ABSCISSAE ',/
     *       ' VALUES RESPECTIVELY YOU NEED (CURRENT MAXIMUM 200) ')
      READ(*,*) N
      WRITE(*,3)
3     FORMAT(' PLEASE SPECIFY NOW THE A LEVEL OF PRECISION, 3.0D-17 ',/
     *       ' WILL GIVE YOU ROUGHLY 16 DIGITS, I.E. DOUBLE PRECISION.')
      READ(*,*)
      CALL GAULEG(X1,X2,X,W,N,EPS)
      WRITE(*,4)
4     FORMAT(' HALF THE WEIGHTS AND ABSCISSAE STARTING WITH THE LAST ',/
     *       ' WEIGHT(I)=WEIGHT(N-I);-ABSCISSAE(I)=ABSCISSAE(N-I);',/,
     *       ' I=0,1,...N/2 IS THE RECURRENCE RELATIONSHIP. ',/)
      DO 10 I = 1,(N/2)
      WRITE(*,5) X(N/2-I+1),W(N/2-I+1)
5     FORMAT(2F20.16)
10    CONTINUE
      STOP
      END
C
      SUBROUTINE GAULEG(X1,X2,X,W,N,EPS)
C
C     GIVEN THE LOWER AND UPPER LIMITS OF INTEGRATION X1 AND X2,AND
C     GIVEN N, THIS SUBROUTINE RETURNS ARRAYS X AND W OF LENGTH N,
C     CONTAINING THE ABSCISSAE AND WEIGHTS OF THE GAUSS-LEGENDRE N-POINT
C     QUADRATURE FORMULA. THE ALGORITHM IS GIVEN IN "NUMERICAL RECIPES"
C     PAGES 125-126 AND IS REWRITTEN HERE INTO DOUBLE PRECISION.
C
      DOUBLE PRECISION X1,X2,X(N),W(N),XM,XL,Z,P1,P2,P3,PP,Z1,PI,EPS
      INTEGER N,I,J,M
      PI = 2.0D0*DASIN(1.0D0)
      M = (N+1)/2
      XM = 0.5D0*(X2+X1)
      XL = 0.5D0*(X2-X1)
      DO 12 I = 1,M
      Z = DCOS(PI*(I-0.25D0)/(N+0.5D0))
C     STARTING WITH THE ABOVE APPROXIMATION TO THE I-TH ROOT, WE ENTER
C     THE MAIN LOOP OF REFINEMENT BY NEWTON'S METHOD.
1     CONTINUE
      P1 = 1.0D0
      P2 = 0.0D0
      DO 11 J = 1,N
C     LOOP UP THE RECURRENCE RELATION TO GET THE LEGENDRE POLYNOMIAL
C     EVALUATED AT Z.
      P3 = P2
      P2 = P1
      P1 = ((2.0D0*J-1.0D0)*Z*P2-(J-1.0D0)*P3)/J
11    CONTINUE
C     P1 IS NOW THE DESIRED LEGENDRE POLYNOMIAL. WE NEXT COMPUTE PP,
C     ITS DERIVATIVE, BY A STANDARD RELATION INVOLVING ALSO P2, THE
C     POLYNOMIAL OF ONE LOWER ORDER.
      PP = N*(Z*P1-P2)/(Z*Z-1.0D0)
      Z1 = Z
      Z = Z1-P1/PP
      IF (DABS(Z-Z1).GT.EPS) GO TO 1
      X((N/2)+1-I) = XM+XL*Z
C      X(N+1-I)=XM+XL*Z     TAKE ONLY HALF THE POINTS AND WEIGHTS
      W((N/2)+1-I) = 2.0D0*XL/((1.0D0-Z*Z)*PP*PP)
C      W(N+1-I) = W(I)
12    CONTINUE
      RETURN
      END