!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Program file name: carlo.f90                                           !
!                                                                         !
!  © Tao Pang 2006                                                        !
!                                                                         !
!  Last modified: November 13, 2009                                       !
!                                                                         !
!  (1) This F90 program is created for the book, "An Introduction to      !
!      Computational Physics, 2nd Edition," written by Tao Pang and       !
!      published by Cambridge University Press on January 19, 2006.       !
!                                                                         !
!  (2) No warranties, express or implied, are made for this program.      !
!                                                                         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
MODULE CSEED
  INTEGER :: SEED
END MODULE CSEED
!
MODULE CB
  INTEGER :: IACCEPT
  REAL :: X, H, W
END MODULE CB
!
PROGRAM CARLO
!
! An example of Monte Carlo simulation with the Metropolis scheme
! with integrand f(x) = x*x.
!
  USE CSEED
  USE CB
  IMPLICIT NONE
  INTEGER, PARAMETER :: NSIZE = 10000, NSKIP = 15, NTOTAL = NSIZE*NSKIP, &
                        NEQ = 10000
  INTEGER :: I
  INTEGER, DIMENSION (8) :: T
  REAL :: S0, DS, Z, ACCEPT, F, G, RANF, WEIGHT
!
  H = 0.4
  Z = 0.46265167
!
! Initiate the seed from the current date and time
!
  CALL DATE_AND_TIME(VALUES = T)
  SEED = T(1)+70*(T(2)+12*(T(3)+31*(T(5)+23*(T(6)+59*T(7)))))
  IF (MOD(SEED,2).EQ.0) SEED = SEED-1
!
  X = RANF()
  W = WEIGHT(X)
  IACCEPT = 0
  DO I = 1, NEQ
    CALL METROPOLIS
  END DO
!
  S0 = 0.0
  DS = 0.0
  IACCEPT = 0
  DO I = 1, NTOTAL
    CALL METROPOLIS
    IF (MOD(I,NSKIP).EQ.0) THEN
      F = G(X)
      S0 = S0 + F
      DS = DS + F*F
    ELSE
    END IF
  END DO
!
  S0 = S0/NSIZE
  DS = DS/NSIZE
  DS = SQRT(ABS(DS-S0*S0)/NSIZE)
  S0 = Z*S0
  DS = Z*DS
  ACCEPT = 100*IACCEPT/NTOTAL
!
  WRITE(6, *) "S = ",S0, " +- ", DS
  WRITE(6, *) "The accept rate = ", ACCEPT
END PROGRAM CARLO
!
SUBROUTINE METROPOLIS
!
! Subroutine for one Metropolis step.
!
  USE CB
  IMPLICIT NONE
  REAL :: XOLD, WNEW, RANF, WEIGHT
!
  XOLD = X
!
  X = X + 2.0*H*(RANF()-0.5)
  IF ((X.LT.0).OR.(X.GT.1)) THEN
    X = XOLD
  ELSE
    WNEW = WEIGHT(X)
    IF (WNEW.GT.(W*RANF())) THEN
      W = WNEW
      IACCEPT = IACCEPT + 1
    ELSE
      X = XOLD
    END IF
  END IF
END SUBROUTINE METROPOLIS
!
FUNCTION G(X) RESULT (GX)
  IMPLICIT NONE
  REAL :: X, GX
  GX = X*X/(EXP(X*X)-1.0)
END FUNCTION G
!
FUNCTION WEIGHT(X) RESULT (WX)
  IMPLICIT NONE
  REAL :: X, WX
        WX = EXP(X*X)-1.0
END FUNCTION WEIGHT
!
FUNCTION RANF() RESULT (CR)
!
! Function to generate a uniform random number in [0,1]
! following x(i+1)=a*x(i) mod c with a=7** 5 and
! c=2**31-1.  Here the seed is a global variable.
!
  USE CSEED
  IMPLICIT NONE
  INTEGER :: H, L, T, A, C, Q, R
  DATA A/16807/, C/2147483647/, Q/127773/, R/2836/
  REAL :: CR
!
  H = SEED/Q
  L = MOD(SEED, Q)
  T = A*L - R*H
  IF (T .GT. 0) THEN
    SEED = T
  ELSE
    SEED = C + T
  END IF
  CR = SEED/FLOAT(C)
END FUNCTION RANF