********************************************************************************
      PROGRAM SIMPLE_MC
********************************************************************************
*      Simple example of using Monte Carlo techniques                          *
*      All it does is calculate the value of pi, in a really dumb (but         *
*      therefore easy to follow) way.                                          *
********************************************************************************
      IMPLICIT NONE
      CHARACTER*30 sys_call,N_txt
      REAL*8 myrndm,x,y,r,A,l,A_square,dA,N,N_accept,dN_accept,N_counter
      INTEGER*4 IJKL,IJKL_in,i
      INTEGER istat
*------------------------------------------------------------------------------*
*      by default, initial the random seed based on the UTC time
      sys_call = 'date +%s | cut -c3- > tmp.tmp' ! is an 8-digit integer saved 
                                ! to a new file "tmp.tmp"
      CALL SYSTEM(sys_call,istat) ! execute the shell command
      OPEN(UNIT=99,FILE='tmp.tmp',STATUS='OLD') ! open the tmp.tmp file
      READ(99,*)IJKL            ! read in the 8-digit integer
      CLOSE(99)
      sys_call = 'rm -f tmp.tmp' ! remove the temporary file
      CALL SYSTEM(sys_call,istat)
*      check if user wants to type one in; <cr> will accept the time-based one
      WRITE(*,'(a,i8,a,$)')'Enter a random 8-digit integer:  [<cr> = ',
     &  IJKL,']  '
      READ(*,'(i8)')IJKL_in
      IF( (IJKL_in.GT.0)        ! user didn't just hit <cr>
     &  .AND. (IJKL_in.LT.1000000) )THEN ! it isn't too big a number
        IJKL = IJKL_in          ! set the seed to the user input
      ELSE IF(IJKL_in.GT.0)THEN ! number is either negative or too big
        WRITE(*,*)'Number must be 8 digits! You gave IJKL = "',
     &    IJKL_in,'"'
        STOP
      END IF
*      now that we have a seed, initialize the subroutine which will give 
*      us our pseudo-random numbers
      CALL RM48IN(IJKL,0,0)     ! from CERN.  The 2nd and 3rd values can be 
                                !   used to resume a series of psuedo-random #s
*      define a square of arbitrary length, L
      L = 2.d0                  ! let's make it 2 so that the biggest circle 
                                !   that fits has a radius of 1
*      initialize counters
      N_accept  = 0.d0          ! # of accepted events
      N_counter = 2.d0          ! # for frequency of displaying progress
      N = 1.e8                  ! # of events to simulate

*      again, check if user doesn't want to go with the defaults
      WRITE(*,'(a,1pe6.0,a,$)')'Enter the number of events to throw:'//
     &  '  [<cr> = ',N,']  '
      READ(*,'(a)')N_txt        ! a silly trick to enable <cr> responses

      IF(N_txt.NE.' ')THEN      ! user didn't just hit <cr>
        READ(N_txt,*)N          ! read the # to simulate from the text variable
      END IF

*      set up the table of the output to the screen
      WRITE(*,*)
      WRITE(*,'(a)')'       N           pi = A/r^2'
      WRITE(*,'(a)')'----------------------------------'

      A_square = l*l            ! area of the square

*      save the results to a file as we go along
      OPEN(UNIT=1,FILE='pi-calculated.asc',STATUS='UNKNOWN')
      DO i = 1,INT(N)+1
        x = 1.d0 - 2.d0*myrndm() ! x and y are independent random variables 
        y = 1.d0 - 2.d0*myrndm() !   spanning [-1,1]; uniformly fills a square

        r = sqrt(x*x + y*y)     ! calculate the radius

        IF (r.LE.1.d0) THEN     ! if within a circle of radius 1
          N_accept =  N_accept+1.d0 ! "pass" -- accept the event
        END IF

*        stats go like 1/sqrt(N), so best to show info logarithmically with 
*        the number of trials
        IF (log10(REAL(i))*2.d0.GT.REAL(N_counter)) THEN
          N_counter = N_counter+1.d0
          dN_accept = sqrt(REAL(n_accept))

          A  = A_square*N_accept/REAL(i) ! this is the area of the events that 
                                ! landed within the circle.  We know the area 
                                ! of the square, so it's just that times the 
                                ! ratio of accepted to total events
          dA = A*abs(dN_accept/N_accept) ! given that i events have been 
                                ! simulated, the uncertainty is only in how 
                                ! many passed instead of failed.  So dA/A = 
                                ! (dN/N)_accepted

          WRITE(1,'(i12,2(1pe14.6))')i,A,dA ! write out the cumulative result
          WRITE(*,'(1pe12.4,0pf10.4,a,f6.4)')REAL(i),
     &      A,'+/-',dA          ! also display it on the screen in the table
        END IF
      END DO
