	parameter (n=10000)
	dimension r(n),ind(n)
	kkk=0
	pi=4.*atan(1.)
	do i=1,n
	x=2.*ran(kkk)-1.
	y=2.*ran(kkk)-1.
	z=2.*ran(kkk)-1.
	r(i)=sqrt(x**2+y**2+z**2)
	enddo
	call indexx(n,r,ind)

	rhob=float(n)/(2.**3)
	open(1,file='rand_walk.dat')
	do i=1,n
	if(r(ind(i)).ge.1)goto10
	vol=4./3.*pi*r(ind(i))**3
	rho=float(i)/vol
	delta=rho/rhob-1.
	write(1,*)r(ind(i)),delta
	enddo	
 10	close(1)
	end


      SUBROUTINE INDEXX(N,ARRIN,INDX)
      DIMENSION ARRIN(N),INDX(N)
      DO 11 J=1,N
        INDX(J)=J
11    CONTINUE
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          INDXT=INDX(L)
          Q=ARRIN(INDXT)
        ELSE
          INDXT=INDX(IR)
          Q=ARRIN(INDXT)
          INDX(IR)=INDX(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            INDX(1)=INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
          ENDIF
          IF(Q.LT.ARRIN(INDX(J)))THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        INDX(I)=INDXT
      GO TO 10
      END
