def shiftRGB(redF,greenF,blueF,blueshiftr=0,blueshiftc=0,greenshiftr=0,greenshiftc=0,redshiftr=0,redshiftc=0,ext=None):
    """
    this code shift the pixels of three r, g, b images. Using g image as reference and shift the other two images. It will return the shifted r,g,b images.
    each row goes along ra direction
    each col goes along dec direction
    CRVAL1 ; ra direction
    CRVAL2: dec direction
    """
    blueHdr = pf.getheader(blueF,ext)
    greenHdr = pf.getheader(greenF,ext)
    redHdr = pf.getheader(redF,ext)
    bluerow = blueHdr['crval1']*3600./0.27
    bluecol = blueHdr['crval2']*3600./0.27
    greenrow = greenHdr['crval1']*3600./0.27
    greencol = greenHdr['crval2']*3600./0.27
    redrow = redHdr['crval1']*3600./0.27
    redcol = redHdr['crval2']*3600./0.27
    """
    col0=int(blueHdr['datasec'].split('[')[1].split(']')[0].split(',')[0].split(':')[0])-1
    col1=int(blueHdr['datasec'].split('[')[1].split(']')[0].split(',')[0].split(':')[1]) 
    row0=int(blueHdr['datasec'].split('[')[1].split(']')[0].split(',')[1].split(':')[0])-1
    row1=int(blueHdr['datasec'].split('[')[1].split(']')[0].split(',')[1].split(':')[1]) 
    """
    blue = pf.getdata(blueF,ext)
    green = pf.getdata(greenF,ext)
    red = pf.getdata(redF,ext)

    ctgreenrow = (bluerow+greenrow+redrow)/3.
    ctgreencol = (bluecol+greencol+redcol)/3.
    blue = nd.shift(blue,[bluerow - ctgreenrow+blueshiftr,bluecol-ctgreencol+blueshiftc],mode='nearest',order=1)
    green = nd.shift(green,[greenrow - ctgreenrow+greenshiftr,greencol-ctgreencol+greenshiftc],mode='nearest',order=1)
    red = nd.shift(red,[redrow - ctgreenrow+redshiftr, redcol-ctgreencol+redshiftc],mode='nearest',order=1)
    return red,green,blue

def imgray2fits(infile, fitsfile='', overwrite=False, headerfile=None, flip=False):
    if fitsfile == '':
        fitsfile = decapfile(infile) + '.fits'

    if exists(fitsfile):
        if overwrite:
            delfile(fitsfile)
        else:
            print fitsfile, 'EXISTS'
            sys.exit(1)
    
    data = loadgray(infile)  # coeim.py
    
    #hdu = pyfits.PrimaryHDU()
    header = headerfile and pyfits.getheader(headerfile)
    hdu = pyfits.PrimaryHDU(None, header)
    hdulist = pyfits.HDUList([hdu])
    hdulist.writeto(fitsfile)

    try:  # If there's a 'SCI' extension, then that's where the WCS is
        header = pyfits.getheader(headerfile, 'SCI')
    except:
        pass
    
    if header <> None:
        if 'EXTNAME' in header.keys():
            del(header['EXTNAME'])
    
    if flip:
        data = flipud(data)
    
    pyfits.append(fitsfile, data, header)
    
    print fitsfile, 'PRODUCED'