*      we're done looping over all events
*      make sure the area calculation and it's uncertainty are updated
      A  = A_square*N_accept/REAL(i)
      dA = A/sqrt(N_accept)
*      print result to screen
      WRITE(*,'(a)')'----------------------------------'
      WRITE(*,'(a,f6.4,a,f6.4)')'    =>   pi = ',A,'+/-',dA
      WRITE(*,'(a)')'----------------------------------'
*      and also update the output file
      WRITE(1,'(i12,2(1pe14.6))')i,A,dA
      CLOSE(1)

      STOP
      END

********************************************************************************
*
* Id: rm48.F,v 1.2 1996/12/12 16:32:06 cernlib Exp 
*
* Log: rm48.F,v 
* Revision 1.2  1996/12/12 16:32:06  cernlib
* Variables ONE and ZERO added to SAVE statement, courtesy R.Veenhof
*
* Revision 1.1.1.1  1996/04/01 15:02:55  mclareni
* Mathlib gen
*
*
      SUBROUTINE RM48(RVEC,LENV)
C     Double-precision version of
C Universal random number generator proposed by Marsaglia and Zaman
C in report FSU-SCRI-87-50
C        based on RANMAR, modified by F. James, to generate vectors
C        of pseudorandom numbers RVEC of length LENV, where the numbers
C        in RVEC are numbers with at least 48-bit mantissas.
C   Input and output entry points: RM48IN, RM48UT.
C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C!!!  Calling sequences for RM48:                                    ++
C!!!      CALL RM48 (RVEC, LEN)     returns a vector RVEC of LEN     ++
C!!!                   64-bit random floating point numbers between  ++
C!!!                   zero and one.                                 ++
C!!!      CALL RM48IN(I1,N1,N2)   initializes the generator from one ++
C!!!                   64-bit integer I1, and number counts N1,N2    ++
C!!!                  (for initializing, set N1=N2=0, but to restart ++
C!!!                    a previously generated sequence, use values  ++ 
C!!!                    output by RM48UT)                            ++ 
C!!!      CALL RM48UT(I1,N1,N2)   outputs the value of the original  ++
C!!!                  seed and the two number counts, to be used     ++
C!!!                  for restarting by initializing to I1 and       ++  
C!!!                  skipping N2*100000000+N1 numbers.              ++
C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C for 32-bit machines, use IMPLICIT DOUBLE PRECISION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER*4 (I-N)
      DIMENSION RVEC(*)
      COMMON/R48ST1/U(97),C,I97,J97
      PARAMETER (MODCNS=1000000000)
      SAVE CD, CM, TWOM24, NTOT, NTOT2, IJKL,TWOM49, ONE, ZERO
      DATA NTOT,NTOT2,IJKL/-1,0,0/
C
      IF (NTOT .GE. 0)  GO TO 50
C
C        Default initialization. User has called RM48 without RM48IN.
      IJKL = 54217137
      NTOT = 0
      NTOT2 = 0
      KALLED = 0
      GO TO 1
C
      ENTRY      RM48IN(IJKLIN, NTOTIN,NTOT2N)
C         Initializing routine for RM48, may be called before
C         generating pseudorandom numbers with RM48.   The input
C         values should be in the ranges:  0<=IJKLIN<=900 OOO OOO
C                                          0<=NTOTIN<=999 999 999
C                                          0<=NTOT2N<<999 999 999!
C To get the standard values in Marsaglia's paper, IJKLIN=54217137
C                                            NTOTIN,NTOT2N=0
      IJKL = IJKLIN
      NTOT = MAX(NTOTIN,0)
      NTOT2= MAX(NTOT2N,0)
      KALLED = 1
C          always come here to initialize
    1 CONTINUE
      IJ = IJKL/30082
      KL = IJKL - 30082*IJ
      I = MOD(IJ/177, 177) + 2
      J = MOD(IJ, 177)     + 2
      K = MOD(KL/169, 178) + 1
      L = MOD(KL, 169)
