def bootstrapZP(image,cat='zp.cat',refcat='ref.cat',refzp=1.0):

    (ra,dec,mag,star)=np.loadtxt(cat,usecols=(0,1,2,3),unpack=True,
                                     dtype=np.dtype([('ra','<f10'),('dec','<f10'),('mag','<f10'),
                                                     ('star','<f10')]))

    (rra,rdec,rmag,rstar)=np.loadtxt(refcat,usecols=(0,1,2,3),unpack=True,
                                            dtype=np.dtype([('rra','<f10'),('rdec','<f10'),('rmag','<f10'),
                                                            ('rstar','<f10')]))

    #Grab only stars from the reference catalog
    refgood=np.where((rstar >= 0.98) & (rmag != 99.0) & (rmag > 17.0) & (rmag < 22.5))
    refcat=SkyCoord(ra=rra[refgood]*u.degree,dec=rdec[refgood]*u.degree)

    #Sort through and remove anything that is not a star and is not isolated
    #from other sources in the input catalog
    catgood=np.where((star >= 0.98) & (mag != 99.0))
    cat=SkyCoord(ra=ra[catgood]*u.degree,dec=dec[catgood]*u.degree)
    idx,d2d,_=cat.match_to_catalog_sky(refcat)
    _,d2d2,_=cat.match_to_catalog_sky(cat,2)
    final=np.where((d2d.arcsec < 0.5) & (d2d2.arcsec >= 5.0))
    
    diff=rmag[refgood][idx][final]-mag[catgood][final]
    imgZP=np.mean(diff)

    print '\n\tUsing '+str(len(diff))+' stars to calculate ZP...'
    print '\tMean ZP: '+str(round(imgZP,3))+' mag\n'

    scaleFactor=10.0**(0.4*(refzp-imgZP))

    iraf.unlearn('imcalc')
    iraf.imcalc(image,image[:-5]+'_scaled.fits','im1*'+str(scaleFactor))

def tpv2tan_hdr(img, ota):
    image = odi.reprojpath+'reproj_'+ota+'.'+img.stem()
    # change the CTYPENs to be TANs if they aren't already
    tqdm.write('TPV -> TAN in {:s}'.format(image))
    iraf.imutil.hedit.setParam('images',image)
    iraf.imutil.hedit.setParam('fields','CTYPE1')
    iraf.imutil.hedit.setParam('value','RA---TAN')
    iraf.imutil.hedit.setParam('add','yes')
    iraf.imutil.hedit.setParam('addonly','no')
    iraf.imutil.hedit.setParam('verify','no')
    iraf.imutil.hedit.setParam('update','yes')
    iraf.imutil.hedit(show='no', mode='h')

    iraf.imutil.hedit.setParam('images',image)
    iraf.imutil.hedit.setParam('fields','CTYPE2')
    iraf.imutil.hedit.setParam('value','DEC--TAN')
    iraf.imutil.hedit.setParam('add','yes')
    iraf.imutil.hedit.setParam('addonly','no')
    iraf.imutil.hedit.setParam('verify','no')
    iraf.imutil.hedit.setParam('update','yes')
    iraf.imutil.hedit(show='no', mode='h')

    # delete any PV keywords
    # leaving them in will give you trouble with the img wcs
    iraf.unlearn(iraf.imutil.hedit)
    iraf.imutil.hedit.setParam('images',image)
    iraf.imutil.hedit.setParam('fields','PV*')
    iraf.imutil.hedit.setParam('delete','yes')
    iraf.imutil.hedit.setParam('verify','no')
    iraf.imutil.hedit.setParam('update','yes')
    iraf.imutil.hedit(show='no', mode='h')

    def execute(self):
        """Execute pyraf task: gsappwave"""
        log.debug("GsappwaveETI.execute()")

        # Populate object lists
        xcldict = copy(self.clparam_dict)
        for fil in self.file_objs:
            xcldict.update(fil.get_parameter())
        for par in self.param_objs:
            xcldict.update(par.get_parameter())
        iraf.unlearn(iraf.gmos.gsappwave)

        # Use setParam to list the parameters in the logfile
        for par in xcldict:
            # Stderr and Stdout are not recognized by setParam
            if par != "Stderr" and par != "Stdout":
                gemini.gmos.gsappwave.setParam(par, xcldict[par])
        log.fullinfo("\nGSAPPWAVE PARAMETERS:\n")
        iraf.lpar(iraf.gmos.gsappwave, Stderr=xcldict["Stderr"], Stdout=xcldict["Stdout"])

        # Execute the task using the same dict as setParam
        # (but this time with Stderr and Stdout)
        # from pprint import pprint
        # pprint(xcldict)
        gemini.gmos.gsappwave(**xcldict)
        if gemini.gmos.gsappwave.status:
            raise Errors.OutputError("The IRAF task gmos.gsappwave failed")
        else:
            log.fullinfo("The IRAF task gmos.gsappwave completed successfully")

