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 wal(lstfile):
    iraf.noao()
    iraf.twodspec()
    iraf.longslit()
    iraf.identify(images = 'Lamp'
        , section = 'middle column', database = 'database'
        , coordlist = 'linelists$idhenear.dat', units = '', nsum = 10
        , match = -3.0, maxfeatures = 50, zwidth = 100.0
        , ftype = 'emission', fwidth = 20.0, cradius = 5.0
        , threshold = 0.0, minsep = 2.0, function = 'chebyshev'
        , order = 6, sample = '*', niterate = 0
        , low_reject = 3.0, high_reject = 3.0, grow = 0.0
        , autowrite = False, graphics = 'stdgraph', cursor = ''
        , crval = '', cdelt = '')
    iraf.reidentify(reference = 'Lamp'
        , images = 'Lamp', interactive = 'no', section = 'column'
        , newaps = True, override = True, refit = True, trace = False
        , step = 10, nsum = 10, shift = 0.0, search = 0.0, nlost = 5
        , cradius = 7.0, threshold = 0.0, addfeatures = False
        , coordlist = 'linelists$idhenear.dat', match = -3.0, maxfeatures = 50
        , minsep = 2.0, database = 'database', logfiles = 'logfile'
        , plotfile = '', verbose = False, graphics = 'stdgraph', cursor = ''
        , answer = 'yes', crval = '', cdelt = '', mode = 'al')
    iraf.fitcoords(images = 'Lamp'
        , fitname = 'Lamp', interactive = True, combine = False, database = 'database'
        , deletions = 'deletions.db', function = 'chebyshev', xorder = 6
        , yorder = 6, logfiles = 'STDOUT,logfile', plotfile = 'plotfile'
        , graphics = 'stdgraph', cursor = '', mode = 'al')
    iraf.longslit(dispaxis = 2)
    iraf.transform(input = '%ftbo%ftbo%@' + lstfile
        , output = '%wftbo%wftbo%@' + lstfile, minput = '', moutput = ''
        , fitnames = 'LampLamp', database = 'database', interptype = 'spline3'
        , flux = True)

def thar_cal(object_b_fn_ec, object_b_fn_ec_w, colour):
# Import IRAF modules:
  iraf.noao(_doprint=0)
  iraf.onedspec(_doprint=0)
# Check input file and reference extraction exist before proceeding:
  if os.path.isfile(object_b_fn_ec) == True:
# Perform dispersion correction:
    iraf.dispcor(input=object_b_fn_ec, output=object_b_fn_ec_w)
    print ' '
    print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    print colour.capitalize() + ' object spectrum             '
    print '(' + str(object_b_fn_ec) + ')'
    print 'successfully wavelength calibrated in file         '
    print str(object_b_fn_ec_w) + '.'
    print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    print ' '
  else:
    print ' '
    print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    print 'Input frame                                        ' 
    print str(object_b_fn_ec)
    print 'does not exist. Exiting script.                    '
    print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    print ' '
    print ' '
    sys.exit()

    def StandardTask(self, InputFile, OutputFile, FitsFolder, airmass_value, exptime_value):
   
        iraf.noao(_doprint=0)
        iraf.onedspec(_doprint=0)               
        
        #From the fits file determine which is the star being treated
        StarName = self.StarnameFromFileName(InputFile)
        
        #Get the corresponding Calibration file
        CalibrationFile, Bandwidth, Bandsep   = self.getCalibrationFile(StarName)
        
        #Incase no output name is given, we generate one with the provided "preffix" (The defaul format is a_std_wolf.dat)
        if OutputFile == None:
            OutputFile    = 'a_std_' + StarName
        
        #Prepare dictionary with the configuration for the tasks           
        Standard_conf_Comb  = self.StandardAttributes(InputFile, OutputFile, CalibrationFile, FitsFolder, airmass_value, exptime_value, Bandwidth, Bandsep)

        #Display the equivalent command in IRAF
        Command = self.printIrafCommand('standard', Standard_conf_Comb)
        print '--- Using the command'
        print Command
        
        #Run the task        
        iraf.onedspec.standard(**Standard_conf_Comb)
        
        return OutputFile 

