def splitkwaj(fn, outPath, outFN, nStart):
    '''Loads a 16-frame FITS file obtained from Kwajalein MIT/LL telescopes - saves a separate FITS files'''
    
    # Basic idea:
    # Kwaj FITS files have one primary header and 15 extensions, each with an
    # associated 1024x1024 image. We load these and save them individually with the
    # associated headers.
    
    hdu1 = pf.open(fn)
    
    # Write the image in the primary HDU to disk
    d = hdu1[0].data
    h = hdu1[0].header
    h.update('EXTEND', 'F')
    
    pf.writeto(outPath + outFN + "_" + str(nStart) + ".fits", d, h)
    
    # Now write the 15 extensions as individual FITS files
    
    fnX = ["01","02","03","04","05","06","07","08","09","10","11","12","13","14","15"]
    
    for i0 in range(15):
        i = i0 + 1	# skip the first image
        d = hdu1[i].data
        h = hdu1[i].header
        k = h.keys()
        hList = arange(43)+9	# indices of the keywords we want to keep
        hdu0 = pf.PrimaryHDU(d)	# Make a new stub of a header
        for j in hList:
            hdu0.header.update(k[j], h[k[j]])
        pf.writeto(outPath + outFN + "_" + str(nStart+i) + ".fits", d, hdu0.header)
    
    return

 def combine_seg_map(self, filt, out_dir):
     """Combines bright and faint segmentation maps. Regions belonging to
     bright objects are expanded by 5 pixels"""
     cat_name = out_dir + '/' + filt + '_clean.cat'
     bright_name = out_dir + '/' + filt + '_bright_seg_map.fits'
     faint_name = out_dir + '/' + filt + '_faint_seg_map.fits'
     hdu1 = pyfits.open(bright_name)
     hdu2 = pyfits.open(faint_name)
     br = hdu1[0].data
     ft = hdu2[0].data
     hdu2.close()
     hdu1.close()
     cat = Table.read(cat_name, format='ascii.basic')
     new_seg = br
     # Expand bright regions by 5 pixels
     q, = np.where(cat['IS_BRIGHT'] == 1)
     for i in q:
         new_seg = fn.seg_expand(new_seg, buff=5, val=int(i) + 1, set_to=int(i) + 1)
         # +1 to account for renumbering
     q, = np.where(cat['IS_BRIGHT'] == 0)
     s = ft.shape
     for i in q:
         for j in range(s[0]):
             pix, = np.where((ft[j, :] == cat['OLD_NUMBER'][i]) & (new_seg[j, :] == 0))
             new_seg[j][pix] = cat['NUMBER'][i] + 1
     new_seg_name = out_dir + '/' + filt + '_comb_seg_map.fits'
     print "Bright faint combined seg map created at", new_seg_name
     pyfits.writeto(new_seg_name, new_seg, clobber=True)
     os.remove(bright_name)
     os.remove(faint_name)

def imcopy(infile, outfile, dim = None):
    print >> sys.stdout, 'Copying ', infile, ' ----> ', outfile
    
    if len(outfile.split('[')) == 1:
        subprocess.call('cp ' + infile + '  ' + outfile, shell = True)
    else:
        if not dim:
            print >> sys.stderr, 'Error : for image section copying, dim parameter cannot be None. Exiting.'
            sys.exit(-1)
            
        header = pyfits.getheader(infile)
        output = numpy.zeros((dim, dim), dtype = numpy.float32)
        
        try:
            f1 = pyfits.open(infile)
        except:
            print >> sys.stderr, 'Error : Not able to open ', infile, '. Exiting.'
            sys.exit(-1)
    
        x1, x2 = int(outfile.split('[')[1].replace(']', '').split(',')[0].split(':')[0] ), int(outfile.split('[')[1].replace(']', '').split(',')[0].split(':')[1])
        y1, y2 = int(outfile.split('[')[1].replace(']', '').split(',')[1].split(':')[0] ), int(outfile.split('[')[1].replace(']', '').split(',')[1].split(':')[1])
        output[x1:x2, y1:y2] = f1[0].data

        outfile = outfile.split('[')[0]
        subprocess.call('rm -f ' + outfile, shell = True)
        pyfits.writeto(outfile, output, header = header)
        
    return outfile

