Scripting the Prescription

I’ve found it very difficult to keep the photometry steps in my head, and having very limited knowledge of IRAF and its scripting abilities, I’ve generated a Python script. This script (below) does two things: 1. spits out the needed IRAF commands allowing for a quick copy/paste into the IRAF terminal, 2. easy documentation of photometry files for use in time series analysis.

 

——Python to IRAF script—–

author: Jesse Rogerson

date: 9 November 2012

 

import sys
import numpy as np
import pyfits
import math

print ‘This Script is to prompt the user on the steps to calculate’
print ‘the g,r, and i magnitudes from the CFHT data.’
print ‘To use properly, simply follow the instructions and prompts’
print ‘throughout.’
obj=raw_input(‘Name of object:’)
print ”
print ”
print ‘STEP1’
print ‘Unpack and display the 3 fits files to be used (g,r,i)’
print ‘STEP2’
print ‘Use ds9 to determine the coordinates of the quasar to be measured’
print ‘and 5 comparison stars for aperture correction’
print ‘place the comparison coords in ‘+obj+’.coo.1′
print ‘place quasar coords in ‘+obj+’.coo.2′
var=raw_input(‘continue? (y/n)’)
print ”
print ”
print ‘STEP3’
print ‘INPUT FILES:’
gfile=raw_input(‘Name of g filter fits file:’)
rfile=raw_input(‘Name of r filter fits file:’)
ifile=raw_input(‘Name of i filter fits file:’)
ccdNum=raw_input(‘CCDnumber:’)
compCoo=obj+’.coo.1′
objCoo=obj+’.coo.2′
print ‘Name of comparison star coordinatefile:’+compCoo
print ‘Name of quasar coordinate file:’+objCoo
print ”
print ”
print ‘STEP4’
print ‘In IRAF, load noao.digiphot.apphot, set pars files’
print ‘datapars,photpars,fitskypars’
var=raw_input(‘continue? (y/n)’)
print ”
print ”
print ‘STEP5’
#photometer the fits files given the coordinate files
#output files will be *.fits.mag.#
print ‘Photometer commands’
print ‘Run the following commands in IRAF:’
print ‘phot’+’ ‘+gfile+'[ccd’+ccdNum+’] coo=’+compCoo+’ interac-‘
print ‘phot’+’ ‘+gfile+'[ccd’+ccdNum+’] coo=’+objCoo+’ interac-‘
print ‘phot’+’ ‘+rfile+'[ccd’+ccdNum+’] coo=’+compCoo+’ interac-‘
print ‘phot’+’ ‘+rfile+'[ccd’+ccdNum+’] coo=’+objCoo+’ interac-‘
print ‘phot’+’ ‘+ifile+'[ccd’+ccdNum+’] coo=’+compCoo+’ interac-‘
print ‘phot’+’ ‘+ifile+'[ccd’+ccdNum+’] coo=’+objCoo+’ interac-‘
var=raw_input(‘continue? (y/n)’)
print ”
print ”
print ‘STEP6’
#create aperture correction files *.apcor.1
print ‘load noao.digiphot.photcal’
print ‘Run the following commands in IRAF’
print ‘watch for convergence….’
print ”
print ”
print ‘mkapf ‘+gfile+’.mag.1 9 ‘+gfile+’.apcor.1 small=1 large=9′
print ‘mkapf ‘+rfile+’.mag.1 9 ‘+rfile+’.apcor.1 small=1 large=9′
print ‘mkapf ‘+ifile+’.mag.1 9 ‘+ifile+’.apcor.1 small=1 large=9′
var=raw_input(‘continue? (y/n)’)
print ”
print ”
print ‘STEP7’
print ‘Extract quasar magnitude and apply aperture correction’
print ‘Run the following commands in IRAF’
print ‘load noao.digiphot.ptools’
#g:
print ‘txdump ‘+gfile+’.mag.2 mag yes | scan (x)’
print ‘txdump ‘+gfile+’.mag.2 merr yes | scan (y)’
print ‘fields ‘+gfile+’.apcor.1 3 | sort num+ rev+ | scan (z)’
print ‘y=sqrt(y**2.+z**2.)’
print ‘fields ‘+gfile+’.apcor.1 2 | sort num+ rev+ | scan (z)’
print ‘x=x+z’
print ‘printf(“%7.3f %6.3f %8.3f”,x,y,z, >> “‘+gfile+’.mea”)’
var=raw_input(‘continue? (y/n)’)
#r:
print ‘txdump ‘+rfile+’.mag.2 mag yes | scan (x)’
print ‘txdump ‘+rfile+’.mag.2 merr yes | scan (y)’
print ‘fields ‘+rfile+’.apcor.1 3 | sort num+ rev+ | scan (z)’
print ‘y=sqrt(y**2.+z**2.)’
print ‘fields ‘+rfile+’.apcor.1 2 | sort num+ rev+ | scan (z)’
print ‘x=x+z’
print ‘printf(“%7.3f %6.3f %8.3f”,x,y,z, >> “‘+rfile+’.mea”)’
var=raw_input(‘continue? (y/n)’)
#i:
print ‘txdump ‘+ifile+’.mag.2 mag yes | scan (x)’
print ‘txdump ‘+ifile+’.mag.2 merr yes | scan (y)’
print ‘fields ‘+ifile+’.apcor.1 3 | sort num+ rev+ | scan (z)’
print ‘y=sqrt(y**2.+z**2.)’
print ‘fields ‘+ifile+’.apcor.1 2 | sort num+ rev+ | scan (z)’
print ‘x=x+z’
print ‘printf(“%7.3f %6.3f %8.3f”,x,y,z, >> “‘+ifile+’.mea”)’
var=raw_input(‘continue? (y/n)’)
print ”
print ”

