def identify_objects(fnlist,skysiglist,fwhmlist,suffix=".coo"):
    """Runs the IRAF routine 'daofind' to locate objects in a series of images,
    creating coordinate files.
    
    Inputs:
    fnlist -> List of strings, each the path to a fits image.
    skysiglist -> List of floats, each the sky background sigma for an image.
    fwhmlist -> List of floats, each the FWHM of objects in an image.
    suffix -> Suffix for the coordinate files. '.coo' by default.
    
    Outputs:
    coolist -> List of strings, each the path to the coordinate files created.
    
    """
    
    print "Identifying objects in images..."
    
    coolist = []
    
    #Open IRAF packages
    iraf.noao(_doprint=0)
    iraf.digiphot(_doprint=0)
    iraf.apphot(_doprint=0)
    for i in range(len(fnlist)):
        coolist.append(fnlist[i]+suffix)
        iraf.daofind(image=fnlist[i],
                     output=fnlist[i]+suffix,
                     fwhmpsf=fwhmlist[i],
                     sigma=skysiglist[i],
                     threshold=4.0,
                     datamin='INDEF',
                     datamax='INDEF',
                     verify='N')
    return coolist

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 _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 align_combine(fitsdir, myfilter, examine=True):
    from pyraf import iraf 
    
    iraf.noao(_doprint=0)
    iraf.digiphot(_doprint=0)
    iraf.apphot(_doprint=0)
    
    os.chdir(fitsdir)
    listfiles = glob.glob(myfilter)
    listfiles.sort()
    
    if (examine):
        print "Opening ",listfiles[0]," to examine."
        iraf.imexamine(input=listfiles[0], \
                    logfile="coords.dat", \
                    keeplog="yes")
        
        with open("align.list",'w') as f:
            for i in listfiles:
                f.write(i+"\n")
    
    print "Aligning with reference:",listfiles[0]
    iraf.imalign( input   =  "@align.list", referenc= listfiles[0], coords  =  "coords.dat", output  = "[email protected]")  
    
    listfiles = glob.glob("a_"+myfilter)
    listfiles.sort()
    with open("comb.list",'w') as f:
        for i in listfiles:
            f.write(i+"\n")
            
    print "Combining"        
    iraf.imcombine(input = "@comb.list",\
                   output = "out.fits",\
                   combine= "median")

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 set_default(base):
    # Loading necessary IRAF packages
    iraf.digiphot(_doprint=0)
    iraf.daophot(_doprint=0)
    iraf.apphot(_doprint=0)
    iraf.reset(min_lenuserarea='200000')

    iraf.datapars.setParam('datamin','INDEF')
    iraf.datapars.setParam('datamax','60000')
    iraf.datapars.setParam('exposure','EXPTIME')
    iraf.datapars.setParam('airmass', 'AIRMASS')
    iraf.datapars.setParam('filter', 'FILTER')

    iraf.findpars.setParam('threshold', 3.0)
    iraf.findpars.setParam('sharphi',  1)
    iraf.findpars.setParam('roundhi', 1.3)
    iraf.findpars.setParam('roundlo', -1.3)

    iraf.daofind.setParam('verify','no')
    iraf.daofind.setParam('interactive','no')

    iraf.photpars.setParam('zmag',25.0)
    iraf.photpars.setParam('weighti','constant')
    iraf.photpars.setParam('apertur',3.0)

    iraf.phot.setParam('output',base+'default')
    iraf.phot.setParam('coords',base+'default')
    iraf.phot.setParam('verify','no')
    iraf.phot.setParam('interactive','no')

    iraf.fitpsf.setParam('box',10.0)
    iraf.fitpsf.setParam('verify','no')
    iraf.fitpsf.setParam('interactive','no')

    iraf.centerpars.setParam('calgori','none')
    iraf.fitskypars.setParam('salgorithm','mode')

    iraf.daopars.setParam('functio','moffat15')
    iraf.daopars.setParam('varorde','0')
    iraf.daopars.setParam('nclean','0')
    iraf.daopars.setParam('saturat','no')
    iraf.daopars.setParam('fitsky','yes')
    iraf.daopars.setParam('recenter','yes')
    iraf.daopars.setParam('groupsk','yes')
    iraf.daopars.setParam('maxnsta','40000')

    iraf.psf.setParam('photfile',base+'default')
    iraf.psf.setParam('pstfile',base+'default')
    iraf.psf.setParam('psfimage',base+'default')
    iraf.psf.setParam('opstfile',base+'default')
    iraf.psf.setParam('groupfil',base+'default')
    iraf.psf.setParam('interac','no')
    iraf.psf.setParam('matchby','yes')
    iraf.psf.setParam('verify','no')
    iraf.psf.setParam('showplo','no')

    iraf.allstar.setParam('verify','no')

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 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 psf(args):
    """ Calculate the PSF of an image.
    """
    # Load iraf packages: phot, pstselect, psf will be needed
    iraf.noao(_doprint=0)
    iraf.digiphot(_doprint=0)
    iraf.module.daophot(_doprint=0)
    
    # Read the seeing and sigma of the sky from the header
    seeing, sigma = utils.get_from_header(args.input, args.FWHM_key, args.sigma)
    
    # Do photometry on the image
    #print "photometry: \n"
    photfile_name = args.input + ".mag.1"
    utils.if_exists_remove(photfile_name)
    iraf.module.phot(args.input, output=photfile_name, coords=args.stars, 
                     wcsin=args.coords, fwhm=seeing, 
                     sigma=sigma, datamax=args.maxval, datamin=args.minval,
                     ccdread=args.ron_key, gain=args.gain_key, exposure=args.expt_key,
                     airmass=args.airm_key, annulus=6*seeing, dannulus=3*seeing, 
                     apert=2*seeing, verbose="no", verify="no", interac="no")

    # Select stars on the image                 
    #print "pstselect: \n"
    pstfile_name = args.input + ".pst.1"
    utils.if_exists_remove(pstfile_name)
    iraf.module.pstselect(args.input, photfile=photfile_name, pstfile=pstfile_name,
                          maxnpsf=20, fwhm=seeing, sigma=sigma,
                       datamax=args.maxval, ccdread=args.ron_key, gain=args.gain_key,
                       exposure=args.expt_key,  function="auto", nclean=1, 
                       psfrad=2*seeing, fitrad=seeing, maxnstar=20, verbose="no",
                       verify="no")

    # Build psf of the stars
    #print "psf: \n"
    psffile_table = args.input + ".psf.1.fits"  # iraf keeps adding the .fits :(
    psgfile_name = args.input + ".psg.1"
    pstfile_name2 = args.input + ".pst.2"   
    utils.if_exists_remove(psffile_table,psgfile_name, pstfile_name2)
    iraf.module.psf( args.input, photfile=photfile_name, pstfile=pstfile_name,
                     groupfile=psgfile_name, opstfile=pstfile_name2,    
                     psfimage=psffile_table,fwhm=seeing, sigma=sigma, datamax=args.maxval, 
                     datamin=args.minval, ccdread=args.ron_key, gain=args.gain_key, 
                     exposure=args.expt_key, function="moffat25", nclean=1, 
                     psfrad=12, fitrad=seeing, maxnstar=20, interactive="no",
                     varorder=-1, verbose="no",verify="no")
                     
    # Use seepsf to build the image of the psf
    psffile_name = args.input + ".psf.fits" 
    utils.if_exists_remove(psffile_name)
    iraf.module.seepsf(psffile_table, psffile_name)

     
    return psffile_name

