pro vbmv2tg,vbmvfile,vv,bb,lgt,lgg,bcc,err
; This routine attempts to determine the values of logTeff, logg that give
; specified values of V_abs and B-V.
; On input, it accepts the pathname vbmvfile, which is an idl .sav file
; containing vectors logt, logg, and corresponding arrays vabs, bmv, bolometric
; correction bc.
; It also accepts values
;  vv = desired value of vabs
;  bb = desired value of B-V
; On output, it returns estimates of the corresponding log(T), log(g), and BC.
; If any problem occurred during the calculation, then on return err=-1, else 0.

; constants
wtv=0.1          ; range in B-V is about this much times range in Vabs,
                 ; hence this is the relative weighting.
dvmax=0.3        ; insist best point matches vabs better than this
dbmax=0.06       ; ditto for B-V
err=0

restore,vbmvfile
; contains arrays logg(nlg),logt(nlt),vabs(nlt,nlg),bmv(nlt,nlg),
;                 llogt(nlt,nlg),llogg(nlt,nlg),bc(nlt,nlg)
nlg=n_elements(logg)
nlt=n_elements(logt)

; find points closest to desired vabs, B-V.  Weight B-V distance more strongly
dv=vabs-vv
db=bmv-bb
dist=sqrt((dv*wtv)^2+db^2)
md=min(dist,ix)             ; find smallest distance

; test for adequate fit
dv0=abs(dv(ix))
db0=abs(db(ix))
if(dv0 gt dvmax or db0 gt dbmax) then begin
  err=-1
  goto,fini
endif

; fit okay, so where is this point?  boundaries of neighbor region?
ixy=long(ix/nlt)
ixx=ix-ixy*nlt
ixbot=(ixx-1) > 0
ixtop=(ixx+1) < (nlt-1)
iybot=(ixy-1)> 0
iytop=(ixy+1) < (nlg-1)

; retrieve the deltas within this box, fit planes
dv1=dv(ixbot:ixtop,iybot:iytop)
db1=db(ixbot:ixtop,iybot:iytop)
ltb=llogt(ixbot:ixtop,iybot:iytop)
lgb=llogg(ixbot:ixtop,iybot:iytop)
bcb=bc(ixbot:ixtop,iybot:iytop)
npt=n_elements(dv1)
funs=fltarr(npt,3)
funs(*,0)=1.
funs(*,1)=reform(ltb)
funs(*,2)=reform(lgb)
wt=fltarr(npt)+1
cc1=lstsqr(reform(dv1,npt),funs,wt,3,rms,chisq,outpv,1)
cc2=lstsqr(reform(db1,npt),funs,wt,3,rms,chisq,outpb,1)
cc3=lstsqr(reform(bcb,npt),funs,wt,3,rms,chisq,outpc,1)

; now solve the linear equations to force both delta vabs and delta B-V =0 
rhs=[-dv(ix),-db(ix)]
aa=[[cc1(1:2)],[cc2(1:2)]]
aai=invert(aa)
sol=reform(rhs#aai)  ; desired displacement from central logt, logg in $
                   ; units of logt, logg

; check to see that displacements are smaller than one step in each coord.
lt0=llogt(ix)
lg0=llogg(ix)
dlt=abs(logt(1)-llogt(0))
dlg=abs(logg(1)-llogg(0))
if(abs(sol(0)) gt dlt or abs(sol(1)) gt dlg) then begin
  err=-1
  goto,fini
endif

; if so, accept them and report result.
lgt=lt0+sol(0)
lgg=lg0+sol(1)
bcc=cc3(0)+lgt*cc3(1)+lgg*cc3(2)

fini:

end