if !d.name eq 'PS' then begin
  device,xsize=18,ysize=7,yoffset=3
  !p.charthick=2 & !p.thick=2 & !x.thick=2 & !y.thick=2
end
;
;  mv idl.ps ../fig/psolution.ps
;
!p.charsize=2.3
!x.margin=[8.5,.5]
!y.margin=[3.2,.5]
!p.multi=[0,3,1]
;
;  input parameters for the model
;
Teff=16300.
Teff=5200.
Teff=9300.
a=1. & b=0.
Ttop=Teff/2^.25
print,Ttop
;
;  units of BB14
;
rho0=4e-4
T0=38968.
sigmaSB=5.67d-5 ;(erg/cm2/s/K4)
Rgas=8.31434d7  ;erg/mol/K
kap0=1e6/1e8    ;(cm2/g)
grav=274e2      ;(cm/s^2)
Mm=1e8          ;(1 Mm)
mu=.6
;
gam=5./3.
nabad=1.-1./gam
cp=Rgas/mu/nabad
;
P0=(Rgas/mu)*rho0*T0
K0=16.*sigmaSB*T0^3/(3.*kap0*rho0)
Frad=sigmaSB*Teff^4
nabrad0=Frad*P0/(K0*T0*rho0*grav)
;
;  uniform array in lnP
;
nP=1000
lnP=10.+15.*findgen(nP)/(nP-1)
P=exp(lnP)
n=(3.-b)/(1.+a)
RHS=(n+1.)*nabrad0*(P/P0)^(1.+a)+(Ttop/T0)^(4.+a-b)
T=T0*RHS^(1./(4.+a-b))
rho=P/(Rgas*T/mu)
;
print,nabrad0
z=fltarr(nP)
f=-P/(rho*grav) ;(integrand)
for i=1,nP-1 do z[i]=z[i-1]+.5*(f[i]+f[i-1])*(lnP[i]-lnP[i-1])
;
;  determine bottom value zbot
;
good=where(T le T0)
ibot=max(good)
zbot=z(ibot)
znew=(z-zbot)/Mm ;(in Mm)
;
;  entropy in units of cp
;
S=alog(P)/gam-alog(rho)
Sbot=S[ibot]
S=S-Sbot
;
!x.title='!8z!6 [Mm]'
plot,znew,T,xr=[0,5],yr=[0,T0],ytit='!8T!6 [K]'
plot_io,znew,rho,xr=[0,5],ytit='!7q!6 [g/cm!u3!n]'
plot,znew,S,xr=[0,5],ytit='!8S!6/!8c!dP!n!6'
;
!x.title='!6'
!p.multi=0
END