def test_whereis(tmpdir):
    iraf.plot(_doprint=0)
    iraf.images(_doprint=0)
    outfile = str(tmpdir.join('output.txt'))

    cases = ("pw", "lang", "st", "std", "stdg", "stdpl", "star", "star man",
             "vi", "noao", "impl", "ls", "surf", "surf man", "smart man",
             "img", "prot", "pro", "prow", "prows", "prowss", "dis", "no")

    # Test the whereis function
    for arg in cases:
        args = arg.split(" ")  # convert to a list
        iraf.printf("--> whereis " + arg + '\n', StdoutAppend=outfile)
        kw = {'StderrAppend': outfile}
        iraf.whereis(*args, **kw)  # catches stdout+err

    with open(outfile) as f:
        for line in f:
            if 'task not found' in line:
                continue

            row = line.split()
            if line.startswith('-->'):
                cmd = row[2]
                if len(row) > 3:
                    cmd2 = row[3]
            elif row[0] == 'user.man' and cmd2 == 'man':
                pass
            else:
                assert all([s.split('.')[1].startswith(cmd) for s in row]), \
                    '{}'.format(row)

def display_image(img,frame,_z1,_z2,scale,_xcen=0.5,_ycen=0.5,_xsize=1,_ysize=1,_erase='yes'):
    goon='True'
    import glob, subprocess, os, time
    ds9 = subprocess.Popen("ps -U {:d} u | grep -v grep | grep ds9".format(os.getuid()),shell=True,stdout=subprocess.PIPE).stdout.readlines()
    if len(ds9)== 0 :   
       subproc = subprocess.Popen('ds9',shell=True)
       time.sleep(3)

    if glob.glob(img):
       from pyraf import iraf
       iraf.images(_doprint=0)
       iraf.tv(_doprint=0)
       import string,os
       if _z2: 
          try:
              sss=iraf.display(img + '[0]', frame, xcen=_xcen, ycen=_ycen, xsize=_xsize, ysize=_ysize, erase=_erase,\
                                   fill='yes', zscale='no', zrange='no', z1=_z1, z2=_z2,Stdout=1)
          except:
              print ''
              print '### ERROR: PROBLEM OPENING DS9'
              print ''
              goon='False'                 
       else:
        try:  
            sss=iraf.display(img + '[0]', frame, xcen=_xcen, ycen=_ycen, xsize=_xsize, ysize=_ysize, erase=_erase, fill='yes', Stdout=1)
        except:
            print ''
            print '### ERROR: PROBLEM OPENING DS9'
            print ''
            goon=False
 
       if scale and goon:
          answ0 = raw_input('>>> Cuts OK ? [y/n] ? [y] ')
          if not answ0: answ0='y'
          elif answ0=='no' or answ0=='NO': answ0='n' 

          while answ0=='n':
              _z11=float(string.split(string.split(sss[0])[0],'=')[1])
              _z22=float(string.split(string.split(sss[0])[1],'=')[1])
              z11 = raw_input('>>> z1 = ? ['+str(_z11)+'] ? ')
              z22 = raw_input('>>> z2 = ? ['+str(_z22)+'] ? ')
              if not z11: z11=_z11
              else: z11=float(z11)
              if not z22: z22=_z22
              else: z22=float(z22)
              print z11,z22
              sss=iraf.display(img + '[0]',frame,fill='yes', xcen=_xcen, ycen=_ycen, xsize=_xsize, ysize=_ysize, erase=_erase,\
                                   zrange='no', zscale='no', z1=z11, z2=z22, Stdout=1)
              answ0 = raw_input('>>> Cuts OK ? [y/n] ? [y] ')
              if not answ0: answ0='y'
              elif answ0=='no' or answ0=='NO': answ0='n'
       if goon:
          _z1,_z2=string.split(string.split(sss[0])[0],'=')[1],string.split(string.split(sss[0])[1],'=')[1]
    else:
        print 'Warning: image '+str(img)+' not found in the directory '
    return _z1,_z2,goon

def imarith(operand1, op, operand2, result):
    from pyraf import iraf
    iraf.images()

    pars = iraf.imarith.getParList()
    iraf.imcombine.unlearn()

    print "%s %s %s -> %s" % (operand1, op, operand2, result)
    iraf.imarith(operand1=operand1, op=op, operand2=operand2, result=result)

    iraf.imarith.setParList(pars)

