def __init__(self, input, **kwargs):
        """
        Init.

        :param input: FITS association
        :type input: string
        :param kwargs: additional keyword arguments
        :type kwargs: dictionary
        """
        self.input = input
        self.settings = dict(asndir='asn', rawdir='raw',
                             jref='/grp/hst/cdbs/jref/',
                             mtab='/grp/hst/cdbs/mtab/',
                             sourceImage='POL0V_drz.fits',
                             sigma=5.0)
        self.settings.update(kwargs)

        #add env variables to both system and IRAF
        os.putenv('jref', self.settings['jref'])
        os.putenv('mtab', self.settings['mtab'])
        iraf.set(jref=self.settings['jref'])
        iraf.set(mtab=self.settings['mtab'])

        #IRAF has funny True and False
        self.yes = iraf.yes
        self.no = iraf.no

def get_seeing(file, scale, ref):
    coords = file.replace('fits','dao')
    log = file.replace('fits','psf')
    out = file.replace('fits','seeing')

    iraf.noao()
    iraf.noao.obsutil()
    iraf.set(stdgraph="uepsf")

    data = iraf.psfmeasure(file, coords='markall', wcs='logical', display='no', frame=1, level=0.5, size='FWHM', radius=10.0, sbuffer=1.0, swidth=3.0, iterations=1, logfile=log, imagecur=coords, Stdout=1)

    fwhm = data.pop().split().pop()

    fwhm_pix = float(fwhm)
    fwhm = fwhm_pix*scale
    seeing = dimm_seeing(fwhm_pix, ref, scale)
    print "Seeing: ", seeing

    fp = open(out, 'w')
    fp.write("%f %f\n" % (fwhm_pix, seeing))
    fp.close()
    os.system("echo \"image;text 85 500 # text={Spot FWHM = %5.2f pixels}\" | xpaset WFS regions" % fwhm_pix)
    os.system("echo \'image;text 460 500 # text={Seeing = %4.2f\"}\' | xpaset WFS regions" % seeing)
    os.system("echo \"set wfs_seeing %4.2f\" | nc hacksaw 7666" % seeing)
    # New function call to insert seeing value into MySQL database.
    fits2mysql(file, mode, fwhm, seeing)

 def initialize_instrument(self, output_dir):
     self.logger.info("Initializing instrument.")
     inst_file_name = output_dir + '/tmp/cp.dat'
     if self.debug:
         iraf.set(debug=1)
     iraf.set(use_new_imt='no')
     iraf.noao.imred(_doprint=0)
     iraf.noao.imred.ccdred(_doprint=0)
     if not os.path.exists(inst_file_name):
         with open(inst_file_name, 'w+') as inst_file:
             inst_file.write('subset          FILTER\n\n')
             inst_file.write('darktime        EXPTIME\n\n')
             inst_file.write('\'Dark Frame\'    dark\n')
             inst_file.write('\'Bias Frame\'    zero\n')
             inst_file.write('\'Light Frame\'   object\n')
             inst_file.write('\'Flat Field\'    flat\n')
     iraf.noao.imred.ccdred.setinst(
         instrument = 'cp',
         review = 'no',
         mode = 'h',
         dir = output_dir + '/tmp/',
         site = ''
     )
     iraf.noao.digiphot(_doprint=0)
     iraf.noao.digiphot.apphot(_doprint=0)
     iraf.dataio(_doprint=0)
     iraf.noao.digiphot.ptools(_doprint=0)

    def display(self,extn=None):
        """Display SCI extensions from the Astrodata object
           and mark objects and reference stars found in the
           input file tables OBJCAT and REFCAT.

           DS9 must be running.
        """
        try:
            from stsci.numdisplay import display as ndisplay
        except ImportError:
            from numdisplay import display as ndisplay

        #file = self.file
        #ad = AstroData(file)

        ad = self.adout
        nscext = ad.count_exts('SCI')

        if extn != None:
            si,ei = extn,extn+1
        else:
            si,ei = 1,nscext+1
            
        fr = 1
        for xt in range(si,ei):
            #ff = '%s[%d]' % (file,xt)
            ff = ad['SCI',xt].data
            
            iraf.set (stdimage='imt47')
            ndisplay(ff, frame=fr)
    
            xtnad = ad['OBJCAT',xt]
            if xtnad:
                x = xtnad.data.field('x')
                y = xtnad.data.field('y')
                print 'Marking %d objects with red.' % len(x)
                np.savetxt('/tmp/xy.txt', zip(x,y), fmt='%.2f', delimiter=' ')
                tvmark(frame=fr,coords='/tmp/xy.txt',color=204)
            else:
                print "No stars found in extension %d."%xt
        
            xtnad = ad['REFCAT',xt]
            if xtnad:
                x = xtnad.data.field('x')
                y = xtnad.data.field('y')
                print 'Marking %d refstars with green.' % len(x)
                np.savetxt('/tmp/xy.txt', zip(x,y), fmt='%.2f', delimiter=' ')
                tvmark(frame=fr,coords='/tmp/xy.txt',color=205)
            else:
                print "No reference stars found in extension %d."%xt

            fr += 1

    def build(self, fwhm=None, maxlin=40000, zeropoint=None):
        """

        @param fwhm:  initial guess of the FWHM, used for selecting sources that might be stars
        @param thresh: the detection threshold for stars (in sky sigma)
        @param maxlin: maximum linearity of the CCD, brighter sources are saturated.
        @param zeropoint: the cts/exptime mag zeropoint of the image
        @return:
        """
        if not zeropoint:
            zeropoint = storage.get_zeropoint(self.expnum,
                                              self.ccd,
                                              self.prefix,
                                              self.version)
        if not fwhm:
            fwhm = storage.get_fwhm(self.expnum,self.ccd,prefix=self.prefix, version=self.version)

        filename = storage.get_image(self.expnum,
                                     self.ccd,
                                     version=self.version,
                                     prefix=self.prefix)
        basename = os.path.splitext(filename)[0]

        iraf.set(writepars=0)

        iraf.kinky
        iraf.pipeline
        phot_filename = "{}.{}".format(basename, PHOT_EXT)
        psf_filename = "{}.{}".format(basename, PSF_EXT)
        zeropt_filename = "{}.{}".format(basename, ZERO_PT_USED)
        with open(zeropt_filename, 'w') as z:
            z.write("{}\n".format(zeropoint))

        apmax = AP_OUT_SCALE * fwhm
        swidth = fwhm

        # Turn off verbosity and parameter checking in DAOPHOT
        iraf.daophot.verify=iraf.no
        iraf.daophot.verbose=iraf.no
        iraf.daophot.update=iraf.no
        iraf.psf.showplots=iraf.no
        iraf.psf.interactive=iraf.no

        iraf.centerpars.calgori="centroid"
        iraf.centerpars.maxshift=3

        # set the basic 'datapars' for the images.
	    iraf.datapars.exposure=""
        # Note the itime is set to 1 as the SGwyn ZP includes the exposure time.
	    iraf.datapars.itime=1

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 buildtmc(tmcname):
    from pyraf import iraf
    from iraf import stsdas,hst_calib,synphot
    out=open('buildtmc.log','w')
    f=pyfits.open(tmcname)
    flist=f[1].data.field('filename')
    iraf.set(crrefer='./') #work locally

    for k in range(len(flist)):
        oldname=iraf.osfn(flist[k]).split('[')[0]
        newname=fincre(oldname)
        if os.path.exists(newname):
            flist[k]=fincre(flist[k])
        else:
            out.write("%s: no change necessary\n"%oldname)
    f.writeto(tmcname.replace(".fits","_new.fits"))
    out.close()