def findStars(frame,fwhm,sigma,thresh,extra="",overwrite=True,interactive=True):
	if overwrite:
		if os.path.isfile(frame+".fullcoo"):
			os.system("rm "+frame+".fullcoo")
	iraf.noao()
	iraf.digiphot()
	iraf.daophot()
	iraf.datapars.setParam("fwhmpsf",fwhm)
	iraf.datapars.setParam("sigma",sigma)
	iraf.findpars.setParam("thresh",thresh)
	foo = iraf.daofind(image=frame+extra,output=frame+".fullcoo",Stdin=cr,Stdout=1)

def aperture_photometry(image,photfilesuffix,threshold=4.,wcs='logical',mode='small', telescope='Pairitel'):
  '''perform aperture photometry on a given image
  
  This procedure perfrom IRAF/DAOPHOT aperture phtometry on (all extensions of) a fits image.
  It sets some global daophot parameters, then perfroms daophot.daofind and daophot.phot
  The output .mag file is filtered for valid magnitudes and a mag.sex file computed, where
  all coordinates are transformed to sexagesimal coordinates.
  
  input: image: filename
  keywords:  threshold=4 :daofind detection threshold in sigma_sky
             wcs='logical' :working wcs system, can be  'logical','physical' or 'tv'
                        in any case an additional sexagesimal output file will be computed
         mode='small' : 'standard' = 'large' or 'psf' = 'small' choses a pre-defined aperture size for 
                        aperture photometry of standards or just as starting point for psf photometry
  '''
  print 'Performing aperture photometry on '+image
  if telescope == 'Pairitel':
      set_Pairitel_params()
  elif telescope == 'FLWO':
      set_FLWO_params()
  else:
      print 'No parameters found for this telescope!'
  
  iraf.mscred(_doprint=0)
  iraf.digiphot(_doprint=0)
  iraf.daophot(_doprint=0)
  iraf.daophot.verify = False
  iraf.daophot.wcsin = wcs
  iraf.daophot.wcsout = wcs
  iraf.daophot.wcspsf = wcs
  iraf.centerpars.cbox=2.*iraf.datapars.fwhmpsf
  iraf.fitskypars.annulus=20 # usually one would use 5.*iraf.datapars.fwhmpsf
  #annulus must be wide enough to have good sky statistics
  iraf.fitskypars.dannulus=10. # this is standard
  if 'standard'.find(mode.lower()) != -1 or 'large'.find(mode.lower()) != -1:
    iraf.photpars.aperture = 4.*iraf.datapars.fwhmpsf
    iraf.centerpars.calgorithm = 'centroid'
  elif 'small'.find(mode.lower()) != -1 or 'psf'.find(mode.lower()) != -1:
    iraf.photpars.aperture = 4. # close in for crowded regions
    iraf.centerpars.calgorithm = 'none'
  else:
    print '''use: aperture_photometry(image,threshold=4.,wcs="logical",mode="small")
    
    The following mode keywords are implemented: standard, large, psf, small'''
    raise NotImplementedError
  #for imagebi in fits_ext(image):
  sky,skydev=get_sky(image)
  iraf.datapars.datamin=sky-5*skydev
  iraf.datapars.sigma=skydev      
  #this could be changed to use Stdin and Stdout so that fewer files are written (saves space and time), but for debug purposes better keep all that
  iraf.daofind(image, output='default', verbose=False, verify=False, threshold=threshold) 
  iraf.daophot.phot(image, coords='default', output=image+photfilesuffix, verbose=False, verify=False, interactive=False, wcsout=wcs)