def imcombine(filelist, out, options, bpmask=None, reject="none", nlow=None,
        nhigh=None):

    '''Convenience wrapper around IRAF task imcombine

    Args:
        filelist: The list of files to imcombine
        out: The full path to the output file
        options: Options dictionary
        bpmask: The full path to the bad pixel mask
        reject: none, minmax, sigclip, avsigclip, pclip
        nlow,nhigh: Parameters for minmax rejection, see iraf docs
    
    Returns:
        None

    Side effects:
        Creates the imcombined file at location `out'
    '''

    #TODO: REMOVE Iraf and use python instead. STSCI Python has
    # A builtin routine.
    from pyraf import iraf
    iraf.images()


    filelist = [("%s[0]" % f) for f in filelist]
    pars = iraf.imcombine.getParList()
    iraf.imcombine.unlearn()

    path = "flatcombine.lst"
    f = open(path, "w")
    for file in filelist:
        f.write(file + "\n")
    f.close()

    s = ("%s," * len(filelist))[0:-1]
    s = s % tuple(filelist)

    f = open("flatcombinelog.txt", "w")
    if reject == 'minmax':
        t = iraf.imcombine("@%s" % path, out, Stdout=f,
            reject=reject, nlow=nlow, nhigh=nhigh)
    else:
        t = iraf.imcombine(s, out, Stdin=filelist, Stdout=f,
            reject=reject)
    f.close()
    f=open("flatcombinelog.txt")
    for line in f:
        info(line.rstrip("\n"))
    f.close()

    iraf.imcombine.setParList(pars)

def get_ao_performance_data():
    # Retrieve information from image headers
    outfile = workdir + 'data/ao_perform.txt'

    ir.images()
    ir.imutil()
    ir.hselect('*_img.fits', 'DATAFILE,DATE-OBS,UTC,AOLBFWHM,LGRMSWF', 'yes', Stdout=outfile)

    # Retrieve MASS/DIMM data files
    # get_mass_dimm.go('20081021')   # this one doesn't exist
    get_mass_dimm.go('20100815')
    get_mass_dimm.go('20100828')
    get_mass_dimm.go('20100829')

def imcombine(filelist, out, listfile=None, bpmask=None, reject="none", 
                nlow=None, nhigh=None):

    """Convenience wrapper around IRAF task imcombine

    Args:
        filelist (list of str): The list of files to imcombine
        out (str): The full path to the output file
        bpmask (str): The full path to the bad pixel mask
        reject (str): none, minmax, sigclip, avsigclip, pclip
        nlow,nhigh (int,int): Parameters for minmax rejection, see iraf docs
    
    Returns:
        None

    Side effects:
        Creates the imcombined file at location `out`

    """

    #TODO: REMOVE Iraf and use python instead. STSCI Python has
    # A builtin routine.
    from pyraf import iraf
    iraf.images()

    filelist = [("%s[0]" % f) for f in filelist]
    pars = iraf.imcombine.getParList()
    iraf.imcombine.unlearn()

    if listfile is None:
        path = "flatcombine.lst"
    else: path = listfile
    f = open(path, "w")
    for file in filelist:
        f.write(file + "\n")
    f.close()

    s = ("%s," * len(filelist))[0:-1]
    s = s % tuple(filelist)

    if reject == 'minmax':
        t = iraf.imcombine("@%s" % path, out, combine="average",
            reject=reject, nlow=nlow, nhigh=nhigh)
    elif reject == 'sigclip':
        t = iraf.imcombine("@%s" % path, out, combine="average",
            reject=reject, lsigma=nlow, hsigma=nhigh)
    else:
        t = iraf.imcombine(s, out, Stdin=filelist, Stdout=1, combine="average",
            reject=reject)

    iraf.imcombine.setParList(pars)

def sub_bias(image, combined_bias, image_b):
# Import IRAF modules:
  iraf.images(_doprint=0)
  iraf.imutil(_doprint=0)
  parList = "bias_subtraction_imarith.par"
# Check input file and combined_bias frame exists before proceeding:
  if os.path.isfile(image) == True:
    if os.path.isfile(combined_bias) == True:
      if os.path.isfile(parList) == True:
