pro find_prop1,colvec,errvec,rmag,gb,$
    teff,logz,logg,logl,ext,props,properr,chmin,rats,xform=xform,flg=flg,$
    bias=bias

; This routine searches the color vs stellar parameter table colblock for
; the item that best matches (in a chi^2 sense) the input color vector
; colvec, with estimated errors errvec.
; Also included in the maximum-likelihood calculation is allowance for the
; known distributions of stars with teff, luminosity, metallicity, and
; distance from the galactic plane.
; Inputs for the single input star are the vector of colors colvec(9),
; the corresponding vector of estimated errors errvec(9),
; the stellar r magnitude rmag, and the galactic latitude gb (degrees).
; Also input are the array of model colors colblock and its coordinate vectors
; teffu, loggu, logzu,
; and the HR diagram log probability array hrb, with its coordinate vectors
; logtn, logln.
; Estimated properties are returned in the vector props, in the order
; {teff, log(Z), log(g), A_V}.  Estimated errors on these quantities
; are returned in properr.
; If keyword bias is set, giants are preferred over dwarfs a priori by an amount
; proportional to the bias value.

if(keyword_set(flg)) then sflg=1 else sflg=0      ; stop at checkpoints if set

; constants
wfact0=1.0               ; weighting of ln prior prob relative to chi^2
common tables,logln,logtn,hrb,$
              logtp,loggp,bcp,bmvp,loglump,logmassp,rabsp,vabsp,vmrp,$
              teffu,loggu,logzu,colblock