def collapse_cube(w1, w2):
    """ Collapse a MUSE data cube.

    Arguments

    cube : MUSE data cube name containing both data and stat extensions.
    iext : Initial extension to be used. Default is one for combined cubes.

    """
    fits = "slice_w{0}_{1}.fits".format(w1, w2)
    outfits = "collapsed_w{0}_{1}.fits".format(w1, w2)
    data = pf.getdata(fits, 0)
    error = pf.getdata(fits, 1)
    h = pf.getheader(fits, 0)
    h2 = pf.getheader(fits, 1)
    h["NAXIS"] = 2
    del h["NAXIS3"]
    h2["NAXIS"] = 2
    del h2["NAXIS3"]
    print "Starting collapsing process..."
    start = time.time()
    w = wavelength_array(fits)
    # newdata = np.trapz(data, dx=np.diff(w)[0], axis=0)
    # newdata = np.nansum(data, axis=0) * np.diff(w)[0]
    newdata = np.nanmedian(data, axis=0)
    noise = 1.482602 / np.sqrt(6.) * np.nanmedian(np.abs(2.* data - \
           np.roll(data, 2, axis=0) - np.roll(data, -2, axis=0)), \
           axis=0)
    end = time.time()
    print "Collapsing lasted {0} minutes.".format((end - start)/60.)
    hdu = pf.PrimaryHDU(newdata, h)
    hdu2 = pf.ImageHDU(noise, h2)
    hdulist = pf.HDUList([hdu, hdu2])
    hdulist.writeto(outfits, clobber=True)
    return

    def do_test_3_fqpm_tilt(self):
        """ Test FQPM tilting (no FQPM yet), no field mask. Verify proper behavior in Lyot plane"""

        osys = poppy.OpticalSystem("test", oversample=self.oversample)
        osys.addPupil('Circle', radius=6.5/2)
        osys.addPupil('FQPM_FFT_aligner')
        osys.addImage()  # perfect image plane
        osys.addPupil('FQPM_FFT_aligner', direction='backward')
        osys.addPupil('Circle', radius=6.5/2)
        osys.addDetector(pixelscale=self.pixelscale, fov_arcsec=3.0)
            #TODO testing of odd and even focal plane sizes?
        
        plt.clf()
        poppy._FLUXCHECK=True
        psf = osys.calcPSF(wavelength=self.wavelength, save_intermediates=True, display_intermediates=True)
        psf.writeto('test3a_psf.fits', clobber=True)

        # after the Lyot plane, the wavefront should be all real. 
        check_wavefront('wavefront_plane_004.fits', test='is_real', comment='(Lyot Plane)')

        cen = webbpsf.measure_centroid('wavefront_plane_002.fits', boxsize=50)
        head = pyfits.getheader('wavefront_plane_002.fits')
        desired_pos = (head['NAXIS1']-1)/2.0
        self.assertAlmostEqual( cen[0], desired_pos, delta=0.025) #within 1/50th of a pixel of desired pos?
        self.assertAlmostEqual( cen[1], desired_pos, delta=0.025) #within 1/50th of a pixel of desired pos?
                # This is likely dominated by uncertainties in the simple center measuring algorithm...

        _log.info("FQPM FFT half-pixel tilting is working properly in intermediate image plane")

        cen2 = webbpsf.measure_centroid('wavefront_plane_005.fits', boxsize=50)
        head2 = pyfits.getheader('wavefront_plane_005.fits')
        desired_pos2 = (head2['NAXIS1']-1)/2.0
        self.assertAlmostEqual( cen2[0], desired_pos2, delta=0.05) #within 1/20th of a pixel of desired pos?
                                    
        _log.info("FQPM FFT half-pixel tilting is working properly in final image plane")

