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 fit_elli(clus_id, line_s):
    #try:
    size = c.size
    values = line_s.split()
    mask_file = 'ell_mask_' + str(imagefile)[:6] + '_'  + str(clus_id) + '.fits'
    image_file = 'image_' + str(imagefile)[:6] + '_'  + str(clus_id) + '.fits' 
    print image_file
    xcntr_o  = float(values[1]) #x center of the object
    ycntr_o  = float(values[2]) #y center of the object
    xcntr = (size/2.0) + 1.0 + xcntr_o - int(xcntr_o)
    ycntr = (size/2.0) + 1.0 + ycntr_o - int(ycntr_o)
    mag    = float(values[7]) #Magnitude
    radius = float(values[9]) #Half light radius
    mag_zero = 25.256 #magnitude zero point
    sky	 = float(values[10]) #sky
    if(float(values[11])>=0 and float(values[11])<=180.0): 
        pos_ang = float(values[11]) - 90.0 #position angle
    if(float(values[11])<0 and float(values[11])>=-180.0):
        pos_ang = 90.0 - abs(float(values[11]))  #position angle
    if(float(values[11])>180 and float(values[11])<=360.0):
        pos_ang = float(values[11]) - 360.0 + 90.0 #position angle
    if(float(values[11])>=-360 and float(values[11])<-180.0):
        pos_ang = float(values[11]) + 360.0 - 90.0 #position angle	
    axis_rat = 1.0 / float(values[12]) #axis ration b/a
    eg = 1 - axis_rat
    if(eg<=0.05):
        eg = 0.07
    major_axis = float(values[14])#major axis of the object
    iraf.imcopy(mask_file, 'image'+str(mask_file)[8:]+'.pl')
    run_elli(image_file, xcntr, ycntr, eg, pos_ang, major_axis)

 def test_imcopy(self):
     iraf.imcopy('dev$pix', 'image.short', verbose=False)
     with fits.open('image.short.fits') as f:
         assert len(f) == 1
         assert f[0].header['BITPIX'] == 16
         assert (f[0].header['ORIGIN'] ==
                 'NOAO-IRAF FITS Image Kernel July 2003')
         assert f[0].data.shape == (512, 512)

def custom1(filename): # for NACO timing mode cubes - removes horizontal banding
    #iraf.imarith(filename,'-','dark','temp')
    iraf.imarith(filename,'/','flatK','temp')
    im = pyfits.getdata('temp.fits')
    med = median(im.transpose())
    out = ((im).transpose()-med).transpose()
    (pyfits.ImageHDU(out)).writeto("temp2.fits",clobber=True)
    iraf.imdel('temp')
    iraf.imcopy('temp2[1]','temp')

def trim_img(img):
    x1,x2 = 2508,15798
    y1,y2 = 2216,15506
    input = img[:-5]+'['+repr(x1)+':'+repr(x2)+','+repr(y1)+':'+repr(y2)+']'
    output =  img[:-5]+'.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 shiftTrim(files, xi, xf, yi, yf, dx=None, dy=None, offset_x=0, offset_y=0, 
    fnew='', newdir=TRIMDIR, test=False):
    """
    Given initial and final x and y values of shifted stars, will compute shift in x and y and 
    trim files to compensate for shifting image.

    Warning:
    - Assumes that movement across CCD is uniform and predictable in x and y (can be zero)
    - Discard bad images after trim is complete (requires complete series to trim accurately)
    - Only use reduced science images 
    - Pixels are discretely counted

    Parameters:
    ----------------------
    parameter: (dtype) [default (if optional)], information

    xi: (int), initial x value(s) for star(s)
    xf: (int), final x values(s) for star(s)
    yi: (int), initial y value(s) for star(s)
    yx: (int), final y value(s) for star(s)
    dx: (int) [None], overwrite x shift value (will ignore xi,xf)
    dy: (int) [None], overwrite y shift value (will ignore yi,yf)
    fnew: (string) [None], add string to new file name
    newdir: (string) [TRIMDIR], trimmed image directory
    test: (boolean) [False], if True will only print pixel output and not trim files
    ----------------------
    """
    xi = np.array(xi)
    xf = np.array(xf)
    yi = np.array(yi)
    yf = np.array(yf)
    
    num = len(files) - 1
    
    if dx == None:
        dx = _dCalc(xi,xf,num,'x')
        
    if dy == None:
        dy = _dCalc(yi,yf,num,'y')
    
    print 'dx: %s, dy: %s' % (dx,dy)
    for i,f in enumerate(files):

        x = _pixelFinder(dx, i, abs(dx*num), 1, CCDx, offset_x)
        y = _pixelFinder(dy, i, abs(dy*num), 1, CCDy, offset_y)
        
        print "x: %s, y: %s, x*y: %s" % (x, y, ((x[1]-x[0])*(y[1]-y[0])))
        
        f_trim = f + '[%s:%s,%s:%s]' % (x[0],x[1],y[0],y[1])

        if fnew != None:
            f = f + '.' + fnew

        if test == False:
            iraf.imcopy(f_trim, newdir + f)