def imshift(filename,shifts,center,refFile,name_ext='.al',clobber=False):
    f = pyfits.open(filename)

    header = f[0].header
    header['REF_FILE'] = (os.path.basename(refFile),'Reference file')
    header['PRE_FILE'] = (os.path.basename(filename),'Filename before shift')
    header['XSHIFT'] = (shifts[0],'X shift from ref_file')
    header['YSHIFT'] = (shifts[1],'Y shift from ref_file')
    header['XCEN'] = (center[0],'X shift from ref_file')
    header['YCEN'] = (center[1],'Y shift from ref_file')
    header['PALIGN'] = (True,'Aligned')

    newName = os.path.splitext(filename)
    newName = ''.join([newName[0],name_ext,newName[1]])

    if shifts[0] != 0 and shifts[1] != 0:
        newDat = shift(f[0].data,(shifts[0],shifts[1]))
    else:
        newDat = f[0].data
        
    print filename
    print '\tShifting (%.2f,%.2f) pixels' % (shifts[0],shifts[1])
    print '\tWriting to %s' % newName
    pyfits.writeto(newName,newDat,header=header,clobber=clobber)

    return newName

def test_exampleimage():
    """Test application of model compared to an independent implementation that was run on the
    example image.
    """
    shiftcoeff = 1.e-7

    #n, r0, t0, rx, tx, r, t, alpha
    cd = galsim.cdmodel.PowerLawCD(
        5, 2. * shiftcoeff, shiftcoeff, 1.25 * shiftcoeff, 1.25 * shiftcoeff, 0.75 * shiftcoeff,
        0.5 * shiftcoeff, 0.3)
    # model used externally to bring cdtest1 to cdtest2
    image_orig  = galsim.fits.read("fits_files/cdtest1.fits") # unprocessed image
    image_proc  = galsim.fits.read("fits_files/cdtest2.fits") # image with cd model applied with
                                                              # other library
    # Calculate the test image
    image_plcd  = cd.applyForward(image_orig)
    # For debugging: make if True in block below to output difference image.
    # Compare to fits_files/cdtest[1-2].fits above
    if False:
        import pyfits
        pyfits.writeto(
            "junk_test_cdmodel_exampleimage_difference.fits", (image_proc - image_plcd).array,
            clobber=True)
    # These images have a large flux per pixel, so make the typical flux per pixel in each image
    # closer to O(1) for a more transparently meaningful decimal order in the test
    norm = 2.64 / np.std(image_orig.array)
    image_proc *= norm
    image_plcd *= norm
    # Compare
    np.testing.assert_array_almost_equal(
        image_proc.array, image_plcd.array, 4, "Externally and internally processed image unequal")

    def copySub(self, other_mapper, dataset, data_id, use_cols, new_template=None):
        """Copy a subset of a particular file in the directory structure defined by this mapper to
        the same location in a directory structure defined by some other mapper.

        @param[in] other_mapper    The mapper defining the directory structure to which the file
                                   should be copied.
        @param[in] dataset         Type of dataset to get; must be one of the keys in self.mappings.
        @param[in] data_id         A dict of values with which to expand the path template
                                   (the first value in self.mappings).
        @param[in] use_cols        A list of columns to copy over (i.e., neglect the others).
        @param[in] new_template    Naming template to use for output catalog, if different from
                                   previous.
        @param[out] outfile        The new file name.
        """
        import numpy
        import pyfits
        # Read in the catalog.
        template, reader, writer = self.mappings[dataset]
        infile = os.path.join(self.full_dir, template % data_id)
        incat = readCatalog(infile)

        # Choose the subset of data to save.
        outcat = numpy.zeros(len(incat),
                             dtype=numpy.dtype(use_cols))
        for col in use_cols:
            outcat[col[0]] = incat[col[0]]

        # Write to output file.
        if new_template is None:
            new_template = template
        outfile = os.path.join(other_mapper.full_dir, new_template % data_id)
        pyfits.writeto(outfile + ".fits", outcat, clobber = True)
        return outfile+'.fits'

