function get_allobs,ra,dec,ids,filtnams,cosmic,$ xglob,ism,isx,extparms,chpparms,jd,mag,cflg,zdb=zdb ; this routine makes a time series of all observations of a given star, corrected ; for extinction. On input, ; ra = RA of desired star, in decimal degrees (used only to get correct tile) ; dec = Dec of desired star, ditto ; ids = starid of desired star (used to identify star) ; filtnams = name of desired filter ; cosmic,xglob = global cosmic color and zero airmass parms ; ism(nblk) = min sequence number for each extinction block ; isx(nblk) = max sequence number for each extinction block ; extparms(4,7,nblk) = extinction a,b,c,k values for each filter ; chpparms(4,7,nblk) = chip parms for each filter (uniformly zero as this is written). ; On return: ; jd(nobs) = dates of observations ; mag(nobs) = extinction-corrected magnitudes. ; cflg(nobs) = cloud flag for obs (0 = okay, not zero = cloudy) ; Returned value is the number of time samples found. ; constants zdbase=getenv('ZDBASE') if(keyword_set(zdb)) then zdbase=zdb path=zdbase+'/survey/' imdbase=path+'images/images' blockdbase=path+'/blocks/blocks' starpath=path+'stars/' obspath=path+'observations/' filtdir=['u','g','r','i','z','gr','d5'] colstar=[['u','g'],['g','r'],['r','i'],['i','z'],['i','z'],['g','r'],$ ['g','r']] filtnam=['u ','g ','r ','i ','z ',$ 'Gred ','D51 '] ;astrolib nfilt=n_elements(filtdir) rad=180.d00/!pi ; degrees per radian ; get tile required, make db pathnames tilenames,0,1,fix(15.*ra),fix(dec),stile tilenames,0,1,fix(15.*ra),fix(dec),otile,/obj stardb=starpath+stile for i=0,6 do begin if(strtrim(filtnams,2) eq strtrim(filtnam(i),2)) then begin ifilt=i endif endfor obsdbase=obspath+filtdir(ifilt)+'/'+otile(0) ; open stars dbase, find star, make colors. name=findfile(stardb+'.dbh',count=exists) if(not exists) then begin nobs=0 jd=[-1.d00] mag=[99.99] cflg=[1] goto,fini endif dbopen,stardb ll=dbcircle(ra,dec,.05,dis,/silent) ; see if this is fast because of indexed search dbext,ll,'starid,u,g,r,i,z,gred,d51',sid,uu,gg,rr,ii,zz,gred,d51 s=where(sid eq ids,ns) if(ns eq 1) then begin colors=[uu(s)-gg(s),gg(s)-rr(s),rr(s)-ii(s),ii(s)-zz(s),ii(s)-zz(s),gg(s)-rr(s),$ gg(s)-rr(s)] endif else begin nobs=0 jd=[-1.d00] mag=[99.99] cflg=[1] goto,fini ; no unambiguous match, so bail out. endelse dbclose ; look to see if observations dbase exists. If not, set output parms and bail. name=findfile(obsdbase+'.dbh',count=exists) if(not exists) then begin nobs=0 jd=[-1.d00] mag=[99.99] cflg=[1] goto,fini endif ; open observations dbase, get all obs for star dbopen,obsdbase ll=dbcircle(ra,dec,.05,dis,/silent) ; see if this is fast because of indexed search dbext,ll,'starid,iseq,mag,x,chip',sido,iseqo,mago,xo,chipo dbclose s1=where(sido eq ids,ns1) if(ns1 gt 0) then begin ll=ll(s1) sido=sido(s1) iseqo=iseqo(s1) mago=mago(s1) xo=xo(s1) chipo=chipo(s1) nobs=ns1 endif else begin nobs=0 jd=[-1.d00] mag=[99.99] cflg=[1] goto,fini endelse ; open images database, get jd values to go with sequence numbers dbopen,imdbase dbext,iseqo,'jd,cldflg',jdo,cfgo dbclose jd=jdo cflg=cfgo ; make block number for each observation blkno=lonarr(nobs) for i=0,nobs-1 do begin sb=where(ism le iseqo(i) and isx ge iseqo(i),nsb) if(nsb ge 0) then begin blkno(i)=sb(nsb-1) endif else begin print,'sequence number with no associated block!', iseqo(i) endelse endfor ; make arrays of a,b,c,k,chipcor for each observation aa=reform(extparms(0,ifilt,blkno)) bb=reform(extparms(1,ifilt,blkno)) cc=reform(extparms(2,ifilt,blkno)) kk=reform(extparms(3,ifilt,blkno)) chipcor=reform(chpparms(chipo,ifilt,blkno)) ; make extinction-corrected magnitudes mag=mago-aa-bb*colors(ifilt)-kk*(xo-xglob)-cc*(colors(ifilt)-cosmic(ifilt))*$ (xo-xglob)-chipcor fini: return,nobs end pro pointing_eval1,filin,pointid,npoint,nsamp,rmspnt,rmsall,good,zdb=zdb ; This routine attempts to assess the quality of data for each pointing ; described in the ascii input file filin. ; *** This version of the routine, modified at Dave Latham's behest, tries ; to compute variability of the r magnitude and of colors g-r, r-i, ; i-z, and g-d51, not merely the individual magnitudes. ; *** ; For each pointing, it computes the extinction-corrected magnitudes for ; a sample of 20 bright (but not too bright) stars near the center of the ; pointing. For each filter it then computes the number of distinct ; times (separated by 1 hour or more) that the pointing has been observed, ; and the typical scatter among pointings (if more than 1) for each star. ; Because of long-short exposure pairs, for some filters there are ; more exposures than pointings. Statistics for the number and scatter ; of all exposures are also returned. ; On return, ; pointid(npt) = ID number of pointint, encoding RA,Dec eg 39019408 ; npoint(nfilt,npt) = number of distinct visits (> 1 hour separation) ; nsamp(nfilt,npt) = number of different images at this pointing ; rmsall(nfilt,npt) = typical scatter among all images taken at this pointing ; rmspnt(nfilt,npt) = typical scatter when results for each visit are averaged. ; good(npt) = a quality index for the pointing, defined as follows: ; -1 = no data for this pointing ; 0 = data for 1 visit but not including all of g,r,i,z,D51 filters ; 1 = all filters present for visit; no 2nd visit in any filter ; 2 = at least 2 visits, but some filters have only 1 visit ; 3 = all filters present for 2 or more visits, but typical scatter ; values for some or all filters exceed norms. ; 4 = all filters present for 2 or more visits, typical scatter values ; are within norms ; constants zdbase=getenv('ZDBASE') if(keyword_set(zdb)) then zdbase=zdb path=zdbase+'/survey/' imdbase=path+'images/images' blockdbase=path+'/blocks/blocks' starpath=path+'stars/' obspath=path+'observations/' filtdir=['u','g','r','i','z','gr','d5'] colstar=[['u','g'],['g','r'],['r','i'],['i','z'],['i','z'],['g','r'],$ ['g','r']] filtnam=['u ','g ','r ','i ','z ',$ 'Gred ','D51 '] astrolib nfilt=n_elements(filtdir) rad=180.d00/!pi ; degrees per radian dw=0.30 ; edge length of box centered on pointings to search nwant=20 ; want this many stars per pointing for robust estimates jdnow=2454000.d00 ; JD to get global parms dt=1./24. ; images at least 1 hour apart constitute separate visits rmsthrsh=.03 ; threshold between good/bad pointing repeatability ; get global parameters. choose a JD that gives parms for KeplerCam, since the cosmic ; colors and xglob never change, ; chip corrections are all zero, and the extinction parms all come out of the ; blocks database anyway. rd_global,jdnow,gnfilt,gfilt_names,gparms,geparms,chipz,cosmic,xglob,zdb=zdb ; get extinction parms for all blocks (not just ones in 'blocks' input parm) dbopen,blockdbase dbext,-1,'iseqmin,iseqmax',ism,isx nblka=n_elements(ism) extparms=fltarr(4,nfilt,nblka) chpparms=fltarr(4,nfilt,nblka) for i=0,nblka-1 do begin dbext,i+1,'a,b,c,k,chip2,chip3,chip4',a,b,c,k,c2,c3,c4 extparms(0,*,i)=reform(a,1,nfilt) extparms(1,*,i)=reform(b,1,nfilt) extparms(2,*,i)=reform(c,1,nfilt) extparms(3,*,i)=reform(k,1,nfilt) chpparms(1,*,i)=reform( c2,1,nfilt) chpparms(2,*,i)=reform( c3,1,nfilt) chpparms(3,*,i)=reform( c4,1,nfilt) endfor ; open input file, read the list of pointings and ra, dec values openr,iun,filin,/get_lun ss='' nline=0 while(not eof(iun)) do begin readf,iun,ss nline=nline+1 endwhile point_lun,iun,0 pid=strarr(nline) rapt=dblarr(nline) decpt=dblarr(nline) f1='(i4,1x,a8,6x,a10,3x,a8)' s1='' s2='' s3='' for i=0,nline-1 do begin readf,iun,n1,s1,s2,s3,format=f1 pid(i)=s1 hh=double(strmid(s2,0,2)) hm=double(strmid(s2,3,2)) hs=double(strmid(s2,6,4)) dd=double(strmid(s3,0,2)) dm=double(strmid(s3,3,2)) ds=double(strmid(s3,6,2)) rapt(i)=15.*(hh+hm/60.+hs/3600.) ; RA in decimal degrees decpt(i)=dd+dm/60.+ds/3600. ; Dec in decimal degrees endfor close,iun free_lun,iun ; loop over pointings npoint=lonarr(nfilt,nline) nsamp=lonarr(nfilt,nline) rmsall=fltarr(nfilt,nline) rmspnt=fltarr(nfilt,nline) good=intarr(nline) for ipt=0,nline-1 do begin ; determine tiles that are included in this pointing decptp=decpt(ipt)+dw/2. idecptp=fix(decptp) decptm=decpt(ipt)-dw/2. idecptm=fix(decptm) raptp=rapt(ipt)+dw/(2.*cos(decpt(ipt)/rad)) iraptp=fix(raptp) raptm=rapt(ipt)-dw/(2.*cos(decpt(ipt)/rad)) iraptm=fix(raptm) tilecr=[iraptm,iraptm,iraptp,iraptp] tilecd=[idecptm,idecptp,idecptm,idecptp] tilenames,0,4,tilecr,tilecd,tnames tilenames,0,4,tilecr,tilecd,tonames,/obj so=sort(tnames) tnames=tnames(so) tonames=tonames(so) tu=uniq(tnames) tunames=tnames(tu) tuonames=tonames(tu) ntile=n_elements(tu) ; find all stars in each of the target tiles that are in right brightness ; range. Concatenate lists for different tiles. idlist=[''] rmag=[0.] ras=[0.d00] decs=[0.d00] for i=0,ntile-1 do begin dbname=starpath+tunames(i) name=findfile(dbname+'.dbh',count=exists) if(not exists) then begin print,'notile' goto,notile endif dbopen,dbname ll=dbfind('10.5 < r < 15.',/silent) if(ll(0) ge 0) then begin dbext,ll,'starid,r,ra,dec',sid,r,ra,dec idlist=[idlist,sid] rmag=[rmag,r] ras=[ras,ra] decs=[decs,dec] endif dbclose notile: endfor ; extract only the stars within ra, dec ranges sc=where(15.*ras ge raptm and 15.*ras le raptp and decs ge decptm and $ decs le decptp,nsc) if(nsc gt 0) then begin rmag=rmag(sc) idlist=idlist(sc) ras=ras(sc) decs=decs(sc) endif else begin ; stop goto,nostars endelse ; sort this list according to brightness. Choose top nwant, or fail. nstar=n_elements(rmag)-1 ;********testing********** if(nstar lt nwant) then print,'nostars! ',rapt(ipt),decpt(ipt) ;************************** if(nstar lt nwant) then goto,nostars ; bail if too few stars sob=sort(rmag) rmag=rmag(sob(0:nwant-1)) idlist=idlist(sob(0:nwant-1)) ras=ras(sob(0:nwant-1)) decs=decs(sob(0:nwant-1)) ; loop over stars nobsi=lonarr(nwant,nfilt) nvisi=lonarr(nwant,nfilt) rmsi=fltarr(nwant,nfilt) rmsp=fltarr(nwant,nfilt) for istar=0,nwant-1 do begin ; loop over filters, making sequential lists of ifilt, jd, mag jdv=[0.d0] magv=[0.] ifiltv=[0] for ifilt=0,nfilt-1 do begin ; get a time series for each star. rass=ras(istar) decss=decs(istar) idlists=idlist(istar) filtnams=filtnam(ifilt) nobs=get_allobs(rass,decss,idlists,filtnams,cosmic,xglob,$ ism,isx,extparms,chpparms,jd,mag,cfg) ; reject obs for which cloud flag is set, redefine nobs scf=where(cfg eq 0,nobs) if(nobs gt 0) then begin jd=jd(scf) mag=mag(scf) cfg=cfg(scf) endif if(nobs gt 0) then begin jdv=[jdv,jd] magv=[magv,mag] ifiltv=[ifiltv,intarr(nobs)+ifilt] endif nobsi(istar,ifilt)=nobs endfor nts=n_elements(jdv)-1 if(nts gt 0) then begin jdv=jdv(1:*) magv=magv(1:*) ifiltv=ifiltv(1:*) endif else begin goto,nodat endelse ; soj=sort(jd) ; sort into time order ; jd=jd(soj) ; mag=mag(soj) ; if(nobs eq 0) then goto,endfilt ; first tag each observation with a visit number soj=sort(jdv) magv=magv(soj) ifiltv=ifiltv(soj) jdv=jdv(soj) nptt=n_elements(jdv) ivi=lonarr(nptt)-1 ii=0 jd0=jdv(0) for iv=0,nptt-1 do begin delt=jdv(iv)-jd0 if(abs(delt) le dt) then begin ivi(iv)=ii endif else begin ii=ii+1 jd0=jdv(iv) ivi(iv)=ii endelse endfor for ifilt=0,nfilt-1 do begin sf=where(ifiltv eq ifilt,nsf) if(nsf gt 0) then begin uvisi=uniq(ivi(sf)) nvisi(istar,ifilt)=n_elements(uvisi) ; num of distinct visits per star,filter endif else begin nvisi(istar,ifilt)=0 endelse endfor nvismax=max(ivi)+1 ; num of visits to current star in all filters together ; make averages of magnitudes in separate visits magavg=fltarr(nfilt,nvismax)+99.9 for ifilt=0,nfilt-1 do begin for ivis=0,nvismax-1 do begin s=where(ivi eq ivis and ifiltv eq ifilt,ns) if(ns gt 0) then magavg(ifilt,ivis)=total(magv(s))/ns endfor endfor ; make photometric indices for each visit. Indices in order r, g-r, r-i, i-z, g-d51 pindex=fltarr(nfilt,nvismax)+99.9 pindex(0,*)=magavg(2,*) ; r mag difs=[[1,2],[2,3],[3,4],[1,6]] ; filter indices of desired differences for iiv=0,nvismax-1 do begin for idd=0,3 do begin if(magavg(difs(0,idd),iiv) lt 40 and magavg(difs(1,idd),iiv) lt 40) then begin pindex(idd+1,iiv)=magavg(difs(0,idd),iiv)-magavg(difs(1,idd),iiv) endif else begin pindex(idd+1,iiv)=99.9 endelse endfor endfor ; estimate data dispersion among observations. Measure used depends on nobs. nobs=0 => 0. for ifilt=0,nfilt-1 do begin sobs=where(ifiltv eq ifilt and magv lt 40,nobs) if(nobs gt 0) then mag=magv(sobs) if(nobs eq 1) then rmsi(istar,ifilt)=0. if(nobs eq 2) then rmsi(istar,ifilt)=abs(mag(1)-mag(0))/sqrt(2.) if(nobs eq 3) then begin som=sort(mag) magso=mag(som) rmsi(istar,ifilt)=min([magso(1)-magso(0),magso(2)-magso(1)]) endif if(nobs ge 3 and nobs le 6) then rmsi(istar,ifilt)=stddev(mag) if(nobs gt 6) then begin quartile,mag,med,q,dq rmsi(istar,ifilt)=dq/1.349 endif endfor ; estimate data dispersion among visits, as above for ifilt=0,4 do begin sii=where(pindex(ifilt,*) lt 40.,nii) if(nii gt 0) then magt=reform(pindex(ifilt,sii)) if(nii eq 1) then rmsp(istar,ifilt)=0. if(nii eq 2) then rmsp(istar,ifilt)=abs(magt(1)-magt(0))/sqrt(2.) if(nii eq 3) then begin sov=sort(magt) magsv=magt(sov) rmsp(istar,ifilt)=min([magsv(1)-magsv(0),magsv(2)-magsv(1)]) endif if(nii ge 4 and ii le 6) then rmsp(istar,ifilt)=stddev(magt) if(nii gt 6) then begin quartile,magt,med,q,dq rmsp(istar,ifilt)=dq/1.349 endif endfor nodat: endfor ; stars loop ; pointid(npt) = ID number of pointint, encoding RA,Dec eg 39019408 ; npoint(nfilt,npt) = number of distinct visits (> 1 hour separation) ; nsamp(nfilt,npt) = number of different images at this pointing ; rmsall(nfilt,npt) = typical scatter among all images taken at this pointing ; rmspnt(nfilt,npt) = typical scatter when results for each visit are averaged. ; good(npt) = a quality index for the pointing, defined as follows: ; -1 = no data for this pointing ; 0 = data for 1 visit but not including all of g,r,i,z,D51 filters ; 1 = all filters present for visit; no 2nd visit in any filter ; 2 = at least 2 visits, but some filters have only 1 visit ; 3 = all filters present for 2 or more visits, but typical scatter ; values for some or all filters exceed norms. ; 4 = all filters present for 2 or more visits, typical scatter values ; are within norms ; compute stats for this pointing for j=0,nfilt-1 do begin npoint(j,ipt)=median(nvisi(*,j)) nsamp(j,ipt)=median(nobsi(*,j)) s2=where(rmsi(*,j) gt 0.,ns2) if(ns2 ge 3) then rmsall(j,ipt)=median(rmsi(s2,j)) s3=where(rmsp(*,j) gt 0.,ns3) if(ns3 ge 3) then rmspnt(j,ipt)=median(rmsp(s3,j)) endfor ; messy logic to set values in good array if(max(npoint(*,ipt)) le 0) then good(ipt)=-1 sgf=[1,2,3,4,6] ; corresp [g,r,i,z,d51] filters sgf1=[0,1,2,3,4] ; corresp [r,g-r,r-i,i-z,g-d51] if(max(npoint(sgf,ipt)) eq 1 and min(npoint(sgf,ipt)) le 0) then good(ipt)=0 if(max(npoint(sgf,ipt)) eq 1 and min(npoint(sgf,ipt)) eq 1) then good(ipt)=1 if(max(npoint(sgf,ipt)) ge 2 and min(npoint(sgf,ipt)) eq 1) then good(ipt)=2 if(max(npoint(sgf,ipt)) ge 2 and min(npoint(sgf,ipt)) ge 2) then begin if(max(rmspnt(sgf1,ipt)) ge rmsthrsh and min(rmspnt(sgf1,ipt)) gt 0.) $ then good(ipt)=3 else good(ipt)=4 endif nostars: if(max(npoint(*,ipt)) le 0) then good(ipt)=-1 if((ipt mod 10) eq 0) then print,ipt endfor ; pointings loop pointid=pid end