def stdsensfunc(fs=None):
    iraf.cd('work')
    if fs is None:
        fs = glob('x1d/sci*x1d*c?.fits')
    if len(fs) == 0:
        print "WARNING: No extracted spectra to create sensfuncs from."
        iraf.cd('..')
        return

    if not os.path.exists('std'):
        os.mkdir('std')
    for f in fs:
        # Put the file in the std directory, but last 3 letters of sens
        outfile = 'std/' + f.split('/')[1]
        outfile = outfile.replace('x1d', 'sens').replace('sci', 'std')
        outfile = outfile.replace('.fits', '.dat')
        # if the object name is in the list of standard stars from pysalt
        if isstdstar(f):
            # We use pysalt here because standard requires a
            # dispersion correction which was already taken care of above
            # Write out an ascii file that pysalt.specsens can read
            asciispec = 'std/std.ascii.dat'
            spectoascii(f, asciispec)
            # run specsens
            stdfile = pysaltpath + '/data/standards/spectroscopic/m%s.dat' % pyfits.getval(f, 'OBJECT').lower().replace('-','_')
            extfile = pysaltpath + '/data/site/suth_extinct.dat'
            iraf.unlearn(iraf.specsens)
            iraf.specsens(asciispec, outfile, stdfile, extfile,
                          airmass=pyfits.getval(f, 'AIRMASS'),
                          exptime=pyfits.getval(f, 'EXPTIME'), function='poly',
                          order=11, clobber=True, mode='h', thresh=1e10)
            # delete the ascii file
            os.remove(asciispec)
    iraf.cd('..')

def mark_boundaries(flist, slit_size, combined_img, rot_img, ax1, img_angle, ref_angle, fudge_factor, target_font_size):
    flist_r = flist[::-1]

    angle =  (ref_angle - img_angle) - 90
    angle_rad = angle/180.0 * math.pi
    ydim, xdim = np.shape(combined_img)
    x = np.arange( (ydim) * math.cos(angle_rad))
    
    slope = -math.tan(angle_rad)
    for indx, i in enumerate(np.arange(xdim/int(slit_size)+1)):                             
        intercept = (ydim) * math.sin(angle_rad) + ((i)*slit_size/math.cos(angle_rad))
        xnew = x+((i)*math.sin(angle_rad)*slit_size)
        y = slope*xnew+intercept
        ax1.plot(xnew, y, 'k')
        try:
            if np.mod(i, 2) == 1:
                ax1.text(math.ceil(xnew[-1]+1), math.ceil(y[-1]), '%s' %(pyfits.getval(flist_r[indx].replace('collapse_flt', 'combined_img'), 'targname', 0)) , fontsize = target_font_size, ha = 'left', va = 'bottom')
                ax1.plot(np.arange(50) +  math.ceil(xnew[-1]), np.ones(50) * math.ceil(y[-1]) - 1, 'k')
                ax1.plot((-np.arange(50)) +  math.ceil(xnew[0]), np.ones(50) * math.ceil(y[0]) - 1, 'k')
            else:
                ax1.text(math.ceil(xnew[0]-1), math.ceil(y[0]), '%s' %(pyfits.getval(flist_r[indx].replace('collapse_flt', 'combined_img'), 'targname', 0)) , fontsize = target_font_size, ha = 'right', va = 'bottom')
                ax1.plot((-np.arange(50)) +  math.ceil(xnew[0]) - 1, np.ones(50) * math.ceil(y[0]) - 1, 'k')
                ax1.plot(np.arange(50) +  math.ceil(xnew[-1]) - 1, np.ones(50) * math.ceil(y[-1]) - 1, 'k')
        except:
            pass
    ax1 = mark_pix_loc(xdim, ydim, angle_rad, ax1)
    #pdb.set_trace()
    return ax1

    def inputs(cls, gainmap):
        """ return the pulse height image name from the input gainmap

        """

        gainmap_path, gainmap_name = os.path.split(gainmap)
        segment = pyfits.getval(gainmap, 'SEGMENT')
        dethv = int(pyfits.getval(gainmap, 'DETHV'))

        both_inputs = [gainmap]

        if segment == 'FUVA':
            other_root = gainmap.replace(FUVA_string, FUVB_string).replace('_{}_'.format(dethv), '_???_')
        elif segment == 'FUVB':
            other_root = gainmap.replace(FUVB_string, FUVA_string).replace('_{}_'.format(dethv), '_???_')

        other_gainmap = glob.glob(other_root)
        if len(other_gainmap) != 1:
            raise IOError("too many gainmaps found {}".format(other_gainmap))
        else:
            other_gainmap = other_gainmap[0] 

        both_inputs.append(other_gainmap)
        both_inputs.sort()

        return tuple(both_inputs)

