pro phase_space

!P.MULTI=[0,2,1]
  !X.MARGIN=[6.5,0.5]
   readmatrix,aa,'plummer_perl.dat',maxcolumns=10000
   ss=size(aa)

   ntime=ss(1)
   npoints=(ss(2)-1)/8

   print,'Time bins:',ntime
   print,'Particles: ',npoints

   FOR t=0,ntime-1 DO BEGIN
      plot,[0],[0],xr=[-2,2],yr=[-2,2],/nodata,xst=1,yst=1,xtit='x',ytit='y',tit='time='+string(aa(t,0),form='(f5.2)')
      FOR i=0,npoints-1 DO BEGIN
         oplot,[aa(t,1+i*8)],[aa(t,1+i*8+1)],psym=4
      END

      plot,[0],[0],xr=[0,2],yr=[-10,10],/nodata,xst=1,yst=1,xtit='r',ytit='!Bv_r!N'
      FOR i=0,npoints-1 DO BEGIN
         rr=sqrt(aa(t,1+i*8+0)^2 + aa(t,1+i*8+1)^2 + aa(t,1+i*8+2)^2)
         vv=(aa(t,1+i*8+0)*aa(t,1+i*8+3)+aa(t,1+i*8+1)*aa(t,1+i*8+4)+aa(t,1+i*8+2)*aa(t,1+i*8+5))/rr
         oplot,[rr],[vv],psym=4
      END
      name='phase_'+string(t,form='(i04)')+'.png'
      write_png,name,tvrd(/true)
   END

end