def test_pv_invert_many(find=True):
    import es_util
    from glob import glob
    sep = '-'*70

    dir=os.path.expanduser("~/data/astrometry")
    imdir=os.path.join(dir,'image')
    plotdir=os.path.join(dir,'plots')
    testdir=os.path.join(dir,'test')

    x = numpy.array([21.34, 1000.0, 1500.17], dtype='f8')
    y = numpy.array([11.21, 1000.0, 1113.92], dtype='f8')

    pattern = os.path.join(imdir,'*')
    images=glob(pattern)


    hdr=pyfits.getheader(images[0])
    n1=hdr['naxis1']
    n2=hdr['naxis2']
    xrang=numpy.array([1.0, n1], dtype='f8')
    yrang=numpy.array([1.0, n2], dtype='f8')
    #n=100
    n=10
    x,y = es_util.make_xy_grid(n, xrang, yrang)

    rms = []
    xdiffs = []
    ydiffs = []
    for imname in images:
        sys.stdout.write('image: %s\n' % imname)
        sys.stdout.write(sep+'\n')
        hdr=pyfits.getheader(imname)
        wcs=wcsutil.WCS(hdr)

        sys.stdout.write('doing xforms\n')
        ra,dec = wcs.image2sky(x,y)
        xn,yn = wcs.sky2image(ra,dec, find=find)

        xdiff = xn-x
        ydiff = yn-y

        t=(xdiff)**2 + (ydiff)**2
        trms = numpy.sqrt( t.sum()/t.size )
        rms.append(trms)
        xdiffs.append(xdiff)
        ydiffs.append(ydiff)

        sys.stdout.write('rms: %s\n' % trms)

    rms = numpy.array(rms,dtype='f8')
    out={}
    out['rms'] = rms
    out['x'] = x
    out['y'] = y
    out['xdifflist'] = xdiffs
    out['ydifflist'] = ydiffs
    return out

 def use_correlation(self):
     """
     Use correlation data-cube.
     """
     import numpy
     
     from pyfits import getdata, getheader, writeto
     from glob import glob
     from os.path import splitext
     from sys import stdout
     
     self.print("\n A correlation cube will be used.")
     self.print(" Looking for an existing correlation data-cube in the current folder.")
     candidates = glob("*.fits")
                     
     corr_cube = None
     for candidate in candidates:
         if 'CORRFROM' in getheader(candidate):
             if getheader(candidate)['CORRFROM'] == self.input_file:
                 self.print(" Correlation cube to be used: %s" % candidate)
                 return candidate
         
     if corr_cube == None:
         self.print(" Correlation cube not found. Creating a new one.")
         data = getdata(self.input_file)
         corr_cube = numpy.empty_like(data)
         
         x = numpy.arange(self.width)
         y = numpy.arange(self.height)
         X, Y = numpy.meshgrid(x, y)
         x, y = numpy.ravel(X), numpy.ravel(Y)
             
         for i in range(x.size):
             s = data[:,y[i],x[i]]
             s = s / s.max()  # Normalize
             s = s - s.mean() # Remove mean to avoid triangular shape
             s = numpy.correlate(s, self.ref_s, mode='same')
             corr_cube[:,y[i],x[i]] = s
                 
             temp = (((i + 1) * 100.00 / X.size))
             stdout.write('\r %2d%% ' % temp)
             stdout.write(self.loading[int(temp * 10 % 5)])
             stdout.flush()
         
         self.print(" Done.")
         corr_name = splitext(self.input_file)[0] + '--corrcube.fits'
         self.print(" Saving correlation cube to %s" % corr_name)
             
         corr_hdr = self.header.copy()
         corr_hdr.set('CORRFROM', self.input_file,'Cube used for corrcube.')
         corr_hdr.set('', '', before='CORRFROM')
         corr_hdr.set('', '--- Correlation cube ---', before='CORRFROM')
         
         writeto(corr_name, corr_cube, corr_hdr, clobber=True)
         del corr_hdr
         del corr_cube
         
         return corr_name

def time_steps(path, num_exp):
    import pyfits as pf

    label = [0] * num_exp
    label[0] = '0 minutes'
    for i in range(num_exp - 1):
        label[i + 1] = str(round((pf.getheader(path, 4 + i + 1)['tai-beg'] -
                                  pf.getheader(path, 4)['tai-beg']) / 60)) + ' minutes'
    return label

def data_collect(file_directory, file_prefix, file_postfix, first_image, num_images, num_parameters):

    # Extracts data and header from first image block, gets dimensions image block
    # Determines the number of frames total for the file range

    x = p.getdata(file_directory + file_prefix + str(first_image) + file_postfix)

    h = p.getheader(file_directory + file_prefix + str(first_image) + file_postfix)

    cube_shape = x.shape

    num_frames = num_images * cube_shape[0]

   
    # Preallocates the 3d data cube
 
    data_cube = np.zeros((num_frames, cube_shape[1], cube_shape[2]))


    # Preallocates the 2D parameter array: (number of frames,  number of 
    # parameters)

    photometry = [None]*((num_images) * cube_shape[0])  

    for i in np.arange(num_frames):
        photometry[i] = [None] * num_parameters

 
    # Loop for each image cube, extracts frame data to data cube and updates
    # photometry parameters.

    image_num = 0

    for image in np.arange(num_images):

        data = p.getdata(file_directory + file_prefix + str(first_image + image) + file_postfix)

        hdr = p.getheader(file_directory + file_prefix + str(first_image + image) + file_postfix)
    
        # Reads in time data from header
        time     = hdr['UTCS_OBS']
        interval = hdr['FRAMTIME'] 
    
        for i in np.arange(cube_shape[0]):
         
            # Reads in each frame in the image to the data cube
            data_cube[image_num + i] = data[i]
     
            # Writes image number to photometry parameter 0 for each frame
            photometry[image_num + i][0] = first_image + image_num
       
            # Writes time the frame was taken to photometry parameter 1
            photometry[image_num + i][1] = (interval * i) + time
        
        image_num += cube_shape[0]

    return(photometry, data_cube, num_frames, cube_shape)

