CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      subroutine metropolis(x0,seed,dx,ntherm,nsamp,nskip,
     & avg1,davg1,avg0,davg0)
C
C  metropolis subroutine
C
C  INPUT:
C   x0: starting point
C   seed: random number seed
C   dx : scale of random steps
C   ntherm: number of thermalization steps
C   nsamp: # of MC samples wanted
C   nskip: for decorrelation
C  OUTPUT:
C     avg1,avg0: MC integrals in numerator, denominator
C     davg1,davg0 : error in MC integral
C
C  CALLS
C    function FUNC1,FUNC0
C    defined elsewhere (in metatron.f)
C
C   given a WEIGHT function used in metrostep,
C   computes avg1 = < func1 > and avg0 = < func0 >
C   as well as statistical averages
C

      implicit none

C....... INPUT...............

      real x0  ! starting point
      integer seed ! random seed
      real dx       ! scale of random steps
      integer ntherm,nsamp ! # of thermalization steps, # of MC samples
      integer nskip

C........OUTPUT...............

      real avg1,avg0  ! MC integral
      real davg1,davg0  ! error in MC integral

C..........INTERMEDIATE..........

      real x   ! location
      logical accept ! flag for acceptance

      integer naccept, ntries ! # of acceptances and ntries
      real acceptrate         ! acceptance rate

      integer i
      integer isamp
      integer iskip

      real FUNC1,FUNC0        ! defined elsewhere



      x = x0
C.............. THERMALIZATON LOOP

      naccept = 0

      do i = 1,ntherm
         call metrostep(x,dx,seed,accept)
         if(accept)naccept=naccept+1
      enddo
      acceptrate = float(naccept)/float(ntherm)
      print*,' during thermalization, acceptance rate = ',
     & acceptrate       

C........... OPEN FILES FOR OUTPUT TO BE POST-PROCESSED.......

      open(unit=10,file='walker.dat',status ='unknown')
      open(unit=11,file='func1.dat',status='unknown')
      open(unit=12,file='func0.dat',status='unknown')

C............. COMPUTE INTEGRAL

      avg1 = 0.
      davg1 = 0.
      avg0 = 0.0
      davg0 = 0.
      isamp = 0
      iskip = 0
      
      naccept = 0
      ntries = 0
      do while(isamp .lt. nsamp)
         call metrostep(x,dx,seed,accept)
         ntries  = ntries +1
         if(accept)then
             naccept = naccept + 1
             iskip = iskip + 1
             write(11,*)func1(x)
             write(12,*)func0(x)
             if(iskip .eq. nskip)then

                isamp = isamp + 1
                avg1 = avg1 + FUNC1(x)
                davg1 = davg1 + FUNC1(x)**2
                avg0 = avg0 + FUNC0(x)
                davg0 = davg0 + FUNC0(x)**2

                write(10,*)x
                iskip = 0
             endif
         endif

      enddo
      print*,' acceptance rate = ',float(naccept)/float(ntries)

      avg0 = avg0 /float(nsamp)
      davg0 = davg0/float(nsamp) - avg0*avg0
      davg0 = sqrt(davg0/float(nsamp))

      avg1 = avg1 /float(nsamp)
      davg1 = davg1/float(nsamp) - avg1*avg1
      davg1 = sqrt(davg1/float(nsamp))
      return
      end



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc
C
      subroutine metrostep(x,dx,seed,accept)
C
C takes a step and then accepts or rejects it 
C
C input:
C    x:  position
C    dx : scale of random step (from -dx to dx)
C    seed: random number seed
C OUTPUT:
C    accept: logical
C
C needs externally defined weight function Weight(x)
C
      implicit none

C..............INPUT......................
      real x   ! starting point
      real dx    ! range of random step
      integer seed  ! random seed

C.............OUTPUT......................

      logical accept   ! flag for accept/reject

C............INTERMEDIATE..................

      real ran3   ! random number generator
      real xstep   ! actual random step
      real xtmp    ! temp position

      real Weight   ! weight function
      real ratio    ! ratio of weight functions
      real die      ! random number to compare

      xstep = 2*dx*ran3(seed)-dx   ! generates a uniform random # between -dx and dx
      xtmp = x+xstep

      ratio = weight(xtmp)/weight(x)   ! note, for real applications, if weight is
                                        ! complicated function, save weight(x) from last time

      die = ran3(seed)       ! random number between 0 and 1

      if(ratio .ge. die)then 

          accept = .true.
          x = xtmp

      else

          accept = .false.  ! reject
      endif
      return
      end