pro linked_props2,logto,loggo,m,b,logl,logd,a_v
; This routine accepts vectors logt,logg, and scalars m, b describing
; stellar models, and returns corresponding parameter vectors logl, logd, a_v.
; On input,
;  logto = alog10(teff) of the models, logto=logto(nt)
;  loggo = alog10(g) (cgs units), loggo = loggo(ng)
;  m = r magnitude observed for star
;  b = galactic latitude (degrees)
; On output,
;  logl = alog10(L/L_sun), logl=logl(nt,ng)
;  logd = alog10(distance), with distance in pc, logd=logd(nt,ng)
;  a_v = total V-band extinction along the line of sight, in stellar magnitudes, 
;         a_v(nt,ng)

; constants
radian=180./!pi
vabssun=4.83        ; V band abs magnitude of sun
tsun=5777.          ; effective temperature of sun
loggsun=4.438        ; log(g) for sun
kv = 0.001          ; V-band extinction (mags per pc)

;kv=0.00044           ; **TEST** for M67 only!!!

kr = 0.88*kv        ; r-band extinction (mags per pc)
; Looks like the first-principles value is about kr = 0.88 * kv
niter=21
ee=.001             ; satisfied with last adjustment to logd less than this
hdust=150.          ; dust scale height in pc
filin='/home/tbrown/d/basel/colors/Corrected/lcb97_p00.cor.UBVRIJHKLM'

; make 2-dimensional versions of logt, logg
;nt=n_elements(logto)
;ng=n_elements(loggo)
;logt2=rebin(logto,nt,ng)
;logg2=rebin(reform(loggo,1,ng),nt,ng)

; look up star parameters on desired grid:  log(L), BC, V-r
;parmblkintrp,logt2,logg2,logl,bcc,vmrc
parmblkintrp,logto,loggo,logl,bcc,vmrc

; compute log(L), assuming mass = solar mass
;ltsun=alog10(tsun)
;tg2vbmv,filin,logto,loggo,vabs2,bmv2,bc2
;logmass=starmass(logt,logg)        ; estimate star mass from posn in HR diag
;logl=4.*(logt-ltsun)-(logg-loggsun) ; +logmass

; compute distance iteratively
;bcc=interpol(bc,logbct,logt)      ; bolometric corr at each input teff
;vmrc=interpol(vmr,logbct,logt)    ; V-R color ditto (note: not really V-r,
                                  ; which is what we actually want here)
; start with maximum distance possible
zp=vabssun-bcc-vmrc-5.            ; stuff that doesn't depend on distance
lcur=0.2*m+0.5*logl-0.2*zp
sinb=sin(abs(b)/radian) > 0.005        ; forces b > 0.28 degree 
c0=hdust*kr/sinb
np=n_elements(lcur)
lbot=fltarr(np)                  ; starting min distance = 1 pc
ltop=lcur
for i=0,niter-1 do begin
  lcur=(lbot+ltop)/2.             ; next guess bisects current low, high limits
  cur=10.^lcur
  tau=c0*(1.-exp(-cur*sinb/hdust))
  elhs=lcur-logl/2.+(tau+zp-m)/5.  ; if positive, distance is too large
  sp=where(elhs ge 0.,nsp)
  if(nsp gt 0) then ltop(sp)=lcur(sp)
  sn=where(elhs lt 0.,nsn)
  if(nsn gt 0) then lbot(sn)=lcur(sn)
  dl=lcur-(ltop+lbot)/2.
  if(max(abs(dl)) le ee) then goto,fini
endfor

fini: logd=lcur
a_v=tau*kv/kr

end