def fixpix(fs=None):
    iraf.cd('work')
    if fs is None:
        fs = glob('nrm/sci*nrm*.fits')
    if len(fs) == 0:
        print "WARNING: No rectified images to fix."
        iraf.cd('..')
        return
    if not os.path.exists('fix'):
        os.mkdir('fix')
    for f in fs:
        outname = f.replace('nrm', 'fix')
        # Copy the file to the fix directory
        shutil.copy(f, outname)
        # Set all of the BPM pixels = 0
        h = pyfits.open(outname, mode='update')
        h['SCI'].data[h['BPM'].data == 1] = 0
        # Grab the CRM extension from the lax file
        laxhdu = pyfits.open(f.replace('nrm', 'lax'))
        h.append(pyfits.ImageHDU(data=laxhdu['CRM'].data.copy(),
                                 header=laxhdu['CRM'].header.copy(),
                                 name='CRM'))
        h.flush()
        h.close()
        laxhdu.close()

        # Run iraf's fixpix on the cosmic rays, not ideal,
        # but better than nothing because apall doesn't take a bad pixel mask
        iraf.unlearn(iraf.fixpix)
        iraf.flpr()
        iraf.fixpix(outname + '[SCI]', outname + '[CRM]', mode='hl')
    iraf.cd('..')

    def execute(self):
        """Execute pyraf task: gmosaic"""
        log.debug("GmosaicETI.execute()")

        # Populate object lists
        xcldict = copy(self.clparam_dict)
        for fil in self.file_objs:
            xcldict.update(fil.get_parameter())
        for par in self.param_objs:
            xcldict.update(par.get_parameter())
        iraf.unlearn(iraf.gmos.gmosaic)

        # Use setParam to list the parameters in the logfile
        for par in xcldict:
            #Stderr and Stdout are not recognized by setParam
            if par != "Stderr" and par != "Stdout":
                gemini.gmos.gmosaic.setParam(par, xcldict[par])
        log.fullinfo("\nGMOSAIC PARAMETERS:\n")
        iraf.lpar(iraf.gmos.gmosaic, Stderr=xcldict["Stderr"], \
            Stdout=xcldict["Stdout"])

        # Execute the task using the same dict as setParam
        # (but this time with Stderr and Stdout)
        #from pprint import pprint
        #pprint(xcldict)
        try:
            gemini.gmos.gmosaic(**xcldict)
        except:
            # catch hard crash
            raise RuntimeError("The IRAF task gmos.gmosaic failed")
        if gemini.gmos.gmosaic.status:
            # catch graceful exit upon error
            raise RuntimeError("The IRAF task gmos.gmosaic failed")
        else:
            log.fullinfo("The IRAF task gmos.gmosaic completed successfully")