def main():
    parser = argparse.ArgumentParser(description='Cross correlate images and return shift necessary to align image2 to image1')
    parser.add_argument('image1',type=str, help='FITS file of image1')
    parser.add_argument('image2',type=str, help='FITS file of image2')
    parser.add_argument('-s',metavar='size',type=int, default=None, help='Specify box size for correlation. Default is the full image, which can be very slow')
    parser.add_argument('-c',metavar=('x_cen', 'y_cen'),type=int,nargs=2, default=None,help="If '-size' specified, optionally include a center for the box region. Default is the center of image1.")
    parser.add_argument('-o',type=str,nargs='?',metavar='outfile',const='-1',default=None,help="If '-o' specified, shift image2 and write to [image2].shft.fits.  If '-o [filename]', shift image2 and write to [filename].")

    args = parser.parse_args()

    print 'Cross-correlating\n\t%s\n\t%s' % (args.image1,args.image2)
    xcorr_im = xcorr(args.image1,args.image2,size=args.s,center=args.c)

    print 'Calculating shift'
    shiftx, shifty = find_shift(xcorr_im)

    print '\t(%i, %i)' % (shiftx, shifty)

    # if outfile specified, perform shift of second image
    if args.o:
        outfile = args.o if args.o != '-1' else \
                  os.path.splitext(args.image2)[0] + '.shft.fits'

        image2, header = pyfits.getdata(args.image2, header=True)
        image2 = shift(image2, (shifty,shiftx), cval=np.nan)
        header['SHFT_REF'] = (args.image1, 'Reference image of shift')
        header['SHFT_X'] = (shiftx, 'X shift pix')
        header['SHFT_Y'] = (shifty, 'Y shift pix')

        print 'Performing shift on %s' % args.image2
        print '\tWriting to %s' % outfile
        pyfits.writeto(outfile,image2,header=header,clobber=True)
    
    return 0

def ccdgap(name):
	fimg = pft.open(name)
	prihdr = fimg[0].header
	scidata = fimg[0].data

	n1 = prihdr['NAXIS1']
	n2 = prihdr['NAXIS2']

	#below are the 4 coordinates of edge of the ccd gap
	a = ccd_locate(scidata)[0]-4;b = ccd_locate(scidata)[1]+2;c = ccd_locate(scidata)[2]-2;d = ccd_locate(scidata)[3]+3
	e = n2 #ccd height

	gap1_part1 = (scidata[:,a-6:a-1].sum(axis=1))/5.0
	gap1_part2 = (scidata[:,b+1:b+6].sum(axis=1))/5.0
	gap2_part1 = (scidata[:,c-6:c-1].sum(axis=1))/5.0
	gap2_part2 = (scidata[:,d+1:d+6].sum(axis=1))/5.0
	
	grad1 = (gap1_part2-gap1_part1)/((b-a)+5.0)
	grad2 = (gap2_part2-gap2_part1)/((d-c)+5.0)

	for i in range(a,b):
		scidata[:,i] = grad1*((i-a)+2)+gap1_part1

	for i in range(c,d):
		scidata[:,i] = grad2*((i-c)+2)+gap2_part1

	namec = "c"+name
	pft.writeto(namec,data=scidata,header=prihdr,clobber=True)
	fimg.close()
	os.system('mv %s history/' % (name))
	return

def main():
    parser = argparse.ArgumentParser(description='Normalize image by exptime. If image already normalized, it is skipped.')
    parser.add_argument('filelist', nargs='+', help='Files to normalize')
    parser.add_argument('-o',metavar='outfile', type=str, help='Specify optional output file, otherwise rewrite file')
    parser.add_argument('-key', type=str, default='EXPTIME', help='Specify exposure time keyword (default=EXPTIME)')
    parser.add_argument('-normkey',type=str, default='NORM', help='Specify normalized keyword (default=NORM)')
    parser.add_argument('--c',action='store_true',help='If specified, force clobber on write')

    args = parser.parse_args()

    for filename in args.filelist:
        data,header = pyfits.getdata(filename, header=True)

        # If already normalized, skip this file
        if is_normalized(header, args.normkey):
            print 'Skipping %s.  Already normalized.' % filename
            continue

        # Else, normalize
        data, header = normalize(data, header, args.key, args.normkey)

        if args.o:
            outname = args.o
        else:
            outname = filename

        print 'Writing to %s' % outname
        pyfits.writeto(outname, data, header=header, clobber=args.c)

def fits_write(sim_dir, filename, image):
    pathname = join(sim_dir, filename)
    header = pyfits.header.Header([('SIMPLE', True), ('NAXIS', 2),
        ('NAXIS1', image.shape[0]), ('NAXIS2', image.shape[1])])
    if (os.path.exists(pathname)):
        os.remove(pathname)
    pyfits.writeto(pathname, image, header)

