!    Integracija Monte Karlo metodom
!     Program
  Program metropolis
  implicit none
!
	 real*8, parameter::a=-2., b=2.
  	  real*8 x,delta,sf,sd,x0,xp ! x- vrednost nezavisno promenljive iz nekog intervala [a,b], y-niz vrednosti podintegralne f-je
      real*8 ran2,r,y(50000),f ,eta,d ! korak integracije
      integer i,idum,neprihvati,n,k ! inicijator generatora sluchajnih brojeva
! Otvaranje fajla (datoteke) u koji ce se smestiti rezultati izrachunavanja
    open(30,file='metro.dat',status='unknown')
!



x0=0.
delta=4.
d=0.2
sf=0.
sd=0.
neprihvati=0
n=0

idum=-123456


do while(neprihvati.lt.4000)
n=n+1
x=ran2(idum)
xp=x0+delta*(2.*x-1.)
r=f(xp)/f(x0)

if(r.lt.1)then
   eta=ran2(idum)
   if(r.lt.eta)then
     neprihvati=neprihvati+1
   else
     x0=xp
     sf=sf+f(xp)
     sd=sd+f(xp)*f(xp)
   endif
else
x0=xp
sf=sf+f(xp)
sd=sd+f(xp)*f(xp)
endif

y(n)=xp

enddo
sf=sf/(n-neprihvati)
sd=sd/(n-neprihvati)


print*,float(n),float(neprihvati),(b-a)*sf,(b-a)*dsqrt(sd-sf*sf)/dsqrt(dfloat(neprihvati)),sqrt(6.284)

do i=1,n
write(30,100)(b-a)/n*i ,y(i)
enddo

100 format(10(e13.7,2x))

stop


end

function f(x)
real*8 f,x
f=dexp(-0.5*x*x)
end

!     *****************************RAN2****************************
      function ran2(idum)
      integer*4 idum,IM1,IM2,INM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      real*8 ran2,AM,EPS,RNMX
      parameter(IM1=2147483563,IM2=2147483399,AM=1.d0/IM1,INM1=IM1-1,IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,NTAB=32,NDIV=1+INM1/NTAB,EPS=1.2d-7,RNMX=1.-EPS)
!     Long period (>2x10^{18}) random number generator of L'Ecuyer
!     with Bays-Durham shuffle and added safeguards. Returns a uniform
!     random deviate between 0.0 and 1.0 (exclusive of the endpoint
!     values). Call with idum a negative integer to initialise;
!     thereafter, do not alter idum between successive deviate in a
!     sequence. RNMX should approximate the largest floating value
!     that is less than 1.
      integer*4 idum2,j,k,iv(NTAB),iy
      save iv,iy,idum2
      data idum2/123456789/,iv/NTAB*0/,iy/0/

      if(idum .le. 0) then             !Initialise
         idum = max(-idum,1)           !Be sure to prevent idum=0
         idum2 = idum
         do 100, j=NTAB+8,1,-1         !Load the shuffle table (after
            k=idum/IQ1                 ! 8 warm-ups)
            idum = IA1*(idum-k*IQ1)-k*IR1
            if(idum .lt. 0) idum = idum + IM1
            if(j .le. NTAB) iv(j) = idum
 100     continue
         iy = iv(1)
      end if
      k = idum/IQ1                     !Start here when not initialising
      idum = IA1*(idum-k*IQ1)-k*IR1    !Compute idum=mod(IA1*idum,IM1)
                                       !without overflow by Schrage's
                                       !method
      if(idum .lt. 0) idum = idum + IM1
      k = idum2/IQ2
      idum2 = IA2*(idum2-k*IQ2)-k*IR2  !Compute idum2=mod(IA2*idum2,IM2)
                                       !likewise
      if(idum2 .lt. 0) idum2 = idum2 + IM2
      j = 1 + iy/NDIV                  !Will be in the range 1:NTAB
      iy = iv(j) - idum2               !Here idum is shuffled, idum and
      iv(j) = idum                     !idum2 are combined to generate
                                       !output
      if(iy .lt. 1) iy = iy + INM1
      ran2 = min(AM*iy,RNMX)           !Because users do not expect
                                       !endpoint values
      return
      end
!     ******************************************************************