def deimos_preproc(image):

    '''Take a MEF DEIMOS image, extract the relevant extensions, trim the
    LVM slit masks appropriately, and rename.'''

    # Needing grating
    graname = get_head(image, "GRATENAM", extn=0)
    
    # "Blue" chip
    iraf.imcopy("%s[%i]" % (image, BEXT), "d%s_B.fits" % image[6:10])
    iraf.ccdproc("d%s_B.fits" % image[6:10], overscan=yes, trim=yes,
                 fixpix=yes, biassec=BBIASSEC, trimsec=BTRIMSEC,
                 fixfile="%s_%i.fits" % (MASK, BEXT))
    iraf.imcopy("d%s_B.fits%s" % (image[6:10], DGRATINGS[graname]["blvmreg"]),
                "td%s_B.fits" % image[6:10])
    iraf.imtranspose("td%s_B.fits" % image[6:10], "rtd%s_B.fits" % image[6:10])

    # "Red" chip
    iraf.imcopy("%s[%i]" % (image, REXT), "d%s_R.fits" % image[6:10])
    iraf.ccdproc("d%s_R.fits" % image[6:10], overscan=yes, trim=yes,
                 fixpix=yes, biassec=RBIASSEC, trimsec=RTRIMSEC,
                 fixfile="%s_%i.fits" % (MASK, BEXT))
    iraf.imcopy("d%s_R.fits%s" % (image[6:10], DGRATINGS[graname]["rlvmreg"]),
                "td%s_R.fits" % image[6:10])
    iraf.imtranspose("td%s_R.fits" % image[6:10], "rtd%s_R.fits" % image[6:10])
    iraf.rotate("rtd%s_R.fits" % image[6:10], "rtd%s_R.fits" % image[6:10], 180.0)

    return

def _iraf_dqbits_init(_data):
    """Initialize common IRAF tasks for dqbits tests
    """
    if not HAS_IRAF:
        return

    # imports & package loading
    iraf.stsdas(_doprint=0)
    iraf.imgtools(_doprint=0)
    iraf.artdata(_doprint=0)
    iraf.mstools(_doprint=0)

    # create two data files as input (dont care if appropriate to mscombine)
    iraf.imcopy('dev$pix', _data['dqbits']['input1'])
    iraf.imcopy('dev$pix', _data['dqbits']['input2'])

def HandleEllipseTask(cutimage, xcntr, ycntr, SizeX, SizeY, sky, out):
    """Running the ellipse task. SizeX, SizeY are the total size"""
    manual_profile = 0
    try:
        raise ImportError() #Temporarily kill this loop as the new flagging does not work with pyraf-functions, yet
        from pyraf import iraf
        from fitellifunc import run_elli
        use_pyraf = 1
    except ImportError:
        use_pyraf = 0
        print 'No pyraf installed!'
        WriteError('Cannot find pyraf installation! Trying manual 1d ' + \
                   'profile finder\n')
    if use_pyraf:
        if out:
            ell_mask_file = 'OEM_' + c.fstring + '.fits'
            ell_out = 'OE_' + c.fstring + '.txt'
        else:
            ell_mask_file = 'EM_' + c.fstring + '.fits'
            ell_out = 'E_' + c.fstring + '.txt'
        plfile = 'GalEllFit.fits.pl'
        CleanEllipse(ell_out, 0)
        try:
            iraf.imcopy(ell_mask_file, plfile, verbose='no')
            iraf.flpr()
        except:
            pass
        try:
            run_elli(cutimage, ell_out, xcntr, ycntr, c.eg, \
                     c.pos_ang, c.major_axis, sky)
            CleanEllipse(ell_out, 1)
            try:
                iraf.flpr()
            except:
                pass
            if exists(ell_out):
                pass
            else:
                manual_profile = 1
        except:
            manual_profile = 1
            WriteError('Error in ellipse task. Trying manual profile finder\n')
            try:
                c.Flag = SetFlag(c.Flag, GetFlag('ELLIPSE_FAIL'))
            except badflag:
                pass
    if use_pyraf == 0 or manual_profile:
        FitEllipseManual(cutimage, xcntr, ycntr, SizeX, SizeY, sky, out)

def zeropadfits(smfits, bigfits, padfits):
    """Pads smfits with zeros to match size of bigfits.
       Result is padfits, centered as was smfits.
       Assumes smfits & bigfits are squares w/ odd # of pixels across.
    """
    NY, NX = fits.getheader(bigfits)["NAXIS2"], fits.getheader(bigfits)["NAXIS1"]
    ny, nx = fits.getheader(smfits)["NAXIS2"], fits.getheader(smfits)["NAXIS1"]
    print "\nPadding 'smfits' at %ix%i to match 'bigfits' at %ix%i\n" % (nx, ny, NX, NY)
    center = (NY + 1) / 2
    border = ny / 2
    lo = center - border
    hi = center + border
    croprange = "[%d:%d,%d:%d]" % (lo, hi, lo, hi)

    imarith(bigfits, "*", 0, padfits)
    imcopy(smfits, padfits + croprange)

	def trimMyself(self, outname="remcut.fits", region="[400:500,400:500]", verbose=False):
		
                self._logger["trimMyself"] = []

		if self._Name == "None":
			self._logger["trimMyself"].append("No filename set, please check.")

		imCopy = imFits()
		imCopy._Name = outname
		iraf.imcopy("%s%s" % (self._Name,region), imCopy._Name)
		
		if verbose:
			self._logger["trimMyself"].append(iraf.imstat(self._Name))	

			for log in self._logger["trimMyself"]:
				print log

		return imCopy

def lris_bpreproc(image):

    '''Take an LRIS blue image, combine into a single extension, subtract overscan,
    trim, and rotate.'''

    iraf.keck()
    iraf.lris()
    
    # Need grating name to start
    grism = get_head("../rawdata/%s" % image, "GRISNAME", extn=0)

    # Convert to single extension fits image
    iraf.multi2simple("../rawdata/%s" % image, "b%s.fits" % image[8:12],
                      overscan=yes, header=yes, trim=yes, verbose=no, debug=no)

    # Trim and rotate
    iraf.imcopy("b%s[%s]" % (image[8:12], LGRISMS[grism]["trimreg"]),
                "tb%s" % image[8:12])
    iraf.imtranspose("tb%s" % image[8:12], "rtb%s" % image[8:12])

    return

def split_frames(cubefile, range_pairs, target_dir):
    dirname, filename = os.path.split(cubefile)
    filebase = filename.rsplit('.', 1)[0]
    if len(filebase) > 8:
        warn("IRAF doesn't like long filenames. "
             "Consider shortening the cube filename ({0})".format(filebase))
    
    outfiles = []
    
    for fromidx, toidx in range_pairs:
        for i in range(fromidx, toidx+1):
            infile = cubefile + "[*,*,{0}]".format(i)
            outfile = '{0}/frame_{1:04}.fit'.format(target_dir, i)
            debug("imcopy", infile, outfile)
            iraf.imcopy( # easier to use imcopy and preserve headers than to use pyfits I think
                input=infile,
                output=outfile
            )
            outfiles.append(outfile)
        
        # f = open(inlst, 'w')
        # f.writelines(infiles)
        # f.write('\n')
        # f.close()
        
        # outfile = '{0}/{1}_{2}-{3}.fit'.format(target_dir, filebase, fromidx, toidx)
        # debug("imcombine input={input} output={output} combine=sum reject=none".format(
        #     input="@{0}".format(inlst), #','.join(infiles),
        #     output=outfile,
        # ))
        # outfiles.append(outfile)
        # iraf.imcombine(
        #     input=','.join(infiles),
        #     output=outfile,
        #     combine="sum",
        #     reject="none",
        #     # project='no', # IRAF wants bools specified / default is nonboolean?
        #     # mclip='no',
        # )
    return outfiles

