CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      program metatron
C
C  sample use of metropolis algorithm for MC quadrature
C
C  computes < x^2 exp(-x^2 > / < exp(-x^2) >
C
C  version 0: uses Gaussian weighting
C
C  CWJ -SDSU - Sept 2006
C
C  SUBROUTINES CALLED
C    metropolis (in metrolib)
C
C  also: need to define real functions:
C  weight(x)   (weight factor)
C  func1(x)   (for numerator)
C  func0(x)   (for denominator)
C
C  then metropolis generates {x_i} according to weight(x)
C
C  and avg1 = < func1(x) >, avg0 = < func0(x) >
C
C  and final result = < func1(x) > / < func0(x) >
C
C  for VERSION 0,  weight = exp(-x^2)
C  so func0 = 1 and func1 = x^2
C

      implicit none

      real avg1, davg1  ! = < x^2 exp(-x^2) >, error
      real avg0, davg0  ! = <  exp(-x^2) >, error

      real dx           ! random step size

      real x0		! starting point
      
      integer seed	! random seed
      
      integer ntherm, nsamp,nskip
			! thermalization steps, # of MC samples, steps to skip
			! in order to decorrelate

      real ratio	! = < x^2 exp(-x^2 > / < exp(-x^2) >
      real dratio	!  error in ratio

     
      print*,' Enter seed '
      read*,seed
      print*,' enter starting point, dx ' 
      read*,x0,dx
 
      print*,' Enter ntherm, nskip, nsamp '
      read*,ntherm,nskip,nsamp

      call metropolis(x0,seed,dx,ntherm,nsamp,nskip,
     & avg1,davg1,avg0,davg0)

C      print*,' integer = ',avg, ' +/- ',davg

C................ NOW COMBINE.......

      ratio = avg1/avg0
C.......... relative errors add quadratically.........
      dratio = sqrt( (davg1/avg1)**2 + (davg0/avg0)**2)*ratio
      print*,' Answer = ',ratio,' +/- ',dratio

      print*,avg1,davg1,avg0,davg0

      end


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


      real function weight(x)
C
C  gaussian weight function
C 
      implicit none
      real x

      weight = exp(-x**2)
      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


      real function FUNC1(x) 

      implicit none
      real x

      func1 = x*x
      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


      real function FUNC0(x) 

      implicit none
      real x

      func0 = 1.
      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC