pro m67_color_fit,names,aa,bb,rms,gmrf,dataf,ngood,zdb=zdb,xform=xform
; This routine uses m67_basel_colors to get synthetic colors, RA, Dec for
; all of Latham's stars in M67.  It then matches the star positions to
; those in the stars database to obtain magnitudes for all matching
; stars in all available filters.  It creates the colors
; u-g (?), g-r, r-i, i-z, z-J, J-H, H-K, g-Gred (?), and g-D51
; for the observed stars, and for each color it performs a robust fit
; to determine coefficients aa,bb such that, as nearly a possible,
;   aa_(g-r) + bb_(g-r)*((g-r_bas)-.5) + (g-r_bas) = (g-r_obs)
;   aa_(r-i) + bb_(r-i)*((g-r_bas)-.5) + (r-i_bas) = (r-i_obs)
;                     etc.
; Thus, if the synthetic colors are an optimal representation of the observed
; colors, then aa = bb = 0 for all colors.
; On return, rms(ncolor) contains the robust rms magnitude of the fit
; residuals.
; Array dataf contains the observed colors for stars that were retained
; in the fit for each color, gmrf contains the basel g-r colors for these stars,
; and ngood(9) contains the number of retained stars.
; The routine makes a plot of the fitted values vs g-r for each color,
; overplotted with the fitted function.  The plot is entitled Color_Fit.ps
; If keyword zdb is set, its value overrides zdbase
; If keyword xform is set, the synthetic colors are transformed using
; the coefficients found in file xform before fits are done.
;
;*******************************
; Hacked to provide only constant offsets, not linear fits, 7 July 2008
; This version used for colblock_cas3 - related reduction, at least.
;*******************************

; constants
zdbase=getenv('ZDBASE')
if(keyword_set(zdb)) then zdbase=zdb
starpath=zdbase+'/survey/stars/'
m67db=['t131+010','t131+011','t131+012','t132+010','t132+011','t132+012',$
       't133+010','t133+011','t133+012']  ; tiles to search
ntiles=9
rcap=10./60.        ; capture radius = 2 arcsec expressed in arcmin
astrolib
common models,colblock,teffu,logzu,loggu
;restore,'~/d/basel/color_block1.sav'
;restore,'~/d/basel/colblock_castelli_smth.sav'
restore,'~/d/basel/colblock_cas3_smth.sav'
cblock=colblock

; get the Latham/Montgomery-derived properties of the M67 stars.
m67_basel_colors,ra,dec,v,bmv,teff,logg,bas_colors,names,/mix
bas_colors=transpose(bas_colors)
if(keyword_set(xform)) then xform_colors,bas_colors,xform,bas_colors
bas_colors=transpose(bas_colors)
nst=n_elements(ra)

; The bas_colors have untraceable provenance, so ignore them and create new ones
; using colors_blk_intrp.  This is more consistent, anyway, because this routine
; will be used to provide colors for star classification.
nst=n_elements(teff)
logzi=fltarr(nst)-.1             ; appropriate  for M67
colors_blk_intrp,cblock,teffu,loggu,logzu,teff,logg,logzi,colsout,$
      xform=xform
red_vector,3.1*.09,dcolor,xform=xform
stop
colsred=colsout+rebin(dcolor,9,nst)
bas_colors1=transpose(colsred)

stop

; find matching stars in the 'stars' database
mago=fltarr(nst,10)-1.
for i=0,ntiles-1 do begin
  path=starpath+m67db(i)
  dbopen,path
  for j=0,nst-1 do begin
    ll=dbcircle(ra(j)/15.,dec(j),rcap,dis)
    if(ll(0) gt 0) then begin
      dbext,ll(0),'u,g,r,i,z,j,h,k,gred,d51',u,g,r,ii,z,j2m,h2m,k2m,gred,d51
      mago(j,*)=reform([u,g,r,ii,z,j2m,h2m,k2m,gred,d51],1,10)
    endif
  endfor
endfor

; strip out unmatched stars
s=where(mago(*,1) gt 0.,ns)
if(ns gt 0) then begin
  mago=mago(s,*)
  ra=ra(s)
  dec=dec(s)
  v=v(s)
  bmv=bmv(s)
  teff=teff(s)
  logg=logg(s)
  bascm=bas_colors1(s,*)
endif
      
; make observed colors corresp to synthetic ones
obsc=fltarr(ns,9)
obsc(*,0)=mago(*,0)-mago(*,1)
obsc(*,1)=mago(*,1)-mago(*,2)
obsc(*,2)=mago(*,2)-mago(*,3)
obsc(*,3)=mago(*,3)-mago(*,4)
obsc(*,4)=mago(*,4)-mago(*,5)
obsc(*,5)=mago(*,5)-mago(*,6)
obsc(*,6)=mago(*,6)-mago(*,7)
obsc(*,7)=mago(*,1)-mago(*,8)
obsc(*,8)=mago(*,1)-mago(*,9)
s1=where(obsc lt (-2.) or obsc gt 40. or obsc eq 0.,ns1)
if(ns1 gt 0) then obsc(s1)=99.9

; do the fits
aa=fltarr(9)
bb=fltarr(9)
rms=fltarr(9)
gmrf=fltarr(300,9)
dataf=fltarr(300,9)
ngood=lonarr(9)
for i=0,8 do begin
  s=where(obsc(*,i) lt 5. and obsc(*,i) gt (-2.),ns)
  if(ns gt 5) then begin
    xx=bascm(s,1)-0.5
    yy=obsc(s,i)-bascm(s,i)
    wt=fltarr(ns)+1.
    funs=fltarr(ns,2)
    funs(*,0)=1.
    funs(*,1)=xx
;   cc=lstsqr(yy,funs,wt,2,rmsf,chisq,resid,1)

    cc=[mean(yy)]
    resid=yy-cc(0)

    quartile,resid,med,q,dq
    s1=where(abs(resid/dq) ge 4./1.349,ns1)
    if(ns1 gt 0) then wt(s1)=0.
;   cc=lstsqr(yy,funs,wt,2,rmsf,chisq,resid,1)

    sgg=where(wt gt 0.,nsgg)
    cc=[mean(yy(sgg))]
    resid=yy-cc(0)
    rmsf=stddev(resid(sgg))

    aa(i)=cc(0)
;   bb(i)=cc(1)
    bb(i)=0.
    rms(i)=rmsf
    s2=where(wt ne 0.,ns2)
    dataf(0:ns2-1,i)=yy(s2)
    gmrf(0:ns2-1,i)=funs(s2,1)+.5
    ngood(i)=ns2
  endif
endfor

; make the plot
psll,name='Color_Fit_cas3.ps'
xtit='Basel g-r color'
xr=[0.2,1.4]
yr=[-0.25,.25]
stop
!p.multi=[0,3,3]
for i=0,8 do begin
  nn=ngood(i)-1
  if(nn gt 10) then begin
    plot,gmrf(0:nn,i),dataf(0:nn,i),psym=1,xtit=xtit,ytit=names(i),xran=xr,$
      yran=yr,/xsty,/ysty,charsiz=1.5,symsiz=0.7
    oplot,gmrf(0:nn,i),aa(i)+bb(i)*(gmrf(0:nn,i)-0.5)
  endif else begin
    plot,gmrf(0:10,1),dataf(0:10,1),xtit=xtit,ytit=names(i),xran=xr,yran=yr,$
      /xsty,/ysty,/nodata,charsiz=1.5
    xyouts,0.60,0.,'No Data',charsiz=1.5
  endelse
endfor
psend
sun
!p.multi=[0,1,1]

stop

end