!     Last change:  DOS  31 Jul 2000    8:01 pm
!                                       *** copyright 2000 ***
! *** filename quad2.f95                *** John F. Monahan **
!                                       **********************
      program quad2
!  demonstration of one dimension numerical integration
!
!  compute posterior mean & variance for Example 10.1
!
      implicit none
      real t,s0,s1,s2,tm,se,ab,wgt,ps,pstar
      integer i,np
!
  20  format(/'  Number of evaluations',i5/23x,'Normalization',5x,    &
     &       'Posterior',7x,'Posterior'/26x,'constant',9x,'mean',     &
     &       11x,'variance')
  21  format(' Mod Gauss-Hermite',3x,e12.4,2f16.8)
  22  format(' Shift Gauss-Legendre',e12.4,2f16.8)
  28  format(4x,2f16.8)
!
!                   put output here
      open( unit=6, file='quad2.out' )
!                   read from these tables
      open( unit=8, file='qgaushm.tab' )
      open( unit=9, file='qgauslg.tab' )
!
!                   mode of posterior for normal approx
      tm = .527167
!                   standard error from second deriv at mode of post
      se = 1./sqrt(42.036)
!
!                   these numbers of points
      do np = 3,10
      write(6,20) np
!                   do the modified gauss-hermite first      
      s0 = 0.
      s1 = 0.
      s2 = 0.
      do i = 1,np
!                   read abscissa and weight from table
      read(8,28) ab,wgt
!                   relocate and rescale by mode and 2nd deriv      
      t = tm + se*ab
!                   divide by normal density
      ps = wgt*pstar(t)*2.506628*exp(ab*ab/2.)
      s0 = s0 + ps
      s1 = s1 + ps*t        
      s2 = s2 + ps*t*t  
      end do  ! loop on i
!                  now convert to get posterior means and variances
      s1 = s1 / s0
      s2 = s2/s0 - s1*s1
!                  rescale s0 by se
      s0 = s0 * se
      write(6,21) s0,s1,s2
!                   now the shifted gauss-legendre      
      s0 = 0.
      s1 = 0.
      s2 = 0.
      do i = 1,np
!                   read abscissa and weight from table
      read(9,28) ab,wgt
!                   gauss-legendre shifted to (0, 1)      
      t = ab
      ps = wgt*pstar(t)
      s0 = s0 + ps
      s1 = s1 + ps*t        
      s2 = s2 + ps*t*t  
      end do  ! loop on i
!                  now convert to get posterior means and variances
      s1 = s1 / s0
      s2 = s2/s0 - s1*s1
      write(6,22) s0,s1,s2
      end do  ! loop on np
      stop
      end program quad2
      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 quad2.f95         *********************