def sn_marshall():
    """

    """
    from pyraf import iraf

    import unicorn.go_3dhst as go
    import threedhst.process_grism as proc
    import unicorn.analysis

    os.chdir(unicorn.GRISM_HOME+'SN-MARSHALL')

    #### Copy necessary files from PREP_FLT to DATA
    os.chdir('PREP_FLT')
    grism_asn  = glob.glob('MARSHALL-2??-G141_asn.fits')
    files=glob.glob('MARSHALL-2*-G141_shifts.txt')
    files.extend(grism_asn)
    for file in files:
        status = os.system('cp '+file+' ../DATA')
        #shutil.copy(file, '../DATA')

    try:
        iraf.imcopy('MARSHALL-F125W_drz.fits[1]', '../DATA/f125w_sci.fits')
    except:
        os.remove('../DATA/f125w_sci.fits')
        iraf.imcopy('MARSHALL-F125W_drz.fits[1]', '../DATA/f125w_sci.fits')
    
    os.chdir('../')

    #### Initialize parameters
    go.set_parameters(direct='F160W', LIMITING_MAGNITUDE=26)
    threedhst.options['PREFAB_DIRECT_IMAGE'] = '../PREP_FLT/MARSHALL-F160W_drz.fits'
    threedhst.options['OTHER_BANDS'] = [['f125w_sci.fits', 'F125W' , 1248.6, 26.25]]

    #### Main loop for reduction
    for i in range(len(grism_asn)):
        asn = grism_asn[i]
        proc.reduction_script(asn_grism_file=asn)
        unicorn.analysis.make_SED_plots(grism_root=asn.split('_asn.fits')[0])
        go.clean_up()

def GetSkyMedian(gain,trim):
	
	# for this and subtracted images there is a new way to get the sky level.
	# this involves measuring the stddev in the subtracted images and then squaring it 
	# and multiplying by the gain to get the flux per pix which can then be used in the
	# sky calculations, see below to do so.
	
	# first imcopy the subtracted images into a new file
	# this is to stopp the truncation problem which is common. 
	
	t=cmd.getoutput('ls /Volumes/DATA/NITES/data/2011-08-12/reduced/CheckSkyInSubtracted/s_M71*.fits').split('\n')
	
	if trim == 1:
		for i in range(0,len(t)):
			image=str(t[i]+"[1:500,1:500]")
			image2=str(t[i])

			# clobber the old image
			iraf.imcopy(input=image,output=image2)
	
	# then define a box in each image (same box for bright and dark time) and 
	# get the standard deviation in it for each frame that night. Then get the average 
	# nightly stddev which is used to get the sky level.
	
	std_list=np.empty(len(t))
	
	for i in range(0,len(t)):
		h=pf.open(t[i])
		d=h[0].data[243:313,274:344]
		
		std_list[i]=np.std(d)
		#print "%s std: %.2f" % (t[i],std_list[i])
		
	av=pow(np.average(std_list),2)*gain
	print "Sky_bg_pp: %.2f e-" % (av)
	
	
	return av

def apply_rectangle_mask(image, mask, threshold, min_width, min_height):
    #finds rectangular middle part with good exposure and trims images to that part.
    hdulist = pyfits.open(image)
    imdata = hdulist[0].data
    #print len(imdata)
    hdulist.close()
    hdulist = pyfits.open(mask)
    immask = hdulist[0].data
    #print len(immask)
    hdulist.close()
    
    good_row = np.zeros(immask.shape[1])
    good_col = np.zeros(immask.shape[0])
    for i in np.arange(0, immask.shape[1]):
        good_col[i] = len(np.where(immask[i,:] > threshold*immask.max())[0])
    
    #print good_row
    for i in np.arange(0, immask.shape[0]):
        good_row[i] = len(np.where(immask[:,i] > threshold*immask.max())[0])
    
    # if sufficient exposure exists (i.e. not a bad frame), trim the image and the mask accordingly:
    if (len(np.where(good_row > min_width)[0]) > min_height) & (len(np.where(good_col > min_height)[0]) > min_width):
        good_width = 0.8 * max(good_row)
        good_height = 0.8 * max(good_col)
	ind_row = [np.where(good_row >= good_width)[0].min(), np.where(good_row >= good_width)[0].max()]
	ind_col = [np.where(good_col >= good_height)[0].min(), np.where(good_col >= good_height)[0].max()]
	
	print ind_row
	print ind_col
	
	outdata = deepcopy(imdata[ind_row[0]:ind_row[1], ind_col[0]:ind_col[1]])
	outmask = deepcopy(immask[ind_row[0]:ind_row[1], ind_col[0]:ind_col[1]])
	
	outimage = image.replace('.fits', '_trimmed.fits')
	outmask = mask.replace('.fits', '_trimmed.fits')
	
	iraf.imcopy(image + '[' + str(ind_row[0]+1) + ':' + str(ind_row[1]+1) + ',' + str(ind_col[0]+1) + ':' + str(ind_col[1]+1) + ']', outimage)
	iraf.imcopy(mask + '[' + str(ind_row[0]+1) + ':' + str(ind_row[1]+1) + ',' + str(ind_col[0]+1) + ':' + str(ind_col[1]+1) + ']', outmask)

