def p60photsub(ofile, inlis, refimage, refstars, ot, filter):
   
    refstars.set_mag(filter)
    coofile=open('ot.coo','w')
    ot.wcs2pix(refimage)
    coofile.write('%10.3f%10.3f\n' % (ot[0].xval, ot[0].yval))
    coofile.close()

    outfile=open(ofile, 'a+')
  
    p60photutils.p60zeropt(inlis, refstars, filter)

    subimage=iraffiles('@%s' % inlis)
    for image in subimage:

        root,ext=image.split('.')
        iraf.phot('%s.sub' % root, coo='ot.coo', output='%s.grb' % root,
                  interac=no)

        stars=Starlist('%s.mag' % root)
        refstars.wcs2pix(image)
        zp,zpu=stars.zeropt(refstars,method='mean',rejout=0)

        otstars=Starlist('%s.grb' % root)    
        
        [mjd,utshut]=get_head(image,['OBSMJD','UTSHUT'])

        if (len(stars)>0):
            outfile.write('%s\t%s\t%.3f\t%.3f\t%.3f\t%.3f\n' % (mjd, utshut, otstars[0].mag+zp, sqrt(pow(zpu,2)+pow(otstars[0].magu,2)), otstars[0].magu, zpu))

    outfile.close()

 def daophot(self, framenum, coords, outfile, apertures, verbose = "no"):
     """
     Aperture photometry of stars in the coords file using IRAF PHOT routine.
     
     Parameters
     ----------
     framenum : int
         Frame number in the image cube to perform aperture photometry on.
         
     coords : string
         Text file with coordinate of the stars
         
     outfile : string
         Text file to which photometry results are written
     
     Returns
     -------
     None
     """
     # load iraf packages
     self.setiraf()
     
     iraf.delete(outfile)
     iraf.phot(image = self.sci_file + "[,," + str(framenum) + "]", coords = coords, output = outfile, fwhmpsf = self.fwhmpsf, sigma = self.sigma, readnoise = self.readnoise, epadu = self.epadu, exposure = self.exposure, calgorithm = self.calgorithm, cbox = self.cbox, maxshift = self.maxshift, salgorithm = self.salgorithm, annulus = self.annulus, dannulus = self.dannulus, apertures = apertures, zmag = self.zmag, interactive = "no", verify = "no", verbose = verbose)        
     
     return

    def _get_photometry(self):
        """ Get the photometry for the target.

        If the target is a standard star, aperture photometry will be performed. For the moment nothing is done with
        the others, but in due time (TODO) photometry.py will be included here. """

        basename = "standards"
        fd, coords_file = tempfile.mkstemp(prefix=basename, suffix=".coords")
        os.write(fd, "{0} {1} \n".format(self.RA, self.DEC))
        os.close(fd)

        if self.objtype == "standard":
            iraf.noao(_doprint=0)
            iraf.digiphot(_doprint=0)
            iraf.apphot(_doprint=0)
            seeing = self.header.hdr[self.header.seeingk]
            photfile_name = self.header.im_name + ".mag.1"
            utilities.if_exists_remove(photfile_name)
            kwargs =  dict(output=photfile_name, coords=coords_file,
                      wcsin='world', fwhm=seeing, gain=self.header.gaink, exposure=self.header.exptimek,
                      airmass=self.header.airmassk, annulus=6*seeing, dannulus=3*seeing,
                      apert=2*seeing, verbose="no", verify="no", interac="no")
            iraf.phot(self.header.im_name, **kwargs)
            [counts] = iraf.module.txdump(photfile_name, 'FLUX', 'yes', Stdout=subprocess.PIPE)
            utilities.if_exists_remove(coords_file)
            return float(counts)