print ‘Extracting Header Info for calculations’

files=(gfile,rfile,ifile)
#files=(‘1572573p.fits’,’1572590p.fits’,’1572602p.fits’)
outfile=open(path+obj+’.tbl’,’w’)
#outfile=open(path+’002845.77+010648.3.tbl’,’w’)

outfile.write(‘#inst_mag merr AIRMASS PHOTO_C’
+’PHOTO_CS PHOTO_X PHOT_K MJD \n’)

for f in files:
hdulist=pyfits.open(f)
airmass=hdulist[0].header[‘AIRMASS’]
pc=hdulist[0].header[‘PHOT_C’]
pcs=hdulist[0].header[‘PHOT_CS’]
px=hdulist[0].header[‘PHOT_X’]
pk=hdulist[0].header[‘PHOT_K’]
mjd=hdulist[0].header[‘MJDATE’]
hdulist.close()
mag=np.genfromtxt(path+f+’.mea’,dtype=str)
outfile.write(mag[0]+’ ‘+mag[1]+’ ‘+str(airmass)+’ ‘+str(pc)+’ ‘+str(pcs)
+’ ‘+str(px)+’ ‘+str(pk)+’ ‘+str(mjd)+’\n’)

outfile.close()

dT=np.genfromtxt(path+obj+’.tbl’,skip_header=1,dtype=float)
g=dT[0,:]
r=dT[1,:]
i=dT[2,:]

#R = r + r_c + r_K*(X_r-1)
#Rerr= sqrt(r_err**2.+r_cs**2.)
R=r[0]+r[3]+(r[6]*(r[2]-1.))
Rerr=math.sqrt(r[1]**2.+r[4]**2.)

#G = [ g + g_c + g_K*(X_g-1) – g_X*R ] / [1-g_X]
#Gerr = sqrt(g_err**2.+g_cs**2.)
G=(g[0]+g[3]+(g[6]*(g[2]-1.)) – (g[5]*R))/(1.-g[5])
Gerr=math.sqrt(g[1]**2.+g[4]**2.)

#I = [ i + i_c + i_K*(X_i-1) + i_X*R ] / [1+i_X)
#Ierr = sqrt(i_err**2.+i_cs**2.)
I=(i[0]+i[3]+(i[6]*(i[2]-1.)) + (i[5]*R))/(1.+i[5])
Ierr=math.sqrt(i[1]**2.+i[4]**2.)

