pro mk_hrblock,gwid,hrb,logtn,logln,hrbout
; This routine makes a smooth, continuous, interpolatable table called hrb
; containing the ln of the relative probability of finding a star within
; a specified cell in log(teff), log(bolometric luminosity) space.  Also
; returned are arrays logtn, logln giving the corresponding teff and luminosity
; (both in solar units).  Results are returned to the calling program, and
; also written as an IDL save file to hrbout.
; Method is to edit and smooth a histogram derived from ~9800 stars for
; which there are good Hipparcos magnitudes, colors, and parallaxes.
; On input, parameter gwid(2) specifies the gaussian smoothing width in 
; logt and logl 

; constants
tran=[3200.,15000.]        ; histogram teff range
dlt=.01                    ; raw histogram bin size in log(teff)
loglran=[-2.,4.]           ; range in log(lum)
dll=0.1                    ; bin size in log(lum)
bgn=3.                     ; scales total number of bright giants
fms=2.5                    ; scales number of added faint MS stars
lteffb=[3.47712,3.51188,3.54407,3.57403,3.60206,3.62839,$
      3.65321,3.67669,3.69897,3.72016,3.74036,3.75967] ; ordinates for BC vs teff
bcc=[-2.61567,-2.38699,-1.65687,-1.31560,-0.986509,-0.756203,$
    -0.589964,-0.409369,-0.290209,-0.229958,-0.172984,-0.124102]
;  The above values come from parmblkintrp run on a grid with 250K steps in Teff
;  assuming logg = 4.5, which supposedly makes sense for lower MS stars
logtsun=alog10(5777.)

; get histogram
prob_hr,tran,dlt,loglran,dll,logt,logl,hist

; make 2d coordinate tables
sz=size(hist)
nlt=sz(1)
nll=sz(2)
logt2=rebin(logt,nlt,nll)
logl2=rebin(reform(logl,1,nll),nlt,nll)

; trim out questionable data
s=where(logt2 lt (-.21) and logl2 lt (1.),ns)
if(ns gt 0) then hist(s)=0
s=where(logl2 lt (-0.6+8.*logt2),ns)
if(ns gt 0) then hist(s)=0
s=where(hist lt 2,ns)
if(ns gt 0) then hist(s)=0
hist(5,1)=0
hist(9,4)=0      ; squash 2 points that make trouble later.

uu1=hist
;stop

; add a branch for bright giants
hist0=hist            ; save it for later comparison
hist=float(hist)
sx=where(logt le (-.1082),nsx)
ltx=logt(sx)
xx=ltx+.1082
ff=1./(1.+(xx/.07)^2)
yy=-0.2-21.97*ltx-28.33*ltx^2   ; fit to M67 isochrone, brighter by .1 in logL
;yy=1.96-8.61*xx+19.3*xx^2
yy1=yy+.15-2.*xx
sy=(yy-loglran(0))/dll
sy1=(yy1-loglran(0))/dll
for i=0,nsx-1 do begin
  hist(sx(i),sy(i))=hist(sx(i),sy(i))+bgn*ff(i)
  if(sy1(i) lt nll) then hist(sx(i),sy1(i))=hist(sx(i),sy1(i))+bgn*ff(i)
endfor

uu2=hist
;stop

; scale results to allow for limited survey volume for faint stars
bcc2=interpol(bcc,lteffb,logt2+logtsun)
logrvol=1.5*(logl2 + 1.51 + 0.6*bcc2) < 0.
factor=10.^logrvol
hist=hist/factor
uu3=hist

;stop

; and another for faint MS stars
sx=where(logt le (-.07),nsx)
ltx=logt(sx)
xxa=ltx+.07
ff2=1.-atan(xxa/.03)*!pi*(abs(xxa)/.04)^1.5
;yy2=-0.053+8.02*ltx+5.90*ltx^2
yy2=-0.003+8.02*ltx+5.90*ltx^2
sy2=(yy2-loglran(0))/dll
for i=0,nsx-1 do begin
  if(sy2(i) ge 0) then hist(sx(i),sy2(i))=hist(sx(i),sy2(i))+fms*ff2(i)
endfor

uu4=hist
hist(5,1)=0
hist(9,4)=0      ; squash 2 points that make trouble later.
;stop

; make smoothed version, brute force (to get the far wings continuous and >0)
vhrb=dblarr(nlt,nll)
for i=0,nlt-1 do begin
  dt=abs(logt2-logt(i))
  argt=(dt/gwid(0))^2 < (dt/gwid(0))*10.
  for j=0,nll-1 do begin
    dl=abs(logl2-logl(j))
    argl=(dl/gwid(1))^2 < (dl/gwid(1))*10.
    wts=exp(-((argt+argl) < 60.))
    vhrb(i,j)=total(hist*wts)
  endfor
endfor

hrb=float(alog(vhrb))
logtn=logt
logln=logl

stop

save,hrb,logtn,logln,file=hrbout

end