def getSkyMeanSDinAnnulus(ann,delta=5):
	iraf.noao()
	iraf.digiphot()
	iraf.apphot()
	iraf.photpars.setParam('apertures','1')
	iraf.phot.setParam('interactive','no')
	iraf.phot.setParam('image',fitsDir+fitsFile)
	iraf.phot.setParam('coords',fitsDir+fitsFile+".coo")
	outmag=".maglim"
	try:	
		os.remove(fitsDir+fitsFile+outmag)
	except:
		print "File does not exist BEFORE running phot, so no need to delete."
	iraf.phot.setParam('output',fitsDir+fitsFile+outmag)
	iraf.phot.setParam('interac','no')	
	iraf.fitskypars.setParam('annulus',str(ann))	
	iraf.fitskypars.setParam('dannulus',str(delta))
	## NO SIGMA CLIPPING! JUST TO BE SAFE: (6/2013)
	iraf.fitskypars.setParam('sloclip',"0")
	iraf.fitskypars.setParam('shiclip',"0")

	iraf.phot(fitsDir+fitsFile,Stdin=cr)	
	aa, nn, xx, ss = ao.readPhotMagFile(fitsDir,fitsFile,outmag)
	try:	
		os.remove(fitsDir+fitsFile+outmag)
	except:
		print "File not found to delete AFTER running phot; that's odd."
	return xx

def daofind(filelist):
	iraf.datapars.setParam("exposure", "EXPTIME")
	iraf.datapars.setParam("airmass", "AIRMASS")
	iraf.datapars.setParam("filter", "INSFILTE")
	iraf.datapars.setParam("ccdread", "RDNOISE")
	iraf.datapars.setParam("gain", "GAIN")

	iraf.centerpars.setParam("calgorithm", "centroid")

	iraf.photpars.setParam("apertures", "17.56")
	
	iraf.fitskypars.setParam("annulus", "20.56")
	iraf.fitskypars.setParam("dannulus", "5")
	
	

	with open(filelist) as file1:
		for line in file1:
			filename = line.strip()
			print("Filename %s"  % filename)
			y,z = np.loadtxt(filename + "-daoedit", usecols=[3,4], unpack=True)
			sigma = np.mean(y) 
			fwhm = np.mean(z)
			iraf.datapars.setParam("fwhmpsf", "%s" % str(fwhm)  )
			iraf.datapars.setParam("sigma", "%s" % str(sigma)  )
			iraf.daofind.setParam("mode", "h")
			iraf.phot.setParam("mode", "h")
			iraf.datapars.setParam("mode", "h")
			iraf.centerpars.setParam("mode", "h")
			iraf.fitskypars.setParam("mode", "h")
			iraf.photpars.setParam("mode", "h")
	
			iraf.daofind(filename, filename + "-daofind", mode='h')
			
			iraf.phot(filename, coords=(filename + "-daofind"), output=(filename+".mag"), mode='h')

def psffit2(img, fwhm, psfstars, hdr, _datamax, psffun='gauss',  fixaperture=False):
    ''' 
    giving an image, a psffile,  calculate the magnitudes of strs in the file _psf.coo
    '''
    import lsc
    _ron = lsc.util.readkey3(hdr, 'ron')
    _gain = lsc.util.readkey3(hdr, 'gain')
    if not _ron:
        _ron = 1
        print 'warning ron not defined'
    if not _gain:
        _gain = 1
        print 'warning ron not defined'

    iraf.digiphot(_doprint=0)
    iraf.daophot(_doprint=0)
    zmag = 0.
    varord = 0  # -1 analitic 0 - numeric

    if fixaperture:
        print 'use fix aperture 5 8 10'
        hdr = lsc.util.readhdr(img+'.fits')
        _pixelscale = lsc.util.readkey3(hdr, 'PIXSCALE')
        a1, a2, a3, a4, = float(5. / _pixelscale), float(5. / _pixelscale), float(8. / _pixelscale), float(
                    10. / _pixelscale)
    else:
        a1, a2, a3, a4, = int(fwhm + 0.5), int(fwhm * 2 + 0.5), int(fwhm * 3 + 0.5), int(fwhm * 4 + 0.5)

    iraf.fitskypars.annulus = a4
    iraf.fitskypars.salgori = 'mean'  #mode,mean,gaussian
    iraf.photpars.apertures = '%d,%d,%d' % (a2, a3, a4)
    iraf.datapars.datamin = -100
    iraf.datapars.datamax = _datamax
    iraf.datapars.readnoise = _ron
    iraf.datapars.epadu = _gain
    iraf.datapars.exposure = 'EXPTIME'  
    iraf.datapars.airmass = ''
    iraf.datapars.filter = ''
    iraf.photpars.zmag = zmag
    iraf.centerpars.calgori = 'centroid'
    iraf.centerpars.cbox = a2
    iraf.daopars.recenter = 'yes'
    iraf.delete('_psf2.ma*', verify=False)

    iraf.phot(img+'[0]', '_psf2.coo', '_psf2.mag', interac=False, verify=False, verbose=False)

    iraf.daopars.psfrad = a4
    iraf.daopars.functio = psffun
    iraf.daopars.fitrad = a1
    iraf.daopars.fitsky = 'yes'
    iraf.daopars.sannulus = a4
    iraf.daopars.recenter = 'yes'
    iraf.daopars.varorder = varord
    iraf.delete("_als2,_psf.grp,_psf.nrj", verify=False)
    iraf.group(img + '[0]', '_psf2.mag', img + '.psf', '_psf.grp', verify=False, verbose=False)
    iraf.nstar(img + '[0]', '_psf.grp', img + '.psf', '_als2', '_psf.nrj', verify=False, verbose=False)
    photmag = iraf.txdump("_psf2.mag", 'xcenter,ycenter,id,mag,merr', expr='yes', Stdout=1)
    fitmag = iraf.txdump("_als2", 'xcenter,ycenter,id,mag,merr', expr='yes', Stdout=1)
    return photmag, fitmag

