def flatCorrection():
	print "FlatCorrection start"
	#put all others params to no , they may be set by previous actions
	#WITH ccdproc
#	iraf.ccdproc.setParam('flatcor', 'yes')
#	iraf.ccdproc.setParam('fixpix', 'no')
#	iraf.ccdproc.setParam('darkcor', 'no')
#	iraf.ccdproc.setParam('illumcor', 'no')
#	iraf.ccdproc.setParam('trim', 'no')
#	iraf.ccdproc.setParam('overscan', 'no')
#	iraf.ccdproc.setParam('zerocor', 'no')
#	iraf.ccdproc.setParam('trimsec', '')
#	iraf.ccdproc.setParam('biassec', '')
#	#online
#	iraf.ccdproc.setParam('output', '')
	for f in FILTERS:	
		flatfilename = os.path.join(OUTPUTDIR, "flat", f, "FlatNorm.fits")
		if os.path.isfile(flatfilename):
#			#WITH ccdproc
#			iraf.ccdproc.setParam('flat', flatfilename)
#			iraf.ccdproc.setParam("images", os.path.join(OUTPUTDIR, "object", f) + "/*.fits")
#			iraf.ccdproc()
			with open(os.path.join(OUTPUTDIR, "object", f, "list")) as file1:
				for line in file1:
					objfilename = line.strip()
					if objfilename!="":
						objFullfilename = os.path.join(OUTPUTDIR, "object", f,	line.strip())
						iraf.imarith(objFullfilename, '/', flatfilename, objFullfilename)
		else:
			print "Flat file %s not present" % flatfilename
	print "FlatCorrection end"

 def sub_darks(self,data):
     for exp in data.keys():
         inputlist = data[exp]['raw']
         outputlist = [s[s.rfind('/')+1:s.find('.FIT')]+\
                           '_ds.fits'\
                           for s in inputlist]
         
         if self.darks[exp]['combined'] == None: self.gen_dark(exp)
         if len(inputlist) > 10:
             while len(inputlist) > 0:
                 inputstring = ','.join(inputlist[:10])
                 outputstring = ','.join(outputlist[:10])
                 iraf.imarith(inputstring,'-',\
                                  self.darks[exp]['combined'],\
                                  outputstring)
                 data[exp]['ds'] += outputlist[:10]
                 inputlist = inputlist[10:]
                 outputlist = outputlist[10:]
             
         else:
             iraf.imarith(','.join(inputlist),\
                              '-',\
                              self.darks[exp]['combined'],\
                              ','.join([s[s.rfind('/')+1:s.find('.FIT')]+\
                                            '_ds.fits' for s in inputlist]))
             data[exp]['ds'] += outputlist

def skySub(dfView):
    """
	given the dfview, and using the skyflat calculated in makeSkyFlat,
	subtract the skyflat from each image
	the output will be one sky-subtracted image for each image in the dfView
	the names of the sky-subtracted fits images have 'sky' infront of them
	"""
    # grab the paths to images to be sky-subtracted
    images = dfView.file.values.tolist()
    # adding a 's' infront of each file name, save to a seperate list
    skySubImages = ["s-" + i[5:] for i in images]
    # grab the YYMMDD observation date from this view
    date = str(dfView.Date.values[0])

    # organize the input, skyflat, and output in iraf-friendly ways
    inputFiles = joinStrList(images)
    skyFlat = "scratch/" + date + "sky.fits"
    outputSkySub = joinStrList(skySubImages, scratch=True)
    iraf.imarith(
        inputFiles,
        "-",
        skyFlat,
        outputSkySub,
        divzero=0.0,
        hparams="",
        pixtype="",
        calctype="",
        verbose="yes",
        noact="no",
    )
    return

def createFlatFiles():
	import re
	print "CreateFlatFiles start"
	for f in FILTERS:
		if(os.listdir(os.path.join(OUTPUTDIR, "flat", f))):
			iraf.imcombine.setParam("input", os.path.join(OUTPUTDIR, "flat", f) + "/*.fits")
			flatFile = os.path.join(OUTPUTDIR, "flat", f , "Flat.fits")
			if os.path.exists(flatFile):
				print("flatFile %s alreday exists deleting"  % flatFile)
				os.remove(flatFile)
			iraf.imcombine.setParam("output", flatFile)
			#from doc:	
			#http://www.iac.es/sieinvens/siepedia/pmwiki.php?n=HOWTOs.PythonianIRAF
			#--> iraf.listpix(mode='ql')     # confirms parameter
			#--> iraf.listpix(mode='h')     # doesn't ask for [email protected]@
			iraf.imcombine(mode="h")
			#NORMALIZE
			#imstat
			res = iraf.imstat(flatFile, Stdout=1)
			print(res[0].strip()) 
			print(res[1].strip()) 
			resArray = re.split("\s+", res[1].strip())
			#max value
			#toDivValue = float(resArray[5])
			#meanValue
			toDivValue = float(resArray[2])
			flatNormFile = os.path.join(OUTPUTDIR, "flat", f , "FlatNorm.fits")
			if os.path.exists(flatNormFile):
				print("flatNormFile %s alreday exists deleting"  % flatNormFile)
				os.remove(flatNormFile)
			#divide by max value
			iraf.imarith(flatFile, '/', toDivValue, flatNormFile)
		else:
			print("NO FLAT FILES for filter %s PRESENT" %f)
	print "CreateFlatFiles end"

