program trapezoidal_demo
C
C   this program illustrates integration of the exponential function
C   using the trapezoidal rule
C
C   this version uses error traps
C
C   CWJ SDSU 2/3/2005
C

      implicit none

C....... define lattice................

      real a,b	! endpoints
      real temp  ! to be used in error checking
      real dx    ! 
      integer N  ! # of points

C....... define array.................

      integer size
      parameter (size=500)
      real f(0:size)

C........ dummy indices..............

      integer i
      real x    
C......... integral................
 
      real sum
      real exactsum   ! for comparison

C............... NOW BEGIN..........


      print*,' Demonstration of trapezoidal rule '
      print*,' on function exp(x) on interval (a,b) '
      print*,' '
    
C...............ENTER PARAMETERS.......

      print*,' enter a,b of interval '
      read*,a,b

C................ CHECK ORDER.............

      if(a.gt.b)then
	temp = b
        b = a
        a = temp
      endif

C.................. READ IN NUMBER OF POINTS
    
      print*,' Enter # of points to be used '
      read*,N
 
C............ ERROR TRAP....................

      if(N.gt.size)then
	print*,' Sorry, must be less than ',size
        print*,' You might have to change SIZE and recompile '
        stop
      endif

      dx = (b-a)/float(n)

C.............. NOW FILL OUT ARRAY...........

      do i = 0,N
  	x = a + i*dx
        f(i) = exp(x)
      enddo
      
C..............INITIALIZE INTEGRAL .........

      sum = 0.0    

C..............apply trapezoidal rule......

      do i = 1,N
	sum = sum + 0.5*dx*(f(i)+f(i-1))
      enddo

C..............print results

      exactsum = exp(b)-exp(a)
      print*,' exact = ',exactsum,' approx = ',sum,
     & ' % error = ',100.*abs(exactsum-sum)/exactsum


      end