;$Id: pdt.pro,v 1.6 2019/01/03 18:30:06 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)
ssigma=dz*total_half(rho)
print,'Sigma=',ssigma
;
!p.multi=0
!p.multi=[0,2,2]
!p.multi=[0,1,3]
!p.multi=[0,3,1]
plot,z,tt,yst=0,ytit='!8T!6 [K]'
;
circ_sym,1.3,1
Sigma=integr(rho,x=z)
tau=integr(var.kapparho,x=z,/rev)
tau2=integr(var.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
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]

;           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))
;
nt=n_elements(ts.dt)
dt=ts.dt[nt-1]
kNy=!pi/dz
Ccool=4.2
dtcrit=Ccool*ell/cgam
dtcrit2=.03*dz^2/chi
dtcrit2=.35*dz^2/chi
dtcritg=1.0*dz/cgam
dtcritc=1.0*dz/cs
ytit='!7d!8t!6 [ks]'
plot_io,z,dtcrit+dtcrit2,ytit=ytit,yr=minmax([dtcrit,dtcrit2])
oplot,z,dtcrit,col=55,li=2
oplot,z,dtcrit2,col=122
oplot,z,dtcritg,col=155
oplot,z,dtcritc,col=188
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
;
;plot_io,z,ell,ytit='ell'
;
ytit='!8c!d!7c!n!6 [km/s]'
plot_io,z,cgam,ytit=ytit,yr=minmax([cgam,cs]),yst=0
oplot,z,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,dtcritc,itau1,dt,tt,cs,cgam
;save,file='rhoTz.sav',rho,TT,z,tau,chi,cp,ss
;
print,'chi[0]=',chi[0]
print,'cs[0]=',cs[0]
print,'dt*chi[0]/dz^2=',dt*chi[0]/dz^2
print,'dt*cs[0]/dz=',dt*cs[0]/dz
END