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

Posted in Research and tagged as , , , ,

Comments are closed.