implicit none
      integer i,j,d,L,MCS,ii,N,k
      parameter(L=35,N=850,MCS=23000)
      real x(N),y(N),Nr(100000)
      real r,ran1,jj,q,w,T,xn,yn,dx,dy,dr2,U0,Ut,dU,ro,dr,r1,licznik

      d=-1
      open(1,file='zad2.txt')

      jj=-0.5
      do i=1,N
       x(i)=mod(i-1,L)+0.5
       if(mod(i-1,L).eq.0) then
        jj=jj+1
       endif
       y(i)=jj
      enddo

      T=0.7
      ro=N/(L*L*1.0)
      Nr=0

      do j=1,MCS

       do i=1,N
        U0=0
        Ut=0

          xn=x(i)+(ran1(d)-0.5)*0.1
          if(xn.ge.L) xn=xn-L
          if(xn.lt.0) xn=xn+L
          yn=y(i)+(ran1(d)-0.5)*0.1
          if(yn.ge.L) yn=yn-L
          if(yn.lt.0) yn=yn+L

          do ii=1,N
           if(ii.ne.i) then

             dx=xn-x(ii)
             if(dx.ge.(0.5*L)) dx=L-dx
             dy=yn-y(ii)
             if(dy.ge.(0.5*L)) dy=L-dy
             dr2=dx*dx+dy*dy

             if(sqrt(dr2).le.2.5) then
              Ut=Ut+4*(dr2**(-6)-dr2**(-3))
             endif

             dx=x(i)-x(ii)
             if(dx.ge.(0.5*L)) dx=L-dx
             dy=y(i)-y(ii)
             if(dy.ge.(0.5*L)) dy=L-dy
             dr2=dx*dx+dy*dy

             if(sqrt(dr2).le.2.5) then
              U0=U0+4*(dr2**(-6)-dr2**(-3))
             endif

             if(j.gt.30) then
              if(mod(j,10).eq.0) then
               licznik=licznik+1
                do k=1, 1000
                  dr=k*0.01
                  r1=sqrt(dr2)
                  if(r1.le.dr) then
                    if(r1.gt.(dr-0.01)) then
                     Nr(k)=Nr(k)+1
                    endif
                  endif
                enddo
               endif
             endif

            endif
           enddo

          dU=Ut-U0
           w=min(1.0,exp(-dU/T))
           q=ran1(d)
            if(q.le.w) then
              x(i)=xn
              y(i)=yn
            endif

       enddo

      enddo

      do k=1,1000
        Nr(k)=Nr(k)/(2.0*3.1415*k*0.01*0.01*N*200*ro)
        write(1,*) 0.01*k, Nr(k)
      enddo

      do i=1, N
       write(1,*) x(i), y(i)
      enddo

      write (*,*) 'end'
      close(1)

      read(*,*)
      end

      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