;$Id: pdt.pro,v 1.1 2019/01/02 07:36:22 brandenb Exp $
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=4,yoffset=3
  !p.charthick=1.6 & !p.thick=1.6 & !x.thick=1.6 & !y.thick=1.6
end
;
!p.charsize=1.6
!x.margin=[8.2,.5]
!y.margin=[3.2,.5]
;
@parameters
!x.title='!8z!6 [Mm]'
;
;  mv idl.ps ~/tex/notes/raddisk/fig/pvar.ps
;
cs0=start.cs0
rho0=start.rho0
lnrho0=alog(rho0)
grav=-start.gravz
gam=start.gamma
cp=start.cp
lnrho=alog(rho0)
cs2=cs0^2*exp((gam-1.)*(var.lnrho-lnrho0)+gam*var.ss/cp)
cs=sqrt(cs2)
Hp=cs2/(gam*grav)
TT=cs2/((gam-1.)*cp)
rho=exp(var.lnrho)
pp=(1.-1./gam)*cp*TT*rho
zzz=var.z
KK=16.*sigmaSB*TT^3/(3.*var.kapparho)
chi=KK/(rho*cp)
cgam=16*sigmaSB*tt^3/(cp*rho)
ell=1./var.kapparho
;
;  compute mass outside the simulated domain
;
TTmin=min(TT)
;
;  compute sigma
;
;dz=z(1)-z(0)
dz=0.
z=[0,1]
;
!p.multi=0
!p.multi=[0,2,2]
!p.multi=[0,1,3]
!p.multi=[0,3,1]
plot,z,[tt,tt],yst=0,ytit='!8T!6 [K]'
;
;  panel 2
;
ss=alog(tt/pp^.4)
yrs=[12,15]
;
nt=n_elements(ts.dt)
dt=ts.dt[nt-1]
kNy=!pi/dz
Ccool=1.
dtcrit=4.2*Ccool*ell/cgam
dtcrit2=.03*dz^2/chi
dtcritg=.8*dz/cgam
ytit='!7d!8t!6 [ks]'
plot_io,z,replicate(dtcrit+dtcrit2,2),ytit=ytit
oplot,z,replicate(dtcrit,2),col=55,li=2
oplot,z,replicate(dtcrit2,2),col=122
oplot,z,replicate(dtcritg,2),col=155
;oplot,[1,1]*z(itau1),[1,1]*dtcrit(itau1),ps=8,col=122
;oplot,[1,1]*z(itau2),[1,1]*dtcrit(itau2),ps=8,col=122
oplot,z,[0,0]+z*0+dt,li=3
;plot,z,deriv(z,ss),yst=0,xtit='!8z!6 [Mm]',ytit='!6d!8s!6/!8c!6!dp!n/d!8z!6'
;oplot,z,z*0,li=3
;
;plot_io,z,ell,ytit='ell'
;
ytit='!8c!d!7c!n!6 [km/s]'
plot_io,z,[0,0]+cgam,ytit=ytit,yr=minmax([cgam,cs]),yst=0
oplot,z,[0,0]+cs,li=2
;
;oplot,z,5.*(6.4)^2*(-z^2)+39e3,col=122
;plot,tau,tt^4,xr=[0,1e-3],yr=[0,8000]
;
heat_uniform=param.heat_uniform
zheat_uniform_range=param.zheat_uniform_range
;flux=var.kr_frad(*,2)/var.kapparho
;Teff=abs(minmax(flux/sigmaSB))^.25
;Teff0=(heat_uniform*zheat_uniform_range/sigmaSB)^.25
;print,'Teff=',Teff0,Teff
;
cwd,run
print,'$mv idl.ps ~/tex/upasana/tstep/fig/pdt_'+run+'.ps'
;
!p.multi=0
stop
save,file='pdt.sav',z,dtcrit,dtcrit2,dtcritg,itau1,dt,tt,cs,cgam
;save,file='rhoTz.sav',rho,TT,z,tau,chi,cp,ss
;
END