pro show_wave
  !P.MULTI=[0,1,4]
  L=20
  gamma = 5./3.
  rho_ic=0.125/L
  P_ic = 1.0/L
  csnd = sqrt(gamma * P_ic/rho_ic)
  
  for i=0,100 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
     
     x0=1+csnd*h.time
     while x0 GT L do begin
        x0-=L
     end
     x1=L-x0

     plot,x[0,*],p[*]/P_ic,psym=3,yr=1+0.03*[-1,1],xst=1,yst=1,xtit='x',ytit='P',charthick=2,charsize=2
     oplot,x0*[1,1],[-10,10],col=64
     oplot,x1*[1,1],[-10,10],col=64

     plot,x[0,*],rho[*]/rho_ic,psym=3,yr=1+0.02*[-1,1],xst=1,yst=1,xtit='x',ytit='rho',charthick=2,charsize=2
     oplot,x0*[1,1],[-10,10],col=64
     oplot,x1*[1,1],[-10,10],col=64
     
     plot,x[0,*],v[0,*],psym=3,xst=1,yst=1,yr=[-0.06e,0.06],xtit='x',ytit='vx',charthick=2,charsize=2
     oplot,x0*[1,1],[-10,10],col=64
     oplot,x1*[1,1],[-10,10],col=64

     plot,x[0,*],v[1,*],psym=3,xst=1,yst=1,yr=[-0.1,0.1],xtit='x',ytit='vy',charthick=2,charsize=2
     oplot,x0*[1,1],[-10,10],col=64
     oplot,x1*[1,1],[-10,10],col=64

     
     save_screen,"frame_"+snr
;     stop
  end
end