def mkflatINTWFC(data, combine='median', reject='avsigclip'):
    """
    Creates a master flat from a given list of flat frames.

    Reads readnoise and gain from the FITS header, as each WFC chip has
    different values.

    :note: ccdproc
    """
    for filter in data:
        input1 = ''
        for file in data[filter]:
            input1 += file + '[1],'

        input2 = input1.replace('[1]', '[2]')
        input3 = input1.replace('[1]', '[3]')
        input4 = input1.replace('[1]', '[4]')

        inputs = [input1, input2, input3, input4]

        for ext, input in enumerate(inputs):
            print input
            iraf.unlearn('flatcombine')
            iraf.flatcombine(input=input[:-1],
                             output='FLAT_{0:>s}_{1:d}.fits'.format(filter, ext + 1),
                             combine=combine,
                             reject=reject,
                             rdnoise='READNOIS',
                             gain='GAIN')

    def execute(self):
        """Execute pyraf task: gemcombine"""
        log.debug("GemcombineETI.execute()")

        # Populate object lists
        xcldict = copy(self.clparam_dict)
        for fil in self.file_objs:
            xcldict.update(fil.get_parameter())
        for par in self.param_objs:
            xcldict.update(par.get_parameter())
        iraf.unlearn(iraf.gemcombine)

        # Use setParam to list the parameters in the logfile 
        for par in xcldict:
            #Stderr and Stdout are not recognized by setParam
            if par != "Stderr" and par !="Stdout":
                gemini.gemcombine.setParam(par,xcldict[par])
        log.fullinfo("\nGEMCOMBINE PARAMETERS:\n")
        iraf.lpar(iraf.gemcombine, Stderr=xcldict["Stderr"], \
            Stdout=xcldict["Stdout"])

        # Execute the task using the same dict as setParam
        # (but this time with Stderr and Stdout)
        try:
            gemini.gemcombine(**xcldict)
        except:
            # catch hard crash of the primitive
            raise Errors.OutputError("The IRAF task gemcombine failed")
        if gemini.gemcombine.status:
            # catch graceful exit on error
            raise Errors.OutputError("The IRAF task gemcombine failed")
        else:
            log.fullinfo("The IRAF task gemcombine completed successfully")

 def subtract_dark(self):
     """Subtracting Dark Frame"""
     iraf.unlearn(iraf.ccdproc)
     iraf.ccdproc(self.data.iraf.modatfile(), 
         ccdtype="", fixpix="no", overscan="no", trim ="no", zerocor="no", darkcor="yes", flatcor ="no", 
         dark=self.dark.iin("Dark"))
     self.data.idone()

def trim_img(img,x1,x2,y1,y2):
    """
    Trim a stacked image based on the coordinates given. The image is trimmed
    using ``imcopy`` through pyraf, so the x and y pixel ranges should be given
    in the correct ``imcopy`` format. ``[x1:x2,y1:y2]``

    Parameters
    ---------
    img : str
        String containing name of the image currently in use
    x1 : int
        Pixel coordinate of x1
    x2 : int
        Pixel coordinate of x2
    y1 : int
        Pixel coordinate of y1
    y2 : int
        Pixel coordinate of y2

    Returns
    -------
    img : str
        The new image is given the extension ``.trim.fits``.

    """
    x1,x2 = x1,x2
    y1,y2 = y1,y2
    input = img.nofits()+'['+repr(x1)+':'+repr(x2)+','+repr(y1)+':'+repr(y2)+']'
    output =  img.nofits()+'.trim.fits'
    if not os.path.isfile(output):
        print 'Trimming image: ' ,img
        iraf.unlearn(iraf.imcopy)
        iraf.imcopy(input = input,output = output,verbose='no',mode='h')

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 make_bpms(images):
  for i in range(len(images)):
    for key in OTA_dictionary:
      mask,masked_array = mask_ota(images[i],OTA_dictionary[key])
      hdu = fits.PrimaryHDU(mask.astype(float))
      # create string for mask fits name
      mask_name = bppath+'mask_'+OTA_dictionary[key]+'.'+str(images[i][-27:49])
      BPM = mask_name.replace('fits','pl')
      if not os.path.isfile(mask_name):
        hdu.writeto(mask_name,clobber=True)
      if not os.path.isfile(mask_name.replace('fits','pl')):
        iraf.unlearn(iraf.imutil.imcopy)
        iraf.imutil.imcopy.setParam('input',mask_name)
        iraf.imutil.imcopy.setParam('output',mask_name.replace('fits','pl'))
        iraf.imutil.imcopy.setParam('verbose','no')
        iraf.imutil.imcopy(mode='h')
      iraf.unlearn(iraf.imutil.hedit)
      iraf.imutil.hedit.setParam('images',images[i]+'['+str(key)+']')
      iraf.imutil.hedit.setParam('fields','BPM')
      iraf.imutil.hedit.setParam('value',BPM)
      iraf.imutil.hedit.setParam('add','no')
      iraf.imutil.hedit.setParam('addonly','no')
      iraf.imutil.hedit.setParam('verify','no')
      iraf.imutil.hedit.setParam('update','yes')
      iraf.imutil.hedit(mode='h')
  return