def coroverbiastrim(lstfile):
    iraf.noao()
    iraf.imred()
    iraf.ccdred()
    x1,x2,y1,y2 = get_trim_sec()
    iraf.ccdproc(images = '@' + lstfile + '//[1]'
        , output = '%bo%bo%@' + lstfile
        , ccdtype = '', max_cache = 0, noproc = False
        , fixpix = False, overscan = True, trim = False
        , zerocor = True, darkcor = False, flatcor = False
        , illumcor = False, fringecor = False, readcor = False
        , scancor = False, readaxis = 'line', fixfile = ''
        , biassec = '[5:45,%s:%s]'%(y1,y2), trimsec = '[%s:%s,%s:%s]'%(x1,x2,y1,y2)
        , zero = 'Zero', dark = '', flat = '', illum = '', fringe = ''
        , minreplace = 1.0, scantype = 'shortscan', nscan = 1
        , interactive = False, function = 'chebyshev', order = 1
        , sample = '*', naverage = 1, niterate = 1
        , low_reject = 3.0, high_reject = 3.0, grow = 1.0)
    iraf.ccdproc(images = '%bo%bo%@' + lstfile
        , output = '%tbo%tbo%@' + lstfile
        , ccdtype = '', max_cache = 0, noproc = False
        , fixpix = False, overscan = False, trim = True
        , zerocor = False, darkcor = False, flatcor = False
        , illumcor = False, fringecor = False, readcor = False
        , scancor = False, readaxis = 'line', fixfile = ''
        , biassec = '[5:45,%s:%s]'%(y1,y2), trimsec = '[%s:%s,%s:%s]'%(x1,x2,y1,y2)
        , zero = 'Zero', dark = '', flat = '', illum = '', fringe = ''
        , minreplace = 1.0, scantype = 'shortscan', nscan = 1
        , interactive = False, function = 'chebyshev', order = 1
        , sample = '*', naverage = 1, niterate = 1
        , low_reject = 3.0, high_reject = 3.0, grow = 1.0)
    iraf.flpr()

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 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 _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 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 standard():
    stdpath = os.path.split(os.path.realpath(__file__))[0] + os.sep + 'standarddir' + os.sep
    print('standard dir is ' + stdpath)
    extpath = os.path.split(os.path.realpath(__file__))[0] + os.sep + 'LJextinct.dat'
    iraf.noao()
    iraf.twodspec()
    iraf.longslit(dispaxis = 2, nsum = 1, observatory = 'Lijiang', 
            extinction = extpath, caldir = stdpath)
    for objname in stdgroup:
        stdname, stdmag, stdmagband = get_std_name(objname)
        print('the standard star is ' + stdname)
        stdmag = float(stdmag)
        outname1 = 'stdawftbo' + stdgroup[objname][0]
        inname   = ''
        for tmpname in stdgroup[objname]:
            inname = inname + 'awftbo' + tmpname + ','
        inname = inname[0:-1]
        iraf.standard(input = inname
                , output = outname1, samestar = True, beam_switch = False
                , apertures = '', bandwidth = 30.0, bandsep = 20.0
                , fnuzero = 3.6800000000000E-20, extinction = extpath
                , caldir = stdpath, observatory = ')_.observatory'
                , interact = True, graphics = 'stdgraph', cursor = ''
                , star_name = stdname, airmass = '', exptime = ''
                , mag = stdmag, magband = stdmagband, teff = '', answer = 'yes')
    for name in stdgroup:
        inpar = 'stdawftbo' + stdgroup[name][0]
        iraf.sensfunc(standards = inpar, sensitivity = 'sensawftbo' + stdgroup[name][0], 
            extinction = extpath, function = 'spline3', order = 9)

