program bordel
      implicit none
      real*8 frandom,x,y,z,d,half,dx,dy,dz,s(100),ave,a2,dd,h,p,pi,o
      integer*4 i,ix,iy,iz,nat,iq,j,n
      character*80 s80
      common /randc/ ix, iy, iz
      ix=1
      iy=1000
      iz=99
      pi=4.0d0*datan(1.0d0)
      half=0.5d0
      s=0.0d0
      ave=0.0d0
      a2=0.0d0
      if(iargc().ne.1)then
       write(6,*)'Usage: bordel <d>'
       stop
      endif
      call getarg(1,s80)
      read(s80,*)d
      dd=d/100.0d0
      o=0.2d0*d
      open(9 ,file='FILE.X')
      open(91,file='B.X')
      read(9,90)s80
 90   format(a80)
      write(91,90)s80
      read(9,*)nat
      write(91,*)nat
      n=0
      do 1 i=1,nat
      read(9,*)iq,x,y,z
      n=n+3
      h=(dble(n)/d*dd)/d/dsqrt(pi)

100   dx=d*(half-frandom(i))
      j=nint((dx+half*d)/dd)
      if(s(j).gt.h*exp(-(dx/o)**2))goto 100
      if(j.gt.0.and.j.lt.100)s(j)=s(j)+1.0d0/d

101   dy=d*(half-frandom(i))
      j=nint((dy+half*d)/dd)
      if(s(j).gt.h*exp(-(dy/o)**2))goto 101
      if(j.gt.0.and.j.lt.100)s(j)=s(j)+1.0d0/d

102   dz=d*(half-frandom(i))
      j=nint((dz+half*d)/dd)
      if(s(j).gt.h*exp(-(dz/o)**2)) goto 102
      if(j.gt.0.and.j.lt.100)s(j)=s(j)+1.0d0/d

      ave=ave+dx+dy+dz
      a2=a2+dx**2+dy**2+dz**2
1     write(91,91)iq,x+dx,y+dy,z+dz
 91   format(i6,3f12.6,' 0 0 0 0 0 0 0 0.0')
      close(9)
      close(91)
      open(9,file='distr.prn')
      do 3 i=1,100
  3   write(9,92)-d*half+dble(i)*dd,s(i)
  92  format(f10.4,f12.2)
      close(9)
      ave=ave/dble(3*nat)
      a2 =a2 /dble(3*nat)
      write(6,60)ave,dsqrt(a2-ave**2)
  60  format(' Average: ',f12.6,', RMS:',f12.6)
      end

      function frandom(i)
      implicit none
      real*8 frandom
      integer*4 i
c
c     Algorithm AS 183 Appl. Statist. (1982) vol.31, no.2
c
c     Returns a pseudo-random number rectangularly distributed
c     between 0 and 1.   The cycle length is 6.95E+12 (See page 123
c     of Applied Statistics (1984) vol.33), not as claimed in the
c     original article.
c
c     IX, IY and IZ should be set to integer values between 1 and
c     30000 before the first entry.
c
      integer*4 ix, iy, iz
      common /randc/ ix, iy, iz
c
      ix = mod(171 * ix, 30269)
      iy = mod(172 * iy, 30307)
      iz = mod(170 * iz, 30323)
      frandom = mod(dble(ix) / 30269.0d0 + dble(iy) / 30307.0d0 +
     1                        dble(iz) / 30323.0d0, 1.0d0)
      return
      end