def atmofile(imgstd, imgout=''):
    # print "LOGX:: Entering `atmofile` method/function in %(__file__)s" %
    # globals()
    from pyraf import iraf
    import os
    import ntt

    iraf.noao(_doprint=0)
    iraf.onedspec(_doprint=0)
    iraf.set(direc=ntt.__path__[0] + '/')
    _cursor = 'direc$standard/ident/cursor_sky_0'
    if not imgout:
        imgout = 'atmo_' + imgstd
    os.system('rm -rf ' + imgout)
    iraf.noao.onedspec.bplot(imgstd, cursor=_cursor,
                             spec2=imgstd, new_ima=imgout, overwri='yes')
    return imgout

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.')

 def setiraf(self):
     """
     Set IRAF global parameters and load DAOPHOT package for aperture
     photometry.
     
     Parameters
     ----------
     
     Returns
     -------
     None
     """
     iraf.prcacheOff()
     iraf.set(writepars=0)
 
     # Load IRAF packages
     iraf.noao(_doprint = 0)
     iraf.noao.digiphot(_doprint = 0)
     iraf.noao.digiphot.daophot(_doprint = 0)
 
     return

def get_seeing(file, scale, ref):

#    log = file.replace('fits','psf')		# never used

# Gets called with full path to the fits file
#    print "daofind get_seeing: ", file

    iraf.noao()
    iraf.noao.obsutil()
    iraf.set(stdgraph="uepsf")

    # name of file with coordinates of star images.
    coords = file.replace('fits','dao')

    data = iraf.psfmeasure(file, coords='markall', wcs='logical', display='no', frame=1, level=0.5, size='FWHM', radius=10.0, sbuffer=1.0, swidth=3.0, iterations=1, logfile=log, imagecur=coords, Stdout=1)

    # print "daofind psfmeasure data: ", data
    # last = data.pop()
    # data.append(last)
    # print "daofind psfmeasure last line of data: ", last
    # psfmeasure returns lots of information (data for each spot), but we just
    # want the last line (the average fwhm), which looks like:
    # Average full width at half maximum (FWHM) of 6.0163

    fwhm = data.pop().split().pop()
    print "daofind psfmeasure fwhm: ", fwhm

    fwhm_pix = float(fwhm)
    fwhm = fwhm_pix*scale
    seeing = dimm_seeing(fwhm_pix, ref, scale)
    print "Seeing: ", seeing

    seeingfile = file.replace('fits','seeing')
    fp = open(seeingfile, 'w')
    fp.write("%f %f\n" % (fwhm_pix, seeing))
    fp.close()

    os.system("echo \"image;text 85 500 # text={Spot FWHM = %5.2f pixels}\" | xpaset WFS regions" % fwhm_pix)
    os.system("echo \'image;text 460 500 # text={Seeing = %4.2f\"}\' | xpaset WFS regions" % seeing)
    os.system("echo \"set wfs_seeing %4.2f\" | nc hacksaw 7666" % seeing)

