;
;  summary
;
cwd,run
dz=z[1]-z[0]
Sigma=dz*total(rho)
csmax=max(cs)
Omega=param.Omega
Heat=param.heat_uniform
zHeat=param.zheat_uniform_range
Height=csmax/Omega
Mdot=3.*!pi*2.*Heat*zHeat/(1.5*Omega)^2
alpSS=Mdot/(3.*!pi*Sigma*Omega*Height^2)
;
;  compute mass outside the simulated domain
;
zmax=5.
zext=grange(-zmax,zmax,1000)
TTsurf=min(TT)
cssurf=min(cs)
rhosurf=min(rho)
Hsurf=cssurf/(sqrt(gam)*Omega)
rhoext=rhosurf*exp(-z^2/(2.*Hsurf^2))
fact=rhosurf/min(rhoext)
plot_io,z,rho,xr=[-1,1]*zmax,yr=[1e-7*min(rho),2*max(rho)]
rhoext=rhosurf*exp(-zext^2/(2.*Hsurf^2))
oplot,zext,fact*rhoext,col=122,li=2
good=where(abs(zext) gt max(z))
dzext=z(1)-z(0)
Sigmaext=dzext*total_half(rhoext(good))
print,'Sigma,Sigmaext=',Sigma,Sigmaext
;
;  print results
;
fo='(7e9.2,2x,a)'
openw,1,'Summary.tmp'
printf,1,'   Sigma    Omega   Height   csmax     Heat      Mdot    alpSS   Run'
printf,1,'[g/cm3 Mm] [ks^-1]   [Mm]    [km/s] [g/cm3/Mm [1e21 g/s]           '
printf,1,'                                    (km/s)^3                       '
printf,1,Sigma,Omega,Height,csmax,heat,Mdot,alpSS,run,fo=fo
close,1
spawn,'cat Summary.tmp '
spawn,'cat Summary.tmp >> ../idl/Summary.dat'
END