pro plot_refvsx,blocks,rf,zdb=zdb,ps=ps,log=log
; This routine makes postscript plots of the observed magnitudes of reference
; stars, uncorrected for atmospheric extinction, as a function of airmass.
; It does this for every filter for which observations exist, for every
; block listed in input array blocks.
; If rf=0, it plots all reference stars that are observed at every opportunity
; in the given block and filter
; If rf=1, it plots only stars in the R19 field
; If rf=2, it plots only stars in the NGC6811 field.
; If keyword ps is set, ps and gif files are created
; If keyword log is set, values relating to the fits are written to a log
; file residing in a standard directory
print, '*****current version as of 071210*****'
; Trunk directories where PS and GIF files are created. Make sure they exist!
; Likewise for logfile directory!
psdir='/mnt/server/tbrown/dbase/diagnostics/plots_ps/'
gifdir='/mnt/server/tbrown/dbase/diagnostics/plots_gif/'
;psdir='/knobos/e/kolinski/kepler_diagnostics/plots_ps/'
;gifdir='/knobos/e/kolinski/kepler_diagnostics/plots_gif/'
logdir='/mnt/server/tbrown/dbase/diagnostics/log/'
psflag = keyword_set(ps)
; constants
zdbase=getenv('ZDBASE')
if(keyword_set(zdb)) then zdbase=zdb
blkdbase=zdbase+'/survey/blocks/blocks'
obspath=zdbase+'/survey/observations/'
stardbase=zdbase+'/survey/stars/tstds'
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
ral=[0.,19.5208,19.6367,8.8561]
decl=[0.,36.1167,46.5667,11.8333]
oname=['/tostds','/tostds','/tostds','/to132+011']
sname=['/tstds','/tstds','/tstds','/t132+011']
nblk=n_elements(blocks)
; loop over blocks, get range of images in block
if(keyword_set(log)) then begin
openw,iunl,logdir+'refvsxlog_test',/get_lun,/append
sst=systime()
printf,iunl,'Created: '+sst
endif
for i=0,nblk-1 do begin
ib=blocks(i)
print,'block = ',ib
dbopen,blkdbase
dbext,[ib],'iseqmin,iseqmax,jdmin',ismin,ismax,jdmin
dbclose
;Form date name and create date-named directory for plot repository
;Dave Latham requested plot date change (subtracted 1 from jdmin)
CALDAT, (jdmin-1.0D), mon, day, yr
yr = STRTRIM(STRING(yr[0]),2)
mon = STRTRIM(STRING(mon[0]),2)
IF (STRLEN(mon) LT 2) THEN mon='0'+mon
day = STRTRIM(STRING(day[0]),2)
IF (STRLEN(day) LT 2) THEN day='0'+day
sdatenam=yr+mon+day ;for log file
IF (psflag) THEN BEGIN
pspath = psdir+yr+mon+day+'/'
file0 = FINDFILE(pspath, count=ct)
IF (ct EQ 0) THEN BEGIN
str0 = 'mkdir '+pspath
SPAWN, str0
ENDIF
gifpath = gifdir+yr+mon+day+'/'
file1 = FINDFILE(gifpath, count=ct)
IF (ct EQ 0) THEN BEGIN
str1 = 'mkdir '+gifpath
SPAWN, str1
ENDIF
psfile = pspath+yr+mon+day+'_rvx_b'+strtrim(string(ib),2)+'_test.ps'
set_plot, 'ps'
DEVICE, /helvetica, bits_per_pixel=8, filename=psfile, /landscape, $
/inches, xsize=10, ysize=8, xoffset=.5, yoffset=10.5
set_plot, 'x'
ENDIF
; loop over filters, find observations in tostds in each filter, in
; right range of iseq numbers, in right RA, Dec range if rf ne 0.
for j=0,nfilt-1 do begin
print,'filtno = ',j
window, 5, retain=2, xsize=800, ysize=650 ;open plot window before loop
obsdbase=obspath+filtdir(j)+oname(rf)
ierr=findfile(obsdbase+'.dbd')
berr=byte(ierr)
if(berr(0) eq 0) then goto,endloop ; bail if target dbase doesn't exist
dbopen,obsdbase
sisr=strtrim(string(ismin),2)+' < iseq < '+strtrim(string(ismax),2)
ll=dbfind(sisr)
if(ll(0) gt 0) then begin
if(rf gt 0) then begin
ll1=dbcircle(ral(rf),decl(rf),15.,dis,ll,/silent)
ll=ll1
endif
if(ll(0) lt 0) then goto,nostars
dbext,ll,'starid,iseq,ra,dec,mag,x,chip',sid,iseq,ra,dec,mag,x,chip
dbclose
; make list of unique stars, images, get standard magnitudes for all stars
su=uniq(sid,sort(sid)) & usid=sid(su)
nsu=n_elements(usid)
sq=uniq(iseq,sort(iseq)) & uiseq=iseq(sq)
nui=n_elements(uiseq)
stardbase=starpath+sname(rf)
dbopen,stardbase
ll=dbget('starid',usid)
dbext,ll,filtstar(j),mag0
dbclose
; get fwhm and aperture corrections for the images
dbopen,imdbase
dbext,uiseq,'apercor,fwhm,ra,jd',apercor,fwhm,raio,jdio
dbclose
; count observations for each unique star
nos=lonarr(nsu)
for kk=0,nsu-1 do begin
s=where(sid eq usid(kk),ns)
nos(kk)=ns
endfor
nosmax=max(nos)
; pick out subset of stars that were observed on every possible image but one
; sg=where(nos ge nosmax-1,nsg)
; if(j eq 3) then sg=where(nos ge nosmax-3,nsg) ; ***temporary hack***
; !p.multi = [0,1,2]
; pick out subset of stars that were observed exactly the same number of
; times as many other stars.
ho=histogram(nos,min=0)
hom=max(ho,ixh)
sh=where(ho ge 0.6*hom,nsh) ; find big peaks in the histogram
msh=max(sh) ; big peak with max no of obs
ssh=where(sh ge 0.75*msh,nssh) ; big peaks with similar number of obs
gno=sh(ssh) ; no of obs for each of these peaks
if (nssh gt 4) then gno=gno(nssh-4:nssh-1) ; take at most 4 biggest
ngno=n_elements(gno) ; how many peaks have we got?
sg=[-1]
for iz=0,ngno-1 do begin
ssg=where(nos eq gno(iz),nssg) ; take stars with exactly this no of obs
if(nssg gt 0) then sg=[sg,ssg]
endfor
sg=sg(1:*)
nsg=n_elements(sg)
if(nsg gt 5) then begin
; make medians of (mags of all stars on each image) and (mags-mags0)
; of all stars on each image, likewise of airmass for these stars on each image
; but only for stars with reasonable mag (or mag, mag0 for difference).
medm=fltarr(nui)
meddif=fltarr(nui)
errdif=fltarr(nui)
medx=fltarr(nui)
tmag=fltarr(nsg,nui)+99.9
tmag0=fltarr(nsg,nui)+99.9
tdif=fltarr(nsg,nui)+99.9
tx=fltarr(nsg,nui)+99.9
; also make aperture corrections, fwhm weighted by number of stars for ea image
apca=fltarr(nui)
fwhma=fltarr(nui)
for mm=0,nui-1 do begin
fwa=0.
apa=0.
ntot=0.
for cc=0,3 do begin
si=where(iseq eq uiseq(mm) and chip eq cc,nsi)
apa=apa+nsi*apercor(cc,mm)
fwa=fwa+nsi*fwhm(cc,mm)
ntot=ntot+nsi
endfor
apca(mm)=apa/ntot
fwhma(mm)=fwa/ntot
for nn=0,nsg-1 do begin
st=where(sid eq usid(sg(nn)) and iseq eq uiseq(mm),nst)
if(nst eq 1) then begin
tmag(nn,mm)=mag(st)
tmag0(nn,mm)=mag0(sg(nn))
tdif(nn,mm)=mag(st)-mag0(sg(nn))
tx(nn,mm)=x(st)
endif
endfor
sug=where(tmag(*,mm) lt 40,nsug)
sue=where(tmag(*,mm) lt 40 and tmag0(*,mm) lt 40,nsue)
if(nsue ge 3) then begin
medm(mm)=median(tmag(sug,mm))
meddif(mm)=median(tmag(sue,mm)-tmag0(sue,mm))
quartile,tmag(sue,mm)-tmag0(sue,mm),medofd,qrt,dqrt
errdif(mm)=dqrt
medx(mm)=median(tx(sug,mm))
endif
endfor
; compute linear fits to meddif vs medx, residuals
if(nui ge 2 and max(abs(medx)) gt 0.) then begin
ss=where(medx ne 0.,nss)
if(nss gt 0) then begin
medx=medx(ss)
meddif=meddif(ss)
medm=medm(ss)
nui=nss
rai=raio[ss]
jdi=jdio[ss]
endif
dd=poly_fit(medx,meddif,1)
azero=dd(0)
slope=dd(1)
; look for and reject outliers in the residuals
yyf=poly(medx,dd)
resid=meddif-yyf
nresid=n_elements(resid)
if(nresid ge 5) then begin
quartile,resid,medrr,qrr,dqrr
dmed=abs(resid-medrr)/dqrr
sgdd=where(dmed le 3.*1.349,nsgdd) ; data within 3-sigma
if(nsgdd lt nresid and nsgdd ge 4) then begin
dd=poly_fit(medx(sgdd),meddif(sgdd),1)
azero=dd(0)
slope=dd(1)
print,nresid-nsgdd,'Bad data rejected, filter=',j
endif
endif
endif else begin
azero=meddif(0)
slope=0.
dd=[azero,slope]
if(max(abs(meddif)) eq 0.) then begin
print,'no nonzero meddif -- bailing'
goto,endloop
endif
endelse
yyf=poly(medx,dd)
resid=meddif-yyf
if(n_elements(resid) ge 2) then medrms=stddev(resid) else medrms=99.
minx=min(medx) ; smallest image median airmass
maxx=max(medx) ; largest
nxpts=n_elements(medx) ; number of points to which fit is made
msiga=medrms/sqrt((nxpts-2) > 1)
mxran=(maxx-minx) > .01
msigk=medrms/(mxran*sqrt((nxpts-2) > 1)) ; rough and ready estimates,
; do better later
f1='(i5,2x,1a8,2x,1a4,7(2x,f7.3),i6)'
if(keyword_set(log)) then begin
fs=filtstar(j)
printf,iunl,ib,sdatenam,fs,azero,msiga,slope,msigk,medrms,$
minx,maxx,nsg,format=f1
endif
xx=.8+.5*findgen(4)
yyp=poly(xx,dd)
stra=string(azero+1.215*slope,format='(f6.3)')
strs=string(slope,format='(f6.3)')
strsiga=string(msiga,format='(f7.4)')
strsigk=string(msigk,format='(f7.4)')
strrms=string(medrms,format='(f7.4)')
; get hourangle
get_hourangle, rai, jdi, hrai
ha = INTARR(nui)
sha = WHERE(hrai GT 0, nsha)
IF (nsha GT 0) THEN ha[sha]=1
; make plots of useful stuff
IF (psflag) THEN BEGIN ;plot ps and gif files first, then to screen
set_plot, 'ps'
!p.multi = [0,1,2]
tit='Test: Block '+strtrim(string(blocks(i)),2)+' filter '+filtstar(j)
xtit='Airmass'
ytit='Magnitudes'
xr=[0.9,2.2]
medmed=total(medm)/nui
difmed=total(meddif)/nui
yr1=[medmed-0.2,medmed+0.3]
yr2=[difmed-.3,difmed+.4]
plot,medx,meddif,xran=xr,yran=yr2,tit=tit,xtit=xtit,ytit=ytit, $
/nodata,/xsty,/ysty
for mm=0,nui-1 do begin
IF (ha[mm] EQ 0) THEN symb=4 ELSE symb=5
oplot,[medx(mm)],[meddif(mm)],psym=symb,thick=3
; oplot,tx(*,mm),tdif(*,mm),psym=3
oploterr,[medx(mm)],[meddif(mm)],[errdif(mm)]
endfor
oplot,xx,yyp,line=2
xyouts,1.,yr2(1)-.05, 'Triangle symbol denotes positive hour angle'
xyouts,1.,yr2(1)-.1,'a = '+stra+' +/- '+strsiga
xyouts,1.,yr2(1)-.15,'k = '+strs+' +/- '+strsigk
xyouts,1.,yr2(1)-.2, 'rms = '+strrms
ytit2='Fit Residual (mag)'
xtit1='FWHM (pix)'
xr1=[0.,12.]
yr1=[-.06,.06]
plot,fwhma,resid,psym=2,tit=tit,xtit=xtit1,ytit=ytit2,xran=xr1,$
yran=yr1,/xsty,/ysty
set_plot, 'x'
ENDIF
!p.multi = [0,1,2]
tit='Test: Block '+strtrim(string(blocks(i)),2)+' filter '+filtstar(j)
xtit='Airmass'
ytit='Magnitudes'
xr=[0.9,2.2]
medmed=total(medm)/nui
difmed=total(meddif)/nui
yr1=[medmed-0.2,medmed+0.3]
yr2=[difmed-.3,difmed+.4]
plot,medx,meddif,xran=xr,yran=yr2,tit=tit,xtit=xtit,ytit=ytit,/nodata,$
/xsty,/ysty
for mm=0,nui-1 do begin
IF (ha[mm] EQ 0) THEN symb=4 ELSE symb=5
oplot,[medx(mm)], [meddif(mm)], psym=symb, thick=3
; oplot,tx(*,mm),tdif(*,mm),psym=3
oploterr,[medx(mm)],[meddif(mm)],[errdif(mm)]
endfor
oplot,xx,yyp,line=2
xyouts,1.,yr2(1)-.05, 'Triangle symbol denotes positive hour angle'
xyouts,1.,yr2(1)-.1,'a = '+stra+' +/- '+strsiga
xyouts,1.,yr2(1)-.15,'k = '+strs+' +/- '+strsigk
xyouts,1.,yr2(1)-.2, 'rms = '+strrms
ytit2='Fit Residual (mag)'
xtit1='FWHM (pix)'
xr1=[0.,12.]
yr1=[-.06,.06]
plot,fwhma,resid,psym=2,tit=tit,xtit=xtit1,ytit=ytit2,xran=xr1,$
yran=yr1,/xsty,/ysty
;Create gif
IF (psflag) THEN BEGIN
giffname=gifpath+yr+mon+day+'_rvx_b'+ $
strtrim(string(ib),2)+'_'+strtrim(filtnam(j),2)+'_test.gif'
tvlct, r, g, b, /get
r[0]=255 & g[0]=255 & b[0]=255
r[255]=0 & g[255]=0 & b[255]=0
gif_image = tvrd()
write_gif, giffname, gif_image, r, g, b
ENDIF
endif ;from: if(nsg gt 5) then begin
nostars:
wdelete, 5
endif
endloop:
endfor ; end of filters loop
if (psflag) then begin
set_plot, 'ps'
device, /close
set_plot, 'x'
endif
endfor ; end of blocks loop
!p.multi=[0,1,1]
if (keyword_set(log)) then begin
close, iunl
free_lun,iunl
endif
;stop
end