def get_seeing(file, scale, ref):
    coords = file.replace("fits", "dao")
    log = file.replace("fits", "psf")
    out = file.replace("fits", "seeing")

    iraf.noao()
    iraf.noao.obsutil()
    iraf.set(stdgraph="uepsf")

    data = iraf.psfmeasure(
        file,
        coords="markall",
        wcs="logical",
        display="no",
        frame=1,
        level=0.5,
        size="FWHM",
        radius=10.0,
        sbuffer=1.0,
        swidth=3.0,
        iterations=1,
        logfile=log,
        imagecur=coords,
        Stdout=1,
    )

    fwhm = data.pop().split().pop()

    fwhm_pix = float(fwhm)
    fwhm = fwhm_pix * scale
    seeing = dimm_seeing(fwhm_pix, ref, scale)
    print seeing
    fp = open(out, "w")
    fp.write("%f %f\n" % (fwhm_pix, seeing))
    fp.close()
    os.system('echo "image;text 85 500 # text={Spot FWHM = %5.2f pixels}" | xpaset WFS regions' % fwhm_pix)
    os.system("echo 'image;text 460 500 # text={Seeing = %4.2f\"}' | xpaset WFS regions" % seeing)
    os.system('echo "set wfs_seeing %4.2f" | nc hacksaw 7666' % seeing)

