Next: BURN3D.F Up: Listing of key Previous: AC3D.F

GAUSS.F

      real x(4000),w(4000),xi(4000),wi(4000)
      open(unit=62,file='gauss20.dat')
      data pi/3.14159 26535 89793 23846 26434/

c***********************************************************************
c   This program (gauss.f) evaluates zeros and weights of Legendre     *
c   polynomials to be used in Gaussian integrals.                      *
c   Pick N to be the order of the Gaussian quadrature that you         *
c   wish to use.  File output will contain Gaussian points in the      *
c   first column, weights in the second column.                        *
c   The limit on N is currently 4000.  If a larger value of N is       *
c   desired (not recommended), then change the dimensions of           *
c   x,w,xi, and wi.                                                    *
c   This program was supplied by Professor W.W. Repko, of the          *
c   Michigan State University Physics and Astronomy Department (1983). *
c***********************************************************************

c  Number N of Gaussian points desired.
      N=20
      M=(N+1)/2

c  Subroutine grule calculates weights and points

      call grule(N,x,w)

      do 10 i=1,M
      xi(i)=-x(i)
      xi(i+M)=x(M+1-i)
      wi(i)=w(i)
      wi(i+M)=w(M+1-i)
10    continue

c  output weights and points

      xc=0.0
      do 20 i=1,n
      write(62,3) xi(i),wi(i)
3     format(' ',2f20.12)

c  test with integral from -1 to 1 of exp(x), should be e - 1/e = 2.3504024

      xc=xc+wi(i)*exp(xi(i))
20    continue
      print *,' Numerical value of integral of exp(x) from -1 to 1 = '
      print *,xc
      print *,' Actual value is 2.3504024'
      
      end
      subroutine grule(N,x,w)
      real x(N),w(N)
      data pi/3.14159 26535 89793 23846 26434/
      data eps/1.e-14/
      dn=float(N)
      M=(N+1)/2
      e1=dn*(dn+1.0)
      do 10 i=1,M
      t=(4.0*float(i)-1.0)*pi/(4.0*dn+2.0)
      x0=(1.0-(1.0-1.0/dn)/(8.0*dn*dn))*cos(t)
20    call legendr(N,x0,pn,pnm1,pnp1)
      den=1.0-x0*x0
      d1=dn*(pnm1-x0*pn)
      dpn=d1/den
      d2pn=(2.0*x0*dpn-e1*pn)/den
      u=pn/dpn
      v=d2pn/dpn
      x1=x0-u*(1.0+0.50*u*v)
      if((x1-x0).lt.eps) go to 30
      x0=x1
      go to 20
30    x0=x1
      call legendr(N,x0,pn,pnm1,pnp1)
      x(i)=x0
10    w(i)=2.0*(1.0-x0**2)/(dn*pnm1)**2
      if(M+M.gt.N) x(M)=0.0
      return
      end
      subroutine legendr(N,x,pn,pnm1,pnp1)
      pkm1=1.0
      pk=x
      do 10 k=2,N
      t1=x*pk
      pkp1=t1-pkm1-(t1-pkm1)/float(k)+t1
      pkm1=pk
10    pk=pkp1
      pn=pk
      pnm1=pkm1
      t1=x*pn
      pnp1=t1-pnm1-(t1-pnm1)/float(k+1)+t1
      return
      end


 

 Next:  BURN3D.F Up:  Listing of key   Previous:  AC3D.F