
C----------------------------------------------------------------------C
C                                                                      C
C  Initialize values for accumulation of measurements.                 C
C                                                                      C
C----------------------------------------------------------------------C

      SUBROUTINE init_measure
      INCLUDE "common.h"

C ----  Initialize sums and sum of squares.
      ibin = 0
      sum_energy = 0.0
      sum_enersq = 0.0
      sum_spins  = 0.0
      sum_spinsq = 0.0
      bin_energy = 0.0
      bin_enersq = 0.0
      bin_spins  = 0.0
      bin_spinsq = 0.0
      avbin_spheat   = 0.0
      avbin_spheatsq = 0.0
      avbin_susp     = 0.0
      avbin_suspsq   = 0.0
      avbin_enersq = 0.0
      avbin_spinsq = 0.0
 
      RETURN
      END



C----------------------------------------------------------------------C
C                                                                      C
C  Take measurements of energy and magnetization.                      C
C                                                                      C
C----------------------------------------------------------------------C

      SUBROUTINE measure
      INCLUDE "common.h"
      INTEGER i,j, ener,magn
      REAL energy,magnetization
      REAL*8 bin_mean, bin_meansq, bin_spheat, bin_susp

      ener = 0
      magn = 0
      DO i=1,side
      DO j=1,side
        ener = ener + spin(i,j) * (spin(pos(i),j) + spin(i,pos(j)))
        magn = magn + spin(i,j)
      ENDDO
      ENDDO
 
C --- Normalize.
      IF (magn.LT.0) THEN
        magn = - magn
      ENDIF
      energy = - FLOAT(ener) / volume
      magnetization = FLOAT(magn) / volume
      sum_accept = sum_accept + FLOAT(accept)/volume
 
C
C  Bin the data in order to calculate errors 
C
      ibin = ibin + 1

      bin_energy = bin_energy + energy
      bin_enersq = bin_enersq + (energy * energy)
      
      bin_spins  = bin_spins + magnetization
      bin_spinsq = bin_spinsq + (magnetization * magnetization)

      IF (ibin.EQ.bin_size) THEN
C
C     Add the bin values to the total value
C
         sum_energy = sum_energy + bin_energy
         sum_enersq = sum_enersq + bin_enersq
C        
         sum_spins  = sum_spins  + bin_spins
         sum_spinsq = sum_spinsq + bin_spinsq
C
C     Calculate bin means and mean squares 
C
         bin_mean = bin_energy / bin_size
	 bin_meansq = bin_mean * bin_mean
	 bin_spheat = (bin_enersq/bin_size) - bin_meansq
	 avbin_spheat = avbin_spheat + bin_spheat
	 avbin_spheatsq = avbin_spheatsq + (bin_spheat * bin_spheat)
	 avbin_enersq = avbin_enersq + bin_meansq
C        
	 bin_mean = bin_spins / bin_size
	 bin_meansq = bin_mean * bin_mean
	 bin_susp = (bin_spinsq/bin_size) - bin_meansq
	 avbin_susp = avbin_susp + bin_susp
	 avbin_suspsq = avbin_suspsq + (bin_susp * bin_susp)
	 avbin_spinsq = avbin_spinsq + bin_meansq
C
C     Reset bin sums to zero for new bin
C
	 ibin = 0
         bin_energy = 0.0
         bin_enersq = 0.0
         bin_spins  = 0.0
         bin_spinsq = 0.0

      ENDIF

      RETURN
      END



C----------------------------------------------------------------------C
C                                                                      C
C   Calculate averages and errors and print out results.               C
C                                                                      C
C----------------------------------------------------------------------C

      SUBROUTINE output
      INCLUDE "common.h"
      REAL*8 av_energy, av_enersq, av_spins, av_spinsq, av_accept
      REAL*8 err_ener, err_spin, err_spheat, err_susp
      REAL*8 specific_heat, susceptibility, variance, scale
C
C  Take averages.
C
      av_energy = sum_energy / measurements
      av_enersq = sum_enersq / measurements
      av_spins  = sum_spins  / measurements
      av_spinsq = sum_spinsq / measurements
      av_accept = sum_accept / measurements

      avbin_enersq   = avbin_enersq   / num_bins
      avbin_spheat   = avbin_spheat   / num_bins
      avbin_spheatsq = avbin_spheatsq / num_bins
      avbin_spinsq   = avbin_spinsq   / num_bins
      avbin_susp     = avbin_susp     / num_bins
      avbin_suspsq   = avbin_suspsq   / num_bins
C
C  Work out variances and errors.
C
C --- Energy
C
      variance = avbin_enersq - (av_energy * av_energy)
      err_ener = SQRT( variance / num_bins )
C
C --- Magnetization
C
      variance = avbin_spinsq - (av_spins * av_spins)
      err_spin = SQRT( variance / num_bins )
C
C --- Specific heat 
C
      specific_heat  = av_enersq - (av_energy * av_energy)
      scale = volume * beta * beta
      specific_heat  = specific_heat * scale
      variance = avbin_spheatsq - (avbin_spheat * avbin_spheat)
      err_spheat = SQRT( variance / num_bins )
      avbin_spheat = avbin_spheat * scale
      err_spheat = err_spheat * scale
C
C --- Susceptibility 
C
      susceptibility = av_spinsq - (av_spins * av_spins)
      susceptibility = susceptibility * scale
      variance = avbin_suspsq - (avbin_susp * avbin_susp)
      err_susp = SQRT( variance / num_bins )
      avbin_susp = avbin_susp * scale
      err_susp = err_susp * scale

C
C  Print out results.
C
      PRINT *, ' '
      PRINT 10, ' Energy         is ', av_energy,' +/-', err_ener
      PRINT *, ' '
      PRINT 10, ' Magnetization  is ', av_spins, ' +/-', err_spin
      PRINT *, ' '
      PRINT 30, ' Specific heat  is ', specific_heat, ' +/-', err_spheat
      PRINT *, ' '
      PRINT 30, ' Susceptibility is ', susceptibility, ' +/-', err_susp
      PRINT *, ' '
      PRINT 20, ' Acceptance     is ', av_accept
  10  FORMAT(A,F10.6,A,F9.6)
  20  FORMAT(A,F12.5)
  30  FORMAT(A,F12.5,A,F11.5)

      RETURN
      END