def makedark(files, output):
    """
    Make dark image for NIRC2 data. Makes a calib/ directory
    and stores all output there. All output and temporary files
    will be created in a darks/ subdirectory. 

    files: integer list of the files. Does not require padded zeros.
    output: output file name. Include the .fits extension.
    """
    redDir = os.getcwd() + '/'  # Reduce directory.
    curDir = redDir + 'calib/'
    darkDir = util.trimdir(curDir + 'darks/')
    rawDir = util.trimdir(os.path.abspath(redDir + '../raw') + '/')

    util.mkdir(curDir)
    util.mkdir(darkDir)
    
    _out = darkDir + output
    _outlis = darkDir + 'dark.lis'
    util.rmall([_out, _outlis])

    darks = [rawDir + 'n' + str(i).zfill(4) + '.fits' for i in files]

    f_on = open(_outlis, 'w')
    f_on.write('\n'.join(darks) + '\n')
    f_on.close()
    
    ir.unlearn('imcombine')
    ir.imcombine.combine = 'median'
    ir.imcombine.reject = 'sigclip'
    ir.imcombine.nlow = 1
    ir.imcombine.nhigh = 1
    ir.imcombine('@' + _outlis, _out)

def reproject_ota(img, ota, rad, decd):
    from pyraf import iraf
    image = odi.illcorpath+'illcor_'+ota+'.'+str(img[16:])
    imout = odi.reprojpath+'reproj_'+ota+'.'+str(img[16:])
    iraf.mscred(_doprint=0)
    iraf.clobber='no'
    iraf.unlearn(iraf.mscred.mscimage)
    iraf.mscred.mscimage.format='image'
    iraf.mscred.mscimage.pixmask='yes'
    iraf.mscred.mscimage.verbose='yes'
    iraf.mscred.mscimage.wcssour='parameters'
    iraf.mscred.mscimage.ref=''
    iraf.mscred.mscimage.ra=rad
    iraf.mscred.mscimage.dec=decd
    iraf.mscred.mscimage.scale=0.11
    iraf.mscred.mscimage.rotation=0.0
    iraf.mscred.mscimage.blank=-999
    iraf.mscred.mscimage.interpo='poly5'
    iraf.mscred.mscimage.minterp='poly5'
    iraf.mscred.mscimage.nxbl=4096
    iraf.mscred.mscimage.nybl=4096
    iraf.mscred.mscimage.fluxcon='yes'
    iraf.mscred.mscimage(image,imout)
    
    return

def stack_images(refimg):
    from astropy.io import fits
    from pyraf import iraf
    print refimg
    fitsref = fits.open(refimg)
    hduref = fitsref[0]
    objname = hduref.header['object']
    filter_name = hduref.header['filter']
    sky_med = odi.find_new_bg(refimg, filter_name)
    # sky_med = hduref.header['skybg']
    output = objname+'_'+filter_name+'.fits'
    output_bpm = objname+'_'+filter_name+'_bpm.pl'
    iraf.unlearn(iraf.immatch.imcombine, iraf.imutil.imarith)
    iraf.immatch.imcombine(odi.scaledpath+'*'+filter_name+'*.fits', 'temp', combine='average', reject='none', offsets='wcs', masktype='goodvalue', maskval=0, blank=-999, scale='none', zero='none', lthresh=-900, hthresh=60000)
    # iraf.imutil.imarith.setParam('operand1','temp')
    # iraf.imutil.imarith.setParam('op','+')
    # iraf.imutil.imarith.setParam('operand2',sky_med)
    # iraf.imutil.imarith.setParam('result',output)
    # iraf.imutil.imarith.setParam('verbose','yes')
    # iraf.imutil.imarith(mode='h')
    iraf.imutil.imexpr('(a != -999) ? a + b : -999',output,'temp.fits',sky_med)
    iraf.imutil.imexpr('a < 0',output_bpm, output)
    iraf.imutil.imdelete('temp', verify='no')
    iraf.unlearn(iraf.imutil.hedit)
    iraf.imutil.hedit.setParam('images',output)
    iraf.imutil.hedit.setParam('fields','BPM')
    iraf.imutil.hedit.setParam('value',output_bpm)
    iraf.imutil.hedit.setParam('add','yes')
    iraf.imutil.hedit.setParam('addonly','no')
    iraf.imutil.hedit.setParam('verify','no')
    iraf.imutil.hedit.setParam('update','yes')
    iraf.imutil.hedit(show='no', mode='h')
    # iraf.immatch.imcombine(reprojpath+'*.fits', 'test', expm='exp.pl', combine='average', reject='none', offsets='wcs', masktype='goodvalue', maskval=0, blank=-999, scale='none', zero='none', lthresh=-900, hthresh=60000)  

    return output

