***INTEGRATION OF UNEQUALLY SPACED DATA*********************************
*     PROGRAM FOR NUMERICAL INTEGRATION FOR UNEQUALLY SPACED DATA      *
*                USING TRAPEZOIDAL AND SIMPSON'S RULE                  *
************************************************************************
*
      program UNEQSIMPSON
      implicit doubleprecision (a-h,o-z)
      parameter (nd = 100)
      doubleprecision X(nd),Y(nd),H(nd),integral
*
      write(6,*)
      print*, 'Enter number of data points, n:'
      read(*,*) n
      write(6,*)
      print*, 'Enter values of independent variable:'
      read(*,*) (X(i), i = 1,n)
      write(6,*)
      print*, 'Enter values of dependent variable:'
      read(*,*) (Y(i), i = 1,n)
c
c     evaluation of segment width
      m = n - 1
      do i = 1, m
         H(i) = X(i + 1) - X(i)
      enddo
      nn = n + 1
      do i = n, nn
         H(i) = 0
      enddo 
*
      sum1 = 0
      sum2 = 0
      sum3 = 0
      i = 1
c
c     loop for integration
10    if (i < n) then
c        Simpson's 3/8 rule
         if ((H(i) == H(i+1)).and.(H(i) == H(i+2))) then
            sum1 = sum1 + (3*H(i) *(Y(i) + 3*(Y(i+1)+Y(i+2)) + Y(i+3))) 
     $                                                              /8d0
            i = i + 3
c        Simpson's 1/3 rule
         elseif (H(i) == H(i+1)) then
            sum2 = sum2 + (2*H(i) * (Y(i)+ 4*Y(i+1)+ Y(i+2))) /6d0
            i = i + 2
c        trapezoidal rule
         elseif (H(i).ne.H(i+1)) then
            sum3 = sum3 + H(i) * (Y(i)+ Y(i+1)) /2d0
            i = i + 1
         endif
         go to 10
      endif
*
      integral = sum1 + sum2 + sum3
*
      write(6,*) '-----------------------------------------------------'
      write(6,101) X(1),X(n),integral
101   format('',5x,'The value of integral with integration limits',f7.2,
     $2x,'and',f7.2,3x,'is:'/30x,f13.4)
      write(6,*) '-----------------------------------------------------'
*
      stop
      end