c A simple composite quadrature program written by Peter Jimack.
c Version 1: March 1st 2007.
c
c The program divides the domain into s subdomains (where s is provided by
c the user at run time. The integral on each subdomain is computed using
c three point Guass quadrature.
c
      program quad
      implicit none

      double precision xl,xr,h,total,findint
      integer s,m

      write(6,*)'Input the two end points for the integration:'
      read(5,*)xl,xr
      write(6,*)'Input the number of sub-intervals for the quadrature:'
      read(5,*)s
      h = (xr-xl)/s
      xr = xl + h
      total = 0.0d0
      do 10 m = 1,s
         total = total + findint(xl,xr)
         xl = xr
         xr = xl + h
   10 continue

      write(6,*)'The estimate of the integral is: ',total

      stop
      end

cccccccccccccccccccccccccc
c
c This function defines the particular quadrature rule that us being used.
c In this case it is the three point Gauss rule.
c
      double precision function findint(xl,xr)
      implicit none

      double precision xl,xr,h,x1,x2,x3,f

      h = xr-xl
      x1 = xl + h*(1.0d0-dsqrt(0.6d0))/2.0d0
      x2 = xl + h/2.0d0
      x3 = xl + h*(1.0d0+dsqrt(0.6d0))/2.0d0
      findint = h*(5.0d0*f(x1)+8.0d0*f(x2)+5.0d0*f(x3))/18.0d0

      return
      end

cccccccccccccccccccccccccc
c
c This function defines the integrand
c
      double precision function f(x)
      implicit none

      double precision x,pi

      pi = 3.1415926535897932d0
      f = dsin(pi*x)

      return
      end