def rotate_2005_lgsao():
    # Need to rotate an image.
    ir.unlearn('rotate')
    ir.rotate.boundary = 'constant'
    ir.rotate.constant = 0
    ir.rotate.interpolant = 'spline3'
    ir.rotate.ncols = 1040
    ir.rotate.nlines = 1040
    
    ir.rotate.xin = 528
    ir.rotate.xout = 540
    ir.rotate.yin = 645
    ir.rotate.yout = 410
    ir.rotate(workdir + 'mag05jullgs_kp.fits', workdir + 'mag05jullgs_kp_rot.fits', 190)

    ir.rotate.xin = 535
    ir.rotate.xout = 540
    ir.rotate.yin = 655
    ir.rotate.yout = 410
    ir.rotate(workdir + 'mag05jullgs_h.fits', workdir + 'mag05jullgs_h_rot.fits', 190)

    ir.rotate.xin = 572
    ir.rotate.xout = 540
    ir.rotate.yin = 694
    ir.rotate.yout = 410
    ir.rotate(workdir + 'mag05jullgs_lp.fits', workdir + 'mag05jullgs_lp_rot.fits', 190)

    ir.rotate.xin = 483
    ir.rotate.xout = 483
    ir.rotate.yin = 612
    ir.rotate.yout = 612
    ir.rotate(workdir + 'mag04jul.fits', workdir + 'mag04jul_rot.fits', 1)

def setup_drizzle(imgsize):
    """Setup drizzle parameters for NIRC2 data.
    @param imgsize: The size (in pixels) of the final drizzle image.
    This assumes that the image will be square.
    @type imgsize: int
    @param type: str
    """
    print 'Setting up drizzle'
    # Setup the drizzle parameters we will use
    ir.module.load('stsdas', doprint=0, hush=1)
    ir.module.load('analysis', doprint=0, hush=1)
    ir.module.load('dither', doprint=0, hush=1)
    ir.unlearn('drizzle')
    ir.drizzle.outweig = ''
    ir.drizzle.in_mask = ''
    ir.drizzle.wt_scl = 1
    ir.drizzle.outnx = imgsize
    ir.drizzle.outny = imgsize
    ir.drizzle.pixfrac = 1
    ir.drizzle.kernel = 'lanczos3'
    ir.drizzle.scale = 1
    ir.drizzle.coeffs = distCoef
    ir.drizzle.xgeoim = distXgeoim
    ir.drizzle.ygeoim = distYgeoim
    ir.drizzle.shft_un = 'input'
    ir.drizzle.shft_fr = 'output'
    ir.drizzle.align = 'center'
    ir.drizzle.expkey = 'ITIME'
    ir.drizzle.in_un = 'counts'
    ir.drizzle.out_un = 'counts'

 def subtract_bias(self):
     """Subtracting Bias Frame"""
     iraf.unlearn(iraf.ccdproc)
     iraf.ccdproc(self.data.iraf.modatfile(), 
         ccdtype="", fixpix="no", overscan="no", trim ="no", zerocor="yes", darkcor="no", flatcor ="no", 
         zero=self.bias.iin("Bias"))
     self.data.idone()