def dark_sky_flat(filter):
    med_otalist = []
    print 'making dark sky flats for',filter
    for key in tqdm(odi.OTA_dictionary):
        image_list = odi.OTA_dictionary[key]+'.'+filter+'.lis'
        med_out = image_list.replace('.lis','.med.fits')
        med_otalist.append(med_out)
        iraf.unlearn(iraf.immatch.imcombine)
        iraf.immatch.imcombine.setParam('input','@'+str(image_list))
        iraf.immatch.imcombine.setParam('output',odi.skyflatpath+med_out)
        iraf.immatch.imcombine.setParam('combine','median')
        iraf.immatch.imcombine.setParam('masktype','goodvalue')
        iraf.immatch.imcombine.setParam('maskvalue',0)
        iraf.immatch.imcombine.setParam('scale','median')
        # iraf.immatch.imcombine.setParam('zero','none')
        iraf.immatch.imcombine.setParam('zero','median')
        iraf.immatch.imcombine(logfile='imcombine.log.txt', mode='h')
        if key == 1:
            data,header = odi.fits.getdata(odi.skyflatpath+med_out,header=True)
            mean, median, std = odi.sigma_clipped_stats(data, sigma=3.0)
            normalization_factor = median
    iraf.set(clobber = 'yes')
    print 'smoothing dark sky flats for',filter
    for i in tqdm(range(len(med_otalist))):
        iraf.unlearn(iraf.imutil.imarith,iraf.imfilter.median)
        iraf.imutil.imarith.setParam('operand1',odi.skyflatpath+med_otalist[i])
        iraf.imutil.imarith.setParam('op','/')
        iraf.imutil.imarith.setParam('operand2',normalization_factor)
        iraf.imutil.imarith.setParam('result',odi.skyflatpath+med_otalist[i])
        iraf.imutil.imarith(verbose='no', mode='h')
        iraf.imfilter.median.setParam('input',odi.skyflatpath+med_otalist[i])
        iraf.imfilter.median.setParam('output',odi.skyflatpath+med_otalist[i])
        iraf.imfilter.median.setParam('xwindow',3)
        iraf.imfilter.median.setParam('ywindow',3)
        iraf.imfilter.median(verbose='no', mode='h')
    iraf.set(clobber = 'no')

    return normalization_factor

def get_star_data(asteroid_id, mag, expnum, header):
    """
    From ossos psf fitted image, calculate mean of the flux of each row of the rotated PSF
    """

    # calculate mean psf
    uri = storage.get_uri(expnum.strip('p'), header[_CCD].split('d')[1])
    ossos_psf = '{}.psf.fits'.format(uri.strip('.fits'))
    local_psf = '{}{}.psf.fits'.format(expnum, header[_CCD].split('d')[1])
    local_file_path = '{}/{}'.format(_STAMPS_DIR, local_psf)
    storage.copy(ossos_psf, local_file_path)

    # pvwcs = wcs.WCS(header)
    # x, y = pvwcs.sky2xy(asteroid_id['ra'].values, asteroid_id['dec'].values)
    x = asteroid_id[_XMID_HEADER].values[0]
    y = asteroid_id[_YMID_HEADER].values[0]

    # run seepsf on the mean psf image
    iraf.set(uparm="./")
    iraf.digiphot(_doprint=0)
    iraf.apphot(_doprint=0)
    iraf.daophot(_doprint=0)
    iraf.seepsf(local_file_path, local_psf, xpsf=x, ypsf=y, magnitude=mag)

    with fits.open(local_psf) as hdulist:
        data = hdulist[0].data

    th = math.degrees(asteroid_id[_THETA_HEADER].values[0])
    data_rot = rotate(data, th)
    data_rot = np.ma.masked_where(data_rot == 0, data_rot)

    data_mean = np.ma.mean(data_rot, axis=1)

    os.unlink(local_psf)
    os.unlink(local_file_path)

    return data_mean[np.nonzero(np.ma.fix_invalid(data_mean, fill_value=0))[0]]