def update_psfdat():
    import os
    from astropy.io import ascii
    from astropy.table import Column
    import pyfits
    # read in the psf mags computed by Dan
    psfdat = ascii.read( 'psfphot.dat')

    # get exposure time and number of exposures from the header
    etimelist, nexplist, etcSNRlist, etcOSNRlist, psfSNRlist = [], [], [], [], []
    for row in psfdat :
        imdir = '/store/snreproc/%s.090mas/'%(row['sn'])
        suffix = '-e00_sub_masked' if row['image']=='sub' else '_reg_drz_sci'
        imfile = os.path.join( imdir,
            '%s_%s_%s%s.fits'%(
                row['sn'], row['band'], row['epoch'], suffix ) )
        etime = pyfits.getval( imfile, 'EXPTIME' )
        nexp = pyfits.getval( imfile, 'NCOMBINE' )
        etimelist.append( etime )
        nexplist.append( nexp )
        etcSNR, etcSNRopt = synphotSNR( row['band'], row['stack'], etime=etime, nexp=nexp)
        etcSNRlist.append( etcSNR )
        etcOSNRlist.append( etcSNRopt )
        psfSNR = 1.08574 / row['err']
        psfSNRlist.append( psfSNR )

    nexpCol = Column( nexplist, 'nexp', dtype=int, format='%i')
    etimeCol = Column( etimelist, 'etime', dtype=float, format='%8.2f' )
    etcSNRCol = Column( etcSNRlist, 'SNRetc', dtype=float, format='%6.2f' )
    etcOSNRCol = Column( etcOSNRlist, 'optSNRetc', dtype=float, format='%6.2f' )
    psfSNRCol = Column( psfSNRlist, 'SNRpsf', dtype=float, format='%6.2f' )

    psfdat.add_columns( [ nexpCol, etimeCol, etcSNRCol, etcOSNRCol, psfSNRCol], indexes=[4,4,9,9,9])
    psfdat.write( 'psfphot2.dat', format='ascii.commented_header' )
    psfdat.write( 'psfphot3.dat', format='ascii.fixed_width' )

def combine_medium(file_list , num=70):
    '''(str, int) -> ndarray, asciidata

    Puts list of fits files into correct order for use with
    median function. Adds to fits header histrory.
    '''
    #does medium assumes all have same fits header
    #high ram but quick
    if len(file_list) < num:
        shape = (fits.getval(file_list[0],'NAXIS1'),
                 fits.getval(file_list[0],'NAXIS2'),
                 len(file_list))
        temp = np.zeros(shape)
        for i,j in enumerate(file_list):
            temp[:,:,i],hdr = util.fromfits(j,verbose=False)
        out = np.median(temp,2)
    else:
        #low ram but takes longer
        shape = (fits.getval(file_list[0],'NAXIS1'),
                 fits.getval(file_list[0],'NAXIS2'))
        out = np.zeros(shape)
        temp = np.zeros((shape[0],len(file_list)))
        for i in xrange(shape[1]):
            for j,k in enumerate(file_list):
                temp[:,j] = util.fromfits(k, verbose=False)[0][:,i]
            out[:,i] = np.median(temp,1)
        temp,hdr = util.fromfits(k, verbose=False)
    hdr.add_history('Medium combine')
    return out, hdr