def scireduce(scifiles, rawpath):
    for f in scifiles:
        setupname = getsetupname(f)
        # gsreduce subtracts bias and mosaics detectors
        iraf.unlearn(iraf.gsreduce)
        iraf.gsreduce('@' + f, outimages=f[:-4]+'.mef', rawpath=rawpath, bias="bias",
                      fl_over=dooverscan, fl_fixpix='no', fl_flat=False,
                      fl_gmosaic=False, fl_cut=False, fl_gsappwave=False, fl_oversize=False)

        if is_GS:
            # Renormalize the chips to remove the discrete jump in the
            # sensitivity due to differences in the QE for different chips
            iraf.unlearn(iraf.gqecorr)
            iraf.gqecorr(f[:-4]+'.mef', outimages=f[:-4]+'.qe.fits', fl_keep=True, fl_correct=True,
                         refimages=setupname + '.arc.arc.fits',
                         corrimages=setupname +'.qe.fits', verbose=True)

            iraf.unlearn(iraf.gmosaic)
            iraf.gmosaic(f[:-4]+'.qe.fits', outimages=f[:-4] +'.fits')
        else:
            iraf.unlearn(iraf.gmosaic)
            iraf.gmosaic(f[:-4]+'.mef.fits', outimages=f[:-4] +'.fits')

        # Flat field the image
        hdu = pyfits.open(f[:-4]+'.fits', mode='update')
        hdu['SCI'].data /= pyfits.getdata(setupname+'.flat.fits', extname='SCI')
        hdu.flush()
        hdu.close()

        # Transform the data based on the arc  wavelength solution
        iraf.unlearn(iraf.gstransform)
        iraf.gstransform(f[:-4], wavtran=setupname + '.arc')

def mkbiasINTWFC(filelist, type='median'):
    """
    Creates a master bias from given list of bias frames.
    Saves each extension to a different FITS file.

    Reads readnoise and gain from the FITS header, as each WFC chip has
    different values.

    :note: This function has been written specifically for INT WFC.

    :param filelist:
    :type filelist:
    :param type: combining parameter (default = median)
    :type type: string

    :return: None
    """
    input1 = ''
    for line in filelist:
        input1 += line[0] + '[1],'
    input2 = input1.replace('[1]', '[2]')
    input3 = input1.replace('[1]', '[3]')
    input4 = input1.replace('[1]', '[4]')
    inputs = [input1, input2, input3, input4]

    #note there are four SCI extensions
    for ext, input in enumerate(inputs):
        iraf.unlearn('zerocombine')
        iraf.zerocombine(input=input[:-1],
                         output='BIAS%i.fits' % (ext+1),
                         combine=type,
                         rdnoise='READNOIS',
                         gain='GAIN')

def make_bpms(img, ota):
    # for i in range(len(images)):
    #   for key in OTA_dictionary:
    # create string for mask fits name
    mask_name = odi.bppath+'mask_'+ota+'.'+str(img[16:])
    BPM = mask_name.replace('fits','pl')
    if not os.path.isfile(BPM):
        mask,gaps = odi.mask_ota(img,ota)
        hdu = odi.fits.PrimaryHDU(mask.astype(float))
        if not os.path.isfile(mask_name):
            hdu.writeto(mask_name,clobber=True)
        #if not os.path.isfile(mask_name.replace('fits','pl')):
            iraf.unlearn(iraf.imutil.imcopy)
            iraf.imutil.imcopy.setParam('input',mask_name)
            iraf.imutil.imcopy.setParam('output',mask_name.replace('fits','pl'))
            iraf.imutil.imcopy.setParam('verbose','no')
            iraf.imutil.imcopy(mode='h')
    iraf.unlearn(iraf.imutil.hedit)
    iraf.imutil.hedit.setParam('images',img+'['+ota+']')
    iraf.imutil.hedit.setParam('fields','BPM')
    iraf.imutil.hedit.setParam('value',BPM)
    iraf.imutil.hedit.setParam('add','yes')
    iraf.imutil.hedit.setParam('addonly','no')
    iraf.imutil.hedit.setParam('verify','no')
    iraf.imutil.hedit.setParam('update','yes')
    iraf.imutil.hedit(show='no', mode='h')
    if os.path.isfile(mask_name):
        os.remove(mask_name)
    return