def mask_create(images, no_unlearn=False):
    output_files = []
    # If the output file list doesn't exsist, create it.
    if len(output_files) == 0:
        output_files = [os.path.join(os.path.dirname(i), "mask_" + os.path.basename(i)) for i in images]
    # Create a temporary file to hold the input image list
    tfile = tempfile.mkstemp()
    os.write(tfile[0], "\n".join(images) + "\n")
    os.close(tfile[0])

    # Create a temporary file to hold the input image list
    tofile = tempfile.mkstemp()
    os.write(tfile[0], "\n".join(output_files) + "\n")
    os.close(tfile[0])

    try:

        # Set up to call the task
        iraf.set(clobber="yes")
        iraf.nproto()

        if no_unlearn == False:
            print "Objmasks is unlearning!!"
            iraf.objmasks.unlearn()

        # set parameters
        iraf.objmasks.images = "@" + tfile[1]
        iraf.objmasks.objmasks = "@" + tofile[1]
        iraf.objmasks.omtype = "boolean"
        iraf.objmasks()

    except iraf.IrafError, e:
        print "objmasks failed"
        print "error #" + str(e.errno)
        print "Msg: " + e.errmsg
        print "Task: " + e.errtask
        raise e

def dark_sky_flat():
  med_otalist = []
  for key in OTA_dictionary:
    image_list = OTA_dictionary[key]+'.lis'
    med_out = image_list.replace('.lis','.med.fits')
    med_otalist.append(med_out)
    iraf.unlearn(iraf.immatch.imcombine)
    iraf.immatch.imcombine.setParam('input','@'+str(image_list))
    iraf.immatch.imcombine.setParam('output',med_out)
    iraf.immatch.imcombine.setParam('combine','median')
    iraf.immatch.imcombine.setParam('masktype','goodvalue')
    iraf.immatch.imcombine.setParam('maskvalue',0)
    iraf.immatch.imcombine.setParam('scale','median')
    iraf.immatch.imcombine.setParam('zero','median')
    iraf.immatch.imcombine(mode='h')
    if key == 1:
      data,header = fits.getdata(med_out,header=True)
      mean, median, std = sigma_clipped_stats(data, sigma=3.0)
      normalization_factor = median
  iraf.set(clobber = 'yes')
  for i in range(len(med_otalist)):
    iraf.unlearn(iraf.imutil.imarith,iraf.imfilter.median)
    iraf.imutil.imarith.setParam('operand1',med_otalist[i])
    iraf.imutil.imarith.setParam('op','/')
    iraf.imutil.imarith.setParam('operand2',normalization_factor)
    iraf.imutil.imarith.setParam('result',med_otalist[i])
    iraf.imutil.imarith(mode='h')
    iraf.imfilter.median.setParam('input',med_otalist[i])
    iraf.imfilter.median.setParam('output',med_otalist[i])
    iraf.imfilter.median.setParam('xwindow',3)
    iraf.imfilter.median.setParam('ywindow',3)
    iraf.imfilter.median()
  iraf.set(clobber = 'no')
  
  
  return normalization_factor

#   within ds9. They show the background image with masked regions superposed. 
#   These are useful for testing whether stellar outskirts are sufficiently masked

#Author: Ryan Muther, Union College, June-July 2013
#This program was written as a replacement for the IDL program of the same 
#name, originally written by Cullen Blake 1999.
# Updated May 8, 2014 by Rebecca Koopmann:
#   - added documentation

