```

SUBROUTINE QGAUSN( M, GMU, GWT ) 2,3

c       Compute weights and abscissae for ordinary Gaussian quadrature
c       on the interval (0,1);  that is, such that

c           sum(i=1 to M) ( GWT(i) f(GMU(i)) )

c       is a good approximation to

c           integral(0 to 1) ( f(x) dx )

c   INPUT :    M       order of quadrature rule

c   OUTPUT :  GMU(I)   array of abscissae (I = 1 TO M)
c             GWT(I)   array of weights (I = 1 TO M)

c   REFERENCE:  Davis, P.J. and P. Rabinowitz, Methods of Numerical
c                   Integration, Academic Press, New York, pp. 87, 1975

c   METHOD:  Compute the abscissae as roots of the Legendre
c            polynomial P-sub-M using a cubically convergent
c            refinement of Newton's method.  Compute the
c            weights from EQ. 2.7.3.8 of Davis/Rabinowitz.  Note
c            that Newton's method can very easily diverge; only a
c            very good initial guess can guarantee convergence.
c            The initial guess used here has never led to divergence
c            even for M up to 1000.

c   ACCURACY:  relative error no better than TOL or computer
c              precision (machine epsilon), whichever is larger

c   INTERNAL VARIABLES:

c    ITER      : number of Newton Method iterations
c    MAXIT     : maximum allowed iterations of Newton Method
c    PM2,PM1,P : 3 successive Legendre polynomials
c    PPR       : derivative of Legendre polynomial
c    P2PRI     : 2nd derivative of Legendre polynomial
c    TOL       : convergence criterion for Legendre poly root iteration
c    X,XI      : successive iterates in cubically-convergent version
c                of Newtons Method (seeking roots of Legendre poly.)

c   Called by- SETDIS, SURFAC
c   Calls- D1MACH, ERRMSG
c +-------------------------------------------------------------------+

c     .. Scalar Arguments ..

INTEGER   M
c     ..
c     .. Array Arguments ..

REAL      GMU( M ), GWT( M )
c     ..
c     .. Local Scalars ..

INTEGER   ITER, K, LIM, MAXIT, NN, NP1
REAL      CONA, PI, T
DOUBLE PRECISION EN, NNP1, ONE, P, P2PRI, PM1, PM2, PPR, PROD,
&                 TMP, TOL, TWO, X, XI
c     ..
c     .. External Functions ..

DOUBLE PRECISION D1MACH
EXTERNAL  D1MACH
c     ..
c     .. External Subroutines ..

EXTERNAL  ERRMSG
c     ..
c     .. Intrinsic Functions ..

INTRINSIC ABS, ASIN, COS, FLOAT, MOD, TAN
c     ..
SAVE      PI, TOL

DATA      PI / 0.0 / , MAXIT / 1000 / , ONE / 1.D0 / ,
&          TWO / 2.D0 /

IF( PI.EQ.0.0 ) THEN

PI   = 2.*ASIN( 1.0 )
TOL  = 10.*D1MACH( 4 )

END IF

IF( M.LT.1 ) CALL ERRMSG( 'QGAUSN--Bad value of M',.True.)

IF( M.EQ.1 ) THEN

GMU( 1 ) = 0.5
GWT( 1 ) = 1.0
RETURN

END IF

EN   = M
NP1  = M + 1
NNP1 = M*NP1
CONA = FLOAT( M - 1 ) / ( 8*M**3 )

LIM  = M / 2

DO 30 K = 1, LIM
c                                        ** Initial guess for k-th root
c                                        ** of Legendre polynomial, from
c                                        ** Davis/Rabinowitz (2.7.3.3a)
T  = ( 4*K - 1 )*PI / ( 4*M + 2 )
X  = COS( T + CONA / TAN( T ) )
ITER = 0
c                                        ** Upward recurrence for
c                                        ** Legendre polynomials
10    CONTINUE
ITER   = ITER + 1
PM2    = ONE
PM1    = X

DO 20 NN = 2, M
P    = (            PM2  = PM1
PM1  = P
20    CONTINUE
c                                              ** Newton Method
TMP    = ONE / ( ONE - X**2 )
PPR    = EN*( PM2 - X*P )*TMP
P2PRI  = ( TWO*X*PPR - NNP1*P )*TMP
XI     = X - ( P / PPR )*( ONE +
&            ( P / PPR )*P2PRI / ( TWO*PPR ) )

c                                              ** Check for convergence
IF( ABS( XI - X ).GT.TOL ) THEN

IF( ITER.GT.MAXIT )
&          CALL ERRMSG( 'QGAUSN--max iteration count',.True.)

X  = XI
GO TO  10

END IF
c                             ** Iteration finished--calculate weights,
c                             ** abscissae for (-1,1)
GMU( K ) = -X
GWT( K ) = TWO / ( TMP*( EN*PM2 )**2 )
GMU( NP1 - K ) = -GMU( K )
GWT( NP1 - K ) = GWT( K )
30 CONTINUE
c                                    ** Set middle abscissa and weight
c                                    ** for rules of odd order
IF( MOD( M,2 ).NE.0 ) THEN

GMU( LIM + 1 ) = 0.0
PROD   = ONE

DO 40 K = 3, M, 2
PROD   = PROD * K / ( K - 1 )
40    CONTINUE

GWT( LIM + 1 ) = TWO / PROD**2
END IF

c                                        ** Convert from (-1,1) to (0,1)
DO 50 K = 1, M
GMU( K ) = 0.5*GMU( K ) + 0.5
GWT( K ) = 0.5*GWT( K )
50 CONTINUE

RETURN
END
```