def run_fitparams(weight=True):
    if weight: 
        weighting = 'photometric'
    else:
        weighting = 'uniform'

    ir.delete('standards.transParams')
    ir.digiphot()
    ir.photcal()
    ir.unlearn('fitparams')
    ir.fitparams('standards.instMag', 'photcal_ukirt_faint.dat',
                 'standards.config', 'standards.transParams',
                 weighting=weighting)

def set_FLWO_params():
  '''Set some daophot parameters, which are typical for my FLWO images
  '''
  iraf.digiphot(_doprint=0)
  iraf.daophot(_doprint=0)
  iraf.datapars.fwhmpsf = 4.5
  iraf.datapars.datamax=40000.
  iraf.datapars.ccdread = "rdnoise"
  iraf.datapars.gain = "gain"
  iraf.datapars.exposure = "exptime"
  iraf.datapars.airmass = "air"
  iraf.datapars.filter = "filter"
  iraf.datapars.obstime="date-obs"

def _dao_setup(fwhmpsf, threshold, sigma):
    iraf.digiphot(_doprint=0)
    iraf.digiphot.apphot(_doprint=0)
    datapars = iraf.noao.digiphot.apphot.datapars
    findpars = iraf.noao.digiphot.apphot.findpars
    datapars.scale = 1.0
    datapars.sigma = sigma * (3.0 / 4.0)
    datapars.fwhmpsf = fwhmpsf
    
    # not important for source finding
    datapars.readnoise = 0.0
    datapars.epadu = 1.0
    
    findpars.threshold = threshold # only care about brightest

