***SIMPSON'S RULE*******************************************************
*      PROGRAM FOR NUMERICAL INTEGRATION BY SIMPSON'S RULE FOR         *
*                  TABULATED DATA / KNOWN FUNCTION                     *
************************************************************************
*
      program SIMPSON
      implicit doubleprecision(a-h,o-z)
      doubleprecision Y(1000),integral,integral1,integral2
      integer select,p
*
*********************FUNCTION TO BE INTEGRATED**************************
      F(x) = (0.383d0/10**4*x**3 + 1.0d0*10**4) /
     $       (8.131d0/10**6*(81.0d0- x**4))
************************************************************************
*
      print*, ' Nature of input:'
      write(6,101) '[1]   TABULATED DATA','[2]   KNOWN FUNCTION','Enter
     $Your Selection < 1 - 2 >:'
101   format(''/a30/a30////a34)
      read(*,*) select
      if (select == 1) then         ! tabulated data
         write(6,*)'----------------------------------------------------
     $------------------'
         write(6,*) ' * WARNING: Simpson''s rule can only be used for ev
     $enly spaced data *'
         write(6,*)'----------------------------------------------------
     $------------------'
         print*, 'The number of data points, n:'
         read(*,*) n
         m = n - 1
         write(6,*)
         print*, 'Enter lower limit of integration, a:'
         read(*,*) a
         print*, 'Enter upper limit of integration, b:'
         read(*,*) b
         h = (b - a)/m
         write(6,*)
         print*, 'Enter values of dependent variable (data):'
         read(*,*) (Y(i), i = 1,n)
      elseif(select == 2) then      ! known function
         write(6,*)
         print*, 'Enter lower limit of integration, a:'
         read(*,*) a
         print*, 'Enter upper limit of integration, b:'
         read(*,*) b
         print*, 'Number of segments required:'
         read(*,*) m
         h = (b - a)/m
*
         n = m + 1
         do i = 1, n
            Y(i) = F(a + (i - 1)*h)
         enddo
      endif
c
c     Simpson's rule
      p    = (m/2)*2
      sum1 = 0
      sum2 = 0
*
      if (p == m) then        ! even number of segments
         do i = 2, m, 2
            sum1 = sum1 + Y(i)
         enddo
         mm = m - 1
         do i = 3, mm, 2
            sum2 = sum2 + Y(i)
         enddo
         hight = (Y(1) + 4*sum1 + 2*sum2 + Y(n)) /(3*m)
         integral = (b - a)*hight
*
      elseif (p.ne.m) then    ! odd number of segments
         if (m.ge.5) then
            mm = m - 3
            do i = 2, mm, 2
               sum1 = sum1 + Y(i)
            enddo
            mmm = mm - 1
            do i = 3, mmm, 2
               sum2 = sum2 + Y(i)
            enddo
            hight1 = (Y(1) + 4*sum1 + 2*sum2 + Y(n-3)) /(3*(m-3))
            integral1 = ((b - 3*h) -a) *hight1
         endif
*
         integral2 = 0
         hight2    = (Y(n-3) + 3*(Y(n-2) + Y(n-1)) + Y(n)) /8d0
         integral2 = 3*h *hight2
         integral  = integral1 + integral2
      endif
*
      write(6,*) '-----------------------------------------------------'
      write(6,102) h,a,b,integral
102   format('',21x'The width of segment is:',f12.6//5x,'The value of in
     $tegral with integration limits',f7.2,2x,'and',f7.2,3x,'is:'/30x,f1
     $3.4)
      write(6,*) '-----------------------------------------------------'
*
      stop
      end