1 | include(dom.inc)
   2 | 
   3 |        SUBROUTINE gauleg(x1,x2,x,w,n)
   4 | 
   5 |        IMPLICIT NONE
   6 | 
   7 | !      IN
   8 |        DOM_INT   :: n
   9 |        DOM_REAL  :: x1,x2
  10 | 
  11 | !      OUT
  12 |        DOM_REAL  :: x(n),w(n)
  13 | 
  14 | !      LOCAL
  15 |        DOM_INT             :: i,j,m
  16 |        DOM_REAL, PARAMETER :: EPS=3.e-14
  17 |        DOM_REAL            :: p1,p2,p3,pp,xl,xm,z,z1
  18 | 
  19 |        m=(n+1)/2
  20 |        xm=0.5*(x2+x1)
  21 |        xl=0.5*(x2-x1)
  22 | 
  23 |        DO i=1,m
  24 |           z=cos(3.141592654*(i-.25)/(n+.5))
  25 | 1         continue
  26 |           p1=1.
  27 |           p2=0.
  28 |           DO j=1,n
  29 |              p3=p2
  30 |              p2=p1
  31 |              p1=((2.*j-1.)*z*p2-(j-1.)*p3)/j
  32 |           ENDDO
  33 |           pp=n*(z*p1-p2)/(z*z-1.)
  34 |           z1=z
  35 |            
  36 |           z=z1-p1/pp
  37 |           if(abs(z-z1).gt.EPS)goto 1
  38 |           x(i)=xm-xl*z
  39 |           x(n+1-i)=xm+xl*z
  40 |           w(i)=2.*xl/((1.-z*z)*pp*pp)
  41 |           w(n+1-i)=w(i)
  42 |        ENDDO
  43 | 
  44 | !        x(1)=0.00000
  45 | !        w(1)=0.04500
  46 | !        x(2)=0.15541
  47 | !        w(2)=0.24500
  48 | !        x(3)=0.45000
  49 | !        w(3)=0.32000
  50 | !        x(4)=0.74459
  51 | !        w(4)=0.24500
  52 | !        x(5)=0.90000
  53 | !        w(5)=0.05611
  54 | !        x(6)=0.93551
  55 | !        w(6)=0.05125
  56 | !        x(7)=0.98449
  57 | !        w(7)=0.03764
  58 | 
  59 |        END SUBROUTINE gauleg