def wavelength_array(spec):
    """ Produces array for wavelenght of a given array. """
    w0 = pf.getval(spec, "CRVAL1")
    deltaw = pf.getval(spec, "CDELT1")
    pix0 = pf.getval(spec, "CRPIX1")
    npix = pf.getval(spec, "NAXIS1")
    return w0 + deltaw * (np.arange(npix) + 1 - pix0)

def update_wcs_from_pointing(files,pixscale):
   for f in files:
      if f[-5:]!=".fits":
         continue
      print "%s ....Updating values" % f
      fullfits = pyfits.open(f,mode='update')
      sci = fullfits[0]
      shape = sci.data.shape
      ra = wcs.ra2deg(sci.header['RA'].strip())
      dec = wcs.dec2deg(sci.header['DEC'].strip())
      sci.header.update('CRPIX1',shape[1]/2.)
      sci.header.update('CRVAL1',ra)
      sci.header.update('CRPIX2',shape[0]/2.)
      sci.header.update('CRVAL2',dec)
      sci.header.update('CD1_1',-pixscale/3600.)
      sci.header.update('CD1_2',0.)
      sci.header.update('CD2_1',0.)
      sci.header.update('CD2_2',pixscale/3600.)
      sci.header.update('CTYPE1','RA---TAN')
      sci.header.update('CTYPE2','DEC--TAN')
      sci.header.update('EQUINOX',2000.0)
      sci.header.update('RADESYS','FK5')
      del sci.header['CDELT1']
      del sci.header['CDELT2']
      fullfits.flush()
      print "  n   CTYPEn  CRPIXn    CRVALn       CDn_1        CDn_2"
      print " --- -------- ------- ---------- ------------- -------------"
      print "  1  %s %7.2f %10.6f %13.6e %13.6e" \
       % (getval(f,'CTYPE1'),getval(f,'CRPIX1'),getval(f,'CRVAL1'), \
       getval(f,'CD1_1'),getval(f,'CD1_2'))
      print "  2  %s %7.2f %10.6f %13.6e %13.6e" \
       % (getval(f,'CTYPE2'),getval(f,'CRPIX2'),getval(f,'CRVAL2'), \
       getval(f,'CD2_1'),getval(f,'CD2_2'))

def wavelength_array(spec, axis=1, extension=0):
    """ Produces array for wavelenght of a given array. """
    w0 = pf.getval(spec, "CRVAL{0}".format(axis), extension)
    deltaw = pf.getval(spec, "CD{0}_{0}".format(axis), extension)
    pix0 = pf.getval(spec, "CRPIX{0}".format(axis), extension)
    npix = pf.getval(spec, "NAXIS{0}".format(axis), extension)
    return w0 + deltaw * (np.arange(npix) + 1 - pix0)