# Subtract combined bias frame from input frame (object or flat)
# using IRAF task imarithmetic:
        iraf.imarith.setParList(ParList="bias_subtraction_imarith.par")
        iraf.imarith(operand1=image, operand2=combined_bias, result=image_b)
        print ' '
        print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        print 'Bias frame ' + str(combined_bias) 
        print 'subtracted from input ' + str(image)
        print 'to create ' + str(image_b)
        print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        print ' '
      else:
        print ' '
        print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        print 'Bias frame subtraction IRAF .par file              ' 
        print str(parList)
        print 'does not exist. Exiting script.                    '
        print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        print ' '
        print ' '
        sys.exit()
    else:
      print ' '
      print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      print 'Combined bias frame                                ' 
      print str(combined_bias)
      print 'does not exist. Exiting script.                    '
      print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      print ' '
      print ' '
      sys.exit()
  else:
    print ' '
    print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    print 'Input frame                                        ' 
    print str(image)
    print 'does not exist. Exiting script.                    '
    print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    print ' '
    print ' '
    sys.exit()

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 Geoxy(path):
    print "Starting Geoxy..."
    filename = path
    os.chdir(filename)   
    corrected = glob.glob("*.xy")
                   
    #Removing existing files    
    filelist = glob.glob("*.done")
    for f in filelist:
        os.remove(f)
    
    #Transforming coordinates
    print "Saving files with the new coordinates."
    for i in corrected:
        iraf.images(_doprint=0)
        iraf.immatch(_doprint=0)
        iraf.geoxy(input=i, output=i+'.done', database='data3', transforms='trans3')     

def test_which(tmpdir):
    iraf.plot(_doprint=0)
    iraf.images(_doprint=0)
    outfile = str(tmpdir.join('output.txt'))

    # To Test: normal case, disambiguation, ambiguous, not found, multiple
    #          inputs for a single command
    cases = ("pw", "lang", "stdg", "stdp", "star", "star man", "vi", "noao",
             "impl", "ls", "surf", "surface", "img", "pro", "prot", "prow",
             "prows", "prowss", "dis")

    # Test the which function
    for arg in cases:
        args = arg.split(" ")  # convert to a list
        iraf.printf("--> which " + arg + '\n', StdoutAppend=outfile)
        kw = {"StderrAppend": outfile}
        iraf.which(*args, **kw)  # catches stdout+err

    diff_outputs(outfile, 'data/cli_which_output.ref')

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 rvsao(bin_sci, task, template_spectra, rvsao_file, rvsao_bin_list, interactive="no", linesig=1.5,
          czguess=0., st_lambda="INDEF", end_lambda="INDEF", badlines=None):
    """
    Use the rvsao task emsao or xcsao to measure relative velocity from emission or absoption spectra
    Note: make sure you select the parameters for your desired use. Refer to the IRAF rvsao package help for description
    of what each parameter is/does
    """

    assert task == 'xcsao' or task == 'emsao', "task is not either 'xcsao' or 'emsao'"

    if os.path.exists(rvsao_file):
        print('File {} already exists'.format(rvsao_file))
        return

    fixbad = "no"
    if badlines:
        fixbad = "yes"

    bin_list = []
    for i in range(len(glob.glob(bin_sci.format('*')))):  # to ensure order is 0-61 (not 0, 1, 10, 11, etc)
        bin_list.append(bin_sci.format(i))
    assert len(bin_list) > 0, 'Absorption/emission(?) bin spectra do not exist: {}'.format(em_bin_sci.format('*'))
    np.array(bin_list).tofile(rvsao_bin_list, sep='\n')

    from pyraf import iraf

    iraf.images()
    iraf.rvsao()

    if task == 'xcsao':
        iraf.xcsao('@{}'.format(rvsao_bin_list), templates=bin_sci.format(template_spectra), report_mode=2,
                   logfiles=rvsao_file, displot=interactive, low_bin=10, top_low=20, top_nrun=80, nrun=211,
                   zeropad="yes", nzpass=1, curmode="no", pkmode=2, s_emchop="no", vel_init="guess", czguess=czguess,
                   st_lambda=st_lambda, end_lambda=end_lambda, fixbad=fixbad, badlines=badlines)

    elif task == 'emsao':

        # Run this interactively as the fit can fail often. Depending on S/N you may have to decrease the linesig
        # to detect the lines, or if lines are matched inconsistently play with the st_lambda, end_lambda parameters
        iraf.emsao('@{}'.format(rvsao_bin_list), logfiles=rvsao_file, displot=interactive, report_mode=2,
                   contsub_plot="no", st_lambda=st_lambda, end_lambda=end_lambda,  # vel_init="guess", czguess=czguess,
                   linesig=linesig, fixbad=fixbad, badlines=badlines)