from pyraf import iraf
from astropy.io import fits
import numpy as np
import sys

# set stdimage to 2048
iraf.set(stdimage='imt2048')

# Load the bad pixel mask and reverse its values so the mask exists in a sensible manner.
# Currently, the zeroes represent good pixels and the ones represent masked pixels.  
# This is the IRAF default, but less useful than the reverse when manipulating the masks
# and images. This reversal process allows for the creation of a masked image by simple 
# elementwise multiplication of the  mask and the image to be masked

def reverseMask(filename):
  mask = fits.getdata(filename+"mask.fits")
  
  manyZeros = np.zeros_like(mask)
  manyOnes = np.ones_like(mask)
  
  maskReversed = np.where((mask!=0),manyZeros,manyOnes)
  

    def _setup_iraf(self, datamin=0, datamax=60000, epadu=1., fitrad_fwhm=1.,
                    fitsky='no', function='moffat25', fwhmpsf=3., itime=10.,
                    maxiter=50, maxnpsf=60, mergerad_fwhm=2., nclean=10,
                    psfrad_fwhm=10., ratio=1., readnoi=0, recenter='yes',
                    roundlo=-1.0, roundhi=1.0,
                    saturated='no', sharplo=0.2, sharphi=1.0,
                    sigma=5., theta=0., threshold=3., varorder=1, zmag=20.,
                    sannulus_fwhm=2., wsannulus_fwhm=2.,
                    annulus_fwhm=10., dannulus_fwhm=2.):
        """Sets the IRAF/DAOPHOT configuration parameters.

        Parameters
        ----------
        fitrad_fwhm : float
            The PSF fitting radius in units `fwhmpsf`. Only pixels within the
            fitting radius of the center of a star will contribute to the
            fits computed by the `allstar` task. For most images the fitting
            radius should be approximately equal to the FWHM of the PSF
            (i.e. `fitrad_fwhm = 1`). Under severely crowded conditions a
            somewhat smaller value may be used in order to improve the fit.
            If the PSF is variable, the FWHM is very small, or sky fitting
            is enabled on the other hand, it may be necessary to increase the
            fitting radius to achieve a good fit. (default: 1.)

        fitsky : str (optional)
            Compute new sky values for the stars in the input list (allstar)?
            If fitsky = "no", the `allstar` task compute a group sky value by
            averaging the sky values of the stars in the group. 
            If fitsky = "yes", the `allstar` task computes new sky values for
            each star every third iteration by subtracting off the best
            current fit for the star and and estimating the median of the
            pixels in the annulus defined by sannulus and wsannulus.
            The new group sky value is the average of the new individual
            values. (default: 'no')

        function : str (optional)
            PSF model function. One of "auto", "gauss", "moffat15", "moffat25",
            "lorentz", "penny1", or "penny2". (default: 'moffat25')

        maxiter : int (optional)
            The maximum number of times that the `allstar` task will iterate
            on the PSF fit before giving up. (default: 50)

        maxnpsf : int (optional)
            The maximum number of candidate psf stars to be selected. (default: 60)

        mergerad_fwhm : float (optional)
            The critical separation in units `psffwhm` between two objects
            for an object merger to be considered by the `allstar` task.
            Objects with separations > mergerad will not be merged; faint
            objects with separations <= mergerad will be considered for merging.
            The default value of mergerad is sqrt (2 *(PAR1**2 + PAR2**2)),
            where PAR1 and PAR2 are the half-width at half-maximum along the
            major and minor axes of the psf model. Merging can be turned off
            altogether by setting mergerad to 0.0. (default: 2.0)

        nclean: int
            The number of additional iterations the `psf` task performs to
            compute the PSF look-up tables. If nclean is > 0, stars which
            contribute deviant residuals to the PSF look-up tables in the
            first iteration, will be down-weighted in succeeding iterations. 
            The DAOPHOT manual recommends 5 passes. (default: 10)

        psfrad_fwhm : float
            The radius of the circle in units `fwhmpsf` within which the PSF
            model is defined (i.e. the PSF model radius in pixels is obtained
            by multiplying `psfrad_fwhm` with `fwhmpsf`). 
            Psfrad_fwhm should be slightly larger than the radius at
            which the intensity of the brightest star of interest fades into
            the noise. Large values are computationally expensive, however.
            Must be larger than `fitrad_fwhm` in any case. (default: 10.)

        ratio : float (optional)
            Ratio of minor to major axis of Gaussian kernel for object
            detection (default: 1.0)

        recenter : str (optional)
            One of 'yes' or 'no'.

        roundlo : float (optional)
            Lower bound on roundness for object detection (default: -1.0)

        roundhi : float (optional)
            Upper bound on roundness for object detection (default: 1.0)

        sannulus_fwhm : float (optional)
            The inner radius of the sky annulus in units `fwhmpsf` used 
            by `allstar` to recompute the sky values. (default: 2.)

        saturated : str (optional)
            Use saturated stars to improve the signal-to-noise in the wings
            of the PSF model computed by the PSF task? This parameter should
            only be set to "yes" where there are too few high signal-to-noise
            unsaturated stars in the image to compute a reasonable model
            for the stellar profile wings. (default: "no")

        sharplo : float (optional)
            Lower bound on sharpness for object detection (default: 0.2)

        sharphi : float (optional)
#.........这里部分代码省略.........

# Importing pyraf -- note that this requires IRAF to be installed
# and the "iraf" and "IRAFARCH" environment variables to be set
with warnings.catch_warnings():
    # Ignore the ImportWarning we often get when importing pyraf, because pyraf
    # likes to create working dirs with its name in the local dir.
    warnings.filterwarnings("ignore", message="(.*)Not importing dir(.*)")
    # Disable pyraf graphics:
    from stsci.tools import capable
    capable.OF_GRAPHICS = False
    from pyraf import iraf
    from iraf import noao, digiphot, daophot
    # Turn cache off to enable multiprocessing to be used;
    # cf. http://www.stsci.edu/institute/software_hardware/pyraf/pyraf_faq#3.5
    iraf.prcacheOff()
    iraf.set(writepars=0)

from .utils import timed


class DaophotError(Exception):
    """Raised if Daophot failed to perform an operation."""
    pass


class Daophot(object):
    """DAOPHOT wrapper class.

    The constructor of this class expects to receive
    ALL the non-default configuration parameters because IRAF is very stateful.

