pro cull_refstar,rac,decc
; This routine searches the tstds database for stars that are within one
; field width (12 arcmin) of the coords rac, decc (both in decimal degrees).
; All observations of such stars are retrieved from the tostds database,
; and statistics are accumulated giving number of observations in each
; color, along with fitted average magnitudes, extinction slopes, and rms
; around fits.  Stars with too few observations in 2 or more colors, or
; with discrepant values of slope or rms, are marked for removal from the
; standard star list.  Statistics are printed, and diagnostic plots are
; displayed.  The user is then queried as to whether the flagged stars
; should be removed.  If so, these stars are removed from the tstds database
; and placed instead in their appropriate target stars tiles.  Also, the
; standard magnitudes of the remaining stars are adjusted to match the
; fitted values at X = 1.215, and these new magnitudes replace the old
; ones in the tstds database.  No magnitude adjustment is done for the
; bad stars.

;zdbase=getenv('ZDBASE')
zdbase='/home/tbrown/dbase'
path=zdbase+'/survey/'
stdbase=path+'stars/tstds'
blockdbase=path+'blocks/blocks'
obspath=path+'/observations/'
stop
astrolib
!priv=3
radian=180./!pi
fhwid=0.2               ; half-width of field in degrees
cols=['u','g','r','i','z','gr','d5']
ncol=n_elements(cols)
x0=1.215                ; standard airmass
gfrac=0.67              ; min ok fraction of obs present
rmsthrsh=[.08,.05,.035,.035,.045,.08,.04]
;rmsthrsh=rmsthrsh*2.0

; get list of standard stars in field of interest
dbopen,stdbase
dra=fhwid/(15.*(cos(decc/radian)))
sdra=string(dra,format='(f7.5)')
sddec=string(fhwid,format='(f6.4)')
sra=string(rac/15.,format='(f8.5)')
sra='ra='+sra+'('+sdra+')'
sdec=string(decc,format='(f8.4)')
sdec='dec='+sdec+'('+sddec+')'
ll=dbfind(sra+','+sdec, /silent)
dbext,ll,'ra,dec,starid,u,g,r,i,z,gred,d51',ra,dec,sid,u,g,r,ii,z,gr,d5 
nst=n_elements(ra)
dbclose

; put magnitudes into an array for easy reference, make stats tables
magtab=fltarr(nst,ncol)
magtab(*,0)=u
magtab(*,1)=g
magtab(*,2)=r
magtab(*,3)=ii
magtab(*,4)=z
magtab(*,5)=gr
magtab(*,6)=d5
nobst=intarr(nst,ncol)
zptt=fltarr(nst,ncol)
slopet=fltarr(nst,ncol)
rmst=fltarr(nst,ncol)

; get night zero points for all blocks, make array of zero points vs seq number
; If multiple blocks point to the same sequence number, you get the last one
dbopen,blockdbase
dbext,-1,'iseqmin,iseqmax,a',ismin,ismax,aa
dbclose
msno=max(ismax)
sqno=lindgen(msno)+1
aaseq=fltarr(msno+1,ncol)
nblock=n_elements(ismin)
for i=0,nblock-1 do begin
  s=where(sqno ge ismin(i) and sqno le ismax(i),ns)
  if(ns gt 0) then aaseq(s,*)=rebin(reform(aa(*,i),1,ncol),ns,ncol)
endfor

; loop over colors, get all obs for each star in list
for i=0,ncol-1 do begin
  obsdbase=path+'observations/'+cols(i)+'/tostds'
  dbopen,obsdbase
  dbext,-1,'starid,iseq,mag,x',sida,isa,maga,xxa
  for j=0,nst-1 do begin
    ll=where(sida eq sid(j),nll)
    if(nll gt 0) then begin
      iss=isa(ll)
      mag=maga(ll)
      xx=xxa(ll)
      nobs=nll
      nobst(j,i)=nobs
      mago=mag-aaseq(iss,i)    ; obsd mag corrected for block zero pt.
      xxo=xx-x0                ; airmass relative to std. 
      if(nobs ge 5) then begin   ; do a fit if enough obs.
        cc=poly_fit(xxo,mago,1)
        yy=poly(xxo,cc)
        dif=mago-yy
        rms=stddev(dif)
        rmso=rms
        quartile,dif,med,q,dq  ; reject 5-sigma outliers
        so=where(abs(dif-med) le 5.*dq/1.349,nso)
        if(nso gt 5) then begin
          cc=poly_fit(xxo(so),mago(so),1)
          yy=poly(xxo(so),cc)
          dif=mago(so)-yy
          rms=stddev(dif)
        endif
        if(rms gt 0.) then begin
          zptt(j,i)=cc(0)-magtab(j,i)  ; correction to star's magnitude
          slopet(j,i)=cc(1)
          rmst(j,i)=rms
        endif else begin
          rmst(j,i)=-1.
        endelse
      endif else begin
        rmst(j,i)=-1.                ; rms=-1 signals insufficient data
      endelse
    endif
  endfor
  dbclose
 print,i
