;+ ; NAME: ; onchip ; PURPOSE: (one line) ; Extract and optionally plot differential on-chip photometry. ; DESCRIPTION: ; This program will read an altlog photometry file as generated by ; basphote and extract one object from the file, compute differential ; photometry for the program object (Serial #0000) against the field ; stars. ; CATEGORY: ; Photometry ; CALLING SEQUENCE: ; onchip,logname,objname,filter,object,primary,others,jd,time,mag,err ; INPUTS: ; logname - Name of the altlog format photometry file to read. ; objname - String with the exact name of the object to extract. ; filter - String containing the filter code to extract. ; object - Serial number of the object to extract. ; primary - Serial number of the on-chip reference star to use as ; the primary reference. ; others - Serial number(s) of the other on-chip stars to use as a ; check of the primary reference. Scalar or vector. ; If set to -1, no secondary references are used. ; OPTIONAL INPUT PARAMETERS: ; KEYWORD PARAMETERS: ; BARS - Flag, if true, extend major ticks across plot (default=false) ; MAGPRI - Magnitude of primary reference star (default=0). ; NOPLOT - When set will inhibit generating a plot to the current ; plot device. (same as /NOPOS and /NOPHOT together). ; NOPOS - No plot of the positions. ; NOPHOT - No plot of the photometry. ; COLORS - Colors to use for the plot. ; colors(0) = axis, program object, (primary standard) ; default = !d.n_colors-1 ; colors(1) = primary standard, default = 1 ; colors(2) = secondary standard 1, default = 2 ; colors(3) = secondary standard 2, default = 2 ; and so on. COLORS should have as many elements as ; needed for the object, primary standard, and all ; secondary standards. ; NOERRBAR - Flag, if set suppresses error bars on star check mag plot ; SHOWALL - Flag, if set, prints out all the relevant mangitudes to ; allow seeing which mags are being used. ; TIMEUNIT - Select time units to plot on x-axis, ; 0 - UT hours (default) ; 1 - JD (decimal days) ; UPPERCASE - Flag, if true, forces object names and filter names to upper case. ; OUTPUTS: ; jd - Julian date of observation. ; time - UT of observation, decimal hours. (or JD, depending on TIMEUNIT) ; mag - reduced magnitude of program object. ; err - Uncertainty of the magnitude. ; COMMON BLOCKS: ; SIDE EFFECTS: ; RESTRICTIONS: ; PROCEDURE: ; MODIFICATION HISTORY: ; 2/4/93, Written by Marc W. Buie, Lowell Observatory. ; 6/29/93, MWB, added filter argument. ; 96 Mar, MWB, added MAGPRI keyword ; 96/06/15, MWB, improved the plotting output for the photometry plot. ; 96/08/24, MWB, added TIMEUNIT and BARS keywords ; 96/10/31, MWB, added "bad" flag support from rdphalt ; 97/09/20, MWB, added UPPERCASE flag. ;- pro onchip,logname,objname,filter,object,primary,in_others,jd,time,mag,err, $ NOPLOT=noplot, COLORS=colors, NOPOS=nopos, NOERRBAR=noerrbar, SHOWALL=showall, $ SKIPLINE=skipline,NUMLINE=numline,MAGPRI=magpri,TIMEUNIT=timeunit, $ BARS=bars,UPPERCASE=uppercase if n_params() lt 5 then begin print,'onchip,logname,objname,filter,object,primary,others,[jd,time,mag,err]' return endif if badpar(logname,7,0,caller='ONCHIP: (logname) ') then return if badpar(objname,7,0,caller='ONCHIP: (objname) ') then return if badpar(filter,7,0,caller='ONCHIP: (filter) ') then return if badpar(object,[2,3],0,caller='ONCHIP: (object) ') then return if badpar(primary,[2,3],0,caller='ONCHIP: (primary) ') then return if badpar(in_others,[2,3],[0,1],caller='ONCHIP: (others) ') then return if badpar(skipline,[0,2,3],0,caller='ONCHIP: (skipline) ',default=0) then return if badpar(numline,[0,2,3],0,caller='ONCHIP: (numline) ',default=-1) then return if badpar(magpri,[0,2,3,4,5],0,caller='ONCHIP: (magpri) ',default=0.0) then return if badpar(timeunit,[0,2,3],0,caller='ONCHIP: (timeunit) ',default=0) then return if badpar(bars,[0,2,3],0,caller='ONCHIP: (bars) ',default=0) then return others = in_others sz=size(others) if sz[0] eq 0 then others=[others] if bars then ticklen=1.0 else ticklen=!p.ticklen if !d.n_colors le 256 then begin ncolors=!d.n_colors eightbit = 1 endif else begin ncolors=256 eightbit = 0 endelse if not keyword_set(colors) then begin if !d.name eq 'PS' then begin loadct,0,/silent colors=[0,!d.n_colors*0.9, $ replicate(!d.n_colors*0.75,n_elements(others))] endif else begin if eightbit then begin colors=[!d.n_colors-1,1,replicate(2,n_elements(others))] endif else begin colors=['ffffff'xl,'0000ff'xl,replicate('00ff00'xl,n_elements(others))] endelse endelse endif title=objname+' and stars from '+logname+' for filter '+filter print,title ; Read the data in the log file. rdphalt,logname,filename,obj,fil,rawjd,exptime,gain,rad,sky1,sky2, $ serial,xpos,ypos,fwhm,maxcnt,sky,skyerr,rawmag,rawerr,bad,skipline=skipline,numline=numline if keyword_set(uppercase) then begin obj = strupcase(obj) fil = strupcase(fil) endif ; Locate all measurements relevant to the requested object and filter. idx = where(obj eq objname and fil eq filter and bad eq 0, count) if (count le 0) then begin print,'ONCHIP: Object ',objname,' and filter ',filter,' not found.' allobj = obj[uniq(obj,sort(obj))] print,'Objects in file are:' print,allobj return endif tmp = serial[idx] allobjid = tmp[uniq(tmp,sort(tmp))] others = others[uniq(others,sort(others))] print,' ' print,'Object list',allobjid,format='(a,1x,25(1x,i2))' print,'Object serial number ',object print,'Primary reference star ',primary if others[0] eq -1 then $ print,'No secondary reference stars used' $ else $ print,'Secondary stars ',others,format='(a,1x,25(1x,i2))' print,' ' ;Construct a single vector that records original serial number for later id if others[0] ne -1 then $ serid=[object,primary,others] $ else $ serid=[object,primary] nred = n_elements(serid) ; Validate and make sure requested objects are really in data for i=0,nred-1 do begin z=where(serid[i] eq allobjid) if z[0] eq -1 then begin if i eq 0 then $ print,'ERROR: Object serial number ',serid[i],' is not in the photometry file.' $ else if i eq 1 then $ print,'ERROR: Primary star serial number ',serid[i],' is not in the photometry file.' $ else $ print,'ERROR: Secondary star serial number ',serid[i],' is not in the photometry file.' return endif endfor ; Build the data array to be reduced. ; First, collect information on the program object. This defines ; the total number of observations and the object information. z = where(serial[idx] eq object, nobs) if nobs eq 0 then begin print,'ERROR: Object ',objname,' not found in the photometry file.' return endif jd = rawjd[idx[z]] ; start building the data array with the first row. l_fname = filename[idx[z]] l_mag = rawmag[idx[z]] l_err = rawerr[idx[z]] l_fwhm = fwhm[idx[z]] l_xpos = xpos[idx[z]] l_ypos = ypos[idx[z]] if timeunit eq 0 then begin jd0=long(jd[0]-0.5)+0.5d0 time=float(jd-jd0)*24.0 xlabel = 'UT time (hours)' endif else begin jd0=long(jd[0]-0.5)+0.5d0 jd0=long(jd0/10000.0d0)*10000.0d0 xlabel = string(long(jd0),format='("Julian Date (",i7,"+)")') time=float(jd-jd0) endelse ; Second, grab the data from the primary reference star. At the ; moment, it is REQUIRED that the primary reference star be ; measured on all frames that the object is measured on. ; Check for proper length z = where(serial[idx] eq primary, nobsp) if nobsp ne nobs then begin print,'ERROR: the object and primary star must have the same number of' print,' observations. object=',nobs,' primary star=',nobsp for i=0,min([nobs,nobsp])-1 do begin if filename[idx[z[i]]] ne l_fname[i] then begin print,'First name mismatch is ',i,' ',l_fname[i],' ',filename[idx[z[i]]] i=nobs+1 endif endfor print,'Last star frame ',filename[idx[z[nobsp-1]]] print,'First frame and serial number ',filename[0],serial[0] print,'First frame and serial number ',filename[n_elements(serial)-1],serial[n_elements(serial)-1] return endif ; Check to see if file names match in the two lists. for i=0,nobs-1 do begin if filename[idx[z[i]]] ne l_fname[i] then begin print,'ERROR: object/primary star filename mismatch' print,' point #',i,' object file ',l_fname[i],' star file ',filename[idx[z[i]]] return endif endfor ; add a new row for the primary reference star. l_mag = [[l_mag], [rawmag[idx[z]]]] l_err = [[l_err], [rawerr[idx[z]]]] l_fwhm = [[l_fwhm],[fwhm[idx[z]] ]] l_xpos = [[l_xpos],[xpos[idx[z]] ]] l_ypos = [[l_ypos],[ypos[idx[z]] ]] ; Now add in the measurements for additional stars. for i=2,nred-1 do begin ; locate the observations for this star z = where(serial[idx] eq serid[i], nobsp) ; add new row to arrays. l_mag = [[l_mag], [replicate(199.999,nobs)]] l_err = [[l_err], [replicate( 0.0 ,nobs)]] l_fwhm = [[l_fwhm],[replicate( 0.0 ,nobs)]] l_xpos = [[l_xpos],[replicate( 0.0 ,nobs)]] l_ypos = [[l_ypos],[replicate( 0.0 ,nobs)]] for j=0,nobsp-1 do begin zz = where(filename[idx[z[j]]] eq l_fname) if zz[0] ne -1 then begin l_mag[zz[0],i] = rawmag[idx[z[j]]] l_err[zz[0],i] = rawerr[idx[z[j]]] l_fwhm[zz[0],i] = fwhm[idx[z[j]]] l_xpos[zz[0],i] = xpos[idx[z[j]]] l_ypos[zz[0],i] = ypos[idx[z[j]]] endif endfor endfor ;Time to compute the composite star magnitude. ;Declare array variables for later use. master_star = fltarr(nobs) master_star_err = fltarr(nobs) ; Array to contain mean difference dmeans=fltarr(nred) derrs=fltarr(nred) relfwhm=fltarr(nred) ; Array to contain the objects scaled to primary diffs=fltarr(nobs,nred) scaled = diffs ; Position vectors of the primary. xref=l_xpos[*,1] yref=l_ypos[*,1] ;Compute the mean magnitude difference between each object and ; the primary reference star. nmeas=replicate(nobs,nred) ;Loop over all objects. for i=0,nred-1 do begin ;Compute the difference between object and the primary reference star. diffs[*,i] = l_mag[*,i]-l_mag[*,1] fwhmrat = l_fwhm[*,i]/l_fwhm[*,1] ;Compute a robust mean of the differenced magnitude, skip primary ; since it will be identically zero. if i ne 1 then begin z = where(l_mag[*,i] lt 80,count) if keyword_set(showall) then begin print,i,count,format='(i3,1x,i5)' print,l_mag[*,i],format='(10(1x,f6.3))' endif robomean,diffs[z,i],3.0,0.5,avg,avgdev,stddev,var,skew,kurt,nfinal dmeans[i] = avg if nfinal ne 1 then $ derrs[i] = stddev / sqrt(float(nfinal-1)) $ else $ derrs[i] = stddev nmeas[i]=nfinal robomean,fwhmrat[z],3.0,0.5,avg relfwhm[i] = avg endif else begin dmeans[i] = 0.0 relfwhm[i] = 1.0 endelse ;Compute the positions relative to primary. l_xpos[*,i] = l_xpos[*,i] - xref l_ypos[*,i] = l_ypos[*,i] - yref endfor ;Report results to screen, just for fun. print,'Serial # Mag Diff Sigma Npts used relFWHM' for i=0,nred-1 do begin if serid[i] ne primary then $ print,serid[i],' - ',primary,dmeans[i],' +- ',derrs[i], $ nmeas[i],relfwhm[i], $ format='(1x,i2,a,i2,3x,f7.4,a,f6.4,5x,i5,5x,f5.2)' endfor ;Recompute star magnitudes, scaled to primary mean (note primary is a NOP). for i=1,nred-1 do begin scaled[*,i] = l_mag[*,i] - dmeans[i] diffs[*,i] = diffs[*,i] - dmeans[i] endfor ;Compute the mean star magnitude plus its std dev at each point. This will ; use all the scaled stars. The uncertainty is derived from the scatter. if nred gt 2 then begin for i=0,nobs-1 do begin toavg = scaled[i,1:nred-1] toerr = l_err[i,1:nred-1] orig = l_mag[i,1:nred-1] z = where(orig lt 80,count) if count gt 1 then begin robomean,toavg[z],3.0,0.5, $ avg,avgdev,stddev,var,skew,kurt,nfinal master_star[i] = avg if nfinal ne 1 then $ master_star_err[i] = stddev / sqrt(float(nfinal-1)) $ else $ master_star_err[i] = stddev endif else if count eq 1 then begin master_star[i] = toavg[z] master_star_err[i] = toerr[z] endif else begin print,'ONCHIP: Warning, observation',i,' has all points undefined.' master_star[i] = 99.999 master_star_err[i] = 0.000 endelse endfor endif else begin master_star = l_mag[*,1] master_star_err = l_err[*,1] endelse ; These are the final object magnitudes and errors. mag=l_mag[*,0]-master_star+magpri err=sqrt(l_err[*,0]^2+master_star_err^2) mag=reform(mag,n_elements(mag)) err=reform(err,n_elements(err)) for i=1,nred-1 do begin diffs[*,i] = l_mag[*,i] - dmeans[i] - master_star endfor plx1=0.12 plx2=0.95 ply1=0.10 ply2=0.95 pad=0.005 cs=1.0 blanks=[' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' '] if not keyword_set(noplot) then begin p_multi=!p.multi if not keyword_set(nophot) then begin ; if others(0) eq -1 then !p.multi=[0,1,3] else !p.multi=[0,1,4] if others[0] eq -1 then yplsz=(ply2-ply1)/3.0 else yplsz=(ply2-ply1)/4.0 plotnum=0 ;Everything's ready for the plot. setwin,0 xr=[min(time),max(time)] good = where(l_mag[*,1] lt 80.0) yr=[max(l_mag[good,1]),min(l_mag[good,1])] plpos=[plx1,ply2-(plotnum+1)*yplsz+pad,plx2,ply2-plotnum*yplsz] plot,[0,1],[0,1],xr=xr,yr=yr,/nodata,charsize=cs,ytit='Raw Mag', $ title=title,position=plpos,xstyle=7,yticklen=ticklen axis,xaxis=1,xstyle=3,xtickn=blanks,xticklen=ticklen,xminor=4 axis,xaxis=0,xstyle=3,xtickn=blanks,xticklen=ticklen,xminor=4 for i=nred-1,0,-1 do begin oplot,time,l_mag[*,i]-dmeans[i], $ psym=8,color=colors[i],symsiz=0.5 endfor if nred gt 2 then begin z=where(l_mag[*,1:nred-1] lt 80) good = diffs[*,1:nred-1]+l_err[*,1:nred-1] ymax=max(good[z]) good = diffs[*,1:nred-1]-l_err[*,1:nred-1] ymin=min(good[z]) yr=[ymax, ymin] plotnum=plotnum+1 plpos=[plx1,ply2-(plotnum+1)*yplsz+pad,plx2,ply2-plotnum*yplsz] plot,[0,1],[0,1],xr=xr,yr=yr,/nodata,charsize=cs, $ ytit='Star Check Mags',position=plpos,/noerase,xstyle=7,yticklen=ticklen axis,xaxis=1,xstyle=3,xtickn=blanks,xticklen=ticklen,xminor=4 axis,xaxis=0,xstyle=3,xtickn=blanks,xticklen=ticklen,xminor=4 for i=nred-1,1,-1 do begin good = diffs[*,i] gooderr = l_err[*,i] sel = l_mag[*,i] z=where(sel lt 80) if not keyword_set(noerrbar) then $ oploterror,time[z],good[z],gooderr[z],psym=8,color=colors[i],symsiz=0.5,errcolor=colors[i] $ else $ oplot,time[z],good[z],psym=8,color=colors[i],symsiz=0.5 endfor endif yr=[max(mag+err),min(mag-err)] plotnum=plotnum+1 plpos=[plx1,ply2-(plotnum+1)*yplsz+pad,plx2,ply2-plotnum*yplsz] if magpri eq 0. then ytitle=objname+' - comp' else ytitle=objname plot,[0,1],[0,1],xr=xr,yr=yr,/nodata,charsize=cs,ytit=ytitle, $ position=plpos,/noerase,xstyle=7,yticklen=ticklen axis,xaxis=1,xstyle=3,xtickn=blanks,xticklen=ticklen,xminor=4 axis,xaxis=0,xstyle=3,xtickn=blanks,xticklen=ticklen,xminor=4 oploterror,time,mag,err,psym=8,color=colors[0],symsiz=0.5,errcolor=colors[0] z = where(l_fwhm gt 0.0) good=l_fwhm[z] yr=[min(good),max(good)] plotnum=plotnum+1 plpos=[plx1,ply2-(plotnum+1)*yplsz+pad,plx2,ply2-plotnum*yplsz] plot,[0,1],[0,1],xr=xr,yr=yr,/nodata,charsize=cs, $ ytit='FWHM',position=plpos,/noerase,xstyle=7,yticklen=ticklen axis,xaxis=1,xstyle=3,xtickn=blanks,xticklen=ticklen,xminor=4 axis,xaxis=0,xstyle=3,xticklen=ticklen,xminor=4,xtit=xlabel for i=nred-1,0,-1 do begin oplot,time,l_fwhm[*,i],psym=8,color=colors[i],symsiz=0.5 endfor if !d.name eq 'TEK' and not keyword_set(NOPOS) then begin print,'Hit a key to continue' junk=get_kbrd(1) endif endif ; PHOT block if not keyword_set(NOPOS) then begin setwin,1 if others[0] eq -1 then !p.multi=[0,2,1] else !p.multi=[0,2,2] x_coeff = goodpoly(time,l_xpos[*,0],1,5.) x_rate=x_coeff[1]/60.0 title1=objname+' '+string(x_rate,format='(f7.4)')+' pix/min' plot,time,l_xpos[*,0],psym=3,tit=title1, $ xtit=xlabel,ytit='X position' y_coeff = goodpoly(time,l_ypos[*,0],1,5.) y_rate=y_coeff[1]/60.0 title1=objname+' '+string(y_rate,format='(f7.4)')+' pix/min' plot,time,l_ypos[*,0],psym=3,tit=title1, $ xtit=xlabel,ytit='Y position' if nred gt 2 then begin yr=[min(l_xpos[*,2:nred-1]), $ max(l_xpos[*,2:nred-1])] plot,[0,1],[0,1],xr=xr,yr=yr,/nodata,xtit=xlabel,ytit='X position' for i=2,nred-1 do begin zz = where( l_mag[*,i] lt 80.0 ) oplot,time,l_xpos[zz,i],color=colors[i],psym=3 endfor yr=[min(l_ypos[*,2:nred-1]), $ max(l_ypos[*,2:nred-1])] plot,[0,1],[0,1],xr=xr,yr=yr,/nodata,xtit=xlabel,ytit='Y position' for i=2,nred-1 do begin zz = where( l_mag[*,i] lt 80.0 ) oplot,time,l_ypos[zz,i],color=colors[i],psym=3 endfor endif if !d.name eq 'TEK' then begin print,'Hit a key to continue' junk=get_kbrd(1) endif setwin,2,xsize=400,ysize=400 !p.multi=0 zz = where( l_mag lt 80.0 ) xr=minmax(l_xpos[zz]) yr=minmax(l_ypos[zz]) plot,[0,1],[0,1],xr=xr,yr=yr,/nodata, $ xtit='X position',ytit='Y position',tit=title for i=0,nred-1 do begin zz = where( l_mag[*,i] lt 80.0 ) oplot,l_xpos[zz,i],l_ypos[zz,i], $ color=colors[i],psym=8,symsiz=0.3 endfor endif ; POS block setwin,0 !p.multi=p_multi endif ; PLOT block if timeunit eq 1 then time=jd end