print G, Gerr
print R, Rerr
print I, Ierr

outfile=open(path+’BPV/CFHTprop/PHOT/’+obj+’.CFHTphot.dat’,’w’)
outfile.write(g[7]+’ ‘+str(G)+’ ‘+str(Gerr)+’ ‘+str(R)
+’ ‘+str(Rerr)+’ ‘+str(I)+’ ‘+str(Ierr))
outfile.close()
#G-R; err = sqrt(Gerr**2.+Rerr**2.) #G-R & err
#R-I; err = sqrt(Rerr**2.+Ierr**2.) #R-I & err

Photometry, first steps

As the CFHT pointings (see my queue update) are not designed to hit one, but many targets, I was forced to use the pointing information to determine which *.fits file a given Stripe82 object was located. Further to the problem, the MegaPrime camera as 36 CCDs in its almost 1 square degree field of view. In order to determine on which of these 36 CCD a given Stripe 82 quasar would land, I used the following formulae:

Parameters:

CRVAL1 – RA (in deg) of pointing centre
CRVAL2 – Dec (in deg) of pointing centre

RA/Dec – quasar RA/Dec

CRPIX1, CRPIX2, CD1_1, CD2_2 – constants given by the CCD layout.

Calculation:
x=CRPIX1 + (RA – CRVAL1)/CD1_1
y=CRPIX2 + (Dec- CRVAL2)/CD2_2