endfor

; decide which stars to pitch
scoren=fltarr(nst)
scorer=fltarr(nst)
mnm=intarr(ncol)
for i=0,ncol-1 do begin
  mnm=max(nobst(*,i))*gfrac
  sn=where(nobst(*,i) le mnm,nsn)
  if(nsn gt 0) then scoren(sn)=scoren(sn)+1.      ; bad if not observed enough
  sr=where(rmst(*,i) lt 0. or rmst(*,i) gt rmsthrsh(i),nsr)
  if(nsr gt 0) then scorer(sr)=scorer(sr)+1.      ; bad if scatter is big
endfor

sg=where(scoren le 2 and scorer le 2,nsg)         ; good stars
sb=where(scoren gt 2 or scorer gt 2,nsb)          ; bad stars
sidbad=sid(sb)
lindex=lindgen(nsb)
sidgood=sid(sg)

stop

; get info on good stars from stars dbase, adjust magnitudes, put them back
tstardbase=path+'/stars/tstds'
dbopen,tstardbase
lg=dbget('starid',sidgood, /silent)
ngd=n_elements(lg)
gindex=lindgen(ngd)
dbext,lg,'starid,ra,dec,dra,ddec,u,sig_u,g,sig_g',sids,ras,decs,dras,ddecs,$
       us,sig_us,gs,sig_gs
dbext,lg,'r,sig_r,i,sig_i,z,sig_z,d51,sig_d51,gred,sig_gred',rs,sig_rs,is,$
       sig_is,zs,sig_zs,d51s,sig_d51s,greds,sig_greds
dbext,lg,'J,H,K,id2mass,var,crwd,jd,source',js,hs,ks,id2s,vars,crwds,jds,srcs

; testing
gmrold=gs-rs
rmiold=rs-is
umgold=us-gs

szpt=zptt(sg,*)
sy=where(us lt 40,ny)
if(ny gt 0) then us(sy)=us(sy)+szpt(sy,0)
sy=where(gs lt 40,ny)
if(ny gt 0) then gs(sy)=gs(sy)+szpt(sy,1)
sy=where(rs lt 40,ny)
if(ny gt 0) then rs(sy)=rs(sy)+szpt(sy,2)
sy=where(is lt 40,ny)
if(ny gt 0) then is(sy)=is(sy)+szpt(sy,3)
sy=where(zs lt 40,ny)
if(ny gt 0) then zs(sy)=zs(sy)+szpt(sy,4)
sy=where(greds lt 40,ny)
if(ny gt 0) then greds(sy)=greds(sy)+szpt(sy,5)
sy=where(d51s lt 40,ny)
if(ny gt 0) then d51s(sy)=d51s(sy)+szpt(sy,6)

dbopen,tstardbase,1
dbupdate,lg,'u,g,r,i,z,gred,d51',us,gs,rs,is,zs,greds,d51s
dbclose

; get info on bad stars from stars dbase, put it into target tiles
tstardbase=path+'/stars/tstds'
dbopen,tstardbase
ll=dbget('starid',sidbad, /silent)
nst=n_elements(ll)
lindex=lindgen(nst)
dbext,ll,'starid,ra,dec,dra,ddec,u,sig_u,g,sig_g',sids,ras,decs,dras,ddecs,$
       us,sig_us,gs,sig_gs
dbext,ll,'r,sig_r,i,sig_i,z,sig_z,d51,sig_d51,gred,sig_gred',rs,sig_rs,is,$
       sig_is,zs,sig_zs,d51s,sig_d51s,greds,sig_greds
dbext,ll,'J,H,K,id2mass,var,crwd,jd,source',js,hs,ks,id2s,vars,crwds,jds,srcs
iras=fix(ras*15.)
idecs=fix(decs)
tilenames,0,nsb,iras,idecs,snames
so=sort(snames)
snames=snames(so)
su=uniq(snames)
usnames=snames(su)
nsu=n_elements(su)
if(nsu eq 1) then begin
  nbotu=[0]
endif else begin
  nbotu=[0,su(0:nsu-2)+1]
endelse

; sort all the data so tiles are together
sids=sids(so)
ras=ras(so)
decs=decs(so)
dras=dras(so)
ddecs=ddecs(so)
us=us(so)
sig_us=sig_us(so)
gs=gs(so)
sig_gs=sig_gs(so)
rs=rs(so)
sig_rs=sig_rs(so)
is=is(so)
sig_is=sig_is(so)
zs=zs(so)
sig_zs=sig_zs(so)
d51s=d51s(so)
sig_d51s=sig_d51s(so)
greds=greds(so)
sig_greds=sig_greds(so)
js=js(so)
hs=hs(so)
ks=ks(so)
id2s=id2s(so)
vars=vars(so)
crwds=crwds(so)
jds=jds(so)
srcs=srcs(so)
; (whew)