def combine_sigmaclip(file_list, num=70):
    '''(str, int) -> ndarray, asciidata

    Puts list of fits files into correct order for use with
    sigma clipping function. Adds to fits header histrory.
    '''
    #fast but uses lots of ram
    if len(file_list) < num:
        shape = (fits.getval(file_list[0],'NAXIS1'),
                 fits.getval(file_list[0],'NAXIS2'),
                 len(file_list))
        temp = np.zeros(shape)
        for i,j in enumerate(file_list):
            temp[:,:,i],hdr = fits.fromfits(j,verbose=False)
        out = Sigmaclip(temp,axis=2)
    else:
        #slow ram but takes longer
        shape = (fits.getval(file_list[0],'NAXIS1'),
                 fits.getval(file_list[0],'NAXIS2'))
        out = np.zeros(shape)
        temp = np.zeros((shape[0],len(file_list)))
        for i in xrange(shape[1]):
            for j,k in enumerate(file_list):
                temp[:,j] = util.fromfits(k,verbose=False)[0][:,i]
            out[:,i] = Sigmaclip(temp,axis=1)
        temp, hdr = util.fromfits(k,verbose=False)
    hdr.add_history('Sigmaclip combine')
    return out, hdr

    def __process_fits_object(self, fits_string):
        '''
        Process FITS file object and extract info.
        http://stackoverflow.com/questions/11892623/python-stringio-and-compatibility-with-with-statement-context-manager
        Returns the temporary file name and DataStream object.
        '''
        test = ''
        with self.__tempinput(fits_string) as tempfilename:
            test = tempfilename
            fitsdata = pyfits.getdata(tempfilename)
            bjd_trunci = float(pyfits.getval(tempfilename, 'bjdrefi', ext=1))
            bjd_truncf = float(pyfits.getval(tempfilename, 'bjdreff', ext=1))

            # Note: Times are updated to be in proper reduced barycentric Julian date,
            # RBJD = BJD - 2400000.0
            time = fitsdata['TIME'] + bjd_trunci + bjd_truncf - 2400000.
            pdcflux = fitsdata['PDCSAP_FLUX']
            pdcerror = fitsdata['PDCSAP_FLUX_ERR']
            errorstat = fitsdata['SAP_QUALITY']

            ndx = np.where(errorstat == 0)
            retval = DataStream(arrays=(time[ndx], pdcflux[ndx], pdcerror[ndx]))

            # Fix for windows, returns the filename into main so that os.unlink can be called there
            return test, retval

def get_ims(fs, imtype):
    imtypekeys = {'sci': 'OBJECT', 'arc': 'ARC', 'flat': 'FLAT'}
    ims = []
    grangles = []
    for f in fs:
        if pyfits.getval(f, 'OBSTYPE') == imtypekeys[imtype]:
            ims.append(f)
            grangles.append(pyfits.getval(f, 'GR-ANGLE'))
    return np.array(ims), np.array(grangles)

def make_counts_image(image, outfile='default'):
    '''FUNCTION TO CONVERT CNTS/SEC IMAGE TO COUNTS (IF NECESSARY)'''

    # -- parse output filename & save a copy to file (NOTE: if outfile == input image, data is overwritten).
    if (image != outfile):
        if outfile == 'default': outfile = image.split('.fits')[0] + '_cnts.fits'
        shutil.copy(image,outfile)
    else: print 'OVERWRITING DATA FOR IMAGE: '+image+'.'


    # -- determine if image is flt/flc, crclean, or drz/drc
    prihdr = pyfits.getheader(outfile,ext=0)
    pscale = prihdr.get('D001SCAL',default='NA')
    if pscale != 'NA': imtype = 'drz'
    elif len(image.split('crclean.fits')) > 1: imtype = 'crclean'
    else: imtype = 'flt'


    # -- initialize a few required parameters
    detector = prihdr['DETECTOR']
    exptime = prihdr['EXPTIME']


    # -- multiply by exposure time (only if image is already in cnts/sec)
    #      [notes] -- IR crcleans are actually in cnts, but "BUNIT" still says per second (can't trust bunit for now).
    #              -- We assume drz are cnts/sec, but this does not have to be true.
    #
    if imtype == 'drz':
        # -- save background & pixel scale info
        if prihdr['extend'] == True: back = pyfits.getval(outfile,'MDRIZSKY',ext=1)  
        else: back = pyfits.getval(outfile,'MDRIZSKY',ext=0)
        pscale_nat = pyfits.getval(outfile,'D001ISCL',ext=0)
        pscale_img = pyfits.getval(outfile,'D001SCAL',ext=0)

        # -- assign the number of chips associated with this image
        if (prihdr['detector'] == 'IR'): nchips = 1.0                                          # IR
        elif (prihdr['subarray'] == True) and (len(prihdr['CCDAMP']) == 1): nchips = 1.0      # UVIS sub-array
        elif (prihdr['detector'] == 'UVIS') and (prihdr['subarray'] == False): nchips = 2.0   # UVIS full-frame
        else: raise exception('Image type is not defined.')

        # -- add background and correct for different pixel scale (original backgrd is measured in raw images)
        fdata = pyfits.getdata(outfile,ext=0)
        fdata_cnts = np.copy(fdata) * exptime + np.sum(back)/nchips * (pscale_img/pscale_nat)**2
        hdulist = pyfits.open(outfile,mode='update')
        hdulist[0].data = fdata_cnts
        hdulist.close()

    elif ((detector == 'IR') & (imtype == 'flt')):
        hdulist = pyfits.open(outfile,mode='update')
        for ff in xrange(len(hdulist)):
            if hdulist[ff].name == 'SCI': hdulist[ff].data = hdulist[ff].data * exptime
        hdulist.close()

    else: print 'IMAGE SHOULD ALREADY BE IN UNITS OF COUNTS. RETURNING...'

    return outfile

