pro fitstars,blockid,zdb=zdb ; This routine fits all of the non-standard stars that were observed during ; the block identified by blockid, deriving magnitudes and errors for all ; filters in which each star has ever been observed. If a given star has ; been observed previously, its new observations are combined with its old ; ones in the process of performing the fit. ; Results go into the 'stars' dbase. ; If keyword zdb is set, it overrides the ZDBASE environment variable print, '*****current version as of 051206*****' ;!EXCEPT=2 ; constants zdbase=getenv('ZDBASE') if(keyword_set(zdb)) then zdbase=zdb blkdbase=zdbase+'/survey/blocks/blocks' obspath=zdbase+'/survey/observations/' starpath=zdbase+'/survey/stars' imdbase=zdbase+'/survey/images/images' filtdir=['u','g','r','i','z','gr','d5'] nfilt=n_elements(filtdir) filtnam=['u ','g ','r ','i ','z ',$ 'Gred ','D51 '] filtstar=['u','g','r','i','z','Gred','D51'] ; filter names in stars dbase colstar=[['u','g'],['g','r'],['r','i'],['i','z'],['i','z'],['g','r'],$ ['g','r']] ; star dbase filters needed for colors sigstar=['sig_u','sig_g','sig_r','sig_i','sig_z','sig_Gred','sig_D51'] f1='(f13.5)' f2='(f10.5)' mt=40. ; threshold indicating bad magnitude values et=0.2 ; threshold for bad magnitude errors. qthrsh=5. ; threshold for bad data point rejection, in sigma dp=.22 ; generous half-width of detector FOV, for making pointings deg2rad=57.295 ; degrees per radian errmin=.001 ; min allowable photometric err in mags astrolib deferr=0.3 ; assumed systematic error if default extinction params used !priv=2 ; get block extinction parameters, times for all blocks dbopen,blkdbase dbext,[blockid],'jdmin,jdmax,iseqmin,iseqmax,valid,a,b,c,k',$ jdmin,jdmax,iseqmin,iseqmax,valid,a,b,c,k sjdmin=string(jdmin,format=f1) sjdmax=string(jdmax,format=f1) dbext,-1,'jdmin,jdmax,iseqmin,iseqmax,valid,a,b,c,k,chip2,chip3,chip4',$ jdmina,jdmaxa,iseqmina,iseqmaxa,valida,aa,ba,ca,ka,chip2a,chip3a,chip4a dbclose,blkdbase nblock=n_elements(jdmina) blockin=intarr(nblock) ; will indicate which blocks are repre ; sented in obs. chipa=fltarr(4,nfilt,nblock) chipa(1,*,*)=reform(chip2a,1,nfilt,nblock) ; reformat chip corrections chipa(2,*,*)=reform(chip3a,1,nfilt,nblock) chipa(3,*,*)=reform(chip4a,1,nfilt,nblock) ; make 'fixed' versions of extinction arrays. In these, if any coeff ; for one filter is flagged as bad, all coeffs for that filter take ; default values. Contents of nofit array for this block, filter = 1,else 0 aaf=aa baf=ba caf=ca kaf=ka chipaf=chipa nofit=intarr(nfilt,nblock) for i=0,nblock-1 do begin ; get current global parameters B and C, and cosmic colors rd_global,jdmina(i),gnfilt,gfilt_names,gparms,geparms,chipz,cosmic,xglob,$ zdb=zdb for j=0,nfilt-1 do begin if(aa(j,i) gt mt or ba(j,i) gt mt or ca(j,i) gt mt or ka(j,i) gt mt or $ max(chipa(*,j,i)) gt mt) then begin aaf(j,i)=geparms(0,j) baf(j,i)=gparms(0,j) caf(j,i)=gparms(1,j) kaf(j,i)=geparms(1,j) chipaf(1:3,j,i)=chipz(*,j) nofit(j,i)=1 endif endfor endfor ; get list of corresp images for input dbopen,imdbase st=dbfind(sjdmin+' < jd < '+sjdmax, /silent) dbext,st,'ra,dec',rahi,deci ; nominal center coords of all images dbclose,imdbase nim=n_elements(rahi) rai=rahi*15. ; make RA in decimal degrees ; identify all pointings, make list of unique possibly observed tiles ptra=fltarr(nim,4) ptdec=fltarr(nim,4) for i=0,nim-1 do begin ptra(i,*)=rai(i)+[-dp,-dp,dp,dp] ptdec(i,*)=deci(i)+[-dp,dp,-dp,dp] endfor ptra=fix(reform(ptra,4*nim)) ; a list of the integer RA & Dec of the ptdec=fix(reform(ptdec,4*nim)) ; corners of all the pointing fields tilenames,0,4*nim,ptra,ptdec,tnames tilenames,0,4*nim,ptra,ptdec,tonames,/obj ; make corresp tile names su=uniq(tnames,sort(tnames)) & tnames=tnames(su) ; star tnames suo=uniq(tonames,sort(tonames)) & tonames=tonames(suo) ; obs tnames ntiles=n_elements(tnames) ; loop over tiles for i=0,ntiles-1 do begin ; check to see that needed tile databases exist stardbase=starpath+'/'+tnames(i) fn=findfile(stardbase+'.dbd',count=fncount) if(fncount gt 0) then begin ; loop over filters for j=0,nfilt-1 do begin obsdbase=obspath+filtdir(j)+'/'+tonames(i) fn=findfile(obsdbase+'.dbd',count=obscount) if(obscount eq 1) then begin ; go to images dbase, find all images with this filter in this time range dbopen,imdbase st=dbfind(sjdmin+' < jd < '+sjdmax+',filter='+filtnam(j), /silent) dbclose,imdbase ; go to observations dbase for this tile, filter. Find all obs from the ; above list of images dbopen,obsdbase cto=0 sto=dbget('iseq',st,-1,count=cto, /silent) nobs=n_elements(sto) if(sto(0) gt 0 and nobs gt 0) then begin ; make list of unique stars dbext,sto,'starid',starid sustarid=uniq(starid,sort(starid)) ustarid=starid(sustarid) nstar=n_elements(ustarid) ; go back to observations dbase, find all observations for all of these stars sta=dbget('starid',ustarid,-1, /silent) dbext,sta,'starid,iseq,ra,dec,mag,err,x,chip,xcen,ycen,sky,',$ astarid,iseq,ra,dec,mag,err,x,chip,xcen,ycen,sky dbext,sta,'dra,ddec',dra,ddec dbclose,obsdbase soa=sort(astarid) astarid=astarid(soa) iseq=iseq(soa) ra=ra(soa) dec=dec(soa) mag=mag(soa) err=err(soa) err=err > errmin x=x(soa) chip=chip(soa) xcen=xcen(soa) ycen=ycen(soa) sky=sky(soa) dra=dra(soa) ddec=ddec(soa) uso=uniq(astarid) auso=[-1,uso] ; go back to images dbase, get corresp JD and error flags dbopen,imdbase dbext,iseq,'jd,cldflg,astflg,fitflg',jdobs,cldf,astf,fitf dbclose,imdbase ; go to stars dbase, get needed magnitudes and colors for all these stars dbopen,stardbase sts=dbmatch('starid',ustarid) sstr='entry,ra,dec,'+filtstar(j)+','+sigstar(j)+','+colstar(0,j)$ +','+colstar(1,j) dbext,sts,sstr,$ star_entry,ras,decs,mags,sigs,col1s,col2s dbclose,stardbase stcolor=col1s-col2s ; protect against NaNs in the stars data si=where(finite(mags) eq 0 or finite(stcolor) eq 0,nsi) if(nsi gt 0) then begin mags(si)=99.9 sigs(si)=99.9 stcolor(si)=99.9 endif ; fill these in with defaults where need be sbad=where(col1s gt mt or col2s gt mt,nsbad) if(nsbad gt 0) then stcolor(sbad)=cosmic(j) ssbad=where(sigs ge mt or sigs le 0.,nssbad) if(nssbad gt 0) then sigs(ssbad)=et ; step through stars, estimate magnitude, internal error, scatter for each magout=fltarr(nstar) errout=fltarr(nstar) ; raout=dblarr(nstar) ; decout=dblarr(nstar) draout=fltarr(nstar) ddecout=fltarr(nstar) scatout=fltarr(nstar) varout=intarr(nstar) jdout=dblarr(nstar) badfout=intarr(nstar) for m=0,nstar-1 do begin mm=auso(m)+1 nobs=auso(m+1)-mm+1 ; fit for magnitude, dra, ddec, error. if(nobs eq 1) then begin ; a common special case sb=where(iseqmina le iseq(mm) and iseqmaxa ge iseq(mm),nsb) sb=sb(0) ; if more than one block fits, take first mag0=mag(mm)-aaf(j,sb)-baf(j,sb)*stcolor(m)-$ kaf(j,sb)*(x(mm)-xglob)-$ caf(j,sb)*(stcolor(m)-cosmic(j))*(x(mm)-xglob)-$ chipaf(chip(mm)-1,j,sb) err0=err(mm) if(nofit(j,sb) ne 0) then err0=sqrt(err0^2+deferr^2) jd0=jdobs(mm) dra0=0. ddec0=0. scat0=0. badf0=2*cldf(mm)+astf(mm) ; if (j eq 2 and i eq 21) then stop endif else begin mes=fltarr(nobs) wts=fltarr(nobs) dras=fltarr(nobs) ddecs=fltarr(nobs) jd0s=dblarr(nobs) cld0s=intarr(nobs) ast0s=intarr(nobs) for nn=0,nobs-1 do begin sb=where(iseqmina le iseq(mm+nn) and iseqmaxa ge iseq(mm+nn),$ nsb) sb=sb(0) mes(nn)=mag(mm+nn)-aaf(j,sb)-baf(j,sb)*stcolor(m)-$ kaf(j,sb)*(x(mm+nn)-xglob)-$ caf(j,sb)*(stcolor(m)-cosmic(j))*(x(mm+nn)-xglob)-$ chipaf(chip(mm+nn)-1,j,sb) cld0s(nn)=cldf(mm+nn) ast0s(nn)=astf(mm+nn) wts(nn)=1./(err(mm+nn))^2 if(cld0s(nn) gt 0) then wts(nn)=wts(nn)/1000. if(ast0s(nn) gt 0) then wts(nn)=wts(nn)/100. dras(nn)=dra(mm+nn) ddecs(nn)=ddec(mm+nn) jd0s(nn)=jdobs(mm+nn) endfor twts=total(wts) mag0=total(mes*wts)/twts err0=sqrt(1./twts) if(nofit(j,sb) ne 0) then err0=sqrt(err0^2+deferr^2) dra0=total(dras*wts)/twts ddec0=total(ddecs*wts)/twts scat0=max(mes)-min(mes) resid0=mes-mag0 jd0=total(jd0s*wts)/twts badf0=2*min(cld0s)+min(ast0s) ; if enough obs, test for bad data and, if found, redo fit. if(nobs gt 3) then begin quartile,resid0,med,q,dq if (dq eq 0.) then begin residsd = stddev(resid0) if (residsd gt 0.) then dq=residsd*1.39 else dq=1. endif sq=where(abs(resid0/dq) ge qthrsh/1.39,nsq) sqp=where(abs(resid0/dq) lt qthrsh/1.39,nsqp) if(nsq gt 0) then begin wts(sq)=1.e-20 twts=total(wts) mag0=total(mes*wts)/twts err0=sqrt(1./twts) dra0=total(dras*wts)/twts ddec0=total(ddecs*wts)/twts scat0=max(mes)-min(mes) resid0=mes-mag0 jd0=total(jd0s*wts)/twts if(nsqp gt 0) then begin badf0=2*min(cld0s(sqp))+min(ast0s(sqp)) endif else begin badf0=1 endelse endif endif endelse ; estimate variability index if(scat0 le .07) then var0=0 if(scat0 gt .07 and scat0 le .3) then var0=1 if(scat0 gt .3) then var0=2 ; put fitted data into arrays magout(m)=mag0 errout(m)=err0 ; *** A problem!! this effectively ignores all filters but the last in ; determining the coordinate adjustments *** draout(m)=dra0/(3600.*15.*cos(decs(m)/deg2rad)) ddecout(m)=ddec0/3600. scatout(m)=scat0 varout(m)=var0 jdout(m)=jd0 badfout(m)=badf0 endfor ; end of loop over m=stars ; load results into stars database ; print,i,j,' '+stardbase dbopen,stardbase,1 dbext,star_entry,'badflag',bft bft(j,*)=reform(badfout,1,nstar) dbupdate,star_entry,filtstar(j)+','+sigstar(j)+',dra,ddec,var,jd,' + $ 'badflag',magout,errout,draout,ddecout,varout,jdout,bft,$ /silent dbclose,stardbase endif endif endfor ; end of loop over j=filters endif endfor ; end of loop over i=tiles print, string(7b) ;ring bell when finished end