def ExpNormalize(images, outbase="_n"):
			
	# Build the list of output image names
	out = [os.path.splitext(i) for i in images]	# Split off the extension
	out = [i[0] + outbase + '.fits' for i in out]	# Paste the outbase at the end of the filename 
												# and put the extension back on
	# Get a list of exposure times.
	exp_times = [GetHeaderKeyword(i, 'exptime') for i in images]
	exp_times = [str(e) for e in exp_times]	
	
	# run imarith to do the normalization
	iraf.imarith.unlearn()
	iraf.imarith.op = '/'
	iraf.imarith.mode = 'h'

	for i in range(len(images)):
		iraf.imarith.operand1 = images[i]
		iraf.imarith.operand2 = exp_times[i]
		iraf.imarith.result = out[i]
	
		iraf.imarith()

		# update the exptime keyword		
		iraf.hedit.unlearn()
		iraf.hedit.verify='no'
		iraf.hedit.show='yes'
		iraf.hedit.update='yes'
		iraf.hedit.images=out[i]
		iraf.hedit.fields='exptime'
		iraf.hedit.value=1
		iraf.hedit.mode='h'

		iraf.hedit(Stdout=1)

	return out

 def unit_converter(self):
     '''
     Converts DRZ images from units of counts/sec to counts
     '''
     
     for image, exptime in zip(self.imlist, self.explist):
         iraf.imarith(operand1=image + '[1]', op='*', operand2=exptime, \
                      result=image[:9] + '_counts.fits')

    def test_hedit(self):
        if os.path.exists('image.real.fits'):
            os.remove('image.real.fits')

        iraf.imarith('dev$pix', '/', '1', 'image.real', pixtype='r')
        iraf.hedit('image.real', 'title', 'm51 real', verify=False,
                   Stdout="/dev/null")
        with fits.open('image.real.fits') as f:
            assert f[0].header['OBJECT'] == 'm51 real'

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 normalise(self, out=None):
     """
     Normalise the image.
     """
     if out is None:
         name, ext = os.path.splitext(self.filename)
         out = name + "_norm" + ext
     iraf.imarith(self.filename, "/", self.MEAN, out)
     return out

 def test_imarith(self):
     iraf.imarith('dev$pix', '/', '1', 'image.real', pixtype='r')
     with fits.open('image.real.fits') as f:
         assert f[0].header['BITPIX'] == -32
         assert f[0].data.shape == (512, 512)
     iraf.imarith('dev$pix', '/', '1', 'image.dbl', pixtype='d')
     with fits.open('image.dbl.fits') as f:
         assert f[0].header['BITPIX'] == -64
         assert f[0].data.shape == (512, 512)

def run_scldark(input):
    iraf.image(_doprint=0)
    iraf.image.imutil(_doprint=0)
    scl=float(input[1])/float(input[2])
#    print input[1],scl
    iraf.imarith(input[0],"*",scl,input[3],divzero="0.0",hparams="",pixtype="",calctype="",verbose="no",noact="no")
    iraf.hedit(input[3],"EXPTIME",input[1], add="yes", addonly="yes", delete="no", verify="no", show="no", update="yes")
    iraf.hedit(input[3],"DARKTIME",input[1], add="yes", addonly="yes", delete="no", verify="no", show="no", update="yes")
    iraf.hedit(input[3],"EXPOSURE","", add="no", addonly="no", delete="yes", verify="no", show="no", update="yes")

