;$Id: pdt.pro,v 1.9 2019/01/17 12:57:55 brandenb Exp $
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=11,yoffset=3
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
end
;
!p.charsize=1.3
!x.margin=[8.8,.5]
!y.margin=[3.2,.5]
;
@parameters
@data/pc_constants
;
;  mv idl.ps ~/tex/notes/raddisk/fig/pvar.ps
;
TT=reform(xya.TTmz[*,nt-1])
pp=reform(xya.ppmz[*,nt-1])
rho=reform(xya.rhomz[*,nt-1])
kapparho=reform(xya.kapparhomz[*,nt-1])
Frad=reform(xya.Fradzmz[*,nt-1])
;
;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.*kapparho)
chi=KK/(rho*cp)
cgam=16*sigmaSB*TT^3/(cp*rho)
ell=1./kapparho
;
;  compute mass outside the simulated domain
;
TTmin=min(TT)
;
;  compute sigma
;
dz=z(1)-z(0)
ssigma=dz*total_half(rho)
print,'Sigma=',ssigma
;
;!p.multi=[0,1,2]
;plot,z,TT,yst=0
;plot,tau,TT^4,xr=[0,1e-3],yr=[0,8000]
!p.multi=0
!p.multi=[0,1,3]
!p.multi=[0,2,2]
plot,z,TT,yst=0,ytit='!8T!6 [K]'
;
circ_sym,1.3,1
Sigma=integr(rho,x=z)
tau=integr(kapparho,x=z,/rev)
tau2=integr(kapparho,x=z)
itau1=findex(1.,tau,/rev)
itau2=findex(1.,tau2)
oplot,[1,1]*z(itau1),[1,1]*TT(itau1),ps=8,col=122
oplot,[1,1]*z(itau2),[1,1]*TT(itau2),ps=8,col=122
ztau1=z(itau1)
ztau2=z(itau2)
print,'ztau1,2=',ztau1,ztau2
;tau=tau-mean(tau)
;
;  panel 2
;
ss=alog(TT/pp^.4)
yrs=[12,15]
;plot,z,ss,yr=yrs,ytit='!8s!6/!8c!6!dp!n'

;           if (f(l1-1+l,m,n,ikapparho)**2>dxyz_2(l)) then
;             dt1_rad(l)=4*kappa(l)*sigmaSB*p%TT(l)**3*p%cv1(l)* &
;                 dxyz_2(l)/f(l1-1+l,m,n,ikapparho)**2/cdtrad
;             if (z_cutoff/=impossible .and. cool_wid/=impossible) &
;             dt1_rad(l)=0.5*dt1_rad(l)*(1.-tanh((z(n)-z_cutoff)/cool_wid))
;           else
;             dt1_rad(l)=4*kappa(l)*sigmaSB*p%TT(l)**3*p%cv1(l)/cdtrad
;             if (z_cutoff/=impossible .and. cool_wid/=impossible) &
;             dt1_rad(l)=0.5*dt1_rad(l)*(1.-tanh((z(n)-z_cutoff)/cool_wid))
;
pc_read_ts,obj=ts
ntt=n_elements(ts.dt)
dt=ts.dt[ntt-1]
kNy=!pi/dz
;dtcrit=ell/cgam
;dtcrit2=1./(chi*kNy^2)
Ccool=4.2
dtcrit=Ccool*ell/cgam
dtcrit2=.20*dz^2/chi
plot_io,z,dtcrit+dtcrit2,ytit='dt1+dt2',yr=minmax([dtcrit,dtcrit2])
oplot,z,dtcrit,col=55,li=2
oplot,z,dtcrit2,col=122
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,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
loadct,6
oplot,z,.9*dz/cs,col=122
nu=5e-4
restore,'shock.sav'
Cshock=1.
nu_shockmax=Cshock*shockmax
nu_shockm=Cshock*shockm
oplot,z,z*0+.03*dz^2/nu,col=122,li=2
oplot,z,z*0+.03*dz^2/nu_shockmax,col=122,li=2
oplot,z,z*0+.03*dz^2/nu_shockm,col=122,li=1
loadct,5
;
plot_io,z,ell,ytit='ell'
plot_io,z,cgam,ytit='cgam',yr=minmax([cgam,cs])
oplot,z,cs,li=2
;
save,file='pdt.sav',z,dtcrit,dtcrit2,itau1,itau2,nu_shockmax,nu_shockm,nu,dz,dt,cs
;
;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=Frad
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/pvar_'+run+'.ps'
;
!p.multi=0
stop
save,file='rhoTz.sav',rho,TT,kapparho,z,tau,chi,cp,ss,cs,Frad
;
END