pro prob_hr,tran,dlt,loglran,dll,logt,logl,hist
; This routine reads the hipparcos data for stars with good parallaxes,
; and uses them to create a histogram containing number of stars per
; increment in log(Teff) and per unit log(lum/lum_sun).  The range of
; the histogram is given by tran(2) (units of K, not logK) and loglran(2).
; The bin size in logT is dlt, and in log(lum) is dll.
; Values corresponding to the centers of the bins are returned in logt, logl,
; and the histogram in array hist.  logt is in units of log(teff/teff_sun).

; constants
tsun=5777.                    ; solar teff
bcsun=-.123                   ; solar bolometric correction
vabssun=4.87                  ; solar abs V mag.
vbmvfile='/home/tbrown/d/stardat/vabs_bmv.sav'
ptc=[1.,-.09,0.85]            ; photometric transformation coeffs to get
                              ; from tycho to bessell magnitudes, from
                              ; ESA SP-1200 (1997) via Vizier

; read the data
rd_hip,ra,dec,hip,vmag,bmv,plax,eplax,pmra,pmdec

; ***********************************************************************
; *** do not do this, since it gives really bad comparison between Tycho
; HR diagram and Yale model isochrone for M67
; ***********************************************************************
; transform the magnitudes and colors
;vmag=vmag*ptc(0) + bmv*ptc(1)
;bmv=bmv*ptc(2)

; construct V absolute magnitudes
distpc=1000./plax
dmod=5.*alog10(distpc/10.)
vabs=vmag-dmod

; get logt, logl
nb=n_elements(bmv)
logti=fltarr(nb)
logli=fltarr(nb)
for i=0,nb-1 do begin
  vbmv2tg,vbmvfile,vabs(i),bmv(i),tmpti,tmpgi,tmpbci,err
  logti(i)=tmpti-alog10(tsun)
  logli(i)=-0.4*(vabs(i)-vabssun+tmpbci-bcsun)
;   tl_vs_cv,bmv,vabs,logti,logli
endfor

; make histogram array
ltran=alog10(tran)
nlt=(max(ltran)-min(ltran))/dlt + 1
logtbot=min(ltran)+dlt*findgen(nlt)-alog10(tsun)
logt=logtbot+dlt/2.
nll=(max(loglran)-min(loglran))/dll+1
loglbot=min(loglran)+dll*findgen(nll)
logl=loglbot+dll/2.
hist=lonarr(nlt,nll)

; now fill it up
for i=0,nlt-1 do begin
  for j=0,nll-1 do begin
    s=where(logti ge logtbot(i) and logti lt logtbot(i)+dlt and $
            logli ge loglbot(j) and logli lt loglbot(j)+dll,ns)
    hist(i,j)=ns
  endfor
endfor

end