def MatchNSubtract(TargetImg,Template,OutputImage,fitgeometry="general"):
    """ Creates OutputImage =  TargetImg - Template after scaling and matching Template to TargetImg.
    fitgeometry can be set to 'rotate' when the Template is also TIRSPEC data
    Otherwise if Template is 2MASS or other instrument set it as 'general'  """
    
    AlignedImg = os.path.splitext(TargetImg)[0]+"_"+os.path.basename(Template)
    AlignedImg = os.path.splitext(AlignedImg)[0][:115]+'.fits' # Reduce filename length for iraf geopmap
    TransformDBfile = AlignImage(TargetImg,Template,AlignedImg,fitgeometry=fitgeometry)
    
    # Now get the Good sky region coordinates
    SkyCoordsFile = os.path.splitext(TargetImg)[0]+'_BlankSky.coo'
    if not os.path.isfile(SkyCoordsFile) :
        iraf.display(TargetImg,1)
        print ('For taking coordinates of good sky. Press _x_ over blank sky areas.')
        imx=iraf.imexam(Stdout=1)
        with open(SkyCoordsFile,'w') as foo :    #Creating blank sky coords files
            for line in imx :               
                foo.write(line.split()[0] +'  '+line.split()[1]+'\n')

    # Now get the regions in the image whose brightness has to be cancelled by scaling
    FluxCoordsFile = os.path.splitext(TargetImg)[0]+'_FluxRegions.coo'
    if not os.path.isfile(FluxCoordsFile) :
        iraf.display(TargetImg,1)
        print ('Press _x_ over areas you want to minimise the residual flux after subtraction')
        imx=iraf.imexam(Stdout=1)
        with open(FluxCoordsFile,'w') as foo :    #Creating Flux ares which we have to remove in subtraction
            for line in imx :               
                foo.write(line.split()[0] +'  '+line.split()[1]+'\n')

    #Now we first has to remove background from both the images.
    TargetSkySubtractedFile = os.path.splitext(TargetImg)[0]+'_SkyS.fits'
    if not os.path.isfile(TargetSkySubtractedFile):
        skyvalue = SkySubtractImage(TargetImg,TargetSkySubtractedFile,SkyCoordsFile)
    else:
        print('Warning: Using old {0} file'.format(TargetSkySubtractedFile))

    AlignedSkySubtractedFile = os.path.splitext(AlignedImg)[0]+'_SkyS.fits'
    if not os.path.isfile(AlignedSkySubtractedFile):
        skyvalue = SkySubtractImage(AlignedImg,AlignedSkySubtractedFile,SkyCoordsFile)
    else:
        print('Warning: Using old {0} file'.format(AlignedSkySubtractedFile))

    #We shall now extract the totel Flux in each tiles from both the images
    TargetFluxinTiles = ExtractTiles(TargetSkySubtractedFile,FluxCoordsFile,Summeryfunction=np.sum,hsize=7*1.5)
    TemplateFluxinTiles = ExtractTiles(AlignedSkySubtractedFile,FluxCoordsFile,Summeryfunction=np.sum,hsize=7*1.5)
    
    def DiffSquareSum(x):
        return np.sum([(targetF - x*templateF)**2 for targetF,templateF in zip(TargetFluxinTiles,TemplateFluxinTiles)])
    
    res = scipy.optimize.minimize_scalar(DiffSquareSum)
    Scale = res.x
    print('Scaling to match the fluxes is {0}'.format(Scale))
    iraf.imarith(operand1=AlignedSkySubtractedFile,op="*",operand2=Scale,result=os.path.splitext(AlignedSkySubtractedFile)[0]+'M.fits')

    iraf.imarith(operand1=TargetSkySubtractedFile,op="-",operand2=os.path.splitext(AlignedSkySubtractedFile)[0]+'M.fits',result=OutputImage)

def imarith(operand1, op, operand2, result):
    from pyraf import iraf
    iraf.images()

    pars = iraf.imarith.getParList()
    iraf.imcombine.unlearn()

    print "%s %s %s -> %s" % (operand1, op, operand2, result)
    iraf.imarith(operand1=operand1, op=op, operand2=operand2, result=result)

    iraf.imarith.setParList(pars)

	def get_median(self):
		'''
		Calculates the median subtracted image used to compute the pixel mask
		'''
		iraf.median(input = self.image, output = 'tmp_med.fits', 
					xwindow = 40, ywindow = 40, verbose ='No')

		iraf.imarith(operand1 = self.image, operand2 = 'tmp_med.fits[0]', 
					 op = '-', result = self.output + '_sub.fits')
		
		os.remove('tmp_med.fits')

def mask_norm(mask):
    # normalizes mask.
    # test if trimmed mask exists:
    if os.path.isfile(mask):
        iraf.cd(os.path.dirname(mask))
	hdulist = pyfits.open(mask)
	imdata = hdulist[0].data
	hdulist.close()
	outname = mask.replace('_trimmed.fits', '_normed.fits')
	print imdata.max()
	iraf.imarith(os.path.basename(mask),"/",imdata.max(),os.path.basename(outname))