def get_fake_centroid(filename,x,y,instrument,filt):
    """
    Locate the center of a fake psf

    INPUTS: The fake-SN psf image in filename, the expected x,y position
    of the center of the psf, the instrument and filter being modeled.
    RETURNS: xcentroid, ycentroid, fluxcorr
    """
    from pyraf import iraf

    iraf.digiphot(_doprint=0)
    iraf.apphot(_doprint=0)
    iraf.unlearn(iraf.apphot.phot)
    iraf.unlearn(iraf.datapars)
    iraf.unlearn(iraf.centerpars)
    #Use the centroid algorithm right now as it seems more robust to geometric distortion.
    iraf.centerpars.calgorithm = 'centroid'
    iraf.centerpars.cbox = 5.0

    iraf.unlearn(iraf.fitskypars)
    iraf.unlearn(iraf.photpars)
    photparams = {
        'interac':False,
        'radplot':False,
        }
    iraf.datapars.readnoise = 0.0
    iraf.datapars.itime = 1.0
    iraf.datapars.epadu = 1.0

    # iraf.digiphot.apphot.fitskypars :
    iraf.unlearn(iraf.fitskypars)
    iraf.fitskypars.salgorithm = 'constant'
    iraf.fitskypars.skyvalue = 0.0

    # iraf.digiphot.apphot.photpars :
    iraf.unlearn(iraf.photpars)
    iraf.photpars.weighting = 'constant'
    iraf.photpars.apertures = 20   # TODO : set this more intelligently !
    iraf.photpars.zmag = 25
    iraf.photpars.mkapert = False

    #Write the coordinate file starting as position x and y
    coxyfile = 'centroid.xycoo'
    coxy = open(coxyfile, 'w')
    print >> coxy, "%10.2f  %10.2f" % (x,y)
    coxy.close()
    if os.path.exists('centroid.mag'): os.remove('centroid.mag')
    iraf.phot(image=filename, skyfile='', coords=coxyfile, output='centroid.mag',
              verify=False, verbose=True, Stdout=1, **photparams)
    f = open('centroid.mag', 'r')
    maglines = f.readlines()
    f.close()
    xcentroid = float(maglines[76].split()[0])
    ycentroid = float(maglines[76].split()[1])

    return xcentroid,ycentroid

def photometry(color, images):
	#grab the aligned images
	#images=sorted(glob.glob('reduced/'+color+'align/*fits'))

	for i in images:
		iraf.phot(i, coords='reduced/'+color+'align/coords.lis',
			fwhmpsf=3,datamin=-500.,datamax=10000.,readnoise=20.,epadu=7.2,
			  exposur="EXPTIME",airmass="SECZ",obstime="JD",filter="IRFLTID",annulus=15,
		  	  dannulu=7, apertur="9",verify='no',output=i+".MC")
	return

