pro show_ngp
  window,xsize=1000,ysize=1000

  readnew,'snap_144',h,'HEAD'

  ngrid=20
  vel_grid=fltarr(3,ngrid,ngrid,ngrid)
  n_grid=fltarr(ngrid,ngrid,ngrid)
  
  readnew,'snap_144',x,'POS',part=0
  readnew,'snap_144',v,'VEL',part=0

  FOR i=0L,h.npart[0]-1 DO BEGIN
     ix=FLOOR(x[0,i]/h.BOXSIZE*ngrid)
     iy=FLOOR(x[1,i]/h.BOXSIZE*ngrid)
     iz=FLOOR(x[2,i]/h.BOXSIZE*ngrid)

     vel_grid[*,ix,iy,iz] += v[*,i]
     n_grid[ix,iy,iz] += 1.0    
  END
  FOR k=0,2 DO vel_grid[k,*,*,*] = vel_grid[k,*,*,*]/n_grid[*,*,*]

  jj=WHERE(ABS(x[2,*]-h.BOXSIZE/2) LT 2*h.BOXSIZE/ngrid)
  plot,[1],[1],xst=1,yst=1,xr=[0,h.BOXSIZE],yr=[0,h.BOXSIZE],/nodata
  oplot,x[0,jj],x[1,jj],psym=1,col=64,symsize=0.5

  seed=1234L
  Nlines = 10
  ind_lines = ROUND(randomu(seed,Nlines)*h.npart[0])

  NSTEPS = 40
  dt=5.0
  FOR ix=0,ngrid-1 DO BEGIN
     FOR iy=0,ngrid-1 DO BEGIN
;        FOR iz=0,ngrid-1 DO BEGIN
        FOR iz=ngrid/2-1,ngrid/2-1 DO BEGIN
           x0=([ix,iy,iz]+0.5)*h.BOXSIZE/ngrid
           oplot,[x0[0]],[x0[1]],psym=4,col=255
           FOR k=0L,NSTEPS DO BEGIN
              print,x0

              ix0=FLOOR(x0[0]/h.BOXSIZE*ngrid)
              iy0=FLOOR(x0[1]/h.BOXSIZE*ngrid)
              iz0=FLOOR(x0[2]/h.BOXSIZE*ngrid)

              x1 = x0 + vel_grid[*,ix0,iy0,iz0] * dt
              
              oplot,[x0[0],x1[0]],[x0[1],x1[1]],col=255

              x0 = x1
           END
        END
     END   
  END
  stop
end