def doit(fn,dp,mt):
  """ JHT: this is the module that analyses an ALFOSC image of
           a Vertical Slit, and recommends an alignment offset to the
           current APERTURE WHEEL stepper motor units."""
          

  ## First check if data file exists
  if not os.access(dp+fn+".fits",os.F_OK): 
     messageOut(mt,"File not found:  "+dp+fn+".fits\n")
     return "File not found:  "+dp+fn+".fits"

  messageOut(mt,"\nVertical-Slit analysis of file:  "+dp+fn+".fits\n")

  from pyraf import iraf
  from pyraf import gwm


  ## Read current grism wheel position from image FITS header
  iraf.images.imutil.imgets(dp+fn,"GRISM")
  grismid=iraf.images.imutil.imgets.value
  if grismid.find("pen")==-1: 
     messageOut(mt,"Vertical-Slit mode:  grism wheel should be Open\n")
     return "File  %s:  Vertical-Slit mode:  grism wheel should be Open" % fn




  ## Read current aperture wheel position from image FITS header
  iraf.images.imutil.imgets(dp+fn,"APERTUR")
  slitid=iraf.images.imutil.imgets.value
  if slitid.find("Vert")==-1: 
     messageOut(mt,"Vertical-Slit mode:  no vertical slit in aperture wheel\n")
     return "File  %s:  Vertical-Slit mode:  no vertical slit in aperture wheel" % fn

  slitwidth=float(slitid.split('_',9)[1].strip('\"'))
  messageOut(mt,"Slit width "+str(slitwidth)+"   "+slitid+"\n")

  iraf.noao(_doprint=0)
  iraf.noao.imred(_doprint=0)
  iraf.noao.imred.specred(_doprint=0)

  if not os.access("/tmp/tiasgat/",os.F_OK):
    os.mkdir("/tmp/tiasgat/")
    os.chmod("/tmp/tiasgat/",0777)
  if os.access("/tmp/tiasgat/plot",os.F_OK):     os.remove("/tmp/tiasgat/plot")
  if os.access("/tmp/tiasgat/plot2",os.F_OK):    os.remove("/tmp/tiasgat/plot2")
  if os.access("/tmp/tiasgat/aplast",os.F_OK):   os.remove("/tmp/tiasgat/aplast")
  if os.access("/tmp/tiasgat/tvmarks",os.F_OK):  os.remove("/tmp/tiasgat/tvmarks")
  if os.access("/tmp/tiasgat/logfile",os.F_OK):  os.remove("/tmp/tiasgat/logfile")

  ## Note that this will *not* update any uparm files !! (see pyraf documentation)
  iraf.noao.imred.specred.dispaxis=2
  iraf.noao.imred.specred.database="/tmp/tiasgat/"
  iraf.noao.imred.specred.plotfile="/tmp/tiasgat/plot"
  iraf.noao.imred.specred.logfile="/tmp/tiasgat/logfile"
  width=6*slitwidth    # x6 or else the widest slits go wrong
  if width<4: width=4  # smaller than 4 and apfind will not find it
  messageOut(mt,"Using width of "+str(width)+" pixels \n")
  iraf.noao.imred.specred.apedit.width=width



  ## Display image on ds9
  iraf.set(stdimage="imt512")
  iraf.display(dp+fn,1,fill="no",Stdout="/dev/null")