def marksn2(img,fitstab,frame=1,fitstab2='',verbose=False):
   from pyraf import iraf
   from numpy import array   #,log10
   import lsc
   iraf.noao(_doprint=0)
   iraf.digiphot(_doprint=0)
   iraf.daophot(_doprint=0)
   iraf.images(_doprint=0)
   iraf.imcoords(_doprint=0)
   iraf.proto(_doprint=0)
   iraf.set(stdimage='imt1024')
   hdr=lsc.util.readhdr(fitstab)
   _filter=lsc.util.readkey3(hdr,'filter')
   column=lsc.lscabsphotdef.makecatalogue([fitstab])[_filter][fitstab]

   rasex=array(column['ra0'],float)
   decsex=array(column['dec0'],float)


   if fitstab2:
      hdr=lsc.util.readhdr(fitstab2)
      _filter=lsc.util.readkey3(hdr,'filter')
      _exptime=lsc.util.readkey3(hdr,'exptime')
      column=lsc.lscabsphotdef.makecatalogue([fitstab2])[_filter][fitstab2]
      rasex2=array(column['ra0'],float)
      decsex2=array(column['dec0'],float)

   iraf.set(stdimage='imt1024')
   iraf.display(img + '[0]',frame,fill=True,Stdout=1)
   vector=[]
   for i in range(0,len(rasex)):
      vector.append(str(rasex[i])+' '+str(decsex[i]))

   xy = iraf.wcsctran('STDIN',output="STDOUT",Stdin=vector,Stdout=1,image=img+'[0]',inwcs='world',units='degrees degrees',outwcs='logical',\
                         formats='%10.1f %10.1f',verbose='yes')[3:]
   iraf.tvmark(frame,'STDIN',Stdin=list(xy),mark="circle",number='yes',label='no',radii=10,nxoffse=5,nyoffse=5,color=207,txsize=2)

   if verbose:
 #     print 2.5*log10(_exptime)
      for i in range(0,len(column['ra0'])):
         print xy[i],column['ra0'][i],column['dec0'][i],column['magp3'][i],column['magp4'][i],column['smagf'][i],column['magp2'][i]

   if fitstab2:
      vector2=[]
      for i in range(0,len(rasex2)):
         vector2.append(str(rasex2[i])+' '+str(decsex2[i]))
      xy1 = iraf.wcsctran('STDIN',output="STDOUT",Stdin=vector2,Stdout=1,image=img+'[0]',inwcs='world',units='degrees degrees',outwcs='logical',\
                            formats='%10.1f %10.1f',verbose='yes')[3:]
      iraf.tvmark(frame,'STDIN',Stdin=list(xy1),mark="cross",number='yes',label='no',radii=10,nxoffse=5,nyoffse=5,color=205,txsize=2)