; loop over unique tiles, put data into correct databases
for i=0,nsu-1 do begin
  stardbase=path+'/stars/'+usnames(i)
  dbopen,stardbase,1
  s=lindex(nbotu(i):su(i))
  dbbuild,sids(s),ras(s),decs(s),dras(s),ddecs(s),us(s),sig_us(s),gs(s),$
     sig_gs(s),rs(s),sig_rs(s),is(s),sig_is(s),zs(s),sig_zs(s),$
     d51s(s),sig_d51s(s),greds(s),$
     sig_greds(s),js(s),hs(s),ks(s),id2s(s),vars(s),crwds(s),jds(s),srcs(s)

; sort into correct order
  dbopen,stardbase
  so=dbsort(-1,'ra,dec')
  dbext,so,'starid,ra,dec,dra,ddec,u,sig_u,g,sig_g',v1,v2,v3,v4,v5,v6,v7,v8,v9
  dbext,so,'r,sig_r,i,sig_i,z,sig_z,d51,sig_d51,gred,sig_gred',w1,w2,w3,w4,w5,$
       w6,w7,w8,w9,w10
  dbext,so,'J,H,K,id2mass,var,crwd,jd,source',x1,x2,x3,x4,x5,x6,x7,x8
  dbopen,stardbase,1
  dbupdate,-1,'starid,ra,dec,dra,ddec,u,sig_u,g,sig_g',v1,v2,v3,v4,v5,v6,v7,$
       v8,v9,v10
  dbupdate,-1,'r,sig_r,i,sig_i,z,sig_z,d51,sig_d51,gred,sig_gred',w1,w2,w3,w4,$
       w5,w6,w7,w8,w9,w10
  dbupdate,-1,'J,H,K,id2mass,var,crwd,jd,source',x1,x2,x3,x4,x5,x6,x7,x8
  dbindex
  dbclose                                                                
endfor

; remove data from tstds database
stardbase=path+'/stars/tstds'
dbopen,stardbase,1
dbdelete,ll
dbclose

; Now deal with observations databases
; for each color, get all observed quantities for bad stars
for i=0,ncol-1 do begin
  obsdbase=path+'observations/'+cols(i)+'/tostds'
  dbopen,obsdbase
  ll=dbget('starid',sidbad, /silent)        ; find all entries concerning bad stars
  nst=n_elements(ll)
  lindex=lindgen(nst)
  dbext,ll,'starid,ra,dec,iseq,filter,mag,err,x,chip,xcen',sido,rao,deco,$
      iseqo,filtero,mago,erro,xo,chipo,xceno
  dbext,ll,'ycen,sky,sharp,chi,dra,ddec',yceno,skyo,sharpo,chio,drao,ddeco
  dbclose

; make list of unique tile names, sort data so these are together
  rai=fix(rao*15.)
  deci=fix(deco)
  tilenames,0,nsb,rai,deci,tnames,/obj
  so=sort(tnames)
  tnames=tnames(so)
  rao=rao(so)
  deco=deco(so)
  iseqo=iseqo(so)
  filtero=filtero(so)
  mago=mago(so)
  erro=erro(so)
  xo=xo(so)
  chipo=chipo(so)
  xceno=xceno(so)
  yceno=yceno(so)
  skyo=skyo(so)
  sharpo=sharpo(so)
  chio=chio(so)
  drao=drao(so)
  ddeco=ddeco(so)
  su=uniq(tnames)
  ntu=n_elements(su)
  utnames=tnames(su)
  usnames=snames(su)
  if(ntu eq 1) then begin
    nbotu=[0]
  endif else begin
    nbotu=[0,su(0:ntu-2)+1]   ; indices for tile i run from nbotu(i) to su(i)
  endelse

; index over target tiles.  Open desired tile, add desired data
  for j=0,ntu-1 do begin
    obsdbase=path+'observations/'+cols(i)+'/'+utnames(j)
    dbopen,obsdbase,1
    s=lindex(nbotu(j):su(j))
    dbbuild,sido(s),rao(s),deco(s),iseqo(s),filtero(s),mago(s),erro(s),$
        xo(s),chipo(s),xceno(s),yceno(s),skyo(s),sharpo(s),chio(s),drao(s),$
        ddeco(s),/silent

; sort results back into correct order
    dbopen,obsdbase
    so=dbsort(-1,'ra,dec')
    dbext,so,'starid,ra,dec,iseq,filter,mag,err,x,chip',v1,v2,v3,v4,v5,v6,$
              v7,v8,v9
    dbext,so,'xcen,ycen,sky,sharp,chi,dra,ddec',w1,w2,w3,w4,w5,w6,w7
    dbopen,obsdbase,1
    dbupdate,-1,'starid,ra,dec,iseq,filter,mag,err,x,chip',v1,v2,v3,v4,v5,v6,$
              v7,v8,v9
    dbupdate,-1,'xcen,ycen,sky,sharp,chi,dra,ddec',w1,w2,w3,w4,w5,w6,w7
    dbindex

  endfor

; delete bad stars from tostds databases
  obsdbase=path+'observations/'+cols(i)+'/tostds'
  dbopen,obsdbase,1
  dbdelete,ll
  dbclose

endfor
      
end