pro show_jet
  loadct,5

  !P.MULTI=[0,1,1]
  L1=10
  L2=30
  gamma = 5./3.
  rho_ic=1.0
  P_ic = 1.0
  csnd = sqrt(gamma * P_ic/rho_ic)

  umin=alog10(0.2)
  umax=alog10(10.)
;  umin=10.
;  umax=15.
  cmin=20
  cmax=239

  for i=0,3 do begin
     snr=string(i,form='(i03)')
     name="snap_"+snr
     name1='./'+name
     
     readnew,name1,h1,"HEAD"
     readnew,name1,x1,"POS"
     readnew,name1,id1,'ID'
     readnew,name1,v1,"VEL"
     readnew,name1,rho1,"RHO"
     readnew,name1,u1,"U"

     p1=(gamma-1)*u1*rho1
     
     jj1=where(id1 GE 1)
     jjb1=where(id1 EQ 0)
     plot,[1],[1],/nodata,yr=[0,3],xr=[0,1],xst=1,yst=1,xtit='x',ytit='y',charthick=2,charsize=2,tit="Simple SPH"

     icol1 = cmin + (alog10(u1) - umin)/(umax-umin) * (cmax - cmin)
     ii=where(icol1 LT cmin,nfound)
     if nfound GT 0 THEN icol1[ii]=cmin
     ii=where(icol1 GT cmax,nfound)
     if nfound GT 0 THEN icol1[ii]=cmax
     FOR j=0L,N_ELEMENTS(jj1)-1 DO BEGIN
        sym=3
        oplot,[x1[0,jj1[j]]],[x1[1,jj1[j]]],psym=sym,col=icol1[jj1[j]]
     END
     oplot,x1[0,jjb1],x1[1,jjb1],psym=4

     save_screen,"frame_"+snr

  end
  stop
end