def imsubtract( image1, image2, outfile=None, 
                clobber=False, verbose=False, debug=False):
    """
    Construct a simple subtraction: image2 - image1.  Guards against
    different sized data arrays by assuming that the lower left pixel
    (0,0) is the anchor point.  (i.e. the second image will be
    trimmed or extended if needed.)  
    """
    import os
    import pyfits
    from numpy import ndarray
    import exceptions

    if debug  : import pdb; pdb.set_trace()

    if outfile : 
        if os.path.isfile( outfile ) and not clobber : 
            print("%s exists. Not clobbering."%outfile)
            return( outfile )

    # read in the images
    if not os.path.isfile( image1 ) :
        raise exceptions.RuntimeError(
            "The image file %s is not valid."%image1 )
    im1head = pyfits.getheader( image1 )
    im1data = pyfits.getdata( image1 )

    if not os.path.isfile( image2 ) :
        raise exceptions.RuntimeError(
            "The image file %s is not valid."%image2 )
    im2head = pyfits.getheader( image2 )
    im2data = pyfits.getdata( image2 )

    # sometimes multidrizzle drops a pixel. Unpredictable.
    nx2,ny2 = im2data.shape
    nx1,ny1 = im1data.shape
    if nx2>nx1 or ny2>ny1 : 
        im2data = im2data[:min(nx1,nx2),:min(ny1,ny2)]
        im1data = im1data[:min(nx1,nx2),:min(ny1,ny2)]
    elif nx2<nx1 or ny2<ny1 : 
        im1data = im1data[:min(nx1,nx2),:min(ny1,ny2)]
        im2data = im2data[:min(nx1,nx2),:min(ny1,ny2)]

    diffim =  im2data - im1data
    
    if not outfile :
        return( diffim )
    else : 
        im2head.update("SRCIM1",image1,"First source image = template for subtraction")
        im2head.update("SRCIM2",image2,"Second source image = search epoch image")
        outdir = os.path.split( outfile )[0]
        if outdir and not os.path.isdir(outdir): 
            os.makedirs( outdir )
        pyfits.writeto( outfile, diffim, 
                        header=im2head,
                        clobber=clobber )
        return( outfile )