def get_masked_image(img):

    # Running IRAF
    iraf.noao(_doprint=0)  
    mask = make_mask_cross(img)    
    masked = img.replace(".fits", "_masked.fits")
    
    if (os.path.isfile(masked)):
        os.remove(masked)
    iraf.imarith(img, "*", mask, masked)
    
    return masked

def sub_bias(image, combined_bias, image_b):
# Import IRAF modules:
  iraf.images(_doprint=0)
  iraf.imutil(_doprint=0)
  parList = "bias_subtraction_imarith.par"
# Check input file and combined_bias frame exists before proceeding:
  if os.path.isfile(image) == True:
    if os.path.isfile(combined_bias) == True:
      if os.path.isfile(parList) == True:
# Subtract combined bias frame from input frame (object or flat)
# using IRAF task imarithmetic:
        iraf.imarith.setParList(ParList="bias_subtraction_imarith.par")
        iraf.imarith(operand1=image, operand2=combined_bias, result=image_b)
        print ' '
        print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        print 'Bias frame ' + str(combined_bias) 
        print 'subtracted from input ' + str(image)
        print 'to create ' + str(image_b)
        print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        print ' '
      else:
        print ' '
        print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        print 'Bias frame subtraction IRAF .par file              ' 
        print str(parList)
        print 'does not exist. Exiting script.                    '
        print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        print ' '
        print ' '
        sys.exit()
    else:
      print ' '
      print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      print 'Combined bias frame                                ' 
      print str(combined_bias)
      print 'does not exist. Exiting script.                    '
      print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      print ' '
      print ' '
      sys.exit()
  else:
    print ' '
    print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    print 'Input frame                                        ' 
    print str(image)
    print 'does not exist. Exiting script.                    '
    print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    print ' '
    print ' '
    sys.exit()

def apply_imarith (op1,op,op2):
    if op1.find('fits')!= -1:
       operand = op1
       result = ''
    else:
        operand = '@' + op1
        result = operand
    iraf.imarith.operand1 = operand
    iraf.imarith.op = op
    iraf.imarith.operand2 = op2
    iraf.imarith.result = result
    #print operand,op,op2
    iraf.imarith(mode='h')
    print 'Applying %s to %s' % (op2,op1)

def CombineImages(imglist,output,method='median',zero='none' ,scale='none',norm=False,
                  reject="avsigclip", statsec='[150:900,150:900]'):
    """ Combined the input list of images and write to output fits file. """
    iraf.imcombine.unlearn()
    imglistfname = os.path.splitext(output)[0]+'.comblist'
    with open(imglistfname,'w') as imgs2combinefile:
        imgs2combinefile.write('\n'.join(imglist)+'\n')
    # Now call iraf imcombine with zero scaling
    combineoutputfile = os.path.splitext(output)[0]+'_un.fits' if norm else output

    iraf.imcombine(input='@'+imglistfname, output=combineoutputfile, combine=method, 
                   reject=reject, statsec=statsec, scale=scale, zero=zero)
    if norm:
        mediancounts = np.median(fits.getdata(combineoutputfile))
        iraf.imarith(operand1=combineoutputfile,op='/',operand2=mediancounts,result=output)

def DarkSubtract(images, cal_path, outbase="_d"):
	
	# Get a list of dark images
	darks = [FindDarkFrame(i, cal_path) for i in images]
	out = [os.path.splitext(i)[0]+outbase+'.fits' for i in images]

	iraf.imarith.unlearn() # initial imarith setup
	iraf.imarith.mode='h'

	print "\n******************"
	print "Dark Subtracting: "
	print "******************"
	for i in range(len(images)):

		iexp_time = GetHeaderKeyword(images[i], 'exptime')
		dexp_time = GetHeaderKeyword(darks[i], 'exptime')

		# If exposure times don't match, scale the dark and fill in the operand2 field
		if iexp_time != dexp_time:
			iraf.imarith.operand1 = darks[i]
			iraf.imarith.operand2 = iexp_time/dexp_time	
			iraf.imarith.op = '*'

			# Create temporary dark filename
			tdark = os.path.join(os.path.dirname(darks[i]),'tmpdark')
			iraf.imarith.result = tdark
			iraf.imarith()
		
			iraf.imarith.operand2 = tdark

		else: 
			iraf.imarith.operand2 = darks[i]

		# Perform the dark subtraction
		iraf.imarith.operand1 = images[i]
		iraf.imarith.op = '-'
		iraf.imarith.result = out[i]
		print iraf.imarith.operand1 + ' - ' + iraf.imarith.operand2
		iraf.imarith()

		# Remove the scaled dark
		if iexp_time != dexp_time:
			os.remove(tdark+'.fits')

	# Make sure out files actually exist 
	out = filter(os.path.exists, out)
	
	return out