      program integrate      
      implicit double precision (a-h,o-z)
      parameter (Nmx=10000000)
      dimension x(Nmx), y(Nmx)
      dimension xint(Nmx), yint(Nmx)
      character*80 nomefile, arg
c----------------------------------------------------------------------
c
      x1 = -1.D16
      x2 = -x1
      if (iargc().eq.1) then 
         call getarg(1,nomefile)
      else if (iargc().eq.2) then 
         call getarg(1,arg)
         call getarg(2,nomefile)
         read (arg,*) x1
      else if (iargc().eq.3) then  
         call getarg(1,arg)
         read (arg,*) x1
         call getarg(2,arg)
         read (arg,*) x2 
         call getarg(3,nomefile)
      else
         print*, ' Integrates a given tabulated function, stored     '
         print*, ' in a two-column file of the type x(i), y(i).      '
         print*, ' The function is integrated from x1 to x2.         '
         print*, '                                                   '
         print*, ' Usage: integrate [x1,x2] nomefile                 '
         STOP
      end if
c
c----------------------------------------------------------------------
c
      open (1,file=nomefile,status='unknown')
      ic   = 0
      iint = 0
 10   ic = ic + 1
      read (1,*,end=20) x(ic), y(ic)
      if (x(ic).gt.x1.and.x(ic).le.x2) then 
         iint = iint + 1
         xint(iint) = x(ic)
         yint(iint) = y(ic)
      end if
      go to 10
 20   Npunti = iint
      print*, ' Read  ', ic-1, ' data'
      print*, ' Using ', Npunti, ' for the integration' 
      print*, ' - from x =', xint(1), ' to x =', xint(iint)

      h = (xint(Npunti)-xint(1))/(Npunti-1)
      call simp(Npunti,h,yint,area)
      print*, ' Area =', area

c
c---------------------------------------------------------------------
c



      STOP
      END
c
c**********************************************************************
c
      subroutine simp(n,h,f,area)
      implicit double precision (a-h, o-z)
      parameter (Nmx=10000000)    
      dimension f(Nmx)

      ss=0.0d0
      ss0=0.0d0
      ss1=0.0d0
      ss2=0.0d0
      do 500 i=2, n-1, 2
        ss0=ss0+f(i-1)
        ss1=ss1+f(i)
        ss2=ss2+f(i+1)
  500 continue
      area = h*(ss0+4.0d0*ss1+ss2)/3.0d0

c     ***** if n is even add the last slice separately

      if (mod(n,2).eq.0) then
        area = area +h*(5.0d0*f(n)+8.0d0*f(n-1)-f(n-2))/12.0d0
      endif

      return
      end    
