pro show_nozzle
  loadct,5

  !P.MULTI=[0,1,3]
  L=20
  gamma = 5./3.
  rho_ic=0.125/L
  P_ic = 1.0/L
  csnd = sqrt(gamma * P_ic/rho_ic)

  for i=146,150 do begin
     snr=string(i,form='(i03)')
     name="snap_"+snr
     
     readnew,name,h,"HEAD"
     readnew,name,x,"POS"
     readnew,name,id,'ID'
     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

     jj=where(id GE 1)
;     jj=where(abs(x[2,*]-0.5) LT 0.1 AND id GE 1)
     jjb=where(abs(x[2,*]-0.5) LT 0.05 AND id EQ 0)
     plot,[1],[1],/nodata,yr=[0,1],xr=[10,40],xst=1,yst=1,xtit='x',ytit='y',charthick=2,charsize=3,tit='t='+string(h.time,form='(f5.2)')

     umin=5.
     umax=13.
     cmin=20
     cmax=239
     icol = cmin + (u - umin)/(umax-umin) * (cmax - cmin)
     ii=where(icol LT cmin,nfound)
     if nfound GT 0 THEN icol[ii]=cmin
     ii=where(icol GT cmax,nfound)
     if nfound GT 0 THEN icol[ii]=cmax
     FOR j=0L,N_ELEMENTS(jj)-1 DO BEGIN
        sym=3
;        if v[0,jj[j]] LT -1*csnd THEN sym=4 ELSE sym=3
        oplot,[x[0,jj[j]]],[x[1,jj[j]]],psym=sym,col=icol[jj[j]]
     END
     oplot,x[0,jjb],x[1,jjb],psym=4

     jj=where(abs(x[1,*]-0.5) LT 0.1)
     jjb=where(abs(x[1,*]-0.5) LT 0.1 AND id EQ 0)

;     for j=0,5 do begin
;        x0=13-csnd*(h.TIME-FLOOR(h.TIME)+j)
;        oplot,x0*[1,1],[0,1],col=128
;     end
     
     plot,x[0,jj],p[jj],psym=3,xst=1,yst=1,xr=[10,40],yr=[0.0,0.42],xtit='x',ytit='P',charthick=2,charsize=3
     oplot,x[0,jjb],p[jjb],psym=4

;     plot,x[0,jj],u[jj],psym=3,xst=1,yst=1,xr=[10,40],yr=[umin,umax],xtit='x',ytit='u',charthick=2,charsize=3
;     oplot,x[0,jjb],u[jjb],psym=4
;     for j=0,5 do begin
;        x0=13-csnd*(h.TIME-FLOOR(h.TIME)+j)
;        oplot,x0*[1,1],[0,100],col=128
;     end
    
     plot,x[0,jj],v[0,jj],psym=3,xst=1,yst=1,xr=[10,40],yr=[-5,0.5],xtit='x',ytit='vx',charthick=2,charsize=3
     cc=sqrt(gamma * P/rho)
     oplot,x[0,jj],-1*cc[jj],psym=3,col=128
;     oplot,[0,20],csnd*[-1,-1],col=128
;     oplot,[0,20],1.5*[-1,-1],col=20,l=2
     oplot,x[0,jjb],v[0,jjb],psym=4
;     for j=0,5 do begin
;        x0=13-csnd*(h.TIME-FLOOR(h.TIME)+j)
;        oplot,x0*[1,1],[-100,100],col=128
;     end

     save_screen,"frame_"+snr,/bmp

  end
  stop
end