c------------------------------------------------------------------------------
      PROGRAM AHB1RES
c------------------------------------------------------------------------------
c     This program calculates averages and statistical errors of the binned
c     data in the files 'len.dat' and 'mag.dat' generated by the program AHB1. 
c------------------------------------------------------------------------------
c     Written by A. W. Sandvik (asandvik@abo.fi), Åbo Akademi University.
c     Last modified: Sept. 2, 2002.
c------------------------------------------------------------------------------

      IMPLICIT NONE

      INTEGER maxr
      PARAMETER (maxr=2000)

      INTEGER i,l,d,n,nr
      REAL*8  beta,av1,er1,av2,er2,d1(maxr),d2(maxr),d3(maxr),d4(maxr)

      OPEN(UNIT=10,FILE='read.in',STATUS='old')
      READ(10,*)d,l,beta
      CLOSE(10)
      n=l**d

      nr=0
      OPEN(UNIT=10,FILE='len.dat',STATUS='old')
      DO i=1,maxr
         READ(10,*,END=10)d1(i),d2(i)
         nr=nr+1
      ENDDO
 10   CLOSE(10)      
      WRITE(*,*)'Number of bins ',nr
      WRITE(*,*)
      CALL HEATCAP(d,n,beta,nr,d1,d2,av1,er1,av2,er2)
      WRITE(*,20)' Internal energy             =  ',av1,er1
      WRITE(*,20)' Specific heat               =  ',av2,er2
 20   FORMAT(A,2F14.7)

      nr=0
      OPEN(UNIT=10,FILE='mag.dat',STATUS='old')
      DO i=1,maxr
         READ(10,*,END=30)d1(i),d2(i),d3(i),d4(i)
         nr=nr+1
      ENDDO
 30   CLOSE(10)      
      CALL AVERAGE(d1,av1,er1,nr)
      WRITE(*,20)' Uniform susceptibility      =  ',av1,er1
      CALL AVERAGE(d2,av1,er1,nr)
      WRITE(*,20)' Spin stiffess               =  ',av1,er1
      CALL AVERAGE(d3,av1,er1,nr)
      WRITE(*,20)' Staggered structure factor  =  ',av1,er1
      CALL AVERAGE(d4,av1,er1,nr)
      WRITE(*,20)' Staggered susceptibility    =  ',av1,er1

      END

c-----------------------------------------------------------------------------
      SUBROUTINE AVERAGE(dat,av,er,nr)
c-----------------------------------------------------------------------------

      IMPLICIT NONE

      INTEGER maxr
      PARAMETER (maxr=2000)

      INTEGER i,nr
      REAL*8  dat(maxr),av,er     

      av=0.d0
      er=0.d0
      DO i=1,nr
         av=av+dat(i)
         er=er+dat(i)**2
      ENDDO
      av=av/DFLOAT(nr)
      er=er/DFLOAT(nr)
      er=DSQRT(DABS(er-av**2)/DFLOAT(nr-1))

      RETURN
      END

c-----------------------------------------------------------------------------
      SUBROUTINE HEATCAP(d,n,beta,nr,len1,len2,ave,ere,avc,erc)
c-----------------------------------------------------------------------------
c     Calculates the internal energy and the specific heat per spin based on
c     the formulas E=d/4-<nh>/(beta*N), C=(<nh**2>-<nh>**2-<nh>)/N, where
c     d is the dimensionality, beta the inverse temperature J/T, N the
c     number of spins, and nh the expansion order. The statistical errors
c     are calculated using the bootstrap method. The arrays len1 and len2
c     contain the SSE data; len1=<nh>/(beta*N), len2=<nh**2>/(beta*N)**2.
c-----------------------------------------------------------------------------

      IMPLICIT NONE

      INTEGER maxr
      PARAMETER (maxr=2000)

      INTEGER d,i,k,n,nr,nb,iboot,nboot
      REAL*8  e,c,l1,l2,ave,ere,avc,erc,beta
      REAL*8  RAN,len1(maxr),len2(maxr)    

      nboot=200
      l1=0.d0
      l2=0.d0
      DO i=1,nr
         l1=l1+len1(i)
         l2=l2+len2(i)
      ENDDO
      l1=l1/DFLOAT(nr)
      l2=l2/DFLOAT(nr)
      ave=DFLOAT(d)*0.25d0-l1
      avc=(l2*beta**2-(l1*beta)**2-l1*beta/DFLOAT(n))*DFLOAT(n)
      ere=0.d0
      erc=0.d0
      DO iboot=1,nboot
         l1=0.d0
         l2=0.d0
         DO i=1,nr
            k=MIN(INT(RAN()*nr)+1,nr)
            l1=l1+len1(k)
            l2=l2+len2(k)
         ENDDO
         l1=l1/DFLOAT(nr)
         l2=l2/DFLOAT(nr)
         e=DFLOAT(d)*0.25d0-l1
         c=(l2*beta**2-(l1*beta)**2-l1*beta/DFLOAT(n))*DFLOAT(n)
         ere=ere+(ave-e)**2
         erc=erc+(avc-c)**2
      ENDDO
      ere=DSQRT(ere/DFLOAT(nboot))
      erc=DSQRT(erc/DFLOAT(nboot))

      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE INITRAN(seed)
c------------------------------------------------------------------------------

      IMPLICIT NONE

      INTEGER iir,jjr,kkr,nnr,mzran
      COMMON/BRANDOM/iir,jjr,kkr,nnr,mzran

      INTEGER seed(4)

      OPEN(UNIT=10,FILE='rand.in',STATUS='old')
      READ(10,*)seed(1)
      READ(10,*)seed(2)
      READ(10,*)seed(3)
      READ(10,*)seed(4)
      CLOSE(10)      
      iir=1+iabs(seed(1))
      jjr=1+iabs(seed(2))
      kkr=1+iabs(seed(3))
      nnr=seed(4)

      RETURN
      END

c------------------------------------------------------------------------------
      REAL*8 FUNCTION RAN()
c------------------------------------------------------------------------------

      IMPLICIT NONE

      INTEGER iir,jjr,kkr,nnr,mzran
      COMMON/BRANDOM/iir,jjr,kkr,nnr,mzran

      mzran=iir-kkr
      IF (mzran.LT.0) mzran=mzran+2147483579
      iir=jjr
      jjr=kkr
      kkr=mzran
      nnr=69069*nnr+1013904243
      mzran=mzran+nnr
      RAN=0.5d0+0.23283064d-9*mzran

      RETURN
      END

c------------------------------------------------------------------------------