def psffit2(img, fwhm, psfstars, hdr, _datamax=45000, psffun='gauss', fixaperture=False):
    import agnkey

    iraf.digiphot(_doprint=0)
    iraf.daophot(_doprint=0)
    zmag = 0.
    varord = 0  # -1 analitic 0 - numeric
    if fixaperture:
        print 'use fix aperture 5 8 10'
        hdr = agnkey.util.readhdr(img+'.fits')
        _pixelscale = agnkey.util.readkey3(hdr, 'PIXSCALE')
        a1, a2, a3, a4, = float(5. / _pixelscale), float(5. / _pixelscale), float(8. / _pixelscale), float(
                    10. / _pixelscale)
    else:
        a1, a2, a3, a4, = int(fwhm + 0.5), int(fwhm * 2 + 0.5), int(fwhm * 3 + 0.5), int(fwhm * 4 + 0.5)

    _center='no'
    iraf.fitskypars.annulus = a4
    iraf.fitskypars.dannulus = a4
    iraf.noao.digiphot.daophot.daopars.sannulus = int(a4)
    iraf.noao.digiphot.daophot.daopars.wsannul = int(a4)
    iraf.fitskypars.salgori = 'mean'  #mode,mean,gaussian
    iraf.photpars.apertures = '%d,%d,%d' % (a2, a3, a4)
    #    iraf.photpars.apertures = '%d,%d,%d'%(a2,a3,a4)
    iraf.datapars.datamin = -100
    iraf.datapars.datamax = _datamax
    iraf.datapars.readnoise = agnkey.util.readkey3(hdr, 'ron')
    iraf.datapars.epadu = agnkey.util.readkey3(hdr, 'gain')
    iraf.datapars.exposure = 'exptime'  #agnkey.util.readkey3(hdr,'exptime')
    iraf.datapars.airmass = 'airmass'
    iraf.datapars.filter = 'filter2'
    iraf.centerpars.calgori = 'gauss'
    iraf.centerpars.cbox = 1
    iraf.daopars.recenter = _center
    iraf.photpars.zmag = zmag

    iraf.delete('_psf2.ma*', verify=False)

    iraf.phot(img+'[0]', '_psf2.coo', '_psf2.mag', interac=False, verify=False, verbose=False)

    iraf.daopars.psfrad = a4
    iraf.daopars.functio = psffun
    iraf.daopars.fitrad = a1
    iraf.daopars.fitsky = 'yes'
    iraf.daopars.sannulus = int(a4)
    iraf.daopars.wsannul = int(a4)
    iraf.daopars.recenter = _center
    iraf.daopars.varorder = varord
    iraf.delete("_als,_psf.grp,_psf.nrj", verify=False)
    iraf.group(img+'[0]', '_psf2.mag', img + '.psf', '_psf.grp', verify=False, verbose=False)
    iraf.nstar(img+'[0]', '_psf.grp', img + '.psf', '_als', '_psf.nrj', verify=False, verbose=False)
    photmag = iraf.txdump("_psf2.mag", 'xcenter,ycenter,id,mag,merr', expr='yes', Stdout=1)
    fitmag = iraf.txdump("_als", 'xcenter,ycenter,id,mag,merr', expr='yes', Stdout=1)
    return photmag, fitmag

def recentre(image,refcoordfile):
    """Returns improved shift by centroiding
on the reference star using phot. This can be VERY
sensitive to the parameters in centerpars."""
    xin,yin = load(refcoordfile,unpack=True)
    try:
        iraf.phot(image,refcoordfile,'temp.mag',inter="no",calgorithm='centroid',
              mode='h',verify='no',update='no',verbose='no')
        xout,yout=iraf.pdump('temp.mag','xcen,ycen','yes',Stdout=1)[0].split()
    except:
        print "Recentring failed on", image
        return 0.,0.
    xout,yout = float(xout),float(yout)
    return xout-xin,yout-yin