def wspectext_allap(filepath, apaxis=0):
    """ this can be directly evaluated in terminal

    Parameters
    ----------
    filepath: string
        filepath

    apaxis: integer
        0 or 1.
        0 denotes dispersion along rows (this is the case for Xinglong216HRS)
        1 denotes dispersion along columns

    Examples
    --------
    wspectext_allap('./w20160120001s.fits')

    """
    nrow, ncol = fits.open(filepath)[0].data.shape
    print '@Cham: nrow = %d, ncol=%d' % (nrow, ncol)
    
    iraf.noao()
    iraf.noao.onedspec()
    
    filename = filepath.split(os.path.sep)[-1]
    
    # determine dirname & dirpath
    dirname = filepath.split(os.path.sep)[-1].replace('.fits', '')
    if os.path.dirname(filepath) == '':
        dirpath = dirname
    else:
        dirpath = os.path.dirname(filepath) + os.path.sep + dirname

    if os.path.exists(dirpath):
        # if dirpath exists
        print '@Cham: directory exists ... (%s)' % dirpath
    else:
        # if dirpath doesn't exist, mkdir
        os.mkdir(dirpath)
        print '@Cham: mkdir %s'%dirpath
    
    # execute wspectest
    if not apaxis == 1:
        for apnum in xrange(1, nrow+1):
            _input = '%s[1:%d, %d]' % (filename, ncol, apnum)
            _output = dirpath + os.sep + filename.replace('.fits', '_%04d.dat' % apnum)
            print '@Cham: wspectext running ... [%s]' % _output
            iraf.onedspec.wspectext(input=_input, output=_output, header='no')
    else:
        for apnum in xrange(1, ncol+1):
            _input = '%s[%d, 1:%d]' % (filename, apnum, nrow)
            _output = dirpath + os.sep + filename.replace('.fits', '_%04d.dat' % apnum)
            print '@Cham: wspectext running ... [%s]' % _output
            iraf.onedspec.wspectext(input=_input, output=_output, header='no')

    # print masseges
    print '@Cham: ----------------------'
    print '@Cham: mission complete!'
    print '@Cham: ----------------------'

	def makeObject(self, inObjectFile="mkobject.coo"):
		
		"""IRAF routine: mkobjects. This places a star in the fits image that is the parent object of this class, at the given ra and dec. It takes as input a file of "ra dec magnitude" for a star type object.

		1. Get info from the header files if you want noise
		2. Get FWHM for size of the star for the given image
		3. Generate the star

		radius = FWHM / scale, from statistics, scale = 2*sqrt(2*ln(2))
		"""
	
		# 1.
		self.loadHeader()
		INFODICT = {
			"RON": self.getHeader("RON"), \
			"GAIN": self.getHeader("GAIN"), \
			"EXPTIME": self.getHeader("EXPTIME"), \
			"MAGZP": self._ZPDICT[self.getHeader("FILTER")], \
			"POISSON": "no", \
			"STAR": "gaussian", \
			}	
		# 2.
		scale = 2*scipy.sqrt(2*scipy.log(2))
		self.getMyMedianFWHM()
		INFODICT["RADIUS"] = self._MEDFWHM / scale

		# 3. 
		outputname = self._Name.replace(".fits", "_mko.fits")
		iraf.noao(_doprint=0)
		iraf.artdata(_doprint=0)
		
	        iraf.mkobjects(input=self._Name, \
                        output=outputname, \
                        # When creating new images
                        #title="",\
                        #ncols=self.xdim, \
                        #nlines=self.ydim, \
                        #header=self.hdrfile, \
                        #background=0.0, \
                        # When creating objects
                        objects=inObjectFile, \
                        xoffset=0.0, \
                        yoffset=0.0, \
                        star=INFODICT["STAR"], \
                        radius=INFODICT["RADIUS"], # this is fwhm/scale (pixels), where scale = 2sqrt(2ln2)\
                        #beta=2.5,  # for star=moffat profile\
                        ar=1.0, \
                        pa=0.0, \
                        distance=1.0, \
                        exptime=INFODICT["EXPTIME"], \
                        magzero=INFODICT["MAGZP"], \
                        # Noise parameters
                        gain=INFODICT["GAIN"], \
                        rdnoise=INFODICT["RON"], \
                        poisson=INFODICT["POISSON"], \
                        seed=1, \
                        comments=1)

		return outputname 

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 combinebias(filename):
    iraf.noao()
    iraf.imred()
    iraf.ccdred()
    iraf.zerocombine(input = 'o//@' + filename
	, output = 'Zero', combine = 'average', reject = 'minmax'
	, ccdtype = '', process = False, delete = False
	, clobber = False, scale = 'none', statsec = ''
	, nlow = 0, nhigh = 1, nkeep = 1, mclip = True
	, lsigma = 3.0, hsigma = 3.0, rdnoise = 'rdnoise'
	, gain = 'gain', snoise = 0.0, pclip = -0.5, blank = 0.0)

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 get_masked_image(img):

    # Running IRAF
    iraf.noao(_doprint=0)  
    mask = make_mask_cross(img)    
    masked = img.replace(".fits", "_masked.fits")
    
    if (os.path.isfile(masked)):
        os.remove(masked)
    iraf.imarith(img, "*", mask, masked)
    
    return masked