def set_Pairitel_params():
  '''Set some daophot parameters, which are typical for my Pairitel images
  '''
  iraf.digiphot(_doprint=0)
  iraf.daophot(_doprint=0)
  iraf.apphot(_doprint=0)
  iraf.datapars.fwhmpsf = 4.5
  iraf.datapars.datamax=15000.
  iraf.centerpars.cbox = 9.0
  iraf.centerpars.calgorithm = "none"
  iraf.datapars.datamin = "INDEF"
  iraf.datapars.ccdread="RDNOISE" # 30. in 2MASS?
  iraf.datapars.gain = "GAIN"
  iraf.datapars.exposure = "EXPTIME"
  iraf.datapars.airmass = "AIRMASS"
  iraf.datapars.filter = "FILTER"
  iraf.datapars.obstime="DATE-OBS"

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 make_iraf_psf(fname):

    with fits.open(fname) as stp:  # Read in image info placed in header previously
        header = stp[0].header
        ra = header["RA"]
        dec = header["DEC"]
        for ext in range(len(stp)):  # Look in each header section
            try:
                header = stp[ext].header
                expnum = header["EXPNUM"]
            except:
                continue

    stamp = fits.open(fname)
    for ext in range(
        len(stamp)
    ):  # Establish size of image and location of source in image (source should be centered in cutout)
        try:
            world_coords = wcs.WCS(stamp[ext].header)
            image_coords = world_coords.sky2xy(float(ra), float(dec))
            image_size = stamp[ext].data.shape
        except:
            continue

    # Some cutouts are split over multiple CCDs, meaning the source may not actually be in some of these small image 'pieces'
    if 0 < image_coords[0] < image_size[0] and 0 < image_coords[1] < image_size[1]:

        ccdnum = stamp[ext].header.get("EXTNO", None)

        psf_file = storage.get_file(expnum, ccd=ccdnum, ext="psf.fits")  # Get psf file for image from storage

        iraf.noao()
        iraf.digiphot()
        iraf.daophot(_doprint=0)
        psf_image = psf_file.replace("psf", "seepsf")

        # Seepsf creates psf_image
        iraf.seepsf(psf_file, psf_image, xpsf=image_coords[0], ypsf=image_coords[1])
        os.unlink(psf_file)

        return psf_image

    else:
        return None

def init_iraf():
    """Initializes the pyraf environment. """
    
    # The display of graphics is not used, so skips Pyraf graphics 
    # initialization and run in terminal-only mode to avoid warning messages.
    os.environ['PYRAF_NO_DISPLAY'] = '1'

    # Set PyRAF process caching off to avoid errors if spawning multiple 
    # processes.
    iraf.prcacheOff()

    # Load iraf packages and does not show any output of the tasks.
    iraf.digiphot(_doprint = 0)
    iraf.apphot(_doprint = 0) 
    iraf.images(_doprint = 0)

    # Set the iraf.phot routine to scripting mode.
    iraf.phot.interactive = "no"
    iraf.phot.verify = "no"

def doPhot():
	print 'Running doPhot'
	try:
		mapFile = open('map.dat')
	except:
		print 'map.dat not found. Exiting.'
		sys.exit()
	try:
		changeFile = open('change.dat')
	except:
		print 'change.dat not found. Exiting.'
		sys.exit()

	coordList = list()
	for line in mapFile:
		coords = map(float,line.split())
		coordList.append(coords)
	coordList = np.array(coordList)

	iraf.noao()
	iraf.digiphot()
	iraf.daophot()
	iraf.ptools()
	iraf.set(clobber='yes')
	
	photFile = open('phot.dat', 'w')
	for line in changeFile.readlines():
		elms = line.split()
		imageName = elms[0]
		changeCoords = np.array([float(elms[1]),float(elms[2])])
		newCoords = coordList + changeCoords
		print 'Image: '+imageName
		coordFile = makeCoordFile(newCoords)
		
		iraf.noao.digiphot.daophot.phot(image=imagesDir+imageName, coords=coordFile, output='.temp-phot', skyfile='', verify='no', fwhmpsf=fwhmpsf, sigma=sigma, readnoise=readnoise, epadu=epadu, exposure=exposureHeader, obstime=obstimeHeader, calgorithm=calgorithm, cbox=cbox, apertures=apertures, annulus=annulus, dannulus=dannulus)

		result = iraf.noao.digiphot.ptools.txdump(Stdout=1, textfiles='.temp-phot', fields='mag, merr, otime', expr='yes')
		writeString = result[0].split()[-1] +' '+ ' '.join([' '.join(x.split()[:2]) for x in result])
		photFile.write(writeString+"\n")
		
	photFile.close()
	raw_input('doPhot Done. Hit return key to continue.')