def apply_phot(args):
    # Temporary file
    with tempfile.NamedTemporaryFile(suffix=".mag.1") as fd:
        phot_modstars = fd.name

    with tempfile.NamedTemporaryFile(suffix=".mag.1") as fd2:
        phot_subtstars = fd2.name

    hdr = args.image.header  # for short
    seeing, sigma = utils.get_from_header(args.image.im_name, hdr.seeingk, hdr.sigmak)
    iraf.phot(
        args.image.im_name,
        output=phot_modstars,
        coords=args.model_stars,
        wcsin=args.coords,
        fwhmpsf=seeing,
        sigma=sigma,
        ccdread=hdr.ccdronk,
        gain=hdr.gaink,
        exposure=hdr.exptimek,
        airmass=hdr.airmassk,
        annulus=6 * seeing,
        dannulus=3 * seeing,
        apert=3 * seeing,
        verbose="no",
        verify="no",
        interac="no",
    )

    iraf.phot(
        args.image.im_name,
        output=phot_subtstars,
        coords=args.subt_stars,
        wcsin=args.coords,
        fwhmpsf=seeing,
        sigma=sigma,
        ccdread=hdr.ccdronk,
        gain=hdr.gaink,
        exposure=hdr.exptimek,
        airmass=hdr.airmassk,
        annulus=6 * seeing,
        dannulus=3 * seeing,
        apert=3 * seeing,
        verbose="no",
        verify="no",
        interac="no",
    )
    return phot_modstars, phot_subtstars

def p48coaddphot(inlis, refstars, ot, filter, outfile):

    outf = open(outfile, "w")
    
    for image in inlis:

        root = image.split(".")[0]
        
        # Get the seeing in the image
        iqpkg.iqobjs(image, SIGMA, get_head(image, "SATURATE"), skyval="0.0")
        seepix = get_head(image, "SEEPIX")

        # Calculate zeropoint
        refstars.wcs2pix(image)
        refstars.set_mag(filter)
        xyfile = open("%s.xy" % root, "w")
        for star in refstars:
            xyfile.write("%10.3f%10.3f\n" % (star.xval, star.yval))
        xyfile.close()
        iraf.phot(image, coords="%s.xy" % root, output="%s.mag" % root,
                  aperture=1.2*float(seepix), interac=no)
        stars = Starlist("%s.mag" % root)
        zp, zpu = stars.zeropt(refstars,method="mean",rejout=0)

        # Measure source
        ot.wcs2pix(image)
        coofile = open("%s.coo" % root, "w")
        coofile.write("%10.3f%10.3f\n" % (ot[0].xval, ot[0].yval))
        coofile.close()
        iraf.phot(image, coords="%s.coo" % root, output="%s.ot" % root,
                  aperture=1.2*float(seepix), calgorithm="none", interac=no)
        stars = Starlist("%s.ot" % root)

        # Just return 99 for non-detection
        if len(stars)==1:
            smag = stars[0].mag
            smagu = stars[0].magu
            mag = smag + zp
        else:
            mag = 99.0
            smagu = 99.0
        mjdobs = DateTimeFrom(get_head(image, OBSDATEKEY)).mjd
        exptime = get_head(image, EXPTIMEKEY)
        outf.write("%15.3f%10.1f%10.3f%10.3f%10.3f\n" % (mjdobs, exptime, mag, smagu, zpu))

    outf.close()

def simplepsfphot(image,coords,psf,refstar,centre=True,vary=False):
    """PSF photometry, with a given PSF file in psf used for every image"""
    iraf.dele('temp.mag*')
    iraf.dele('temp.als')
    iraf.dele('temp.sub.fits')
    if centre:
        xsh,ysh = recentre(image,refstar) 
        print "Fine Centring: ", xsh,ysh
    else: xsh,ysh = 0,0
    if vary:
        setaperture(image,refstar)
    shift_file_coords(coords,xsh,ysh,'tempcoords2',sort='als')
    iraf.phot(image,'tempcoords2','temp.mag2',inter="no",calgorithm='none',
                  mode='h',verify='no',update='no',verbose='no')
    iraf.allstar(image,'temp.mag2',psf,'temp.als','temp.mag.arj','temp.sub.fits',
                 mode='h',verify='no',update='no',verbose='no')
    out = iraf.pdump('temp.als','id,mag,merr,msky','yes',Stdout=1)
    return out   

def createMagFile(inputDir, inputFile, useCoordFile, outputFile, fullyAutomatic=False, aperturesString='5,10,15,20', bestSkyAperture='100'):	
	cr=['\n']
	iraf.cd(inputDir)
	iraf.noao()
	iraf.digiphot()
	iraf.apphot()
	iraf.phot.setParam('interactive','no')
	iraf.centerpars.setParam('cbox','5')
	iraf.photpars.setParam('apertures',aperturesString)
	iraf.fitskypars.setParam('annulus',bestSkyAperture)
	iraf.fitskypars.setParam('dannulus','10')
	iraf.phot.setParam('image',inputFile)
	iraf.phot.setParam('coords',useCoordFile)
	iraf.phot.setParam('output',outputFile)
	if fullyAutomatic: ## If you use this option, it will not prompt you for any of the aperture/sky settings
		iraf.phot(mode='h',Stdin=cr)
	else: # Otherwise, you can change them at will, and get to press enter a bunch
		iraf.phot(mode='h')