def imarith(file1, operator, file2, output, verbose=True, keep=1):
    """
    @param file1:       Input fits image name 1.
    @param operator:    Operator ['+','-','*','/'].
    @param file2:       Input fits image name 2.
    @param output:      Output image.
    @keyword verbose:    Turn verbose on? [True]/False.
    @keyword keep:      Choose whether keep header from [file1] or file2.
    """
    from pyfits import getdata, getheader, writeto;
    from time import strftime;
    
    # Reading data
    data1 = get_fits_data(file1);
    data2 = get_fits_data(file2);
    header = getheader(file2) if keep is 2 else getheader(file1);
    
    # Applying operation
    if verbose: print(" %s %s %s = %s" % (file1, operator, file2, output));
    
    
    if (data1.ndim is 3) and (data2.ndim is 2):
        data2 = data2.reshape((1, data2.shape[0], data2.shape[1]));
        data2 = data2.repeat(data1.shape[0], axis=0);
    if (data1.ndim is 3) and (data2.ndim is 1):
        data2 = data2.reshape((data2.shape[0], 1, 1));
        data2 = data2.repeat(data1.shape[1], axis=1);
        data2 = data2.repeat(data1.shape[2], axis=2);
    
    assert data1.shape == data2.shape; 
    if operator is '+': data1 = data1 + data2;
    elif operator is '-': data1 = data1 - data2;
    elif operator is '*': data1 = data1 * data2;
    elif operator is '/': data1 = data1 / data2;
    else: raise (ValueError, "Invalid operator.");
    
    # Updating header
    try:
        header.set('', '');
    except AttributeError:
        from sys import exit;
        from pyfits import __version__
        print(" Sorry, your PyFits version %s is not supported."% __version__);
        print(" Please, download PyFits v3.0 or superior." );
        print(" Leaving now.\n " );
        exit();
        
    header.set('HISTORY', 'imarith applied on %s' % strftime('%Y-%m-%d %H:%M:%S %Z'));
    header.set('HISTORY', file1, 'First file used on imarith.py.');
    header.set('HISTORY', operator, 'Operator used on imarith.py.');
    header.set('HISTORY', file2, 'Second file used on imarith.py.');
    
    # Writing output
    output = safesave(output);
    writeto(output, data1, header);
    
    return None;

 def __init__(self, im, spec, verbose=False, instrument="hydra"):
     self.image = im
     self.id = im.replace('.fits', '')
     self.spec = spec
     self.verbose = verbose        
     self.data = pf.getdata(spec)
     self.header = pf.getheader(spec)
     self.imheader = pf.getheader(im)
     self.shape_info()
     self.check_fields()
     if instrument == "hydra": 
         self.reverse_spec()
     return

def get_spec_parameters(spec_file):
    '''
    This function reads in a spectrum and the primary and first extension headers
    Input:
        spec_file: the name of the file
    Output:
        spec: binary fits table
        hdr0: primary header
        hdr1: first extension header
    '''
    spec = pyfits.getdata(spec_file, 1)
    hdr0 = pyfits.getheader(spec_file, 0)
    hdr1 = pyfits.getheader(spec_file, 1)
    return spec,hdr0, hdr1

def get_sismo_data(ID):
    """
    Retrieve CoRoT timeseries from a local data repository.
    
    The output record array has fields 'HJD', 'flux', 'e_flux', 'flag'.
    
    @param ID: ID of the target: either an integer (CoRoT ID), an SIMBAD-recognised
    target name, or a valid CoRoT FITS file
    @type ID: int or str
    @return: data, header
    @rtype: numpy recarray, dict
    """
    #-- data on one target can be spread over multiple files: collect the
    #   data
    data = []
        
    if isinstance(ID,str) and os.path.isfile(ID):
        header = pyfits.getheader(ID)
        times,flux,error,flags = fits.read_corot(ID)
        data.append([times,flux,error,flags])
    else:
        #-- resolve the target's name: it's either a target name or CoRoT ID.
        try:
            ID = int(ID)
        except ValueError:
            info = sesame.search(ID,db='S')
            IDs = [alias for alias in info['alias'] if 'HD' in alias]
            if len(IDs)!=1:
                logger.error("Data retrieval for %s not possible. Reason: no HD number resolved" % (ID))
                return
            ID = IDs[0]
        #-- collect the files containing data on the target
        catfiles = config.glob((os.sep).join(['catalogs','corot','sismo']),'*.fits')
        for catfile in catfiles:
            try:
                header = pyfits.getheader(catfile)
            except IOError:
                continue
            if header['starname']==ID or header['corotid'].replace(' ','')=='%s'%(ID):
                times,flux,error,flags = fits.read_corot(catfile)
                data.append([times,flux,error,flags])
    #-- now make a record array and sort according to times
    if not data:
        raise ValueError('target {0} not in offline CoRoT data repository'.format(ID))
    data = np.hstack(data)        
    data = np.rec.fromarrays(data,dtype=[('HJD','>f8'),('flux','>f8'),('e_flux','>f8'),('flag','i')])
    sa = np.argsort(data['HJD'])
    return data[sa],header

    def Lee_cubo(spectra,XX,YY):
        global imagen
        imagen=pyfits.getdata(spectra,header=False)
        header=pyfits.getheader(spectra)
        
        #print len(imagen)
        #empty array
        Lambda_t=[]
        Flux_t=[]
    
    
        for i in range (len(imagen)):
            y=imagen[i][XX][YY]
            #        x=i*header['CDELT1']+header['CRVAL1']
            x=i*header['CD3_3']+header['CRVAL3']
            Lambda_t.append(float(imagen[i][XX][YY]))
            #Flux_t.append(float(i*header['CDELT1']+header['CRVAL1']))
            Flux_t.append(float(i*header['CD3_3']+header['CRVAL3']))
            #print x,y

        Flux=np.array(Lambda_t)
        Lambda=np.array(Flux_t)
        x=Lambda
        y=Flux
        return x,y