def info_license(fitsfile):
    """Print license information stored in FITS"""
    try:
        license = pyfits.getval(fitsfile, "LICENSE")
    except KeyError:
        print("License information not found.")
    else:
        licver = pyfits.getval(fitsfile, "LICVER")
        licurl = pyfits.getval(fitsfile, "LICURL")
        print("{lic} {ver} ({url})".format(lic=license, ver=licver, url=licurl))

def determine_correct_reference_files(input_dir, input_list, filetype):
    '''
    This function creates a nested dictionary which assigns a reference file name (from 
    the reffiles folder) to each gain, binaxis1, binaxis2, useafter combination
    Inputs:
        input_dir: directory of files whose headers you want to update
        input_list: names of files whose headers you want to update
        filetype: either drk or bia - the reference file whose name you want to update
    Output:
        mode_dict: a nested dictionary which hich assigns a reference file name (from 
    the reffiles folder) to each gain, binaxis1, binaxis2, useafter combination

    This code assumes that you have set the environment variable myref to point to the 
    reffiles folder
    '''
    dates = []
    for ifile in input_list:
        filename = os.path.join(input_dir, ifile)
        dates.append(pyfits.getval(filename, 'texpstrt', 0))
        dates.append(pyfits.getval(filename, 'texpend', 0))
    data_start = min(dates)
    data_end = max(dates)
    anneal_weeks_4 = divide_anneal_month(data_start, data_end, '/Users/bostroem/science/cte/', 4)
    anneal_weeks_2 = divide_anneal_month(data_start, data_end, '/Users/bostroem/science/cte/', 2)
    #nested dictionary: gain, binaxis1, binaxis2, week
                #gain
    mode_dict = {1:
                    #binaxis1
                    {1:
                        #binaxis2
                        #1x1 gain = 1
                        {1:{round(anneal_weeks_4[0][0], 4):'', round(anneal_weeks_4[1][0], 4):'', round(anneal_weeks_4[2][0], 4):'', round(anneal_weeks_4[3][0], 4):''},   
                        # 1x2 gain = 1
                        2:{round(anneal_weeks_2[0][0], 4):'', round(anneal_weeks_2[1][0], 4):''}},    
                    2: 
                        #2x1 gain = 1
                        {1:{round(anneal_weeks_2[0][0], 4):'', round(anneal_weeks_2[1][0], 4):''},  
                        #2x2 gain = 1
                        2:{round(anneal_weeks_2[0][0], 4):'', round(anneal_weeks_2[1][0], 4):''}}}, 
                4:  
                    #binaxis1
                    {1:
                        #binaxis2
                        #1x1 gain = 4
                        {1:{round(anneal_weeks_2[0][0], 4):'', round(anneal_weeks_2[1][0], 4):''}}}}   
    reffile_list = glob.glob(os.path.join(os.environ['myref'], '*_%s*.fits' %(filetype)))
    for reffile in reffile_list:
        hdr0 = pyfits.getheader(reffile, 0)
        gain = hdr0['ccdgain']
        binaxis1 = hdr0['binaxis1']
        binaxis2 = hdr0['binaxis2']
        useafter = convert_useafter_to_mjd(hdr0['useafter'])
        mode_dict[gain][binaxis1][binaxis2][round(useafter, 4)] = 'myref$%s' %(reffile.split('/')[-1])
                       
    return mode_dict