def apphot(image,coords,refstar=None,centre=False,vary=False):
    """Apperture photometry with centering based on a reference star.
NB: centre refers to shifting the coordinates by centroiding on the
reference star; recentering on the final phot depends on
centerpars.calgorithm ."""
    iraf.dele('temp.mag*')
    if centre:
        xsh,ysh = recentre(image,refstar)
        print "Fine centring: ", xsh,ysh
    else: #no recentreing by reference star (but could still have calgorithm!=none)
        xsh,ysh = 0,0
    if vary:
        setaperture(image,refstar)
    shift_file_coords(coords,xsh,ysh,'tempcoords')
    iraf.phot(image,'tempcoords','temp.mag2',inter="no",
                  mode='h',verify='no',update='no',verbose='no')
    out = iraf.pdump('temp.mag2','id,flux,msky,stdev','yes',Stdout=1)
    return out

def run_daophot(image, outfile='default', coordfile='NA', apertures='5.0,16.66', annulus=17.0, dannulus=3.0, calgorithm='centroid', salgorithm='median', fwhmpsf=2.5, backsigma=-1.0,rdnoise=5.2):

	'''THIS PROCEDURE RUNS DAOPHOT ON INPUT IMAGE'''

	# Parse input parameters
	if outfile == 'default': outfile = image + '0.mag.1'
	if coordfile == 'NA': coordfile = image + '0.coo.1'

	# Read in fits header
	f = pyfits.open(image)
	fheader = f[0].header
	f.close()

	# Extract relevant info from the header
	exptime = fheader['texptime']
	filter = fheader['FILTER2']
	ipxscl = fheader['D001ISCL']
	opxscl = fheader['D001SCAL']
	num_flts = float(fheader['NDRIZIM'])/2.0
	dp_zmag = -2.5 * np.log10(float(fheader['PHOTFLAM'])) + float(fheader['PHOTZPT'])       # zero-pt for each filter
	#df_max = 10000000.      # upper limit for "good" pixels (sometimes used to ID bad pixels). Not using this parameter right now.

	# Perform read noise correction
	rdnoise_corr = np.sqrt(num_flts * (rdnoise * opxscl/ipxscl)**2)

	# Perform background noise calculation
	if backsigma < 0.0:
		backrms=iraf.imstatistics(image+'[0]', fields='stddev', nclip=10, lsigma=3.0, usigma=3.0, cache='yes', format='no',Stdout=1)
		backsigma=float(backrms[0])


	'''Run daophot'''
	# Remove old phot output files
	file_query = os.access(outfile, os.R_OK)      
	if file_query == True: os.remove(outfile)

	# Run phot
	iraf.phot.unlearn()         # reset daophot parameters to default values
	iraf.phot(image=image+'[0]', interactive='no', verify='no', coords=coordfile, output=outfile, fwhmpsf=fwhmpsf, \
      		sigma=backsigma, readnoise=rdnoise_corr, itime=exptime, calgorithm=calgorithm, salgorithm=salgorithm, \
      		annulus=annulus, dannulus=dannulus, apertures=apertures,zmag=dp_zmag)


	return outfile 		# return name of output catalog