def make_wifes_p08_template(ddir, fn, out_dir, star,rv=0.0):
    """From a p08 file, create a template spectrum for future cross-correlation.
    The template is interpolated onto a 0.1 Angstrom grid (to match higher resolution 
    templates.
    
    Parameters
    ----------
    ddir: string
        Data directory for the p08 file
        
    fn: string
        p08 fits filename
        
    out_dir: string
        Output directory
    
    """
    flux_stamp,wave = read_and_find_star_p08(ddir + '/' + fn)
    heliocentric_correction = pyfits.getheader(ddir + '/' + fn)['RADVEL']
    spectrum,sig = weighted_extract_spectrum(flux_stamp)
    dell_template = 0.1
    wave_template=np.arange(90000)*dell_template + 3000
    spectrum_interp = np.interp(wave_template,wave*(1 - (rv - heliocentric_correction)/2.998e5),spectrum)
    outfn = out_dir + '/' + star + ':' + fn
    pyfits.writeto(outfn,spectrum_interp,clobber=True)

def generate_display_id( dataset, version ):
    '''
    Generate an ID from which all similar stackable data will have in common.
    
    @param dataset: Input AstroData or fits filename
    @type dataset: list of AstroData instance
    
    @param version: The version from which to run.
    @type version: string
    
    @return: A display id.
    @rtype: string  
    '''
    if version != version_index['display_id']:
        try:
            # designed to call generateStackableID_
            idFunc = getattr( globals()['IDFactory'], 'generateDisplayID_' + version )
        except:
            raise "Version: '" + version + "' is either invalid or not supported." 
        
        return idFunc( inputf, version )
    
    """
        shaObj = hashlib.sha1()
        phu = pf.getheader( inputf[0], 0 )
        shaObj.update( phu['OBSID'] )
        shaObj.update( phu['OBJECT'] )
    """
    
    if type(dataset) == str:
        phu = pf.getheader(dataset)
        ID = version + "_" + phu['OBSID'] + "_" + phu['OBJECT'] 
    elif type(dataset) == AstroData:
        ID = version + "_" + dataset.phuValue('OBSID') + "_" + dataset.phuValue('OBJECT')
    return ID

def getvsopname(gemfile, instrument, path):
    # Defaults
    validInstruments = ('GMOSN', 'GMOSS')
    instDict = {'N': 'GMOSN',
                'S': 'GMOSS'}

    # Verify instrument
    if instrument == None:
        if matchInst.match(gemfile) == None:
            errmsg = 'Invalid Gemini frame name, \''+gemfile+'\''
            raise IOError, errmsg
        instpref = matchInst.sub(r"\1",gemfile)
        instrument = instDict[instpref]
    elif instrument not in validInstruments:
        print '\nUSAGE ERROR: Invalid value for instrument (%s)' % instrument
        p.print_help()
        raise SystemExit
    
    phu = pyfits.getheader(path+'/'+gemfile, 0)
    dateobs = phu['DATE-OBS']
    timeobs = phu['TIME-OBS']
    
    (year,month,date) = map(string.atoi, dateobs.split('-'))
    d = datetime.date(year,month,date)
    (hour,minutes,seconds) = map(string.atof, timeobs.split(':'))
    microseconds = round((seconds - int(seconds)) * 1e6)
    if int(microseconds) == 0:  microseconds += 1
    t = datetime.time(int(hour),int(minutes),int(seconds),int(microseconds))
    dt = datetime.datetime.combine(d,t)

    return instrument+'.'+dt.isoformat()[:-3]+'_s1d.fits'