If I calculate the (x,y) (in pix) using the constants from all 36 CCDs, only one CCD will satisfy: 0 < x < 2048,0 < y < 4612. That is the CCD on which the quasar is located. Labeled *.fits[ccd##] x y

With chip number and (x,y) coordinates in hand, I can begin photometry. Using Phot, in $/noao/digiphot/apphot, requires the use of 4 parameter files: datapars, centerpars, photpars, fitskypars. Notable parameters: datapars.ccdread, gain, airmass, filter, exptime, photpars.apertur, zmag, fitskypars.annulus,dannulus

photpars.aperture allows you to choose a number of apertures from which to sample the object, showing the flattening of photons gained with ever increasing apertures.

1572573p.fits[ccd08]
1690.088 3079.221 (estimated x=1606, 2963)

aper     sum(counts)     area(pixels)     flux(counts)     mag     merr
2.00     7505.262         12.87262         4447.409         20.171 0.017
3.00     14969.63         28.55213         8187.152         19.508 0.013
4.00     23573.49         50.56833         11561.13         19.133 0.012
5.00     33010.21         78.8496           14279.73         18.904 0.011
6.00     43157.03         113.218           16262.43         18.763 0.011
7.00     54404.07         154.346           17739.63         18.669 0.012
9.00     79897.76         254.7832         19374.78         18.573 0.014
12.00  128163.3          452.8017         20601.61         18.506 0.017
15.00   188872.4         706.9811         20931.22         18.489 0.021

It’s clear the magnitude measured here is levelling off, the best mag_err is at an aperture of 6 pixels. The question now is: how do I convert this magnitude to an an actual magnitude. Looking at the CFHT image header I’m given:

PHOT_C = 26.4850 / Elixir zero point – measured for camera run
PHOT_CS = 0.0098 / Elixir zero point – scatter
PHOT_NS = 56 / Elixir zero point – N stars
PHOT_NM = 18 / Elixir zero point – N images
PHOT_C0 = 26.4600 / Elixir zero point – nominal
PHOT_X = 0.1480 / Elixir zero point – color term
PHOT_K = -0.1500 / Elixir zero point – airmass term
PHOT_C1 = ‘g_SDSS ‘ / Elixir zero point – color 1
PHOT_C2 = ‘r_SDSS ‘ / Elixir zero point – color 2
COMMENT Photometric Analysis is incomplete for this image.
COMMENT MAG_SAT and MAG_LIM cannot be determined.
COMMENT Formula for Photometry, based on keywords given in this header:
COMMENT m = -2.5*log(DN) + 2.5*log(EXPTIME)
COMMENT M = m + PHOT_C + PHOT_K*(AIRMASS – 1) + PHOT_X*(PHOT_C1 – PHOT_C2)

It seems te formulae at the bottom are the key.

 

 

Notes on Spectra Reduction

I’ve examined the flats to determine the image portion of the chip, and the bias portion of the chip.  These correspond to:

TRIMSEC [10:1198, 5:800]

BIASSEC [1205:1230,5:800]

Next step is to combine all bias frames with zerocombine.  Question: what does ‘to clean up’ mean in the observing notes? Perhaps this was a bias image taken to be sure the CCD was cleared?  I will avoid their use.  Bias images correspond to images 0076 through 0127.  Hmm, the standard zerocombine simply uses all with type ‘zero,’ perhaps I should hide all ones that I don’t think are necessary (i.e., the ‘clean up’ zeros).  Removed all ‘clean up’ biases, for fear they would contaminate the average.  can be put back in again if needed.  This included image 0154, which was not labeled as ‘clean-up’ but have been, for it was done directly after dewer refill.

63 zero images were used

ccdproc was run

combine flat-fields, moved aside those with < 25 000 counts.  not useful.

Flats are broken into ones with 30 sec exposure, and those with 40 sec exposure.

had to change ‘rdnoise’ in ‘flatcombine’ package to 3.5, this was not included in the header, so needed to do it manually

flatcombine done: 65 flats used

Normalized the flat field using ‘response’ package.  The flat field seems really bad, needs a order 6 polynomial to fit to it.

ran ccdproc with normalized flat field, however there is still a gradient left in the images.  Need to map the light correction function. hmmm not sure if i did this properly..

 

my initial pass through the above reductions steps seem to be done…okay.  I would like to go through them again, with a little more understanding.

 

Reduction Steps: Spectra

In order to both better familiarize myself with reducing astronomical data (both the basics from a CCD and the additional steps needed for reducing spectra), and create a lasting document from which to quick-reference, I’ve decided to write up a summary of the steps needed.

This assumes that I have calibration exposures:
biases (zero-second exposures)
flat-field exposures
twighlight exposures – (labeled as ‘skyflat’ in my data, I assume) *these are used to create the illumination correction or ‘slit illumination function’
Reductions Steps:

1. Examine flatfield exposure and determine the area of the chip that contains good data. Also determine the area of the chip that contains a flat overscan.
Overscan:  typically the 32 columns at the RIGHT edge of the frames.  We will average the data over all the columns in the overscan region, and fit these values as a function of line number.  This fit will then be subtracted from each column in the image.  The overscan is then chopped off.  Typically chopping also includes the first/last few rows/columns…just in case

2. combine the biases into one image (using zerocombine)

3. before moving forward, remove the overscan region, process the average bias, and trim the first/last few lines/columns
(essentially, subtracting out the bias from the chip)

4. Creating a proper flat field exposure:
a. combine the flat-field exposures
b. fit a function in the dispersion direction to the combined flat-field (normalizing the combined flat-field)
c. process all target frames and skyflats using the combined flat-field (essentially, dividing out the flat field from all science images).
d. if there is some sort of gradient left over, use the the skyflats to create an illumination correction function, then turn on the illumination correction switch specifying these illumination function in ccdproc

The above is the basics for removing the biases and non-uniformity of the CCD.  At this point, it should be considered that ZERO counts in a pixel is equal to zero light.  Now there are additional steps for extracting and calibrating the spectrum you’ve obtained.

1. Find the spectrum
May be done manually or automatically (the latter works only when the desired spectrum has the strongest peak in the image)

2. Define the extraction window and the background window

3. Trace the centre of spatial profile as a function of the dispersion axis

4. Sum the spectrum within the extraction window, subtracting the sky

5. wavelength calibration – using the comparison lamp spectra

(optional) 6. flux calibration /normalization