def get_darks(path=None, combine_type='median', outdir=None,
              Filter=('SET-TEMP','EXPTIME')):
    '''(str,str,str,tuple(str) or str) -> dict(ndarry), dict(asciidata)

    Opens dark directory, combines all *.fits or *.fit files in the 
    directory, sorts by Filter options (can have multiple) and out 
    puts fits file to outdir and also out ndarray of combined fits 
    files and fits header with modified history.  Combine types can 
    include: "sigmaclip","median","sum" and "sigmaclip."
    Known Issues:
    Median and sigmaclip combine give artifacts when use.'''

    comm = comb_Type(combine_type)
    #gui select directory
    if path is None:
        path = gui_getdir(title='Please Select Dark Directory')
        if not path:
            raise ValueError('Must specify directory where files are.')
    if not path.endswith('/'):
        path += '/'
    if outdir is None:
        outdir = gui_getdir(title='Please Select Save Directory')
    #load paths to fits
    fits_path = sorted(glob(path+'*'))
    fits_path = util.get_fits_type(fits_path,'dark')
    #sort by time and temp
    filters = {}
    for i in fits_path:
	if type(Filter) is tuple:
            filt = ''
            for j in Filter:
                filt += str(fits.getval(i,j)) + '_'
        else:
            filt = fits.getval(i,Filter)
        if not filt in filters.keys():
            filters[filt] = [i]
        else:
            filters[filt].append(i)

    out,hdr = {},{}
    for i in filters.keys():
        out[i],hdr[i] = comm(filters[i])
    #save as fits?
    if outdir is None:
        outdir = path
    else:
        # make dir if don't exisits
        if  not os.path.exists(outdir):
            os.mkdir(outdir)
    for i in out.keys():
        basename = os.path.split(filters[i][0])[-1]
        basename = os.path.splitext(basename)[0]
        util.tofits(outdir+basename+'_%s.fits'%combine_type.lower(), out[i],
               hdr[i],verbose=False)
    return out,hdr

def make_fits(spec, outfile):
    hdu = pf.PrimaryHDU(spec)
    miles = [x for x in os.listdir(".") if x.startswith("Mun") and x.endswith(".fits")][0]
    w0 = pf.getval(miles, "CRVAL1")
    deltaw = pf.getval(miles, "CDELT1")
    pix0 = pf.getval(miles, "CRPIX1")
    hdu.header["CRVAL1"] = w0
    hdu.header["CDELT1"] = deltaw
    hdu.header["CRPIX1"] = pix0
    pf.writeto(outfile, hdu.data, hdu.header, clobber=True)
    return

 def getcrpix(self):
     '''
     Get the cr pixel information from the single_sci.fits file.
     '''
     fitsfile = master_images_query = session.query(MasterImages.fits_file).filter(\
         MasterImages.name == os.path.basename(self.filename)).one()[0]
     png_path = os.path.split(self.filename)[0]
     png_path = os.path.split(png_path)[0]
     fitsfile = os.path.join(png_path, fitsfile)
     self.crpix1 = pyfits.getval(fitsfile, 'CRPIX1', 0)
     self.crpix2 = pyfits.getval(fitsfile, 'CRPIX2', 0)