def cutout(fpC, ra, dec, size, cutout, fpMfn=None, psFieldfn=None, invvarfn=None, band=None):

	wcs = Tan(fpC, 0)
	x,y = wcs.radec2pixelxy(ra, dec)
	x,y = int(x),int(y)
	print 'x,y', x,y
	dl = size / 2
	dh = size - dl
	# ??
	xlo,xhi = max(0, x-dl), min(2048-1, x+dh-1)
	ylo,yhi = max(0, y-dl), min(1489-1, y+dh-1)
	os.system('imcopy %s"[%i:%i,%i:%i]" !%s' %
			  (fpC, xlo, xhi, ylo, yhi, cutout))

	if invvarfn is None:
		   return
		
	bandnum = 'ugriz'.index(band)

	fpc = pyfits.open(fpC)[0].data.astype(float)
	fpM = pyfits.open(fpMfn)
	(gain, darkvar, sky, skyerr) = sdss_psfield_noise(psFieldfn, band=bandnum)

	invvar = sdss_noise_invvar(fpc, fpM, xlo, xhi, ylo, yhi,
							   gain, darkvar, sky, skyerr)
	print invvar.shape
	#print 'x', xlo, xhi
	#print 'y', ylo, yhi
	#invvar = invvar[ylo:yhi, xlo:xhi]
	#print invvar.shape
	pyfits.writeto(invvarfn, invvar, clobber=True)

def make_dead(fact=5.):
    """
    To do in the sorted_by_pos directory in the fringe directory
    Make a boolean mask for data. To work on single_masterflat_image.fits.
    """
    frames = gl.glob("*/*/")
    frames.sort()
    
    length_f = len(frames)
    i_f = 1
    for f in frames:
        print "Computing frame : " + str(i_f) + "/" + str(length_f)
        image = gl.glob(f + "single_masterflat_image.fits")
        if len(image) == 1:
            temp_flat = pf.open(image[0])
            d = temp_flat[0].data
            x_step = 143
            y_step = 128
            x_div = np.linspace(0,4004,29).astype(int)
            y_div = np.linspace(0,4096,33).astype(int)
            mask = np.zeros((4004,4096))
            for x in x_div[:-1]:
                for y in y_div[:-1]:
                    median = np.median(d[x:x+x_step,y:y+y_step])
                    std = np.std(d[x:x+x_step,y:y+y_step])
                    mask[x:x+x_step,y:y+y_step][np.fabs(d[x:x+x_step,y:y+y_step] - median) > fact*std] = 1
            pf.writeto("mask.fits", mask, clobber = True)
            os.system("mv mask.fits " + f)
            temp_flat.close()
        i_f +=1

def modify_binning(field):
    if field in ["fieldA", "fieldB"]:
        return
    seg = pf.getdata("sources.fits")
    white = [x for x in os.listdir(".") if "(white)" in x][0]
    ra = wavelength_array(white, axis=1)
    dec = wavelength_array(white, axis=2)
    # Ofset to the center of NGC 3311
    ra -= ra0
    dec -= dec0
    # Convert to radians
    X = D * 1000 * np.deg2rad(ra)
    Y = D * 1000 * np.deg2rad(dec)
    xx, yy = np.meshgrid(X,Y)
    R = np.sqrt(xx**2 + yy**2)
    base = 10
    Rbins = 10 + 35 * np.logspace(0.3,1,4, base=base) / base
    Rbins = np.hstack((10, Rbins))
    for i,rbin in enumerate(Rbins[:-1]):
        deltar = Rbins[i+1] - Rbins[i]
        newbin = seg.max() + 1
        idxbin = np.where((R > rbin) & (R <= rbin + deltar) & (seg==0))
        if i == 3:
            newbin = 0
        seg[idxbin] = newbin
    pf.writeto("sources.fits", seg, clobber=True)

def extract_spec2():
    infile = workdir + '/maps/standards/HD221246_K3III.fits'
    specInfo, hdr = pyfits.getdata(infile, header=True)
    wave = specInfo[0,:] * 1e3  # in nm
    spec = specInfo[1,:]
    print wave[0:10]
    print wave[-10:]
    print spec

    crpix1 = 1
    crval1 = wave[0]
    cdelt1 = wave[1] - wave[0]
    cunit1 = 'nm'

    tmp = np.arange(len(spec), dtype=float)
    tmp = tmp*cdelt1 + crval1
    print tmp[0:10]
    print tmp[-10:]


    hdr.update('CRPIX1', crpix1)
    hdr.update('CRVAL1', crval1)
    hdr.update('CDELT1', cdelt1)
    hdr.update('CUNIT1', cunit1)

    fitsFile = workdir + 'maps/test_spec_standard.fits'
    ir.imdelete(fitsFile)
    pyfits.writeto(fitsFile, spec, header=hdr)

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 calculate_noise_cube(cube=None, velocity_axis=None,
            velocity_noise_range=[-110,-90,90,110], header=None, Tsys=30.,
            filename=None):

    """ Calcuates noise envelopes for each pixel in a cube
    """

    import numpy as np
    import pyfits as pf
    from mycoords import make_velocity_axis

    if velocity_axis is None:
        velocity_axis = make_velocity_axis(header)


    noise_cube = np.zeros(cube.shape)
    for i in range(cube.shape[1]):
        for j in range(cube.shape[2]):
            profile = cube[:,i,j]
            noise = calculate_noise(profile, velocity_axis,
                    velocity_noise_range)
            #noise = 0.1 # measured in line free region
            noise_cube[:,i,j] = calculate_noise_scale(Tsys,
                    profile, noise=noise)

    if filename is not None:
        pf.writeto(filename, noise_cube, header=header)

    return noise_cube

