PRO sstab
;
;  Replace alog(rho*kappa) with S in the lookup table, iv=1
;  21-jan-96:  Changed to replace speudo blanketed Planck-function with entropy.
;              iv=3
;
    common ctable,ttmean,rhm,drhm,eemin,eemax,dee,mtab,itab,tab

    nz=n_elements(rhm)
    ;plot_oi,rhm,/nodata,[-.001,.0005],ytit='S',xtit='rho',yst=1

    dlnr=alog(rhm(nz-1)/rhm(0))/(nz-1)
    for nr=0,nz-1 do begin
      ntab=mtab(nr)
      ltab=3*ntab
      kpg=itab(nr)-1+indgen(ntab)
      ktt=kpg+2*ltab
      kss=ktt+ltab
      tt=tab(ktt)
      dttdee=tab(ktt+ntab)
      pg=exp(tab(kpg))
      dssdee=1./tt
      d2ssdee2=-dttdee/tt^2
      dssdlnr=-pg/rhm(nr)/tt
      ss=fltarr(ntab)
      dss=dee(nr)*(0.5*(dssdee+dssdee(1:*))+(d2ssdee2-d2ssdee2(1:*))/12.)

      if nr gt 0 then begin
        for k=1,ntab-1 do ss(k)=ss(k-1)+dss(k-1)
		if eemin(nr) gt eemin(nr-1) then begin
		  k1=long((eemin(nr)-eemin(nr-1))/dee(nr)+.5)
		  ss=ss+aver(ssp(k1:*)-ss(*)+0.5*dlnr*(dssdlnrp(k1:*)+dssdlnr(*)))
		end else begin
		  k1=long((eemin(nr-1)-eemin(nr))/dee(nr)+.5)
		  ss=ss+aver(ssp(*)-ss(k1:*)+0.5*dlnr*(dssdlnrp(*)+dssdlnr(k1:*)))
		end
      end
      ssp=ss
      dssdlnrp=dssdlnr

      tab(kss)=ss
      dssdee=(ss(2:*)-ss)/2.
      dssdee=[2.*dssdee(0)-dssdee(1),dssdee,2.*dssdee(ntab-3)-dssdee(ntab-4)]
      tab(kss+ntab)=dssdee
      tab(kss+2*ntab)=dssdlnr*dlnr

      ;oplot,replicate(rhm(nr),ntab),ss,psym=3
    end
END