def mosfire_geoxytran(x_kfp, y_kfp, transform="final.pix2mm.4.972.120k",
        database="/platescale/10March2011.4.972.db", direction="forward"):
    '''Conveninece wrapper around IRAF geoxytran'''
    iraf.images()


    path = os.path.join(os.environ["MOSPATH"], "platescale", 
            "10March2011.4.972.db")
    database = path
    pars = iraf.geoxytran.getParList()
    iraf.geoxytran.unlearn()
    t = iraf.geoxytran("STDIN", "STDOUT", Stdin=["%f %f" % (x_kfp, y_kfp)], Stdout=1,
            database=database,
            transform=transform,
            direction=direction)

    iraf.geoxytran.setParList(pars) 
    (x,y) = np.array(t[0].split(), dtype=np.float64)

    return (x,y)

def flat_field(object_b, combined_flat_b_n, object_b_fn):
# Import IRAF modules:
  iraf.images(_doprint=0)
  iraf.imutil(_doprint=0)
# Check input file and combined_flat frame exist before proceeding:
  if os.path.isfile(object_b) == True:
    if os.path.isfile(combined_flat_b_n) == True:
# Divide bias-subtracted object frame by normalized, bias
# subtracted combined flat frame using IRAF task imarithmetic:
      iraf.imarith(operand1=object_b, op="/", operand2=combined_flat_b_n, result=object_b_fn)
      print ' '
      print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      print 'Bias-subtracted object frame                       '
      print str(object_b)
      print 'successfully flat-fielded using division by        '
      print str(combined_flat_b_n)
      print 'to create bias-subtracted, normalized flat-fielded '
      print str(object_b_fn) + ' frame.'
      print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      print ' '
    else:
      print ' '
      print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      print 'Combined flat frame                                ' 
      print str(combined_flat_b_n)
      print 'does not exist. Exiting script.                    '
      print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      print ' '
      print ' '
      sys.exit()
  else:
    print ' '
    print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    print 'Input frame                                        ' 
    print str(object_b)
    print 'does not exist. Exiting script.                    '
    print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    print ' '
    print ' '
    sys.exit()

def reduce_images(progargs):
    """Top level function to perform the reduction of data images. 
    
    The tasks are performed sequentially: generate all masterbias, 
    generate all masterflats and finally reduce data images.
    
    Args:
        progargs: Program arguments.     
    
    """
    
    logging.info("Starting the reduction of images ...")    

    # Load the images package and does not show any output.
    iraf.images(_doprint=0)

    # Generate all the average bias.
    generate_all_masterbias(progargs.target_dir,
                            progargs.bias_directory)

    # Generate all the average dark.
    generate_all_masterdark(progargs.target_dir,
                            progargs.dark_directory,
                            progargs.bias_directory)

    # Generate all the average flat.
    generate_all_masterflats(progargs.target_dir,
                             progargs.flat_directory,
                             progargs.dark_directory,
                             progargs.bias_directory)

    # Reduce all the data images applying the average bias and flats.
    reduce_data_images(progargs.target_dir,
                       progargs.light_directory,
                       progargs.dark_directory,
                       progargs.bias_directory,
                       progargs.flat_directory)
    
    logging.info("Finished the reduction of images.")    

def mosfire_geoxytrans(
    x_kfp, y_kfp, transform="final.pix2mm.4.972.120k", database="ale/10March2011.4.972.db", direction="forward"
):
    """Conveninece wrapper around IRAF geoxytran"""
    iraf.images()
    path = os.path.join(os.environ["MOSPATH"], "platescale", "10March2011.4.972.db")
    database = path
    pars = iraf.geoxytran.getParList()
    iraf.geoxytran.unlearn()
    ins = []
    for i in range(len(x_kfp)):
        ins.append("%f %f" % (x_kfp[i], y_kfp[i]))
    results = iraf.geoxytran(
        "STDIN", "STDOUT", Stdin=ins, Stdout=1, database=database, transform=transform, direction=direction
    )

    iraf.geoxytran.setParList(pars)

    poss = []
    for result in results:
        poss.append(np.array(result.split(), dtype=np.float64))

    return np.array(poss)

