c **** parameters for g(r) histogram bin dr and maximum distance to go
      parameter(dr=0.2,rmax=4.0,nbins=rmax/dr)
      dimension gr(nbins) 


c **** set up the array and normalization constant (2D !!!)
      const = 3.14159 * Npart /bx/by

      do ibin=1,nbins
       gr(ibin)=0
      enddo

c **** This goes Inside the MC loop, where all running  averages are calculated
c **** Loop through pairs

       do i = 1,Npart-1
       do j = i+1,Npart
        dx = x(1,i)-x(1,j)
        dy = x(2,i)-x(2,j)
c **** Here are the Periodic Boundary Conditions
        dx=dx - bx * anint(dx/bx)
        dy=dy - by * anint(dy/by)
        r2ij=dx*dx+dy*dy
        rij=sqrt(r2ij)
        ibin= int(rij/dr)+1
        if (ibin . le. nbins) then
           gr(ibin) = gr(ibin)+2
        endif
       enddo
       enddo

c **** This goes to the very end, after the  MC run is completed
c **** Normalize it, assuming gr has been accumulated after each 
c **** MC cycle, and it is 2D

       do ibin=1,nbins
        rl = (ibin-1)*dr
        ru = rl + dr
        pideal = const * (ru*ru - rl*rl)
        gr(ibin) = gr(ibin) / icycle / Npart / pideal 
       enddo

c **** Here we can put it's output......

