!     Last change:  DOS  31 Jul 2000    7:59 pm
!                                       *** copyright 2000 ***
! *** filename quad1.f95                *** John F. Monahan **
!                                       **********************
      program quad1
!  demonstration of one dimension numerical integration
!
!  compute posterior mean & variance for Example 10.1
!
      implicit none
      real t,trap0,trap1,trap2,midp0,midp1,midp2,simp0,simp1, &
     &     simp2,fkint,ps,pstar
      integer i,kint
!
  20  format(/'  Number of intervals ',i5/18x,'Normalization',5x,     &
     &       'Posterior',7x,'Posterior'/20x,'constant',10x,'mean',    &
     &       10x,'variance')
  21  format(' Trapezoid rule ',e16.6,2f14.8)
  22  format('  Midpoint rule ',e16.6,2f14.8)
  23  format(' Simpson"s rule ',e16.6,2f14.8)
!
      open( unit=6, file='quad1.out' )
!
      do kint = 2,20
      fkint = real(kint)
      write(6,20) kint
!                   kint is number of intervals
!
!                   trapezoid rule first
      trap0 = 0.
      trap1 = 0.
      trap2 = 0.
!                   integration on the interval (0, 1)
!                    half weight on the two endpoints
      t = 0.
      ps = pstar(t)/2.
      trap0 = trap0 + ps
      trap1 = trap1 + ps*t        
      trap2 = trap2 + ps*t*t  
      t = 1.
      ps = pstar(t)/2.
      trap0 = trap0 + ps
      trap1 = trap1 + ps*t        
      trap2 = trap2 + ps*t*t  
!                    full weight on the other points
      do i = 2,kint
      t = float(i-1)/fkint
      ps = pstar(t)
      trap0 = trap0 + ps
      trap1 = trap1 + ps*t        
      trap2 = trap2 + ps*t*t  
      end do  ! loop on i
!                   now the midpoint rule
      midp0 = 0.
      midp1 = 0.
      midp2 = 0.
      do i = 1,kint
      t = real(2*i-1)/(2.*fkint)
      ps = pstar(t)
      midp0 = midp0 + ps
      midp1 = midp1 + ps*t        
      midp2 = midp2 + ps*t*t  
      end do  ! loop on i
!                  Simpson is the weighted average of trap & midp
      simp0 = ( 2.*midp0 + trap0 ) / 3.
      simp1 = ( 2.*midp1 + trap1 ) / 3.
      simp2 = ( 2.*midp2 + trap2 ) / 3.
!
!                  now convert to get posterior means and variances
      trap1 = trap1 / trap0
      trap2 = trap2/trap0 - trap1*trap1
      midp1 = midp1 / midp0 
      midp2 = midp2/midp0 - midp1*midp1
      simp1 = simp1 / simp0 
      simp2 = simp2/simp0 - simp1*simp1
!                   finally rescale weights
      trap0 = trap0 / fkint
      midp0 = midp0 / fkint
      simp0 = simp0 / fkint
      write(6,21) trap0,trap1,trap2
      write(6,22) midp0,midp1,midp2
      write(6,23) simp0,simp1,simp2
      end do  ! loop on kint
      stop
      end program quad1
      real function pstar(t)
!  posterior from logarithmic series example 10.1      
      implicit none
      real, INTENT(IN) :: t
      pstar = 0.
!                   first take care of values out of range
      if( (t .le. 0.) .or. (t .ge. 1.) ) return
!                   evaluate function
      pstar = (1.-t)*( t**16 )/ ( (-log(1.-t))**10 )
      return
      end function pstar
! *** end of filename quad1.f95         *********************