def main():
	iraf.noao()
	iraf.imred()
	iraf.ccdred()
        print '=' * 20, 'Overscan', '=' * 20
	correct_overscan('spec.lst')
	print '=' * 20, 'combine bias', '=' * 20
	combine_bias('bias.lst')
	print '=' * 20, 'correct bias', '=' * 20
	correct_bias('spec_no_bias.lst', 'bias_spec.fits')
        name = os.popen('ls object*.lst').readlines()
        name = [i.split()[0] for i in name]
        for i in name:
            ntrim_flat(i)

def thar_ref(object_b_fn_ec, wave_ref, colour):
# Import IRAF modules:
  iraf.noao(_doprint=0)
  iraf.onedspec(_doprint=0)
#  iraf.echelle(_doprint=0)
  parList = "wavelength_calibration_refspectra.par"
  if os.path.isfile(object_b_fn_ec) == True:
    if os.path.isfile(wave_ref) == True:
      if os.path.isfile(parList) == True:
# Assign reference spectra to input object spectra using IRAF task refspectra:
        iraf.refspectra.setParList(ParList="wavelength_calibration_refspectra.par")
        iraf.refspectra(input=object_b_fn_ec, references=wave_ref)
        print ' '
        print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        print colour.capitalize() + ' object spectrum             '
        print '(' + str(object_b_fn_ec) + ')                      '
        print 'assigned to ' + str(colour) + ' reference spectrum '
        print '(' + str(wave_ref) + ') successfully.              ' 
        print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        print ' '
      else:
        print ' '
        print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        print 'Wavelength reference IRAF .par file                ' 
        print str(parList)
        print 'does not exist. Exiting script.                    '
        print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        print ' '
        print ' '
        sys.exit()
    else:
      print ' '
      print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      print 'ThAr reference frame                               ' 
      print str(wave_ref)
      print 'does not exist. Exiting script.                    '
      print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      print ' '
      print ' '
      sys.exit()
  else:
    print ' '
    print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    print 'Input frame                                        ' 
    print str(object_b_fn_ec)
    print 'does not exist. Exiting script.                    '
    print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    print ' '
    print ' '
    sys.exit()

def gcDAR():
    iraf.noao()
    obs = iraf.noao.observatory

    # Set up Keck observatory info
    obs(command="set", obsid="keck")
    keck = ephem.Observer()
    keck.long = math.radians(-obs.longitude)
    keck.lat = math.radians(obs.latitude)
    keck.elev = obs.altitude
    keck.pressure = 800.0
    keck.temp = 278.0
    
    # Date of Observations: 06maylgs1
    keck.date = '2006/05/03 %d:%d:%f' % (obs.timezone+3, 28, 46.33)

    # Set up the galactic center target
    sgra = ephem.FixedBody()
    sgra._ra = ephem.hours("17:45:40")
    sgra._dec = ephem.degrees("-29:00:10")
    sgra._epoch = 2000
    sgra.compute()

    pos2 = ephem.FixedBody()
    pos2._ra = ephem.hours("17:45:40")
    pos2._dec = ephem.degrees("-29:00:00")
    pos2._epoch = 2000
    pos2.compute()

    print 'Date of Obs: ', keck.date
    sgra.compute(keck)
    pos2.compute(keck)

    print 'Azimuth of Objects:  %s  vs.  %s' % (sgra.az, pos2.az)

    for ii in range(15):
	keck.lat = math.radians(obs.latitude - (ii*2))

	sgra.compute(keck)
	pos2.compute(keck)

	angAbs = ephem.separation((sgra.ra, sgra.dec), (pos2.ra, pos2.dec))
	angAbs *= 206265.0
	angRel = ephem.separation((sgra.az, sgra.alt), (pos2.az, pos2.alt))
	angRel *= 206265.0
	angDiff = angAbs - angRel

	print 'Sgr A*:  %s   vs.  %s  deltaR = %5d (muas)' % \
	    (sgra.alt, pos2.alt, angDiff*1e6)

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)