        program baby_mc
c
c	Monte Carlo for a trivial problem - 10 states with energy E(i)=i
c	so p(i)=wxp(- beta * E(i)) / Z   where  Z=sum exp(- beta * E(i))
c
        dimension n(10)
c
c       warm up random number generator
c
        idum=193728181
        do i=1,1000
                x=ran1(idum)
        enddo
c
c       input parameters
c
	print*, 'enter beta and number of cycles'
        read*, beta,ncycles
c
c       initial state
c
        i=1+10*ran1(idum)

        e=0
        n=0
        do icycle=1,ncycles
                j=1+10*ran1(idum)       ! new state
                if(j.lt.i) then
                        i=j             ! accept if lower energy
                else
                        z=ran1(idum)    ! accept or not by metropolis test
                        if(z.lt.exp(-beta*(j-i))) i=j
                endif
                e=e+i

                n(i)=n(i)+1
        enddo
c
c       output
c
        print*, e/float(ncycles)
        print '(10f8.5)', n/float(ncycles)

        stop
        end

1
      FUNCTION ran1(idum)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      REAL ran1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER j,k,iv(NTAB),iy
      SAVE iv,iy
      DATA iv /NTAB*0/, iy /0/
      if (idum.le.0.or.iy.eq.0) then
        idum=max(-idum,1)
        do 11 j=NTAB+8,1,-1
          k=idum/IQ
          idum=IA*(idum-k*IQ)-IR*k
          if (idum.lt.0) idum=idum+IM
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM
      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)
      return
      END

