c ctrapz1.f c c
      program trapezoid
c
c    Test Trapezoidal Rule of Integration
c
c      John Mahaffy  4/9/95
c
c   Include the contents of file "trapcom.h" in this program
c
      include 'trapcom.h'
      integer i
      do 10 i=1,ntrys
         np=nint1*2**(i-1)+1
         call setcurve
         call integrate
         call compare
  10  continue
      stop
      end
      subroutine setcurve
c
c    Set up a Curve and range for integration
c    Here we will integrate sin(x) from x=0 to x=3.0
c
c    John Mahaffy 4/9/95
c
      include 'trapcom.h'
      integer i
      real dx,xlow,xhigh
      data xlow,xhigh/0.0,3.0/
c
      dx=(xhigh-xlow)/(np-1)
      x(1)=xlow
c
      y(1)=sin(x(1))
c
      do 10 i=2,np
         x(i)=x(i-1)+dx
         y(i)=sin(x(i))
  10     continue
      ansreal=cos(xlow)-cos(xhigh)
      return
      end
      subroutine integrate
c
c    Use Trapezoidal Rule to Integrate the area under the curve
c    specified by the data points in x and y
c
c    John Mahaffy 4/9/95
c
      include 'trapcom.h'
c
c   Local Variables
c
c    esterr   -   estimated error for Trapizoidal integration
c    sum2     -   Trapizoidal integration using every other available point
c
      integer i
      real sum2, esterr
      ansnum=0.0
      do 10 i=1,np-1
   10    ansnum=ansnum+.5*(y(i)+y(i+1))*(x(i+1)-x(i))
c
c   The following only works if np is an odd integer
c
      sum2=0.0
      do 20 i=1,np-1,2
   20    sum2=sum2+.5*(y(i)+y(i+2))*(x(i+2)-x(i))
      esterr=(sum2-ansnum)*.333333
      write(6,2000) np,ansnum,esterr
 2000 format (' For ',i4,' points the integral is ',1p,e12.5,
     & ', Estimated Error=',e10.2)
      return
      end
      subroutine compare
c
c  Compare results of the numerical integration with the actual result
c
c    John Mahaffy 4/9/95
c
      implicit none
c
c  ansreal -  real answer from analytic integration
c  ansnum  -  answer resulting from the numerical integration
c
      include 'trapcom1.h'
      real error
      error=ansnum-ansreal
      write(6,2001) error
 2001 format(' Actual Error = ',1p,e10.2)
      return
      end
c
c c