def shift_master(ref, master):
    """
    Find the shift between two images.
    returns the shift in y direction (along columns)
    """
    #make a narow copy of the reference flat and masterflat 
    iraf.imcopy(input="../"+ref+"[1990:2010,*]", output="tmp/masterflat_ref.fitscut",Stdout="/dev/null")
    iraf.imcopy(input="tmp/masterflat.fits[1990:2010,*]", output="tmp/masterflat.fitscut",Stdout="/dev/null")

    #find the shift with correlation
    iraf.images(_doprint=0,Stdout="/dev/null")
    iraf.immatch(_doprint=0,Stdout="/dev/null")

    shift=iraf.xregister(input="tmp/masterflat_ref.fitscut, tmp/masterflat.fitscut", reference="tmp/masterflat_ref.fitscut", regions='[*,*]', shifts="tmp/mastershift", output="", databasefmt="no", correlation="fourier", xwindow=3, ywindow=51, xcbox=21, ycbox=21, Stdout=1)

    shift=float(shift[-2].split()[-2])

    #Shift the list of apertures
    f=open(master, "r")
    o=open("database/aptmp_masterflat_s", "w")

    for line in f:
        l=line.split()
        if len(l)>0 and l[0]=='begin':
            l[-1]=str(float(l[-1])-shift)
            o.write(l[0]+"\t"+l[1]+"\t"+"tmp/masterflat_s"+"\t"+l[3]+"\t"+l[4]+"\t"+l[5]+"\n")
        elif len(l)>0 and l[0]=='center':
            l[-1]=str(float(l[-1])-shift)
            o.write("\t"+l[0]+"\t"+l[1]+"\t"+l[2]+"\n")
        elif len(l)>0 and l[0]=='image':
            o.write("\t"+"image"+"\t"+"tmp/masterflat_s"+"\n")
        else:
            o.write(line)

    f.close()
    o.close()

def imarith(operand1, op, operand2, result, doAirmass=False):
    from pyraf import iraf
    iraf.images()

    pars = iraf.imarith.getParList()
    iraf.imcombine.unlearn()

    try: os.remove(result)
    except: pass

    print "%s %s %s -> %s" % (operand1, op, operand2, result)
    iraf.imarith(operand1=operand1, op=op, operand2=operand2, result=result)
    iraf.imarith.setParList(pars)   
    if doAirmass:
        # Adjust FITS header
        with pf.open(operand1) as f:
            am1 = f[0].header['airmass']
        with pf.open(operand2) as f:
            am2 = f[0].header['airmass']

        of = pf.open(result)
        of[0].header['airmass1'] = am1
        of[0].header['airmass2'] = am2
        of.writeto(result, clobber=True)

