CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      program quadmaster
C
C  master routine to call, compare quadrature algorithms
C
C  CWJ -SDSU - Aug 2006
C  This version compare trapezoidal and Simpson's, as a function of dx
C
C SUBROUTINES CALLED:
C   trapezoid (in quadraturesubs.f)
C   simpson     (in quadraturesubs.f)
C

      implicit none

      integer maxsize
      parameter(maxsize = 500)
      real array(0:maxsize)

      integer size
      integer i
 
      real dx,x
      real trap_int,simp_int
      real f
      real xmax

      real exact  ! exact answer
      print*,' Enter max x for integral '
      read*,xmax

      exact = 1.-exp(-xmax)*(1+xmax)

      open(unit=2,file='trapezoid.dat',status='unknown')
      open(unit=3,file='simpson.dat',status='unknown')
     
      do size = 4,100,2
         dx = xmax /float(size)
         do i = 0,size
           x = i*dx
           f = x*exp(-x)
           array(i) = f
         enddo
         call trapezoid(maxsize,size,array,dx,trap_int)
         call simpson(maxsize,size,array,dx,simp_int)
         write(2,*)dx,abs(trap_int-exact)
         write(3,*)dx,abs(simp_int-exact)
      enddo

      close(unit=2)  ! not really needed here but good form
      close(unit=3)

      end