def invert_image(image, data, header, prefix=None):
    
    ext = fits_ext(image)
    output = prefix + "_negatives.fits" or image.replace(ext,"_negative.fits")
    newdata = -data
    pyfits.writeto(output, newdata, header, clobber=True)
    return output

def makeband(band='V'):
    
    files = glob.glob('Mantis*[0-9]'+band+'_cal.fit*')
    zsz = len(files)
    reffile = files[zsz/2]
    image0,header0 = readimage(reffile)
    ysz,xsz = np.shape(image0)
    
    refim = h.pyfits.open(reffile)
    refh = h.pyfits.getheader(reffile)
    
    stack = np.zeros((xsz,ysz,zsz))
    for i in range(zsz):
       im = h.pyfits.open(files[i])
       newim = h.hcongrid(im[0].data,im[0].header,refh)
       stack[:,:,i] = newim
       
    final = np.median(stack,axis=2)
    
    if band == 'V':
        tag = 'Blue'
        
    if band == 'R':
        tag = 'Green'
        
    if band == 'ip':
        tag = 'Red'
        
    test = glob.glob(tag+'.fit')
    if test:
        os.remove(tag+'.fit')
    pf.writeto(tag+'.fit',final,header0)

def make_images(model='A', brighten=0, bandsel=['u', 'g', 'r', 'i', 'z', 'Y', 'J', 'H', 'K']):
    if noisetype == 'realistic':
        zp = zp_realistic
        sky = sky_realistic
        exp = exp_realistic
    else:
        zp = zp_flat
        sky = sky_flat
        exp = exp_flat
    print 'Using zeropoints:', zp
    print 'Using sky values:', sky
    gals = glob('model%s.galfit'%model)
    for g in gals:
        print g
        os.system('nice galfit %s > %s.out'%(g,g))
        imgname = g.replace('.galfit', '')
        img = pyfits.open(imgname+'.fits')
        for j, b in enumerate(bands):
            if b in bandsel:
                ext = img['MODEL_'+b]
                print b, j, ext.name
                zp_factor = 10**(0.4*(zp[j]-29-fade))
                ext.data *= zp_factor
                brighten_factor = 10**(0.4*brighten)
                if noisetype == 'simple':
                    sigma = numpy.sqrt(sky[j]/exp[j]*brighten_factor)/brighten_factor
                else:
                    sigma = numpy.sqrt((ext.data+sky[j])/exp[j]*brighten_factor)/brighten_factor
                ext.data += numpy.random.normal(0.0, 1.0, sigma.shape) * sigma
                pyfits.writeto(imgname+'_%i%s_%s%i_sigma.fits'%(j+1, b, noisetype[0], brighten), sigma, clobber=True)
                pyfits.writeto(imgname+'_%i%s_%s%i.fits'%(j+1, b, noisetype[0], brighten), ext.data, clobber=True)

def split(filename, outdir, prefix):
    basename, ext = os.path.splitext(filename)
    basename = os.path.basename(basename)
    basename = os.path.join(outdir,basename)

    f = pyfits.open(filename)

    ydim = f[0].header['NAXIS2']/2

    Adat = f[0].data[ydim:,:]   #top
    Bdat = f[0].data[0:ydim,:]  #bot

    Afile = ''.join([basename,'_A',ext])
    Bfile = ''.join([basename,'_B',ext])

    hdr = f[0].header
    hdr['WOLLY'] = ('A','Image half')
    hdr['FILENUM'] = (get_filenum(filename,prefix),'Observation number')
    
    pyfits.writeto(Afile,Adat,header=hdr,clobber=True)

    hdr['WOLLY'] = ('B','Image half')
    pyfits.writeto(Bfile,Bdat,header=hdr,clobber=True)

    return (Afile, Bfile)