# Suppress IRAF query for number of apertures to find
# This is only necesary for the widest slits: then the call to
# apfind results in an empty database file, as it cannot find an aperture.
# But aptrace works fine anyway (except for the annoying query) !? 
  iraf.noao.imred.specred.apfind.setParam('nfind.p_value', 1)
  iraf.noao.imred.specred.apfind.setParam('nfind.p_mode','h')

  ## 'find' and trace spectrum; this will dump a plot to /tmp/tiasgat/plot
  lines = iraf.noao.imred.specred.apfind(dp+fn,nfind=1,interactive="no", edit="no", nsum=20, Stdout=1)
  for i in range (0,len(lines)):     messageOut(mt,lines[i]+"\n")


  lines = iraf.noao.imred.specred.aptrace(dp+fn,interactive="no",step=5,low_reject=2.5,
                   high_reject=2.5,function="leg",order=2,niterate=5,naverage=1, Stdout=1)
  for i in range (0,len(lines)):     messageOut(mt,lines[i]+"\n")


  ## Start graphics window; select the correct plot; show plot
  gwm.window("Tiasgat!  graphics")
  iraf.plot.gkiextract("/tmp/tiasgat/plot",2,Stdout="/tmp/tiasgat/plot2")
  gwm.getActiveGraphicsWindow().load("/tmp/tiasgat/plot2")


### how to read the aperture file, as output by aptrace ####
### 
### center    line 6        gives zero point
### max,min   lines 24-25   n = (2 * x - (max + min)) / (max - min)
### c1,c2     lines 26-27   
### 
### The polynomial can be expressed as the sum
###     
###             poly = sum from i=1 to order {c_i * z_i} 
###     
### where  the  the  c_i  are  the  coefficients and the z_i are defined
#.........这里部分代码省略.........