c --- program integ: Numerical integration
      common/one/x(500),y(500),s,xl,xh,n
      open(unit=1, file='integ.txt')
      open(unit=8, file='int-out.txt')
c -----------------------------------------
      write(6,5)
      write(8,5)
    5 format('Numerical integration (piecewise quadratic).')
      n=0
   10 n=n+1
      read(1,*,end=20) x(n),y(n)
      goto 10
   20 n= n-1
      if(n .gt. 2) goto 30 
      write(5,*) 'Insufficient data in input file integ.dat'
      stop
   30 write(6,15) n	  
      write(8,15) n	  
   15 format('The',i3,' data points are:')
      do 50 i=1,n
      write(6,25) x(i),y(i)
   50 write(8,25) x(i),y(i)
   25 format(2f12.5)
   60 write(6,65)
   65 format('Enter the limits of integration',
     1       ' (separate with comma or space):')
      read(5,*) xl,xh
      if(xh .gt. xl) goto 70
      write(6,*) 'Limits not ok, try again (lower first).'
      goto 60
   70 continue
      call int
      write(6,95) xl,xh,s 
      write(8,95) xl,xh,s 
   95 format('The integral from ',f10.4,' to',f10.4,' is',f15.7)
      write(6,*) 'Type 1 to repeat with different limits, 0 to stop.'
      read(5,*) k
      if(k .eq. 1) goto 60
      stop
      end
c ----------------------
      subroutine int
c ----------------------
c --- piecewise quadratic
      common/one/x(500),y(500),s,xl,xh,n
      s=0.
      do 500 k=2,n-1
      x1= x(k-1)
      x2= x(k)
      x3= x(k+1)
      y1= y(k-1)
      y2= y(k)
      y3= y(k+1)
      d= (x3-x2)*(x2**2-x1**2)-(x2-x1)*(x3**2-x2**2)
      a= ((x3-x2)*(y2-y1)-(x2-x1)*(y3-y2))/d
      b= ((y2-y1) -a*(x2**2-x1**2))/(x2-x1)
      c= y1 -b*x1 -a*x1**2
      t= 0.5*(x1+x2)
      u= 0.5*(x2+x3)
      if(k .eq. 2) t= xl
      if(k .eq. (n-1)) u= xh
      if(xl .gt. t) t= xl
      if(xh .lt. u) u= xh
      if(t .gt. u) goto 500
      s= s +(a/3.)*(u**3-t**3) +(b/2.)*(u**2-t**2) +c*(u-t)
  500 continue     
      return
      end