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