def guide_star_info():
    """
    Get the guide star information out of the image headers.
    """
    imgDir = ngc1851_data + 'starfinder/'

    imgs = glob.glob(imgDir + '*_G1_1_cr.fits')

    # Keep results in an output file
    f_out = open(ngc1851_data + 'guide_star_info.txt', 'w')
    f_fmt = '{0:8s} {1:15s} '

    for img in imgs:
        hdr = pyfits.getheader(img)

        f_out.write('{0:20s}  '.format(os.path.basename(img)))

        # There are 4 possible guide stars. Try to print them all out
        for ii in range(4):
            try:
                cfg = 'GWFS{0}CFG'.format(ii+1)
                obj = 'GWFS{0}OBJ'.format(ii+1)
                f_out.write(f_fmt.format(hdr[cfg], hdr[obj]))
            except KeyError:
                f_out.write(f_fmt.format('-', '-'))
                
        f_out.write('\n')

    f_out.close()

def image_info():
    """
    Get the exposure time and filter info for each image. Store it in
    an output file called image_info.txt.
    """
    imgDir = ngc1851_data + 'starfinder/'

    imgs = glob.glob(imgDir + '*_G1_1_cr.fits')

    # Keep results in an output file
    f_out = open(ngc1851_data + 'image_info.txt', 'w')
    f_fmt = '{img:15s}  {xoff:7.2f} {yoff:7.2f} {pa:5.1f}  {filt:12s}  '
    f_fmt += '{exp:6.2f} {coad:2d} {tot:7.2f}\n'

    for img in imgs:
        hdr = pyfits.getheader(img)

        xoff = hdr['XOFFSET']
        yoff = hdr['YOFFSET']
        pa = hdr['PA']
        filt = hdr['FILTER1']
        exp = hdr['EXPTIME']
        coadd = hdr['COADDS']

        f_out.write(f_fmt.format(img=os.path.basename(img), xoff=xoff, yoff=yoff,
                                 pa=pa, filt=filt, exp=exp, coad=coadd, tot=exp*coadd))

    f_out.close()

def rdpsfmodel(psfmodelfile):
    """ Read in a psf model from a fits file. Gaussian parameters are in the
    header, and the image array has a lookup table of non-gaussian components
    sub-sampled to a half-pixel grid.
    If the user provides a 2-tuple instead of a filename, then we presume the
    user already has the psf model components, so we just return it back.

    :param psfmodelfile: a fits file containing the psf model
    :return: [gaussparam, lookuptable]
    """
    psfmodelfile = psfmodelfile
    if isinstance(psfmodelfile, str):
        assert os.path.isfile(
            os.path.abspath(os.path.expanduser(psfmodelfile)))
        # read in the psf non-gaussian components array (i.e. the lookup table)
        lookuptable = pyfits.getdata(psfmodelfile)
        # read in the gaussian parameters from the image header
        hdr = pyfits.getheader(psfmodelfile)
        scale = hdr['GAUSS1']  # , 'Gaussian Scale Factor'
        xpsf = hdr['GAUSS2']  # , 'Gaussian X Position'
        ypsf = hdr['GAUSS3']  # , 'Gaussian Y Position'
        xsigma = hdr['GAUSS4']  # , 'Gaussian Sigma: X Direction'
        ysigma = hdr['GAUSS5']  # , 'Gaussian Sigma: Y Direction'
        psfmag = hdr['PSFMAG']  # , 'aperture magnitude of the PSF star
        psfzpt = hdr['PSFZPT']  # , 'zeropoint used to set PSF star mag scaling
        gaussparam = [scale, xpsf, ypsf, xsigma, ysigma]
    elif np.iterable(psfmodelfile):
        assert len(psfmodelfile) == 4
        gaussparam, lookuptable, psfmag, psfzpt = psfmodelfile
    else:
        raise exceptions.RuntimeError(
            "psfmodel must either be a filename or a 4-tuple giving:"
            "[gaussian parameters, look-up table, psf mag, zpt]")
    return gaussparam, lookuptable, psfmag, psfzpt