program NBODY

	integer MM,NN,PX,PY,FX,FY
	parameter (NN=1024)
	parameter (MM=0, PX=1, PY=2, PZ=3, FX=4, FY=5, FZ=6)
C	MM:mass, PX:position in x coordinate, PY:position in y coordinate
C	FX: force in x coordinate, FY: force in y coordinate
	real dx(0:NN-1), dy(0:NN-1), dz(0:NN-1), dist(0:NN-1), sq(0:NN-1)
        real fac(0:NN-1), tx(0:NN-1), ty(0:NN-1), tz(0:NN-1)
	real p(0:6,0:NN-1), q(0:6,0:NN-1) 
	integer i,j,npts

	print *, 'Enter # of data points: '
	read  *, npts
	seed = (thisNode+1)*100000000
	do i = 0, npts-1
C	    p(MM,i) = SRANDOM()
C	    p(PX,i) = SRANDOM()
C	    p(PY,i) = SRANDOM()
	    p(MM,i) = 1.05 + i
	    p(PX,i) = 1000001 - i
	    p(PY,i) = 999999 + i
	    p(FX,i) = 0.0
	    p(FY,i) = 0.0
	end do

	do i = 0, npts-1
	    q(MM,i) = p(MM,i)
	    q(PX,i) = p(PX,i)
	    q(PY,i) = p(PY,i)
	    q(PZ,i) = p(PZ,i)
	    q(FX,i) = 0.0
	    q(FY,i) = 0.0
	    q(FZ,i) = 0.0
	end do
	do i = 0, npts-1
	  do j = i+1, npts-1
  	    dx(i) = p(PX,i) - q(PX,j)
	    dy(i) = p(PY,i) - q(PY,j)
	    dz(i) = p(PZ,i) - q(PZ,j)
	    sq(i) = dx(i)**2+dy(i)**2+dz(i)**2
	    dist(i) = SQRT(sq(i))
	    fac(i) = p(MM,i) * q(MM,j) / (dist(i) * sq(i))
	    tx(i) = fac(i) * dx(i)
	    ty(i) = fac(i) * dy(i)
	    tz(i) = fac(i) * dz(i)
	    p(FX,i) = p(FX,i)-tx(i)
            q(FX,j) = q(FX,j)+tx(i)
	    p(FY,i) = p(FY,i)-ty(i)
	    q(FY,j) = q(FY,j)+ty(i)
	    p(FZ,i) = p(FZ,i)-tz(i)
            q(FZ,j) = q(FZ,j)+tz(i)
	  end do
	end do
CP	print *, 'PE',thisNode,' -1',(p(FX,i),i=0,npts-1)
CP	print *, 'PEq',thisNode,' -1',(q(FX,i),i=0,npts-1)

	do i = 0, npts-1
	  p(FX,i) = p(FX,i) + q(FX,i)
	  p(FY,i) = p(FY,i) + q(FY,i)
	  p(FZ,i) = p(FZ,i) + q(FZ,i)
	end do

CP	print *, 'PE',thisNode,' p=',(p(FX,i),i=0,npts-1)
CP	print *, 'PE',thisNode,' q=',(q(FX,i),i=0,npts-1)
	print *, ' p=',(p(FX,i),i=0,npts-1)
	print *, ' q=',(q(FX,i),i=0,npts-1)
	stop
	end