c      WRITE(*,'(A,I10,2X,2I10)') ' RM48 INITIALIZED:',IJKL,NTOT,NTOT2
CCC      PRINT '(A,4I10)', '   I,J,K,L= ',I,J,K,L
      ONE = 1.d0
      HALF = 0.5d0
      ZERO = 0.d0
      DO 2 II= 1, 97
      S = 0.d0
      T = HALF
      DO 3 JJ= 1, 48
         M = MOD(MOD(I*J,179)*K, 179)
         I = J
         J = K
         K = M
         L = MOD(53*L+1, 169)
         IF (MOD(L*M,64) .GE. 32)  S = S+T
    3    T = HALF*T
    2 U(II) = S
      TWOM49 = T
      TWOM24 = ONE
      DO 4 I24= 1, 24
    4 TWOM24 = HALF*TWOM24
      C  =   362436.d0*TWOM24
      CD =  7654321.d0*TWOM24
      CM = 16777213.d0*TWOM24
      I97 = 97
      J97 = 33
C       Complete initialization by skipping
C            (NTOT2*MODCNS + NTOT) random numbers
      DO 45 LOOP2= 1, NTOT2+1
      NOW = MODCNS
      IF (LOOP2 .EQ. NTOT2+1)  NOW=NTOT
      IF (NOW .GT. 0)  THEN
c      WRITE(6,'(A,I15)') ' RM48IN SKIPPING OVER ',NOW
          DO 40 IDUM = 1, NTOT
          UNI = U(I97)-U(J97)
          IF (UNI .LT. ZERO)  UNI=UNI+ONE
          U(I97) = UNI
          I97 = I97-1
          IF (I97 .EQ. 0)  I97=97
          J97 = J97-1
          IF (J97 .EQ. 0)  J97=97
          C = C - CD
          IF (C .LT. ZERO)  C=C+CM
   40     CONTINUE
      ENDIF
   45 CONTINUE
      IF (KALLED .EQ. 1)  RETURN
C
C          Normal entry to generate LENV random numbers
   50 CONTINUE
      DO 100 IVEC= 1, LENV
      UNI = U(I97)-U(J97)
      IF (UNI .LT. ZERO)  UNI=UNI+ONE
      U(I97) = UNI
      I97 = I97-1
      IF (I97 .EQ. 0)  I97=97
      J97 = J97-1
      IF (J97 .EQ. 0)  J97=97
      C = C - CD
      IF (C .LT. ZERO)  C=C+CM
      UNI = UNI-C
      IF (UNI .LT. ZERO) UNI=UNI+ONE
      RVEC(IVEC) = UNI
C             Replace exact zeros by 2**-49
         IF (UNI .EQ. ZERO)  THEN
            RVEC(IVEC) = TWOM49
         ENDIF
  100 CONTINUE
      NTOT = NTOT + LENV
         IF (NTOT .GE. MODCNS)  THEN
         NTOT2 = NTOT2 + 1
         NTOT = NTOT - MODCNS
         ENDIF
      RETURN
C           Entry to output current status
      ENTRY RM48UT(IJKLUT,NTOTUT,NTOT2T)
      IJKLUT = IJKL
      NTOTUT = NTOT
      NTOT2T = NTOT2
      RETURN
      END





***************************************************************************
      REAL*8 FUNCTION RANGAUSS(AVER,RMS)
***************************************************************************
*     returns normally distributed pseudo-random real number with average
*     value = averag, and  dispersion = rms
*************************************************************
      IMPLICIT NONE
      REAL*8 rms,x1,x2,aver,v1,v2,s,tmp,myrndm
      LOGICAL called            ! number of previous calls
      DATA called/.false./      ! toggle returned value
      
      IF(called) THEN           ! called once with these random pair
        RANGAUSS= rms*x2+aver   ! so, use second number, don't loose it
        called = .false.        ! reset toggle
      ELSE                      ! use 2 new random numbers
 10     v1 = 2.d0*myrndm()-1.d0
        v2 = 2.d0*myrndm()-1.d0
        s  = v1*v1+v2*v2
        IF(s.GT.1.d0)  GOTO 10
        IF(s.EQ.0.d0)  GOTO 10
        tmp = sqrt((-2.d0)*LOG(s)/s)
        x1  = v1*tmp
        x2  = v2*tmp
        RANGAUSS = rms*x1+aver
        called = .true.
      ENDIF
      
      RETURN
      END
      


***************************************************************************
      REAL*8 FUNCTION myrndm()
***************************************************************************
      IMPLICIT NONE
      REAL*8 rndm(1)
      CALL RM48(rndm,1)
      myrndm = rndm(1)
      RETURN
      END