pro show_blast
  loadct,5

  !P.MULTI=[0,3,1]
  gamma = 5./3.
  rho0 = 0.125
  rho_max = (gamma+1)/(gamma-1)*rho0

  for i=80,147 do begin
     snr=string(i,form='(i03)')
     name="snap_"+snr
     
     readnew,name,h,"HEAD"
     readnew,name,x,"POS"
     readnew,name,v,"VEL"
     readnew,name,rho,"RHO"
     readnew,name,u,"U"

     p=(gamma-1)*u*rho

     r1=0.21*h.time^(2./5.)
     u1=0.08*h.time^(-3./5.)
     P1=1e-3*h.time^(-6./5.)
     
     rr=sqrt((x[0,*]-0.5)^2+(x[1,*]-0.5)^2+(x[2,*]-0.5)^2)
     vv=sqrt(v[0,*]^2+v[1,*]^2+v[2,*]^2)
     
     plot,rr,p,psym=3,xtit='r',ytit='p',yr=[1e-20,0.1],charsize=3,charthick=2,yst=1,ylog=1
     oplot,r1*[1,1],[1e-20,1],l=2
     oplot,[0,1],P1*[1,1],l=2
     
     plot,rr,rho,psym=3,xtit='r',ytit='rho',yr=[0,0.5],charsize=3,charthick=2
     oplot,[0,1],rho_max*[1,1],l=2
     oplot,[0,1],rho0*[1,1],l=2
     oplot,r1*[1,1],[0,1],l=2

     plot,rr,vv,psym=3,xtit='r',ytit='v',yr=[0,0.5],charsize=3,charthick=2,yst=1
     oplot,r1*[1,1],[0,1],l=2
     oplot,[0,1],u1*[1,1],l=2
     
     save_screen,"frame_"+snr,/bmp

  end
  stop
end