def get_incl(name, ritter_data = "ritter_v7.19.fits"):

	from math import fabs

	ra = py.getval(name,'RAOBJ',0)
	de = py.getval(name,'DECOBJ',0)

	#print "Finding", ra, de

	array = py.getdata(ritter_data, 0)

	ra2000 = array["_RAJ2000"]
	de2000 = array["_DEJ2000"]

	bool_ra = ( np.fabs(ra2000 -ra) < POSITION_DELTA )
	bool_dec = ( np.fabs(de2000 - de) < POSITION_DELTA)

	#print ra2000

	#print np.sum(bool_ra*bool_dec)

	Found = False

	#Found = (np.sum(bool_ra*bool_dec) != 1)

	for i in range(len(ra2000)):
		if fabs(ra2000[i] - ra) < POSITION_DELTA:
			if fabs(de2000[i] - de) < POSITION_DELTA:
				Found = True
				ipt = i

	if Found:

		inc = array["Incl"]

		types = array["Type1"]

		inclin = inc[ipt]

		cv_type = types[ipt]

		#if inclin == "":
		#	print "Found match,", name, ra, de, "None"
		if inclin < 90.0 and inclin > 0.0:
			print "Found match, and inclination!", name, ra, de, inclin 
		#else:
		#	print "Found match,", name, ra, de, "--" 

		return inclin, cv_type


	else:
		print "no RA/dec match", name, ra, de
		return 0

def getobstypes(fs):
    # get the type of observation for each file
    obstypes = []
    obsclasses = []
    for f in fs: 
        obstypes.append(pyfits.getval(f, 'OBSTYPE', ext=0))
        obsclasses.append(pyfits.getval(f, 'OBSCLASS', ext=0))
        
    obstypes = np.array(obstypes)
    obsclasses = np.array(obsclasses)
    return obstypes, obsclasses

def fluxcal(stdsfolder='./', fs=None):
    iraf.cd('work')
    if fs is None:
        fs = glob('x1d/sci*x1d*c*.fits')
    if len(fs) == 0:
        print "WARNING: No science chip spectra to flux calibrate."
        iraf.cd('..')
        return

    if not os.path.exists('flx'):
        os.mkdir('flx')
    extfile = pysaltpath + '/data/site/suth_extinct.dat'
    stdfiles = glob(stdsfolder + '/std/*sens*c?.dat')
    print(stdfiles)
    for f in fs:
        outfile = f.replace('x1d', 'flx')
        chip = outfile[-6]
        hdu = pyfits.open(f)
        ga = f.split('/')[1][3:8]
        # Get the standard sensfunc with the same grating angle
        stdfile = None
        for stdf in stdfiles:
            if ga in stdf:
                # Get the right chip number
                if chip == stdf[-5]:
                    stdfile = stdf
                    break
        if stdfile is None:
            print('No standard star with grating-angle %s' % ga)
            continue
        # for each extracted aperture
        for i in range(hdu[0].data.shape[1]):
            # create an ascii file that pysalt can read
            asciiname = 'flx/sciflx.dat'
            outtmpname = 'flx/scical.dat'
            spectoascii(f, asciiname, i)
            # Run pysalt.speccal
            iraf.unlearn(iraf.speccal)
            iraf.flpr()
            iraf.speccal(asciiname, outtmpname, stdfile, extfile,
                         airmass=pyfits.getval(f, 'AIRMASS'),
                         exptime=pyfits.getval(f, 'EXPTIME'),
                         clobber=True, mode='h')
            # read in the flux calibrated ascii file and copy its
            # contents into a fits file
            flxcal = np.genfromtxt(outtmpname).transpose()
            hdu[0].data[0, i] = flxcal[1]
            hdu[0].data[2, i] = flxcal[2]
            # delete the ascii file
            os.remove(asciiname)
            os.remove(outtmpname)
        hdu.writeto(outfile, clobber=True)
    iraf.cd('..')