def calc_curve_of_growth(dataDir, image, apertures):
    aper = apertures

    # ==========
    # Make a coordinate file, PSF is in the middle
    # ==========
    # get image size
    data = pyfits.getdata(dataDir + image + '_psf.fits')
    imgSize = data.shape[0]
    
    # Write out coo file
    cooFile = image + '_psf.coo'
    _coo = open(cooFile, 'w')
    _coo.write('%.2f %.2f\n' % (imgSize/2.0, imgSize/2.0))
    _coo.close()
    
    
    # #########
    # Setup phot for running on PSF
    # #########
    ph.setup_phot(dataDir + image + '_psf', apertures=aper)
    ir.fitskypars.salgorithm = 'constant'
    ir.fitskypars.skyvalue = 0.0
    
    # Output into current directory, not data directory
    input = dataDir + image + '_psf.fits'
    output = workDir + image + '_psf.phot.mag'

    ir.phot(input, cooFile, output, verbose="no")
    
    (radius, flux, mag, merr) = ph.get_phot_output(output, silent=True)
    diff = np.zeros(len(radius), dtype=float)
    diff[1:] = mag[1:] - mag[:-1]
    
    # Save output
    pfile = open(workDir + 'cog_' + image + '.dat', 'w')
    pickle.dump(radius, pfile)
    pickle.dump(flux, pfile)
    pickle.dump(mag, pfile)
    pickle.dump(merr, pfile)
    pickle.dump(diff, pfile)
    pfile.close()

    return (radius, flux, mag, merr, diff)

def run_phot(imageRoot, silent=False,
             apertures=[25,50,75,100,125,150,175,200],
             sky_annulus=200, sky_dannulus=50, zmag=0):

    setup_phot(imageRoot, apertures=apertures, zmag=zmag, silent=silent,
               sky_annulus=sky_annulus, sky_dannulus=sky_dannulus)

    image = imageRoot + '.fits'
    coords = imageRoot + '.coo'

    # Output into current directory, not data directory
    rootSplit = imageRoot.split('/')
    output = rootSplit[-1] + '.phot.mag'

    ir.phot(image, coords, output)

    (radius, flux, mag, merr) = get_phot_output(output, silent=silent)

    return (radius, flux, mag, merr)

	def doPhotometryACS(self, file, inputcoords, apertures, zeropoint, bgsigma, skyann, skydann):
		'''
		This function can be used to do photometry from a given image. Input parameters
		can be varied, however, the task assumes that photometry is done from an ACS
		image and the exposure time is found from the header of the file. 
		Object recentering is done with centroid algorithm and shifts up to 6 pixels are
		possible. For skyfitting algorithm mode is adopted.
		'''
		#load packages, should supress the output
		I.digiphot()
		I.apphot()

		#setting up for ACS
		I.datapar(sigma = bgsigma, exposure = 'exptime', gain = 'ccdgain')
		I.centerpars(calgorithm = 'centroid', cbox = 10., maxshift = 6.)
		I.fitskypars(salgorithm = 'mode', annulus = skyann, \
					 dannulus = skydann, skyvalue = 0.)
		I.photpars(apertures = apertures, zmag = zeropoint)

		I.phot(file, coords=inputcoords, verify='no', verbose = 'no')

    def _get_photometry(self, seeing=None, aperture=None):
        """ Get the photometry for the target.

        If the target is a standard star, aperture photometry will be performed. For the moment nothing is done with
        the others, but in due time (TODO) photometry.py will be included here. """

        if self.objtype is not "standards":
            return None

        basename = "standards"
        fd, coords_file = tempfile.mkstemp(prefix=basename, suffix=".coords")
        os.write(fd, "{0} {1} \n".format(self.RA, self.DEC))
        os.close(fd)


        # If aperture was not given by the user, try and use the seeing as a reference for default values
        if not aperture:
            try:
                seeing = self.header.hdr[self.header.seeingk]
                aperture = 3 * seeing
            except ValueError:  # keyword was not correctly guessed
                pass

        # If aperture exists, do the photometry using it
        if aperture:
            annulus = 2 * aperture
            dannulus = max([aperture, 3])  # minimum of 3 pixels thickness for the sky annulus
            fd, photfile_name = tempfile.mkstemp(".mag.1")
            utilities.if_exists_remove(photfile_name)
            kwargs =  dict(output=photfile_name, coords=coords_file, salgorithm='median',
                      wcsin='world', fwhm=seeing, gain=self.header.gaink, exposure=self.header.exptimek,
                      airmass=self.header.airmassk, annulus=annulus, dannulus=dannulus,
                      apertures=aperture, verbose="no", verify="no", interac="no")
            iraf.phot(self.header.im_name, **kwargs)
            [counts] = iraf.txdump(photfile_name, 'FLUX', 'yes', Stdout=subprocess.PIPE)
        else:
            sys.exit("\n \n Sorry, no aperture was passed by you, and a seeing keyword was not "
                     "found in the header. \n\n ")

        utilities.if_exists_remove(coords_file)
        return float(counts)