1 | include(dom.inc)
   2 | 
   3 |       SUBROUTINE gauleg(x1,x2,x,w,n)
   4 | 
   5 | !       ================================================================!
   6 | !       See Gauss-Legendre abscisa and weight calculation for gaussian  !
   7 | !       integration in Numerical Recipes                                !
   8 | !       ================================================================!
   9 | 
  10 |         IMPLICIT NONE
  11 | 
  12 | !       PARAMETERS
  13 |         DOM_REAL, PARAMETER :: EPS=3.e-14
  14 | 
  15 | !       IN
  16 |         DOM_INT n
  17 |         DOM_REAL :: w(n)
  18 | 
  19 | !       OUT
  20 |         DOM_REAL :: x1,x2,x(n)
  21 | 
  22 | !       LOCAL
  23 |         DOM_INT  :: i,j,m
  24 |         DOM_REAL :: p1,p2,p3,pp,xl,xm,z,z1
  25 | 
  26 |         m = (n+1) / 2
  27 |         xm = 0.5 * (x2+x1)
  28 |         xl = 0.5 * (x2-x1)
  29 | 
  30 |         DO i=1,m
  31 |            z=cos(3.141592654*(i-.25)/(n+.5))
  32 | 1          continue
  33 |            p1=1.
  34 |            p2=0.
  35 |            DO j=1,n
  36 |               p3=p2
  37 |               p2=p1
  38 |               p1=((2.*j-1.)*z*p2-(j-1.)*p3)/j
  39 |            ENDDO
  40 |            pp=n*(z*p1-p2)/(z*z-1.)
  41 |            z1=z
  42 |            
  43 |            z=z1-p1/pp
  44 |            if(abs(z-z1).gt.EPS)goto 1
  45 |            x(i)=xm-xl*z
  46 |            x(n+1-i)=xm+xl*z
  47 |            w(i)=2.*xl/((1.-z*z)*pp*pp)
  48 |            w(n+1-i)=w(i)
  49 |         ENDDO
  50 | 
  51 |       END SUBROUTINE gauleg