def telluric_atmo(imgstd):
    # print "LOGX:: Entering `telluric_atmo` method/function in %(__file__)s"
    # % globals()
    import numpy as np
    import ntt
    from pyraf import iraf

    try:        import pyfits
    except:     from astropy.io import fits as pyfits

    iraf.images(_doprint=0)
    iraf.noao(_doprint=0)
    iraf.twodspec(_doprint=0)
    iraf.longslit(_doprint=0)
    iraf.onedspec(_doprint=0)
    toforget = ['imfilter.gauss', 'specred.apall', 'longslit.identify', 'longslit.reidentify', 'specred.standard',
                'onedspec.wspectext']
    for t in toforget:
        iraf.unlearn(t)

    _grism = ntt.util.readkey3(ntt.util.readhdr(imgstd), 'grism')
    imgout = 'invers_atmo_' + imgstd
    ntt.util.delete(imgout)
    iraf.set(direc=ntt.__path__[0] + '/')
    _cursor = 'direc$standard/ident/cursor_sky_0'
    iraf.noao.onedspec.bplot(imgstd, cursor=_cursor,
                             spec2=imgstd, new_ima=imgout, overwri='yes')
    xxstd, ffstd = ntt.util.readspectrum(imgout)
    if _grism in ['Gr13', 'Gr16']:
        llo2 = np.compress((np.array(xxstd) >= 7550) & (
            np.array(xxstd) <= 7750), np.array(xxstd))
        llh2o = np.compress((np.array(xxstd) >= 7100) & (
            np.array(xxstd) <= 7500), np.array(xxstd))
        ffo2 = np.compress((np.array(xxstd) >= 7550) & (
            np.array(xxstd) <= 7750), np.array(ffstd))
        ffh2o = np.compress((np.array(xxstd) >= 7100) & (
            np.array(xxstd) <= 7500), np.array(ffstd))
    elif _grism in ['Gr11']:
        llo2 = np.compress((np.array(xxstd) >= 6830) & (
            np.array(xxstd) <= 7100), np.array(xxstd))
        llh2o = np.compress((np.array(xxstd) >= 7100) & (
            np.array(xxstd) <= 7500), np.array(xxstd))
        ffo2 = np.compress((np.array(xxstd) >= 6830) & (
            np.array(xxstd) <= 7100), np.array(ffstd))
        ffh2o = np.compress((np.array(xxstd) >= 7100) & (
            np.array(xxstd) <= 7500), np.array(ffstd))
    if _grism in ['Gr13', 'Gr16', 'Gr11']:
        _skyfileh2o = 'direc$standard/ident/ATLAS_H2O.fits'
        _skyfileo2 = 'direc$standard/ident/ATLAS_O2.fits'
        atlas_smooto2 = '_atlas_smoot_o2.fits'
        atlas_smooth2o = '_atlas_smoot_h2o.fits'
        _sigma = 200
        ntt.util.delete(atlas_smooto2)
        ntt.util.delete(atlas_smooth2o)
        iraf.imfilter.gauss(_skyfileh2o, output=atlas_smooth2o, sigma=_sigma)
        iraf.imfilter.gauss(_skyfileo2, output=atlas_smooto2, sigma=_sigma)
        llskyh2o, ffskyh2o = ntt.util.readspectrum(atlas_smooth2o)
        llskyo2, ffskyo2 = ntt.util.readspectrum(atlas_smooto2)
        ffskyo2cut = np.interp(llo2, llskyo2, ffskyo2)
        ffskyh2ocut = np.interp(llh2o, llskyh2o, ffskyh2o)
        _scaleh2o = []
        integral_h2o = []
        for i in range(1, 21):
            j = 0.6 + i * 0.04
            _ffskyh2ocut = list((np.array(ffskyh2ocut) * j) + 1 - j)
            diff_h2o = abs(_ffskyh2ocut - ffh2o)
            integraleh2o = np.trapz(diff_h2o, llh2o)
            integral_h2o.append(integraleh2o)
            _scaleh2o.append(j)
        _scaleo2 = []
        integral_o2 = []
        for i in range(1, 21):
            j = 0.6 + i * 0.04
            _ffskyo2cut = list((np.array(ffskyo2cut) * j) + 1 - j)
            diff_o2 = abs(_ffskyo2cut - ffo2)
            integraleo2 = np.trapz(diff_o2, llo2)
            integral_o2.append(integraleo2)
            _scaleo2.append(j)
        sh2o = _scaleh2o[np.argmin(integral_h2o)]
        so2 = _scaleo2[np.argmin(integral_o2)]
        telluric_features = ((np.array(ffskyh2o) * sh2o) +
                             1 - sh2o) + ((np.array(ffskyo2) * so2) + 1 - so2) - 1
        telluric_features = np.array([1] + list(telluric_features) + [1])
        llskyo2 = np.array([1000] + list(llskyo2) + [15000])
        telluric_features_cut = np.interp(xxstd, llskyo2, telluric_features)
        _imgout = 'atmo_' + imgstd

        data1, hdr = pyfits.getdata(imgstd, 0, header=True)
        data1[0] = np.array(telluric_features_cut)
        data1[1] = data1[1] / data1[1]
        data1[2] = data1[2] / data1[2]
        data1[3] = data1[3] / data1[3]
        ntt.util.delete(_imgout)
        pyfits.writeto(_imgout, np.float32(data1), hdr)
        ntt.util.delete(atlas_smooto2)
        ntt.util.delete(atlas_smooth2o)
        ntt.util.delete(imgout)
    else:
        _imgout = ''
        print '### telluric correction with model not possible '
#.........这里部分代码省略.........

import pyraf
from pyraf import iraf
import copy, os, shutil, glob, sys, string, re, math, operator, time
import pyfits
from types import *
from mx.DateTime import *

from iqpkg import *
import ephem

# Necessary packages
iraf.images()
iraf.immatch()
iraf.imfilter()
iraf.noao()
iraf.imred()
iraf.ccdred()
iraf.digiphot()
iraf.apphot()

yes=iraf.yes
no=iraf.no
INDEF=iraf.INDEF
hedit=iraf.hedit
imgets=iraf.imgets
imcombine=iraf.imcombine

pyrafdir="python/pyraf/"
pyrafdir_key='PYRAFPARS'

if os.environ.has_key(pyrafdir_key):