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)

  umin=11.5
  umax=13.
;  umin=10.
;  umax=15.
  cmin=20
  cmax=239

  for i=0,150 do begin
     snr=string(i,form='(i03)')
     name="snap_"+snr
     name1='WithoutConduction/'+name
     name2='WithConduction/'+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"

     readnew,name2,h2,"HEAD"
     readnew,name2,x2,"POS"
     readnew,name2,id2,'ID'
     readnew,name2,v2,"VEL"
     readnew,name2,rho2,"RHO"
     readnew,name2,u2,"U"

     p1=(gamma-1)*u1*rho1
     p2=(gamma-1)*u2*rho2
     
     jj1=where(id1 GE 1)
     jjb1=where(abs(x1[2,*]-0.5) LT 0.1 AND id1 EQ 0)
     plot,[1],[1],/nodata,yr=[0,1],xr=[0,20],xst=1,yst=1,xtit='x',ytit='y',charthick=2,charsize=3,tit="No conduction"

     icol1 = cmin + (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

     
     jj2=where(id2 GE 1)
     jjb2=where(abs(x2[2,*]-0.5) LT 0.1 AND id2 EQ 0)
     plot,[1],[1],/nodata,yr=[0,1],xr=[0,20],xst=1,yst=1,xtit='x',ytit='y',charthick=2,charsize=3,tit="with Conduction"

     icol2 = cmin + (u2 - umin)/(umax-umin) * (cmax - cmin)
     ii=where(icol2 LT cmin,nfound)
     if nfound GT 0 THEN icol2[ii]=cmin
     ii=where(icol2 GT cmax,nfound)
     if nfound GT 0 THEN icol2[ii]=cmax
     FOR j=0L,N_ELEMENTS(jj2)-1 DO BEGIN
        sym=3
        oplot,[x2[0,jj2[j]]],[x2[1,jj2[j]]],psym=sym,col=icol2[jj2[j]]
     END
     oplot,x2[0,jjb2],x2[1,jjb2],psym=4

     

     
     jj1=where(abs(x1[2,*]-0.5) LT 0.1 AND abs(x1[1,*]-0.5) LT 0.1 AND id1 GE 1)
     jjb1=where(abs(x1[2,*]-0.5) LT 0.1 AND abs(x1[1,*]-0.5) LT 0.1 AND id1 EQ 0)
     jj2=where(abs(x2[2,*]-0.5) LT 0.1 AND abs(x2[1,*]-0.5) LT 0.1 AND id2 GE 1)
     jjb2=where(abs(x2[2,*]-0.5) LT 0.1 AND abs(x2[1,*]-0.5) LT 0.1 AND id2 EQ 0)

     plot,x1[0,jj1],u1[jj1],psym=3,xst=1,yst=1,xr=[0,20],yr=[umin,umax],xtit='x',ytit='u',charthick=2,charsize=3
     oplot,x2[0,jj2],u2[jj2],psym=3,col=64
     oplot,[0,20],p_ic/(gamma-1)/rho_ic*[1,1],col=128


     oplot,x1[0,jjb1],u1[jjb1],psym=4

;     plot,x[0,jj],rho[jj],psym=3,xst=1,yst=1,xr=[0,20],yr=[0,0.25],xtit='x',ytit='rho',charthick=2,charsize=2
;     oplot,[0,20],rho_ic*8*[1,1],col=128

   
     save_screen,"frame_"+snr

  end
  stop
end