def slice_rc(img):
    '''
    Slices the Rainbow Camera into 4 different images and adds the 'filter' keyword in the fits file.
    '''
    fname = os.path.basename(img)
    fdir = os.path.dirname(img)    
    
    # Running IRAF
    iraf.noao(_doprint=0)
    
    '''f = pf.open(img)
    h = f[0].header
    i = f[0].data[1000:-1,0:920]
    g = f[0].data[0:910,0:900]
    r = f[0].data[1035:2046, 1050:2046]
    u = f[0].data[0:910,1050:2046]'''
    
    
    corners = {
    "i" : [1, 910, 1, 900],
    "g" : [1, 910, 1060, 2045],
    "r" : [1040, 2045, 1015, 2045],
    "u" : [1030, 2045, 1, 900]
    }
    
    
    filenames = []
        
    for i, b in enumerate(corners.keys()):
        name = fname.replace(".fits", "_%s.fits"%b)
        iraf.imcopy("%s[%d:%d,%d:%d]"%(img, corners[b][0], corners[b][1], corners[b][2], corners[b][3]), name)
         
        fitsutils.update_par(name, 'filter', b)

        filenames.append(name)
    
    return filenames

def sn_primo():
    from pyraf import iraf

    import unicorn.go_3dhst as go
    import threedhst.process_grism as proc
    import unicorn.analysis
    
    os.chdir(unicorn.GRISM_HOME+'SN-PRIMO')
    
    #### Copy necessary files from PREP_FLT to DATA
    os.chdir('PREP_FLT')
    grism_asn  = glob.glob('PRIMO-1???-G141_asn.fits')
    files=glob.glob('PRIMO-1???-G141_shifts.txt')
    files.extend(grism_asn)
    files.append('PRIMO-1026-F160W_tweak.fits')
    for file in files:
        status = os.system('cp '+file+' ../DATA')
        #shutil.copy(file, '../DATA')
    
    try:
        iraf.imcopy('/Users/gbrammer/CANDELS/GOODS-S/PREP_FLT/PRIMO-F125W_drz.fits[1]', '../DATA/f125w.fits')
    except:
        os.remove('../DATA/f125w.fits')
        iraf.imcopy('/Users/gbrammer/CANDELS/GOODS-S/PREP_FLT/PRIMO-F125W_drz.fits[1]', '../DATA/f125w.fits')
    
    os.chdir('../')
    
    #### Initialize parameters
    go.set_parameters(direct='F160W', LIMITING_MAGNITUDE=26)
    
    #### Main loop for reduction
    for i, asn in enumerate(grism_asn):
        threedhst.options['PREFAB_DIRECT_IMAGE'] = '/Users/gbrammer/CANDELS/GOODS-S/PREP_FLT/PRIMO-F160W_drz.fits'
        threedhst.options['OTHER_BANDS'] = [['f125w.fits', 'F125W' , 1248.6, 26.25]]
        proc.reduction_script(asn_grism_file=asn)
        unicorn.analysis.make_SED_plots(grism_root=asn.split('_asn.fits')[0])
        go.clean_up()

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

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

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

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

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

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

    f.close()
    o.close()

def slice_rc(img):
    '''
    Slices the Rainbow Camera into 4 different images and adds the 'filter' keyword in the fits file.
    '''
    fname = os.path.basename(img)
    fdir = os.path.dirname(img)    
    
    # Running IRAF
    iraf.noao(_doprint=0)    
    
    corners = {
    "g" : [1, 910, 1, 900],
    "i" : [1, 910, 1060, 2045],
    "r" : [1040, 2045, 1015, 2045],
    "u" : [1030, 2045, 1, 900]
    }
    
    
    filenames = []
        
    for i, b in enumerate(corners.keys()):
        logger.info( "Slicing for filter %s"% b)
        name = fname.replace(".fits", "_%s.fits"%b)
        
        #Clean first
        if (os.path.isfile(name)):
            os.remove(name)
            
        iraf.imcopy("%s[%d:%d,%d:%d]"%(img, corners[b][0], corners[b][1], corners[b][2], corners[b][3]), name)
         
        fitsutils.update_par(name, 'filter', b)
        is_on_target(name)
        
        filenames.append(name)
    
    return filenames