; indices defining colors from magnitudes
ii=[[0,1],[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[1,8],[1,9]]
sz=size(ii)
nii=sz(2)

; get sizes of things
sz=size(colblock)
ncol=sz(1)
if(ncol ne nii) then stop

; make colors, estimated errors from the magnitudes in colvec
cobs=fltarr(ncol)
sig2=fltarr(ncol)
badc=intarr(ncol)
for i=0,ncol-1 do begin
  cobs(i)=colvec(ii(0,i))-colvec(ii(1,i))
  sig2(i)=errvec(ii(0,i))^2+errvec(ii(1,i))^2
  if(errvec(ii(0,i)) gt 10. or errvec(ii(1,i)) gt 10.) then badc(i)=1
endfor

; make chi^2 and prior probability grids for a slice in teff, logg space
; at constant logz
trang=[3200.,12000.]
ltrang=alog10(trang)
dlt=.02
lgrang=[0.,6.]
dlg=.2
lzrang=[-.2,-.199]
dlz=1.
grid_likeli,cobs,sig2,rmag,gb,$
  ltrang,dlt,lgrang,dlg,lzrang,dlz,chi2,lprior,loglum,lto,lgo,lzo,$
  lto1,lgo1,lzo1,xform=xform,bias=bias

; adjust wfact to reflect min chi^2
wfact=(min(chi2)/10. > 1.)*wfact0

; find minimum of merit fn,and do 2D chunk in T,Z plane centered on it
subg1:
merit=chi2-lprior*wfact

; plotting stuff for IFA talk, sci team talk
;nb=where(sig2 ge .1,nnb)
;if(nnb le 2 and rmag gt 11.922 and rmag lt 11.924 and min(chi2) gt 2.) then begin
;  plots_ifa,lto,lgo,lzo,chi2,lprior,merit,'ifa_fig1.ps',0
;  plots_ifa,lto,lgo,lzo,chi2,lprior,merit,'ifa_fig2.ps',1
;  plots_ifa,lto,lgo,lzo,chi2,lprior,merit,'ifa_fig3.ps',2
;endif
; end plotting stuff

mm=min(merit,ix)
dlt=.02
ltrang=[lto1(ix)-7.*dlt,lto1(ix)+7.01*dlt]
dlg=.2
lgrang=[lgo1(ix),lgo1(ix)+.01*dlg]
dlz=.25
lzrang=[-4.,.99]

if(sflg eq 1) then stop

grid_likeli,cobs,sig2,rmag,gb,$
  ltrang,dlt,lgrang,dlg,lzrang,dlz,chi2,lprior,loglum,lto,lgo,lzo,$
  lto1,lgo1,lzo1,xform=xform,bias=bias

; find minimum again.  This time try a block that spans the range defined
; by  the larger of (range spanned by 10 smallest merit values)
;  or delta(merit) <= 10 plus a bit in T,Z, full range in g on coarse grid.
subg2:
wfact=(min(chi2)/10. > 1.)*wfact0
merit=chi2-lprior*wfact

; more plotting stuff
;nb=where(sig2 ge .1,nnb)
;if(nnb le 2 and rmag gt 11.922 and rmag lt 11.924 and min(chi2) gt 2.) then begin
;  plots_ifa,lto,lgo,lzo,chi2,lprior,merit,'ifa_fig4.ps',3
;  stop
;endif
; end plotting stuff

mm=min(merit,ix)
s=where(merit le mm+10.,ns)
so=sort(merit)
s1=so(0:9)                     ; 10 best merit values

dlt=.01
ltrang=[(min(lto1(s)) < min(lto1(s1)))-.02,(max(lto1(s)) > max(lto1(s1)))+.021]
ltrang=ltrang+.003
dlz=.125
lzrang=[(min(lzo1(s)) < min(lzo1(s1)))-.25,(max(lzo1(s)) > max(lzo1(s1)))+.26]
dlg=.5
lgrang=[0.,6.]+.03

if(sflg eq 1) then stop

grid_likeli,cobs,sig2,rmag,gb,$
  ltrang,dlt,lgrang,dlg,lzrang,dlz,chi2,lprior,loglum,lto,lgo,lzo,$
  lto1,lgo1,lzo1,xform=xform,bias=bias

; do this once more, searching a selected range in g, with slightly improved
; resolution in Z and much improved in g.
subg3:
wfact=(min(chi2)/10. > 1.)*wfact0
merit=chi2-lprior*wfact
mm=min(merit,ix)
s=where(merit le mm+10,ns)
so=sort(merit)
s1=so(0:9)                     ; 10 best merit values

dlt=.01
ltrang=[(min(lto1(s)) < min(lto1(s1)))-.02,(max(lto1(s)) > max(lto1(s1)))+.021]
ltrang=ltrang-.002
dlz=.1
lzrang=[(min(lzo1(s)) < min(lzo1(s1)))-.25,(max(lzo1(s)) > max(lzo1(s1)))+.26]
dlg=.1
lgrang=[(min(lgo1(s)) < min(lgo1(s1)))-.25,(max(lgo1(s)) > max(lgo1(s1)))+.26]
lgrang=lgrang-.02

if(sflg eq 1) then stop

grid_likeli,cobs,sig2,rmag,gb,$
  ltrang,dlt,lgrang,dlg,lzrang,dlz,chi2,lprior,loglum,lto,lgo,lzo,$
  lto1,lgo1,lzo1,xform=xform,bias=bias

; find the absolute minimum on this grid, and fit a general quadratic
; to the surrounding region.  find min of this function, and call that the
; answer.
isc=0                   ; counts scoots of box center point
wfact=(min(chi2)/10. > 1.)*wfact0
merit=chi2-lprior*wfact
mm=min(merit,ix)
subg4:
dlt=.005
ltrang=[lto1(ix)-.015,lto1(ix)+.016]
ltrang=ltrang+.0025/(1.+isc)
dlz=.05
lzrang=[lzo1(ix)-.2,lzo1(ix)+.21]
dlg=.05
lgrang=[lgo1(ix)-.15,lgo1(ix)+.16]
lgrang=lgrang-.025/(1.+isc)

if(sflg eq 1) then stop

grid_likeli,cobs,sig2,rmag,gb,$
  ltrang,dlt,lgrang,dlg,lzrang,dlz,chi2,lprior,loglum,lto,lgo,lzo,$
  lto1,lgo1,lzo1,xform=xform,bias=bias

; find the minimum on this grid, scoot box if min falls on a boundary.
; if no luck after 2 scoots, declare failure.
wfact=(min(chi2)/10. > 1.)*wfact0
merit=chi2-lprior*wfact
mm=min(merit,ix)
lt0=lto1(ix)
lg0=lgo1(ix)
lz0=lzo1(ix)
if(isc ge 4) then begin
  ferr=1
  sol=[lt0,lg0,lz0]
  serr=[9.99,9.99,9.99]            ; values indicate error
  c2m=chi2(ix)
  goto,windup
endif else begin
  if(lt0 eq min(lto1) or lt0 eq max(lto1) or lg0 eq min(lgo1) or $
     lg0 eq max(lgo1) or lz0 eq min(lzo1) or lz0 eq max(lzo1)) then begin
    isc=isc+1

; some debugging
  openw,iune,'~/kdbase/propdebug1.asc',/append,/get_lun
  scin=[0,0,0]
  if(lt0 eq min(lto1)) then scin(0)=1
  if(lt0 eq max(lto1)) then scin(0)=2
  if(lg0 eq min(lgo1)) then scin(1)=1
  if(lg0 eq max(lgo1)) then scin(1)=2
  if(lz0 eq min(lzo1)) then scin(2)=1
  if(lz0 eq max(lzo1)) then scin(2)=2
  printf,iune,isc,scin
  close,iune
  free_lun,iune

; done debugging 
    goto,subg4
  endif
endelse 

npt=n_elements(lto1)
funs=fltarr(npt,10)
x=lto1-lt0
y=lgo1-lg0
z=lzo1-lz0
funs(*,0)=1.
funs(*,1)=x
funs(*,2)=y
funs(*,3)=z
funs(*,4)=x^2
funs(*,5)=y^2
funs(*,6)=z^2
funs(*,7)=x*y
funs(*,8)=x*z
funs(*,9)=y*z
; weight the points near the minimum more strongly
met2=x^2/dlt^2 + y^2/dlg^2 + z^2/dlz^2
wts=1./(1.+met2)
cc=lstsqr(merit,funs,wts,10,rms,ch2,outp,1)  ; coeffs of quadratic fit
; setting derivs to zero gives linear eqns to solve for min
aa=[[2.*cc(4),cc(7),cc(8)],[cc(7),2.*cc(5),cc(9)],[cc(8),cc(9),2.*cc(6)]]
rhs=-[cc(1),cc(2),cc(3)]
ludc,aa,indexl
sol1=lusol(aa,indexl,rhs)
; reject this solution and set error flag if correction is more than 2.5
; grid spacing in any dimension.
if(abs(sol1(0)) le 2.5*dlt and abs(sol1(1)) le 2.5*dlg and abs(sol1(2)) le 2.5*dlz)$
        then begin
  sol=sol1+[lt0,lg0,lz0]
; compute 1-sigma formal errors, ignoring correlated errors
  serr=1./[sqrt(abs(cc(4)*wfact)),sqrt(abs(cc(5)*wfact)),sqrt(abs(cc(6)*wfact))]
  c2m=cc(0) > 0.
  ferr=0
endif else begin
  sol=[lt0,lg0,lz0]
  serr=[9.99,9.99,9.99]
  c2m=chi2(ix)
  ferr=2
  openw,iune,'~/kdbase/propdebug1.asc',/append,/get_lun
  scin=[0,0,0]
  isc=-5
  scin=[0.,0.,0.]
  if(abs(sol1(0) gt 2.5*dlt)) then scin(0)=sol1(0)
  if(abs(sol1(1) gt 2.5*dlg)) then scin(1)=sol1(1)
  if(abs(sol1(2) gt 2.5*dlz)) then scin(2)=sol1(2)
  printf,iune,isc,scin
  close,iune
  free_lun,iune
endelse

windup:
; compute value and errors for extinction
linked_props2,[sol(0)],[sol(1)],rmag,gb,loglf,logd,a_v
errav=sqrt(4.*serr(0)+.25*serr(1))/2.5        ; guess at error in a_v, in mag

; put everything in a form compatable with previous code
teff=10.^sol(0)
errteff=teff*(10.^serr(0)-1.)
logz=sol(2)
logg=sol(1)
logl=loglf
ext=a_v
props=[teff,logz,logg,ext]
properr=[errteff,serr(2),serr(1),errav]
chmin=c2m
rats=[ferr,isc,0.,0.]

if(sflg eq 1) then stop
;stop

end