ccccccccccccccccccccccccc     Program 9.2     cccccccccccccccccccccccccc
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c Please Note:                                                         c
c                                                                      c
c (1) This computer program is part of the book, "An Introduction to   c
c     Computational Physics," written by Tao Pang and published and    c
c     copyrighted by Cambridge University Press in 1997.               c
c                                                                      c
c (2) No warranties, express or implied, are made for this program.    c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      PROGRAM MCMA
C
C Integral evaluated with the Metropolis Monte Carlo scheme.
C The distribution function is W(x)= (exp(x*x)-1)/Z and the
C function sampled is is g(x) = Z*x*x/(exp(x*x)-1).
C
      PARAMETER (ITHM=10000,ISIZ=10000,L=15)
      INTEGER*4 time,STIME,T(9)
      COMMON /CSEED/ ISEED
      COMMON /CPSTN/ X,B,W,ICPT
C
      B = 0.4
      Z = 0.46265167
C
C Initial seed from the system time and and forced to be odd
C
      STIME = time(%REF(0))
      CALL gmtime(STIME,T)
      ISEED = T(6)+70*(T(5)+12*(T(4)
     *            +31*(T(3)+23*(T(2)+59*T(1)))))
      IF (MOD(ISEED,2).EQ.0) ISEED = ISEED-1
C
C Runs to remove the initial configuration dependence
C
      X = RANF()
      W = WFNT(X)
      DO 100 J = 1, ITHM
        CALL MTSP
  100 CONTINUE
C
C Collect data at every L steps
C
      SUM1 = 0.0
      SUM2 = 0.0
      ICPT = 0
      DO 200 J = 1, ISIZ*L
        CALL MTSP
        IF(MOD(J,L).EQ.0) THEN
          SUM1 = SUM1+G(X)
          SUM2 = SUM2+G(X)*G(X)
        ELSE
        ENDIF
  200 CONTINUE
C
      S  = Z*SUM1/ISIZ
      DS = Z*SQRT(ABS(SUM2/ISIZ-(SUM1/ISIZ)**2)/ISIZ)
      ACPT = 100.0*ICPT/(ISIZ*L)
C
      WRITE(6,999) S,DS,ACPT
  999 FORMAT (3F14.8)
      STOP
      END
C
      SUBROUTINE MTSP
C
C  Subroutine for one Metropolis step.
C
      PARAMETER (ITHM=10000,ISIZ=10000,L=15)
      COMMON /CSEED/ ISEED
      COMMON /CPSTN/ X,B,W,ICPT
C
      XSAV = X
C
      RNDX = RANF()
      X    = X + 2.0*B*(RNDX-0.5)
      IF ((X.LT.0).OR.(X.GT.1)) THEN
        X = XSAV
        RETURN
      ELSE
      ENDIF
C
      WTRY = WFNT(X)
      RNDW = RANF()
      IF (WTRY.GT.(W*RNDW)) THEN
        W    = WTRY
        ICPT = ICPT+1
        RETURN
      ELSE
      ENDIF
      X    = XSAV
      RETURN
      END
C
      FUNCTION G(X)
        G = X*X/(EXP(X*X)-1.0)
      RETURN
      END
C
      FUNCTION WFNT(X)
        WFNT = EXP(X*X)-1.0
      RETURN
      END
C
      FUNCTION RANF()
      DATA IA/16807/,IC/2147483647/,IQ/127773/,IR/2836/
      COMMON /CSEED/ ISEED
        IH = ISEED/IQ
        IL = MOD(ISEED,IQ)
        IT = IA*IL-IR*IH
        IF(IT.GT.0) THEN
          ISEED = IT
        ELSE
          ISEED = IC+IT
        END IF
        RANF = ISEED/FLOAT(IC)
      RETURN
      END