import os,numpy as np,copy,sys,glob,logging,scipy.linalg,matplotlib,matplotlib.pyplot as plt,paircarsdata,time
from casatools import msmetadata,table,measures,quanta,agentflagger,image,calibrater,ms
from casatasks import *
from . import access_ms as am
from . import basic_func as B
from . import flagger as fg
from paircarsutils.paircars_casatasks.poltclean import *
from astropy.io import fits
from astropy.wcs import WCS
from . import mwapb as pb
from datetime import datetime 
import scipy.ndimage as ndi
os.system('rm -rf casa*log')
matplotlib.use('Agg')
'''
Code is written by Devojyoti Kansabanik, 17 Jan, 2021
'''
datadir = os.path.dirname(paircarsdata.__file__)
class IntensitySelfcal:
	'''
	Generic class to perform intensity self-calibration

	Parameters
	----------
	msname : str 
		Name of the measurement set
	metafits : str 
		Name of the metafits file
	num_pixel_in_psf : int
		Number of pixels side one point-spread-funtion
	maximum_emission_scale : float 
		Maximum scale of the emission present in the image in arcsec
	largest_scale : float 
		Largest spatial scale in degree used for self calibration
	verbose : bool
		If True keep all the intermediate images, model, residuals, caltables and details of the log to detailed analysis
	use_wsclean : bool
		Use WSClean or not
	includeflag : bool
		Include flag data or not 
	'''
	def __init__(self,msname,metafits,maximum_emission_scale,num_pixel_in_psf=3,largest_scale=20,verbose=False,use_wsclean=True,savelog=True,includeflag=True):
		self.cwd=os.getcwd()
		if msname[-1]=='/':
			self.msname=msname[:-1]
		else:
			self.msname=msname
		self.savelog=savelog
		self.metafits=metafits
		self.mspath=os.path.dirname(os.path.realpath(msname))
		AM=am.AccessMS(self.msname)
		self.bandwidth=round((AM.calc_bandwidth()/10**3),2)
		self.time_width=round(AM.calc_total_time(),2)
		IB=B.ImageBasic(self.msname,includeflag=includeflag)
		self.psf=IB.calc_psf()
		self.max_baseline=IB.calc_max_baseline_lambda()
		self.num_pixel_in_psf=int(num_pixel_in_psf)
		self.cellsize=IB.calc_cellsize(self.num_pixel_in_psf)
		self.imsize=IB.num_pixels(self.num_pixel_in_psf)
		self.max_size=maximum_emission_scale
		self.multiscale_scales=IB.choose_scales(self.num_pixel_in_psf,self.max_size)
		self.uvtaper=IB.calc_uvtaper(includeflag=True)
		calib_uvrange=IB.calc_calib_uvrange(largest_scale,includeflag=includeflag) # Short baselines sensitive to larger than 20 deg are excluded 
		self.calib_minuv_meter=calib_uvrange[1]
		self.calib_maxuv_meter=calib_uvrange[2]
		self.calib_uvrange_lambda=str('>'+str(calib_uvrange[3])+'lambda')
		imaging_uvrange=IB.calc_calib_uvrange(largest_scale,includeflag=False)
		self.imaging_uvrange_lambda=str('>'+str(imaging_uvrange[3])+'lambda')
		self.imaging_minuv_lambda=imaging_uvrange[3]
		self.imaging_maxuv_lambda=imaging_uvrange[4]
		self.rms_box='50,50,'+str(self.imsize-50)+','+str(int(self.imsize/4)) # CASA box to calculate the rms
		boxcenter_y=int(self.imsize/2)-int(3*3600/self.cellsize)
		ywidth=int(3600/self.cellsize)
		boxcenter_x=int(self.imsize/2)
		self.rms_box_nearsun=str(boxcenter_x-int(self.imsize/5))+','+str(boxcenter_y-2*ywidth)+','+str(boxcenter_x+int(self.imsize/5))+','+str(boxcenter_y+ywidth)
																											 # CASA box to calculate the rms near the Sun
		self.verbose=verbose
		self.wsclean=use_wsclean
		formatter = logging.Formatter('%(asctime)s %(levelname)-8s %(message)s',datefmt='%Y-%m-%d %H:%M:%S')
		self.log_verbose = logging.getLogger('selfcal_verbose_log')
		if self.log_verbose.hasHandlers()==False:
			self.log_verbose.setLevel(logging.DEBUG)
			if self.verbose:
				self.console=logging.StreamHandler(sys.stdout)
				self.console.setFormatter(formatter)
				self.log_verbose.addHandler(self.console)
			if self.savelog==True:
				self.filehandle=logging.FileHandler(self.cwd+'/Intensity_Selfcal_verbose.log')
				self.filehandle.setFormatter(formatter)
				self.log_verbose.addHandler(self.filehandle)
			self.log_verbose.propagate = False
		self.log_verbose.info('Initiating Intensity selfcal object.\n')
		wsclean_datadir=os.path.abspath(os.path.dirname(paircarsdata.__file__))
		c=0
		self.wsclean_start=''
		while c<5:
			try:
				wsclean_path=str(np.load(wsclean_datadir+'/wsclean_path.npy',allow_pickle=True))
				os.path.join(wsclean_path)
				a=os.system(wsclean_path+'/wsclean > wsclean_test')
				if a!=0:
					self.log_verbose.info('WSClean is not installed. Using CASA for imaging.\n')
					self.wsclean=False
					self.wsclean_start=None
				else:
					self.wsclean_start=wsclean_path+'/wsclean'
				break
			except:
				time.sleep(1)
				c+=1
		if self.wsclean_start=='':
			self.log_verbose.warning('WSClean is not installed. Using CASA for imaging.\n')
			self.wsclean=False
			self.wsclean_start=None

	def negative_box(self,max_pix,box_width=3):
		'''
		Create a box about the maximum pixel of image to search negative

		Parameters
		----------
		max_pix : list 
			Maximum pixel [xxmax,yymax]
		box_width : float 
			Box width in degree (default : 3 degree)
		Returns
		-------
		str
			CASA box 'xblc,yblc,xrtc,yrtc'
		'''
		max_pix_xx=max_pix[0]
		max_pix_yy=max_pix[1]
		box_length=(float(box_width)*3600.0)/self.cellsize # Taking a box about the max pixel to serach to minimum
		xblc=max_pix_xx-(box_length/2.0)
		xrtc=max_pix_xx+(box_length/2.0)
		yblc=max_pix_yy-(box_length/2.0)
		yrtc=max_pix_yy+(box_length/2.0)
		if xblc<0:
			xblc=0
		if yblc<0:
			yblc=0
		if xrtc>=self.imsize:
			xrtc=self.imsize-1
		if yrtc>=self.imsize:
			yrtc=self.imsize-1
		os.system('rm -rf casa*log')
		return str(int(xblc))+','+str(int(yblc))+','+str(int(xrtc))+','+str(int(yrtc))

	def calc_dyn_range(self,imagename,sigma,box_width=3,stokes_list=['I']):
		'''
		Calculate the dynamic range of the full Stokes image cube

		Parameters
		----------
		imagename : str
			Name of the image
		sigma : float 
			nsigma value to put a mask for calculating total flux
		box_width : float 
			Negative box width around the maximum pixel in degree (default : 3 degree)
		stokes_list : list 
			List of stokes planes in the image
		Returns
		-------
		dict	
			Dynamic range information {'STOKES':[rms dynamic range,rms,total_flux(non-negative)]}
		float
			Negative dynamic range for Stokes I
		'''
		ia=image()
		out_dict={}
		try:
			if os.path.isdir(imagename):
				imageheader=imhead(imagename=imagename,mode='summary')
				if imageheader['ndim']>2 and imageheader['ndim']==4:
					neg_dyn=0
					for stokes in stokes_list:
						if stokes=='I' or stokes=='XX' or stokes =='YY' or stokes=='RR' or stokes =='LL':
							maxpos=imstat(imagename=imagename,stokes=stokes)['maxpos']
							negative_box=self.negative_box(maxpos,box_width=box_width)
							max_pix=imstat(imagename=imagename,stokes=stokes)['max'][0]
							rms=imstat(imagename=imagename,box=self.rms_box,stokes=stokes)['rms'][0]
							min_pix=imstat(imagename=imagename,box=negative_box,stokes=stokes)['min'][0]
							rms_dyn_range=max_pix/rms
							if min_pix!=0:
								neg_dyn+=max_pix/abs(min_pix)
							else:
								neg_dyn=rms_dyn_range
							ia.open(imagename)
							ia.calcmask('\"'+imagename+'\">'+str(sigma*rms),'mymask')
							ia.close()
							try:
								total_flux=imstat(imagename=imagename,stokes=stokes)['flux'][0]
								out_dict[stokes]=[rms_dyn_range,rms,total_flux]							
							except:
								out_dict[stokes]=[rms_dyn_range,rms,np.nan]
							makemask(mode='delete',inpmask=imagename+':mymask')
						else:
							max_pix=imstat(imagename=imagename,stokes=stokes)['max'][0]
							min_pix=imstat(imagename=imagename,stokes=stokes)['min'][0]
							rms=imstat(imagename=imagename,box=self.rms_box,stokes=stokes)['rms'][0]
							if abs(min_pix)>max_pix:
								max_pix=abs(min_pix)
							rms_dyn_range=max_pix/rms
							if os.path.isdir(stokes+'.image'):
								os.system('rm -rf '+stokes+'.image')
							immath(imagename=imagename,outfile=stokes+'.image',mode='evalexpr',expr='abs(IM0)',stokes=stokes)
							ia.open(stokes+'.image')
							ia.calcmask(stokes+'.image'+'>'+str(sigma*rms),'mymask')
							ia.close()
							try:
								total_flux=imstat(imagename=stokes+'.image',stokes=stokes)['flux'][0]
								out_dict[stokes]=[rms_dyn_range,rms,total_flux]
							except:
								out_dict[stokes]=[rms_dyn_range,rms,np.nan]	
							os.system('rm -rf '+stokes+'.image')			
			else:
				self.log_verbose.info('Error in calculating dynamic range : '+imagename+' not found\n')
				out_dict={}
				neg_dyn=np.nan
				out_dict['NAN']=[np.nan,np.nan]
			if stokes_list==['XX','YY'] or stokes_list==['RR','LL']:
				negative_dyn_range=neg_dyn/2
			else:
				negative_dyn_range=neg_dyn
		except Exception as e:
			self.log_verbose.info('Error in calculating dynamic range : '+str(e)+'\n')
			negative_dyn_range=np.nan
			out_dict['NAN']=[np.nan,np.nan]
		os.system('rm -rf casa*log')
		return out_dict,negative_dyn_range

	def get_residual_flux_frac(self,imagename,residualname,sigma,stokes_list=['I']):
		'''
		Function to calculate residual flux fraction
		Parameters
		----------
		imagename : str
			Name of the image
		residualname : str
			Residual image name
		sigma : int
			Sigma value for threshold
		stokes : list
			Stokes list
		Returns
		-------
		float
			Residual flux
		float
			Residual flux fraction
		'''
		rms=imstat(imagename=imagename,box=self.rms_box)['rms'][0]
		maxpos=imstat(imagename=imagename)['maxpos']
		flux_box=self.negative_box(maxpos,box_width=3)
		image_flux=0
		residual_flux=0
		for stokes in stokes_list:
			image_flux+=imstat(imagename=imagename,box=flux_box,stokes=stokes)['flux'][0]
			residual_flux+=imstat(imagename=residualname,box=flux_box,stokes=stokes)['flux'][0]
		image_flux=image_flux/len(stokes_list)
		residual_flux=residual_flux/len(stokes_list)
		return abs(residual_flux),abs(residual_flux/image_flux)

	def change_start_sigma(self,nsigma,sigma_step,min_sigma):
		'''
		Function to calculate start sigma

		Parameters
		----------
		nsigma : float 
			Present sigma value
		sigma_step : float 
			Sigma reduction step
		min_sigma : float 
			Minimum value of sigma to reduce
		Returns
		-------
		float
			Value of the start sigma
		'''
		if nsigma-sigma_step > min_sigma: # If initial CLEAN could not pick up any flux at start sigma, we lower down start sigma upto minsigma
			if sigma_step>1.0:
				print ('WARNING : Choosing sigma step greater than 1 is too risky. Selfcal may diverge\n')
				self.log_verbose.info('Continuing with sigma step :'+str(sigma_step)+'\n')
				if self.verbose==False:
					print ('Continuing with sigma step :'+str(sigma_step)+'\n')
				nsigma-=sigma_step
			else:
				nsigma-=sigma_step
		else:
			self.log_verbose.info('Stop here, since start sigma is lower than '+str(min_sigma)+'\n')
			nsigma=np.nan
		os.system('rm -rf casa*log')
		return nsigma
		
	def reduce_sigma(self,imagename,nsigma,sigma_step,minsigma,pre_residual=0.0,residual_frac=0.01,stokes_list=['I']):
		'''
		Function to determine whether reduce the CLEAN sigma or not

		Parameters
		----------
		imagename : str 
			Prefix of the image name
		nsigma : float 
			Value of the present n-sigma
		sigma_step : float 
			Step to reduce sigma value
		minsigma : float
			Minimum allowed sigma
		pre_residual : float 
			Previous residual fraction to compare (default : 0.0)
		residual_frac : float 
			Residual flux fraction to reduce sigma (default : 0.01)
		stokes_list : list 
			Stokes plane list
		Returns
		-------
		float
			Reduced value of n-sigma if residual flux is more than given percentage (default : 1%) of the total flux in Stokes I or in all Stokes Q,U,V.
		'''
		imagename=imagename
		residual=imagename.split('.image')[0]+'.residual'
		do_reduce_list=[]
		ia=image()
		imagename_path=os.path.dirname(os.path.realpath(imagename))
		cwd=os.getcwd()
		if imagename_path!='':
			os.chdir(imagename_path)
		os.system('rm -rf reduce_sigma_*')
		for stokes in stokes_list:
			self.log_verbose.info('imstat(imagename=\''+imagename+'\',box=\''+self.rms_box+'\',stokes=\''+stokes+'\')[\'rms\'][0]\n')
			rms=imstat(imagename=imagename,box=self.rms_box,stokes=stokes)['rms'][0]
			os.system('rm -rf reduce_sigma_*')
			immath(imagename=imagename,mode='evalexpr',expr='abs(IM0)',stokes=stokes,outfile='reduce_sigma_'+stokes+'.image')
			immath(imagename=residual,mode='evalexpr',expr='abs(IM0)',stokes=stokes,outfile='reduce_sigma_'+stokes+'.residual')
			ia.open('reduce_sigma_'+stokes+'.image')			
			ia.calcmask('\"reduce_sigma_'+stokes+'.image\">'+str((nsigma-sigma_step)*rms),'mymask')
			ia.close()
			makemask(inpimage='reduce_sigma_'+stokes+'.image',inpmask='reduce_sigma_'+stokes+'.image:mymask',output='reduce_sigma_'+stokes+'.residual:mymask',mode='copy')
			try:
				image_pix_sum=imstat(imagename='reduce_sigma_'+stokes+'.image')['sum'][0]
				residual_pix_sum=imstat(imagename='reduce_sigma_'+stokes+'.residual')['sum'][0]
			except:
				image_pix_sum=1
				residual_pix_sum=0
			try:
				maxval=imstat(imagename='reduce_sigma_'+stokes+'.residual')['max'][0]
			except:
				maxval=0
			try:
				minval=imstat(imagename='reduce_sigma_'+stokes+'.residual')['min'][0]
			except:
				minval=0					
			if maxval>(nsigma-sigma_step)*rms:
				max_frac_diff=(maxval-(nsigma-sigma_step)*rms)/maxval
			else:
				max_frac_diff=0
			if ((residual_pix_sum/image_pix_sum>residual_frac or (max_frac_diff>0 and max_frac_diff>residual_frac)) and abs(minval)<(nsigma-sigma_step)*rms):
				if (pre_residual>0 and residual_pix_sum/image_pix_sum<pre_residual) or pre_residual==0:
					do_reduce_list.append(1)
		os.system('rm -rf reduce_sigma_*')
		os.chdir(cwd)
		if int(np.sum(np.array(do_reduce_list)))>=1:
			if sigma_step>1.0:
				self.log_verbose.warning('Choosing sigma step 1 is too risky. Selfcal may diverge\n')
				if self.verbose==False:
					print ('WARNING : Choosing sigma step greater than 1 is too risky. Selfcal may diverge\n')
					os.system('rm -rf casa*log')
					return nsigma-sigma_step
			elif nsigma-sigma_step<minsigma:
				self.log_verbose.warning('Choosing sigma less than '+str(minsigma)+'\n')
				if self.verbose==False:
					print ('WARNING : Choosing sigma less than '+str(minsigma)+'\n')
					os.system('rm -rf casa*log')
					return minsigma
				else:
					self.log_verbose.info('Continuing with sigma step :'+str(sigma_step)+'\n')
					if self.verbose==False:
						print ('Continuing with sigma step :'+str(sigma_step)+'\n')
					os.system('rm -rf casa*log')
					return nsigma
			else:
				self.log_verbose.info('Reducing sigma to:'+str(nsigma-sigma_step)+'\n')
				os.system('rm -rf casa*log')
				return nsigma-sigma_step
		else:
			self.log_verbose.info('Sigma value is not changed. Sigma is at :'+str(nsigma)+'\n')
			os.system('rm -rf casa*log')
			return nsigma
	
	def calc_iter_num(self,safety_factor,quality_factor,scratch=True,bandpass_selfcal=False):
		'''
		Function to calculate minimum number of selfcal iteration based on safety standard and quality factor
		Parameters
		----------
		safety_factor : int 
			Factor to determine the robustness of the selfcal
		quality_factor : int 
			Factor to determine the quality of the images
		scratch : bool
			Whether the selfcal started from scratch or not
		bandpass_selfcal : bool
			Performing bandpass selfcal or not
		Returns
		-------
		int
			Minimum iteration at fixed sigma
		int
			Minimum iteration
		int
			Maximum iteration
		int
			Number of antenna bins
		float
			Residual flux fraction
		float
			Minimum allowed sigma
		float 
			Fractional flux change
		'''
		if bandpass_selfcal==True:
			antenna_bin=1
			min_sigma=6.0
			if quality_factor==0:
				min_num_iter_fixed_sigma=0
				min_iteration=3
				max_iteration=15
				res_flux=0.03
				frac_flux_change=0.08
			elif quality_factor==1:
				min_num_iter_fixed_sigma=1
				min_iteration=5
				max_iteration=20
				res_flux=0.015
				frac_flux_change=0.05
			else:
				min_num_iter_fixed_sigma=3
				min_iteration=7
				max_iteration=30
				res_flux=0.01
				frac_flux_change=0.03
		else:
			if quality_factor==0:     # Low quality (Quick look image making)
				res_flux=0.03
				frac_flux_change=0.08
				if (safety_factor==0):
					min_sigma=9.0
					min_num_iter_fixed_sigma=0
					if (scratch==True):
						min_iteration=3
						max_iteration=50
						antenna_bin=2
					else:
						min_iteration=2
						max_iteration=10
						antenna_bin=1
				elif (safety_factor==1):
					min_sigma=8.0
					min_num_iter_fixed_sigma=3
					if (scratch==True):
						min_iteration=5
						max_iteration=70
						antenna_bin=3
					else:
						min_iteration=3
						max_iteration=30
						antenna_bin=1
				else:
					min_sigma=7.0
					min_num_iter_fixed_sigma=5
					if (scratch==True):
						min_iteration=7
						max_iteration=100
						antenna_bin=5
					else:
						min_iteration=3
						max_iteration=60
						antenna_bin=1
			elif quality_factor==1:  # Medium quality imaging (Computing speed medium)
				res_flux=0.015
				frac_flux_change=0.05
				if (safety_factor==0):
					min_sigma=8.0
					min_num_iter_fixed_sigma=0
					if (scratch==True):
						min_iteration=5
						max_iteration=200
						antenna_bin=3
					else:
						min_iteration=3
						max_iteration=150
						antenna_bin=1
				elif (safety_factor==1):
					min_sigma=7.0
					min_num_iter_fixed_sigma=5
					if (scratch==True):
						min_iteration=7
						max_iteration=400
						antenna_bin=5
					else:
						min_iteration=5
						max_iteration=350
						antenna_bin=1
				else:
					min_sigma=6.0
					min_num_iter_fixed_sigma=7
					if (scratch==True):
						min_iteration=10
						max_iteration=600
						antenna_bin=5
					else:
						min_iteration=5
						max_iteration=550
						antenna_bin=1
			else:  # Best quality imaging (Computing slow)
				res_flux=0.01
				frac_flux_change=0.03
				if (safety_factor==0):
					max_iteration=700
					min_sigma=7.0
					min_num_iter_fixed_sigma=0
					if (scratch==True):
						min_iteration=7
						max_iteration=600
						antenna_bin=5
					else:
						max_iteration=550
						min_iteration=5
						antenna_bin=1
				elif (safety_factor==1):
					min_num_iter_fixed_sigma=5
					min_sigma=6.0
					if (scratch==True):
						max_iteration=650
						min_iteration=10
						antenna_bin=5
					else:
						max_iteration=600
						min_iteration=7
						antenna_bin=1
				else:
					min_sigma=5.0
					min_num_iter_fixed_sigma=10
					if (scratch==True):
						max_iteration=700
						min_iteration=13
						antenna_bin=5
					else:
						max_iteration=650
						min_iteration=7
						antenna_bin=1
		self.log_verbose.info('Quality factor : '+str(quality_factor)+', Safety standard : '+str(safety_factor)+', Scratch : '+str(scratch)+\
				', Minimum number of iteration at fixed sigma : '+str(min_num_iter_fixed_sigma)+', Minimum iteration : '+str(min_iteration)+', Antenna bins : '+str(antenna_bin)+'\n')
		os.system('rm -rf casa*log')
		return min_num_iter_fixed_sigma,min_iteration,max_iteration,antenna_bin,res_flux,min_sigma,frac_flux_change

	def antenna_string(self,antenna_list,antenna_list_index):
		'''
		Function to return antenna string from antenna list

		Parameters
		----------
		antenna_list : list 
			Antenna list
		antenna_list_index : int 
			Bin number of antenna list
		Returns
		-------
		str
			Antenna string
		'''
		antenna_string=''
		for ant in antenna_list[antenna_list_index]:
			antenna_string+=str(ant)+','
		antenna_string=antenna_string[:-1]
		os.system('rm -rf casa*log')
		return antenna_string
		
	def DR_record(self,DR,record_filename,init=True):
		'''
		Function to keep the record of dynamic ranges at different self calibration steps

		Parameters
		----------
		DR : float 
			Dynamic range of the current selfcal iteration to record
		record_filename : str 
			Name of the file to store dynamic ranges
		init : bool
			Initiating a new record from the current selfcal iteration
		Returns
		-------
		numpy.array
			Dynamic range record array
		'''
		if init==True:
			os.system('rm -rf '+record_filename+'.npy')
			DR_array=np.empty([1])
			DR_array[0]=DR
			np.save(record_filename,DR_array)
		else:
			if os.path.isfile(record_filename+'.npy')==False:
				DR_array=np.empty([1])
				DR_array[0]=DR
			else:
				DR_array=np.load(record_filename+'.npy')
				DR_array=np.append(DR_array,DR)
			np.save(record_filename,DR_array)
		os.system('rm -rf casa*log')
		return DR_array
			
	def DR_delta(self,quality_factor,DR_array,frac_increase,safety_factor): # Not tested and not used
		'''
		Function to calculate the dynamic range increse step 
		(This function is not tested in the pipeline and not used. Use it at your own risk)

		Parameters
		----------
		quality_factor : int 
			Factor to determine the image quality
		DR_array : numpy.array 
			Dynamic range records array
		frac_increase : float 
			Fractional increase of the DR step from previous DR steps
		safety_factor : int 
			Factor to determine the robustness of the selfcal
		Returns
		-------
		float
			Dynamic range increase step
		'''
		if quality_factor==0 and safety_factor==0 and len(DR_array)<3:
			DR_delta_step=20
		else:
			DR_delta_step=np.polyfit(DR_array,1)[0]*frac_increase 
		os.system('rm -rf casa*log')
		return DR_delta_step
			
	def calc_selfcal_snr(self,imagename,sigma,num_antenna_to_use):
		'''
		Function to calculate the SNR of the present iteration for selfcal

		Parameters
		----------
		imagename : str 
			Name of the image
		sigma : float 
			Present sigma value used for making the image
		num_antenna_to_use : int 
			Number of antenna used for imaging
		Returns
		-------
		float
			SNR per antenna for self-calibration
		'''
		if os.path.isdir(imagename):
			try:
				if os.path.exists('convolved_selfcal.image'):
					os.system('rm -rf convolved_selfcal.image')
				imsmooth(imagename=imagename,outfile='convolved_selfcal.image',beam={'major':'3640arcsec','minor':'3640arcsec','pa':'0deg'},targetres=True)
				rms=imstat(imagename=imagename,box=self.rms_box)['rms'][0]
				imagename='convolved_selfcal.image'
				ia=image()
				ia.open(imagename)
				ia.calcmask('\''+imagename+'\'>'+str(sigma*rms),'selfcal_snr_mask')
				ia.close()
				flux=imstat(imagename)['max'][0]
				makemask(inpmask=imagename+':selfcal_snr_mask',mode='delete')
				selfcal_snr=flux/(rms*np.sqrt(num_antenna_to_use-3))
				os.system('rm -rf casa*log convolved_selfcal.image')
				return selfcal_snr
			except Exception as e:
				self.log_verbose.info('Error in calculating selfcal SNR : '+str(e)+'\n')
				os.system('rm -rf casa*log')
				return np.nan
		else:
			os.system('rm -rf casa*log')
			return np.nan

	def dirty_image(self,start_sigma,box_width=3,antenna_to_use='',weight='briggs',robust=0.5):
		'''
		Function make dirty image

		Parameters
		----------
		start_sigma : float
			Start sigma
		box_width : float 
			Negative box width around the maximum pixel in degree (default : 3 degree)
		antenna_to_use : list 
			List of antennas
		weight : str
			Visibility weighting during imaging , 'uniform', 'natural' or 'briggs'. Default is 'briggs'
		robust : float
			Robust parameter for briggs weighting, default : 0.5
		Returns
		-------
		float
			rms based dynamic range
		float
			negative based dynamic range
		float
			Per antenna SNR for self-calibration
		int
			Error message code
		int
			Maximum pixel X position
		int
			Maximum pixel Y position
		'''
		imagename=self.msname.split('.ms')[0]+'_dirty_'+str(self.bandwidth)+'kHz_'+str(self.time_width)+'s'
		present_file=glob.glob(imagename+'*')
		robust=robust*2
		if len(present_file)!=0:
			os.system('rm -rf '+imagename+'*')
			self.log_verbose.info('rm -rf '+imagename+'*\n')
		self.log_verbose.info('Making dirty image...........\n')
		if self.wsclean==False or self.wsclean_start==None:
			self.log_verbose.info('tclean(vis=\''+self.msname+'\',imagename=\''+imagename+'\',imsize=['+str(self.imsize)+'],cell=\''+str(self.cellsize)+'arcsec\',niter=0,antenna=\''+\
									antenna_to_use+'\',uvtaper=\''+str(self.uvtaper)+'\',weighting=\''+weight+'\',robust='+str(robust)+')\n')
			tclean(vis=self.msname,imagename=imagename,imsize=[self.imsize],cell=str(self.cellsize)+'arcsec',niter=0,antenna=antenna_to_use,uvtaper=self.uvtaper,\
					weighting=weight,robust=robust)
			casa_imagename=imagename+'.image'
		else:
			AM=am.AccessMS(self.msname)
			nant=AM.get_num_antenna()
			ant_list=antenna_to_use.split(',')
			baseline_flagged=False
			if len(ant_list)<nant:
				baselines=self.ant_to_baselines(ant_list=ant_list,include_baselines=False,total_ant=nant)
				baselines_str=';'.join(baselines)
				flagdata(vis=self.msname,mode='manual',antenna=baselines_str,flagbackup=True)
				baseline_flagged=True
			scales=[str(i) for i in self.multiscale_scales]
			if weight=='briggs':
				weight=weight+' '+str(robust)
			taper_tukey=10
			wsclean_args=['-scale '+str(self.cellsize)+'asec','-size '+str(self.imsize)+' '+str(self.imsize),'-weight '+weight,'-taper-edge-tukey '+str(taper_tukey),'-taper-edge '+str(self.imaging_maxuv_lambda),\
			'-name '+imagename,'-nwlayers 1','-niter 1','-mgain 0.8','-pol i','-minuv-l '+str(self.imaging_minuv_lambda),'-use-weights-as-taper']
			wsclean_args.append('-quiet')
			if self.verbose:
				self.log_verbose.info('wsclean '+' '.join(wsclean_args)+' '+self.msname+'\n')			
			os.system(self.wsclean_start+' '+' '.join(wsclean_args)+' '+self.msname)	
			wsclean_images=glob.glob(imagename+'*image.fits')
			self.log_verbose.info('wsclean_to_casaimage(wsclean_images='+str(wsclean_images)+',casaimage_prefix=\''+imagename+'\',imagetype=\'image\',keep_wsclean_images=False)\n')
			casa_imagename=self.wsclean_to_casaimage(wsclean_images=wsclean_images,casaimage_prefix=imagename,imagetype='image',keep_wsclean_images=False)
			wsclean_models=glob.glob(imagename+'*model.fits')
			self.log_verbose.info('wsclean_to_casaimage(wsclean_images='+str(wsclean_models)+',casaimage_prefix=\''+imagename+'\',imagetype=\'model\',keep_wsclean_images=False)\n')
			casa_modelname=self.wsclean_to_casaimage(wsclean_images=wsclean_models,casaimage_prefix=imagename,imagetype='model',keep_wsclean_images=False)
			flaglist=flagmanager(vis=self.msname,mode='list')
			flaglist_keys=list(flaglist.keys())
			flaglist_keys.remove('MS')
			if len(flaglist_keys)>0:
				last_flag_key=len(flaglist_keys)-1
				last_flagversion=flaglist[last_flag_key]['name']
				if 'flagdata' in last_flagversion:  # Restore flagversion if the last flagversion is from applycal
					self.log_verbose.info('flagmanager(vis=\''+self.msname+'\',mode=\'restore\',versionname=\''+str(last_flagversion)+'\',merge=\'replace\')\n')
					self.log_verbose.info('flagmanager(vis=\''+self.msname+'\',mode=\'delete\',versionname=\''+str(last_flagversion)+'\')\n')
					flagmanager(vis=self.msname,mode='restore',versionname=last_flagversion,merge='replace')
					flagmanager(vis=self.msname,mode='delete',versionname=last_flagversion)	
		out_dict,negative_dyn_range=self.calc_dyn_range(casa_imagename,start_sigma,box_width=box_width,stokes_list=['I']) # Calculating the dynamic range of the image
		selfcal_snr=self.calc_selfcal_snr(imagename+'.image',start_sigma,len(antenna_to_use.split(',')))
		xmax_pos=int(imstat(imagename+'.image')['maxpos'][0])
		ymax_pos=int(imstat(imagename+'.image')['maxpos'][1])
		if self.verbose==False:
			os.system('rm -rf '+imagename+'.image '+imagename+'.model '+imagename+'.residual '+imagename+'.sumwt '+imagename+'.pb '+imagename+'.psf '+imagename+'.mask ')
		out_dict_keys=out_dict.keys()
		if 'NAN' in out_dict_keys:
			os.system('rm -rf casa*log')
			os.chdir(self.cwd)
			os.system('rm -rf casa*log')
			return 2,out_dict,negative_dyn_range,selfcal_snr,np.nan,np.nan	
		else:
			os.system('rm -rf casa*log')
			return 0,out_dict,negative_dyn_range,selfcal_snr,xmax_pos,ymax_pos
		
	def file_remover_and_keeper(self,num_iter,msg_code,do_bandpass=False,ref_time_chan=True):
		'''
		This function keep and remove caltables, ms, imaging related files based on the need

		Parameters
		----------
		num_iter : int 
			Number of self-calibration iteration
		msg_code : int 
			Selfcal message code
		do_bandpass : bool
			Performing bandpass or not
		ref_timechan : bool
			Reference time channel or not
		'''
		msname_str=am.splited_ms_rename(self.msname,ref_time_chan=ref_time_chan,change_msname=False)
		freqstr=os.path.basename(msname_str).split('.ms')[0].split('_freq_')[1].split('_')[0]  # Frequency string in MHz
		datestr_list=os.path.basename(msname_str).split('.ms')[0].split('_freq_')[0].split('time_')[1].split('_')
		datestr_for_file='_'.join(datestr_list[:3])+'_'+'_'.join(datestr_list[3:]) # Datetime string 
		cwd=os.getcwd()
		file_str=os.path.basename(self.msname).split('.ms')[0]+'_'+str(num_iter)+'_'+str(self.bandwidth)+'kHz_'+str(self.time_width)+'s'
		file_str_ms=os.path.basename(self.msname).split('.ms')[0]+'_'+str(num_iter)
		caltable_name=self.msname.split('.ms')[0]+'.cal' # Caltable name
		if do_bandpass==True:
			file_str_prefix='freq_'+freqstr+'_datetime_'+datestr_for_file+'_bp'
		else:
			file_str_prefix='freq_'+freqstr+'_datetime_'+datestr_for_file
		if self.verbose==True and os.path.isdir(self.mspath+'/'+file_str_prefix)==False: 
						# If verbose=True, directory to keep all intermediate images, caltables, models, residuals 
			os.mkdir(self.mspath+'/'+file_str_prefix)
		if self.verbose:
			if os.path.isdir(self.mspath+'/'+file_str_prefix+'/backup_ms')==False:
				os.makedirs(self.mspath+'/'+file_str_prefix+'/backup_ms')
			if os.path.isdir(self.mspath+'/'+file_str_prefix+'/backup_imagemodel')==False:
				os.makedirs(self.mspath+'/'+file_str_prefix+'/backup_imagemodel')
			if os.path.isdir(self.mspath+'/'+file_str_prefix+'/backup_cal')==False:
				os.makedirs(self.mspath+'/'+file_str_prefix+'/backup_cal')
		if num_iter=='dirty':
			if self.verbose and (msg_code==0):
				os.system('cp -r '+file_str+'.model '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'.model') # If Verbose=True, keep all the models
				os.system('cp -r '+file_str+'.image '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'.image') # If Verbose=True, keep all the image
			os.system('rm -rf '+file_str+'.image '+file_str+'.model '+file_str+'.residual '+file_str+'.sumwt '+file_str+'.pb '+file_str+'.psf *psf*') 
			os.chdir(self.cwd)
			os.system('rm -rf casa*log')
			return
		if num_iter==0 and (msg_code==0 or msg_code==8 or msg_code==9):
			os.system('cp -r '+caltable_name+' junk0.cal') # Copying num_iter=0 caltable to junk0.cal
			os.system('cp -r '+self.msname+' junk0.ms') # Copying num_iter=0 ms to junk0.ms
			os.system('cp -r '+file_str+'.image junk0.image') # Copying num_iter=0 image to junk0.image
			os.system('cp -r '+file_str+'.model junk0.model') # Copying num_iter=0 model to junk0.model
			if os.path.isdir(file_str+'.mask'):
				os.system('cp -r '+file_str+'.mask junk0.mask') # Copying num_iter=0 mask to junk0.mask
			os.system('cp -r '+file_str+'.residual junk0.residual') # Copying num_iter=0 residual to junk0.residual
			if self.verbose and (msg_code==0 or msg_code==8 or msg_code==9):	
				os.system('cp -r '+caltable_name+' '+self.mspath+'/'+file_str_prefix+'/backup_cal/'+file_str_ms+'.cal') # Verbose=True, keep all the caltables
				os.system('cp -r '+self.msname+' '+self.mspath+'/'+file_str_prefix+'/backup_ms/'+file_str_ms+'.ms') # If Verbose=True, keep all the ms
				os.system('cp -r '+file_str+'.model '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'.model') # If Verbose=True, keep all the models
				if os.path.isdir(file_str+'.mask'):
					os.system('cp -r '+file_str+'.mask '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'.mask') # If Verbose=True, keep all the masks
				os.system('cp -r '+file_str+'.image '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'.image') # If Verbose=True, keep all the image
				os.system('cp -r '+file_str+'.residual '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'.residual') # If Verbose=True, keep all the residuals
				if os.path.exists(self.mspath+'/mwapbcor.image'):
					os.system('mv '+self.mspath+'/mwapbcor.image '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'_mwapbcor.image')
				if os.path.exists(self.mspath+'/mwapbcor.model'):
					os.system('mv '+self.mspath+'/mwapbcor.model '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'_mwapbcor.model')
				if os.path.exists(self.mspath+'/quvcor.image'):
					os.system('mv '+self.mspath+'/quvcor.image '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'_quvcor.image')
				if os.path.exists(self.mspath+'/quvcor.model'):
					os.system('mv '+self.mspath+'/quvcor.model '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'_quvcor.model')
			os.system('rm -rf '+self.mspath+'/quvcor.model '+self.mspath+'/quvcor.image '+self.mspath+'/mwapbcor.model '+self.mspath+'/mwapbcor.image')
			os.system('rm -rf '+caltable_name)
			os.system('rm -rf '+file_str+'.image '+file_str+'.model '+file_str+'.residual '+file_str+'.sumwt '+file_str+'.pb '+file_str+'.psf '+file_str+'.mask *psf*') 
				# Removing all imaging related files
		elif num_iter==1 and (msg_code==0 or msg_code==8 or msg_code==9):
			os.system('cp -r '+caltable_name+' junk1.cal') # Copying num_iter=1 caltable to junk1.cal
			os.system('cp -r '+self.msname+' junk1.ms') # Copying num_iter=1 ms to junk1.ms
			os.system('cp -r '+file_str+'.image junk1.image') # Copying num_iter=1 image to junk1.image
			os.system('cp -r '+file_str+'.model junk1.model') # Copying num_iter=1 model to junk1.model
			if os.path.isdir(file_str+'.mask'):
				os.system('cp -r '+file_str+'.mask junk1.mask') # Copying num_iter=1 model to junk1.mask
			os.system('cp -r '+file_str+'.residual junk1.residual') # Copying num_iter=1 residual to junk1.residual
			if self.verbose and (msg_code==0 or msg_code==8 or msg_code==9):	
				os.system('cp -r '+caltable_name+' '+self.mspath+'/'+file_str_prefix+'/backup_cal/'+file_str_ms+'.cal') # Verbose=True, keep all the caltables
				os.system('cp -r '+self.msname+' '+self.mspath+'/'+file_str_prefix+'/backup_ms/'+file_str_ms+'.ms') # If Verbose=True, keep all the ms
				os.system('cp -r '+file_str+'.model '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'.model') # If Verbose=True, keep all the models
				if os.path.isdir(file_str+'.mask'):
					os.system('cp -r '+file_str+'.mask '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'.mask') # If Verbose=True, keep all the masks
				os.system('cp -r '+file_str+'.image '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'.image') # If Verbose=True, keep all the image
				os.system('cp -r '+file_str+'.residual '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'.residual') # If Verbose=True, keep all the residuals
				if os.path.exists(self.mspath+'/mwapbcor.image'):
					os.system('mv '+self.mspath+'/mwapbcor.image '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'_mwapbcor.image')
				if os.path.exists(self.mspath+'/mwapbcor.model'):
					os.system('mv '+self.mspath+'/mwapbcor.model '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'_mwapbcor.model')
				if os.path.exists(self.mspath+'/quvcor.image'):
					os.system('mv '+self.mspath+'/quvcor.image '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'_quvcor.image')
				if os.path.exists(self.mspath+'/quvcor.model'):
					os.system('mv '+self.mspath+'/quvcor.model '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'_quvcor.model')
			os.system('rm -rf '+self.mspath+'/quvcor.model '+self.mspath+'/quvcor.image '+self.mspath+'/mwapbcor.model '+self.mspath+'/mwapbcor.image')
			os.system('rm -rf '+caltable_name)
			os.system('rm -rf '+file_str+'.image '+file_str+'.model '+file_str+'.residual '+file_str+'.sumwt '+file_str+'.pb '+file_str+'.psf '+file_str+'.mask *psf*') 
				# Removing all imaging related files
		elif num_iter>1 and (msg_code==0 or msg_code==8 or msg_code==9):
			if os.path.isdir('junk1.cal'):
				os.system('rm -rf junk0.cal')
				os.system('cp -r junk1.cal junk0.cal') # Move the previous round caltable to pre-previous round for num_iter>1
			else:
				os.system('cp -r '+caltable_name+' junk0.cal') # Copying caltable to junk0.cal
			if os.path.isdir('junk1.ms'):
				os.system('rm -rf junk0.ms')
				os.system('cp -r junk1.ms junk0.ms')	# Move the previous round ms to pre-previous round for num_iter>1
			else:
				os.system('cp -r '+self.msname+' junk0.ms') # Copying ms to junk0.ms
			if os.path.isdir('junk1.model'):
				os.system('rm -rf junk0.model')
				os.system('cp -r junk1.model junk0.model') # Move the previous round model to pre-previous round for num_iter>1
			else:
				os.system('cp -r '+file_str+'.model junk0.model') # Copying model to junk0.model
			if os.path.isdir('junk1.image'):
				os.system('rm -rf junk0.image')
				os.system('cp -r junk1.image junk0.image') # Move the previous round image to pre-previous round for num_iter>1
			else:
				os.system('cp -r '+file_str+'.image junk0.image') # Copying image to junk0.image
			if os.path.isdir('junk1.mask'):
				os.system('rm -rf junk0.mask')
				os.system('cp -r junk1.mask junk0.mask') # Move the previous round mask to pre-previous round for num_iter>1
			elif os.path.isdir(file_str+'.mask'):
				os.system('cp -r '+file_str+'.mask junk0.mask') # Copying model to junk0.mask
			if os.path.isdir('junk1.residual'):
				os.system('rm -rf junk0.residual')
				os.system('cp -r junk1.residual junk0.residual') # Move the previous round residual to pre-previous round for num_iter>1
			else:
				os.system('cp -r '+file_str+'.residual junk0.residual') # Copying model to junk0.residual
			os.system('rm -rf junk1.ms junk1.cal junk1.model junk1.mask junk1.image junk1.residual *psf*')
			os.system('cp -r '+caltable_name+' junk1.cal') # Copying caltable to junk1.cal
			os.system('cp -r '+self.msname+' junk1.ms') # Copying ms to junk1.ms
			os.system('cp -r '+file_str+'.model junk1.model') # Copying model to junk1.model
			if os.path.isdir(file_str+'.mask'):
				os.system('cp -r '+file_str+'.mask junk1.mask') # Copying mask to junk1.mask
			os.system('cp -r '+file_str+'.image junk1.image') # Copying image to junk1.image
			os.system('cp -r '+file_str+'.residual junk1.residual') # Copying residual to junk1.resuidual
			if self.verbose and (msg_code==0 or msg_code==8 or msg_code==9):	
				if os.path.exists(self.mspath+'/'+file_str_prefix+'/backup_cal/'+file_str_ms+'.cal'):
					os.system('rm -rf '+self.mspath+'/'+file_str_prefix+'/backup_cal/'+file_str_ms+'.cal')
				os.system('cp -r '+caltable_name+' '+self.mspath+'/'+file_str_prefix+'/backup_cal/'+file_str_ms+'.cal') # Verbose=True, keep all the caltables
				os.system('cp -r '+self.msname+' '+self.mspath+'/'+file_str_prefix+'/backup_ms/'+file_str_ms+'.ms') # If Verbose=True, keep all the ms
				os.system('cp -r '+file_str+'.model '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'.model') # If Verbose=True, keep all the models
				if os.path.isdir(file_str+'.mask'):
					os.system('cp -r '+file_str+'.mask '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'.mask') # If Verbose=True, keep all the masks
				os.system('cp -r '+file_str+'.image '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'.image') # If Verbose=True, keep all the image
				os.system('cp -r '+file_str+'.residual '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'.residual') # If Verbose=True, keep all the residuals
				if os.path.exists(self.mspath+'/mwapbcor.image'):
					os.system('mv '+self.mspath+'/mwapbcor.image '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'_mwapbcor.image')
				if os.path.exists(self.mspath+'/mwapbcor.model'):
					os.system('mv '+self.mspath+'/mwapbcor.model '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'_mwapbcor.model')
				if os.path.exists(self.mspath+'/quvcor.image'):
					os.system('mv '+self.mspath+'/quvcor.image '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'_quvcor.image')
				if os.path.exists(self.mspath+'/quvcor.model'):
					os.system('mv '+self.mspath+'/quvcor.model '+self.mspath+'/'+file_str_prefix+'/backup_imagemodel/'+file_str+'_quvcor.model')
			os.system('rm -rf '+self.mspath+'/quvcor.model '+self.mspath+'/quvcor.image '+self.mspath+'/mwapbcor.model '+self.mspath+'/mwapbcor.image')
			os.system('rm -rf '+file_str+'.image '+file_str+'.model '+file_str+'.residual '+file_str+'.sumwt '+file_str+'.pb '+file_str+'.psf '+file_str+'.mask *psf*') 
				# Removing all imaging related files
		os.chdir(self.cwd)
		os.system('rm -rf casa*log')
		return

	def initial_bpass(self,modelname,num_iter,rms_thresh,ref_ant,minsnr,calmode,calibrator_caltable=[]):
		'''
		Function to perform and apply bandpass calibration

		Parameters
		----------
		modelname : str 
			Name of the initial model
		num_iter : int
			Selfcal iteration number
		rms_thresh : list
			List of n-sigma value for rms based flagging (Do not go below 6)
		ref_ant : int 
			Reference antenna number
		calmode : str
			'p' for phase-only or 'ap' for amplitude-phase calibration
		calibrator_caltable : list 
			List of any previous caltables
		Returns
		-------
		str
			Name of the initial bandpass table
		'''
		AM=am.AccessMS(self.msname)
		bandwidth=AM.calc_bandwidth()
		cent_freq=AM.calc_meanfreq()
		header=imhead(imagename=modelname,mode='list')
		caltable_name=self.msname.split('.ms')[0]+'.cal' # Caltable name
		if header['ctype1']=='Frequency':
			header_key='crval1'
			header_key_1='cdelt1'
		elif header['ctype2']=='Frequency':
			header_key='crval2'
			header_key_1='cdelt2'
		elif header['ctype3']=='Frequency':
			header_key='crval3'
			header_key_1='cdelt3'
		elif header['ctype4']=='Frequency':
			header_key='crval4'
			header_key_1='cdelt4'
		self.log_verbose.info('imhead(imagename=\''+modelname+'\',mode=\'put\',hdkey=\''+header_key+'\',hdvalue={\'value\':'+str(cent_freq)+',\'unit\':\'Hz\'},verbose=False)\n')
		imhead(imagename=modelname,mode='put',hdkey=header_key,hdvalue={'value':cent_freq,'unit':'Hz'},verbose=False)
		self.log_verbose.info('imhead(imagename=\''+modelname+'\',mode=\'put\',hdkey=\''+header_key_1+'\',hdvalue={\'value\':'+str(bandwidth)+',\'unit\':\'Hz\'},verbose=False)\n')
		imhead(imagename=modelname,mode='put',hdkey=header_key_1,hdvalue={'value':bandwidth,'unit':'Hz'},verbose=False)
		self.log_verbose.info('delmod(vis=\''+self.msname+'\',scr=True)\n')
		delmod(vis=self.msname,scr=True)
		self.log_verbose.info('ft(vis=\''+self.msname+'\',model=\''+modelname+'\',usescratch=True)\n')
		ft(vis=self.msname,model=modelname,usescratch=True)
		self.log_verbose.info('bpass_solver(\''+caltable_name+'\',spw=\'\',timerange=\'\',calmode=\''+calmode+'\',uvrange=\''+self.calib_uvrange_lambda+'\',solnorm=True,normtype=\'median\',refant=\''+\
							str(ref_ant)+'\',minsnr='+str(minsnr)+',solmode=\'R\',rmsthresh=[],gaintable='+str(calibrator_caltable)+')\n')
		caltable_name=self.bpass_solver(caltable_name,spw='',timerange='',calmode=calmode,uvrange=self.imaging_uvrange,solnorm=True,normtype='median',refant=str(ref_ant),minsnr=minsnr,\
						solmode='R',rmsthresh=[],gaintable=calibrator_caltable) # Perform bandpass calibration
		calibrator_caltable.append(caltable_name)
		self.log_verbose.info('applycal(vis=\''+self.msname+'\',gaintable='+str(calibrator_caltable)+',applymode=\'calflag\',calwt=[True],flagbackup=True)\n')
		applycal(vis=self.msname,gaintable=calibrator_caltable,applymode='calflag',calwt=[True],flagbackup=True)
		os.system('rm -rf casa*log')
		return caltable_name

	def bpass_solver(self,caltable,spw='',timerange='',calmode='ap',uvrange='',solnorm=True,normtype='mean',refant='1',minsnr=3,solint='inf',solmode='',rmsthresh=[],gaintable=[]):
		'''
		Function to solve bandpass phase-only or amplitude-phase

		Parameters
		----------
		caltable : str 
			Name of the caltable
		spw : str 
			Spectral window range
		timerange : str
			CASA timerange
		calmode : str 
			'p' for phase-only or 'ap' for amplitude-phase calibration
		uvrange : str 
			UV-range for calibration
		solnorm : bool
			Normalise solution or not)
		normtype : str
			Normalisation type (mean or median)
		refant : str 
			Reference antenna
		minsnr : float 
			Minimum gain SNR
		solint : str
			Solution interval
		solmode : str 
			'R' for robust calibration
		rmsthresh : list 
			List of n-sigma values for rms based flagging (Do not put below 6)
		gaintable : list
			List of any previous caltables
		Returns 
		-------- 
		str
			Name of the caltable
		'''
		cb=calibrater()
		if solmode!='R' or len(rmsthresh)==0:
			cb.open(self.msname)
			cb.selectvis(spw=spw,uvrange=uvrange,time=timerange)
			if len(gaintable)!=0:
				for gt in gaintable:
					cb.setapply(table=gt)
			cb.setsolve(type='B',t=solint,refant=str(refant),apmode=calmode,table=caltable,append=False,minsnr=minsnr,solnorm=solnorm,normtype=normtype)
			cb.solve()
			cb.close()
		elif solmode=='R' and len(rmsthresh)!=0:
			os.system('cp -r '+self.msname+' '+self.msname.split('.ms')[0]+'_temp.ms')
			visname=self.msname.split('.ms')[0]+'_temp.ms'
			for rms in rmsthresh:
				c=0
				self.log_verbose.info('Calibrating and flagging on threshold :'+str(rms)+' sigma\n')
				while c==0:
					cb.open(visname)
					cb.selectvis(spw=spw,uvrange=uvrange,time=timerange)
					if len(gaintable)!=0:
						for gt in gaintable:
							cb.setapply(table=gt)
					cb.setsolve(type='B',t=solint,refant=str(refant),apmode=calmode,table=caltable,append=False,minsnr=minsnr,solnorm=solnorm,normtype=normtype)
					cb.solve()
					cb.close()
					applycal(vis=visname,gaintable=gaintable+[caltable],applymode='calflag',calwt=[True],flagbackup=False)
					num_flag=fg.flagger(visname,float(rms))
					if num_flag==0:
						c=1
			os.system('rm -rf '+visname)
		os.system('rm -rf casa*log')
		return	caltable

	def cal_solar_phaseshift(self,imagename,sigma=10):
		'''
		Calculate the different between solar center and phase center of the image

		Parameters
		----------
		imagename : str 
			Name of the image
		sigma : float 
			Threshold for model fitting (default =10)
		Returns
		-------
		float
			New RA of the phasecenter in degree
		float
			New DEC of the phasecenter in degree
		bool
			Whther phase shift required or not. NOt required if less than image pixel size
		'''
		AM=am.AccessMS(self.msname)
		radec_str,radeg,decdeg=AM.get_phasecenter()
		cellsize=np.rad2deg(np.abs(imhead(imagename=imagename,mode='list')['cdelt1'])) # In degree
		image_path=os.path.dirname(os.path.abspath(imagename))
		os.system('rm -rf '+image_path+'/phaseshift*')	
		os.system('cp -r '+imagename+' '+image_path+'/phaseshift.image')
		modelname=imagename.split('.image')[0]+'.model'
		rms=imstat(imagename=image_path+'/phaseshift.image',box=self.rms_box)['rms'][0]
		if os.path.isfile(image_path+'/phaseshift.fits')==True:
			os.system('rm -rf '+image_path+'/phaseshift.fits')
		exportfits(imagename=image_path+'/phaseshift.image',fitsimage=image_path+'/phaseshift.fits',dropdeg=True,dropstokes=True)
		os.system('rm -rf '+image_path+'/phaseshift.image')
		f=fits.open(image_path+'/phaseshift.fits')
		data=fits.getdata(image_path+'/phaseshift.fits')
		data[data<=sigma*rms]=0
		data[data>sigma*rms]=1
		ndim=data.ndim
		if ndim>2:
			if ndim==3:
				data=data[0,:,:]
			elif ndim==4:
				data=data[0,0,:,:]
		imsize=data.shape[0]
		circular_mask=self.create_circular_mask(imsize,imsize,center=(int(imsize/2),int(imsize/2)), radius=int(3/cellsize))	
		data[~circular_mask]=0				
		cy, cx = ndi.center_of_mass(data)
		w=WCS(image_path+'/phaseshift.fits')
		try:
			result=w.array_index_to_world(0,int(cy),int(cx))
			x_cen=result[0].ra.deg
			y_cen=result[0].dec.deg
		except:
			result=w.array_index_to_world(int(cy),int(cx))
			x_cen=result.ra.deg
			y_cen=result.dec.deg
		ra=float(x_cen)
		dec=float(y_cen)
		fits.writeto(image_path+'/phaseshift.fits',data=data,header=f[0].header,overwrite=True)
		os.system('rm -rf '+image_path+'/phaseshift*')
		if np.sqrt((ra-radeg)**2+(dec-decdeg)**2)<((float(self.cellsize))/3600.0):
			os.system('rm -rf casa*log')
			return radeg,decdeg,False
		else:
			os.system('rm -rf casa*log')
			return ra,dec,True

	def shift_phasecenter(self,imagename,ra,dec):
		'''
		Function to shift image reference RA DEC to a certain value (Only CASA image)

		Parameters
		----------
		imagename : str 
			Name of the image
		ra : float 
			RA in degree
		dec : float
			DEC in degree
		Returns
		-------
		int
			Success code (0 or 1)
		'''
		cwd=os.getcwd()
		AM=am.AccessMS(self.msname)
		radec_str,radeg,decdeg=AM.get_phasecenter()
		IB=B.ImageBasic(self.msname)
		image_path=os.path.dirname(os.path.abspath(imagename))
		if os.path.isfile(image_path+'/wcs_model.fits')==True:
			os.system('rm -rf '+image_path+'/wcs_model.fits')
		if os.path.isdir(image_path+'/I.model'):
			os.system('rm -rf '+image_path+'/I.model')
		if imagename[-1]=='/':
			imagename=imagename[:-1]
		try:
			imsubimage(imagename=imagename,outfile=image_path+'/I.model',stokes='I',dropdeg=False)
		except:
			imsubimage(imagename=imagename,outfile=image_path+'/I.model',stokes='XX',dropdeg=False)
		exportfits(imagename=image_path+'/I.model',fitsimage=image_path+'/wcs_model.fits',dropdeg=True,dropstokes=True)
		w=WCS(image_path+'/wcs_model.fits')
		pix=np.nanmean(w.all_world2pix(np.array([[ra,dec],[ra,dec]]),0),axis=0)
		ra_pix=int(pix[0])
		dec_pix=int(pix[1])
		imhead(imagename=imagename,mode='put',hdkey='CRPIX1',hdvalue=str(ra_pix))
		imhead(imagename=imagename,mode='put',hdkey='CRPIX2',hdvalue=str(dec_pix))
		self.log_verbose.info('Image phase center shifted to , RA : '+str(radec_str[0])+', DEC : '+str(radec_str[1])+'\n')
		os.system('rm -rf casa*log '+image_path+'/I.model '+image_path+'/wcs_model.fits')
		os.chdir(cwd)
		return 0

	def image_source_true_loc(self,outdir,do_bandpass=False,weight='briggs',robust=0.5):
		'''
		Function to make quick image for finding source true location with respect to reference time and channel
	
		Parameters
		----------
		Outdir : str 
			Output directory
		do_bandpass : bool 
			Source true location for bandpass self-calibration or not
		weight : str
			Visibility weighting during imaging , 'uniform', 'natural' or 'briggs'. Default is 'briggs'
		robust : float
			Robust parameter for briggs weighting, default : 0.5 
		Returns
		-------
		str
			Output imagename
		'''
		os.system('rm -rf test_loc*')
		robust=robust*2
		poltclean(vis=self.msname,imagename='test_loc',selectdata=True,startmodel='',stokes='I',antenna='',imsize=[self.imsize],cell=str(self.cellsize)+'arcsec',\
				niter=0,gain=0.5,threshold=['10Jy'],deconvolver='multiscale',scales=self.multiscale_scales,uvtaper=self.uvtaper,weighting=weight,robust=robust,\
				interactive=False,mask='')
		maxpix=imstat(imagename='test_loc.image')['max'][0]
		thresh=str(maxpix/5)+'Jy'
		os.system('rm -rf test_loc*')
		if do_bandpass==True:
			poltclean(vis=self.msname,imagename='test_loc',selectdata=True,startmodel='',stokes='I',antenna='',imsize=[self.imsize],cell=str(self.cellsize)+'arcsec',\
					niter=100000,gain=0.5,threshold=[thresh],specmode='cube',deconvolver='multiscale',scales=self.multiscale_scales,uvtaper=self.uvtaper,weighting=weight,\
					robust=robust,interactive=False,mask='')
			self.log_verbose.info('poltclean(vis=\''+self.msname+'\',imagename=\'test_loc\',selectdata=True,startmodel=\'\',stokes=\'I\',antenna=\'\',imsize=['\
					+str(self.imsize)+'],cell=\''+str(self.cellsize)+'arcsec\',niter=100000,gain=0.5,threshold=[\''+\
		str(thresh)+'\'],specmode=\'cube\',deconvolver=\'multiscale\',scales='+str(self.multiscale_scales)+',uvtaper=\''+str(self.uvtaper)+'\',weighting=\''+weight+'\',robust='+\
					str(robust)+',interactive=False,mask=\'\')\n')
		else:
			poltclean(vis=self.msname,imagename='test_loc',selectdata=True,startmodel='',stokes='I',antenna='',imsize=[self.imsize],cell=str(self.cellsize)+'arcsec',\
					niter=100000,gain=0.5,threshold=[thresh],deconvolver='multiscale',scales=self.multiscale_scales,uvtaper=self.uvtaper,weighting=weight,\
					robust=robust,interactive=False,mask='')
			self.log_verbose.info('poltclean(vis=\''+self.msname+'\',imagename=\'test_loc\',selectdata=True,startmodel=\'\',stokes=\'I\',antenna=\'\',imsize=['\
					+str(self.imsize)+'],cell=\''+str(self.cellsize)+'arcsec\',niter=100000,gain=0.5,threshold=[\''+str(thresh)+'\'],deconvolver=\'multiscale\',scales='+\
					str(self.multiscale_scales)+',uvtaper=\''+str(self.uvtaper)+'\',weighting=\''+weight+'\',robust='+str(robust)+',interactive=False,mask=\'\')\n')
		if os.path.isdir(outdir)==False:
			os.makedirs(outdir)
		if outdir[-1]=='/':
			outdir=outdir[:-1]
		imagename=outdir+'/'+os.path.basename(self.msname).split('.ms')[0]+'_true_loc.fits' # Imagename prefix
		if os.path.isfile(imagename):
			os.system('rm -rf '+imagename)
		exportfits(imagename='test_loc.image',fitsimage=imagename)
		os.system('rm -rf test_loc.* casa*log')
		return imagename

	def remove_model_negative(self,imagename,modelname,sigma=10,overwrite=False):
		'''
		Function to remove negatives from model image

		Parameters
		----------
		imagename : str 
			Name of the image
		modelname : str 
			Name of the model
		sigma : float 
			Sigma value for thresholding
		overwrite : bool 
			Overwrite the model image or not
		Returns
		-------
		str					
			Model imagename without negatives
		'''
		ia=image()
		if overwrite==False:
			if os.path.isdir(modelname+'.nonegative')==True:
				os.system('rm -rf '+modelname+'.nonegative')
			os.system('cp -r '+modelname+' '+modelname+'.nonegative')
			modelname=modelname+'.nonegative'
		try:
			rmsi=imstat(imagename=imagename,box=self.rms_box,stokes='I')['rms'][0]
			stokes='I'
		except:
			try:
				rmsxx=imstat(imagename=imagename,box=self.rms_box,stokes='XX')['rms'][0]
				rmsyy=imstat(imagename=imagename,box=self.rms_box,stokes='YY')['rms'][0]
				stokes='XXYY'
			except Exception as e:
				self.log_verbose.error('Either Stokes \'I\' or \'XX,YY\' not present. Error message : '+str(e)+'\n')
				return None
		if stokes=='I':
			ia.open(imagename)
			data=ia.getchunk()
			ia.close()
			datai=data[:,:,0,:]
			sigma_pos=np.where(datai>=(sigma*rmsi))
			sigma_pos1=np.where(datai<(sigma*rmsi))
			ia.open(modelname)
			data=ia.getchunk()
			data_copy=copy.deepcopy(data)
			datai_copy=data_copy[:,:,0,:]
			datai[sigma_pos]=0
			pos=np.where(datai<0)
			datai_copy[pos]=0
			data_copy[:,:,0,:]=datai_copy
			ia.putchunk(data_copy)
			ia.done()
		else:
			ia.open(imagename)
			data=ia.getchunk()
			ia.close()
			dataxx=data[:,:,0,:]
			datayy=data[:,:,1,:]
			sigma_posxx=np.where(dataxx>=(sigma*rmsxx))
			sigma_posxx1=np.where(dataxx<(sigma*rmsxx))
			sigma_posyy=np.where(datayy>=(sigma*rmsyy))
			sigma_posyy1=np.where(datayy<(sigma*rmsyy))
			ia.open(modelname)
			data=ia.getchunk()
			data_copy=copy.deepcopy(data)
			dataxx_copy=data_copy[:,:,0,:]
			datayy_copy=data_copy[:,:,1,:]
			dataxx[sigma_posxx]=0
			datayy[sigma_posyy]=0
			posxx=np.where(dataxx<0)
			dataxx_copy[posxx]=0
			data_copy[:,:,0,:]=dataxx_copy
			posyy=np.where(datayy<0)
			datayy_copy[posyy]=0
			data_copy[:,:,1,:]=datayy_copy
			ia.putchunk(data_copy)
			ia.done()
		return modelname

	def remove_outside_model(self,imagename,modelname,box_width=5,overwrite=False):
		'''
		Function to remove negatives from model image

		Parameters
		----------
		imagename : str 
			Name of the image
		modelname : str 
			Name of the model
		box_width : float
			Box width in degree outside which no emissions are picked up
		overwrite : bool 
			Overwrite the model image or not
		Returns
		-------
		str					
			Model imagename
		'''
		ia=image()
		if modelname[-1]=='/':
			modelname=modelname[:-1]
		if imagename[-1]=='/':
			imagename=imagename[:-1]		
		if overwrite==False:
			if os.path.isdir(modelname+'.nooutside')==True:
				os.system('rm -rf '+modelname+'.nooutside')
			os.system('cp -r '+modelname+' '+modelname+'.nooutside')
			modelname=modelname+'.nooutside'
		imsize=imhead(imagename,mode='list')['shape'][0]
		cellsize=np.rad2deg(np.abs(imhead(imagename=imagename,mode='list')['cdelt1'])) # In degree
		circular_mask=self.create_circular_mask(imsize,imsize, center=(int(imsize/2),int(imsize/2)), radius=int(box_width/cellsize))
		ia.open(modelname)
		data=ia.getchunk()
		data_copy=copy.deepcopy(data)
		for i in range(data_copy.shape[-1]):
			for j in range(data_copy.shape[-2]):
				data_copy[:,:,j,i][~circular_mask]=0
		ia.putchunk(data_copy)
		ia.done()
		return modelname

	def wsclean_to_casaimage(self,wsclean_images=[],casaimage_prefix='CASA',imagetype='image',keep_wsclean_images=True): #TODI : Include other stokes mode
		'''
		Function to convert WSClean image in CASA image (Stokes modes : 'IQUV', 'XXYY', 'I', 'QU', 'IV','IQ')

		Parameters
		----------
		wsclean_images : list 
			List of WSClean images
		casaimage_prefix : str 
			Output CASA image name prefix (default : 'CASA\_')
		imagetype : str 
			'image', 'model', 'residual' or 'dirty' (default : 'image')
		keep_wsclean_images : bool 
			Keep the WSClean images or not
		Returns
		-------
		str
			Output CASA imagename
		'''
		stokes=[]
		wsclean_images=sorted(wsclean_images)
		for i in wsclean_images:
			name_split=os.path.basename(i).split('.fits')[0].split('-')
			if len(name_split)>=3:
				if name_split[-2] not in stokes:
					stokes.append(name_split[-2])
			else:
				if 'I' not in stokes:
					stokes.append('I')
		stokes=sorted(stokes)
		if stokes!=['I','Q','U','V'] and stokes!=['XX','YY'] and stokes!=['I','V'] and stokes!=['Q','U'] and stokes!=['I'] and stokes!=['I','Q']:
			if self.verbose:
				self.log_verbose.info('Stokes axes are not in \'IQUV\',\'I\',\'QU\',\'IV\',\'IQ\' or \'XX,YY\'\n')
			else:
				print('Stokes axes are not in \'IQUV\',\'I\',\'QU\',\'IV\',\'IQ\' or \'XX,YY\'\n')
		elif stokes==['I']:
			imagename=casaimage_prefix+'.'+imagetype
			if os.path.isdir(imagename):
				os.system('rm -rf '+imagename)
			importfits(fitsimage=wsclean_images[0],imagename=imagename,defaultaxes=True,defaultaxesvalues=['ra','dec','stokes','freq'])
			if keep_wsclean_images==False:
				for i in wsclean_images:
					os.system('rm -rf '+i)
			return imagename
		elif stokes==['I','V']:
			for i in wsclean_images:
				if '-I-' in i:
					data=fits.getdata(i)
					header=fits.getheader(i)
				else:
					data=np.append(data,fits.getdata(i),axis=0)
			header['NAXIS4']=2.
			header['CRVAL4']=1.
			header['CDELT4']=3.
			fits.writeto(casaimage_prefix+'_IV.fits',data=data,header=header,overwrite=True)
			imagename=casaimage_prefix+'.'+imagetype
			if os.path.isdir(imagename):
				os.system('rm -rf '+imagename)
			importfits(fitsimage=casaimage_prefix+'_IV.fits',imagename=imagename,defaultaxes=True,defaultaxesvalues=['ra','dec','stokes','freq'])
			os.system('rm -rf '+casaimage_prefix+'_IV.fits')
			if keep_wsclean_images==False:
				for i in wsclean_images:
					os.system('rm -rf '+i)
			return imagename
		elif stokes==['I','Q','U','V']:
			for i in wsclean_images:
				if '-I-' in i:
					data=fits.getdata(i)
					header=fits.getheader(i)
				else:
					data=np.append(data,fits.getdata(i),axis=0)
			header['NAXIS4']=4.
			header['CRVAL4']=1.
			header['CDELT4']=1.
			fits.writeto(casaimage_prefix+'_IQUV.fits',data=data,header=header,overwrite=True)
			imagename=casaimage_prefix+'.'+imagetype
			if os.path.isdir(imagename):
				os.system('rm -rf '+imagename)
			importfits(fitsimage=casaimage_prefix+'_IQUV.fits',imagename=imagename,defaultaxes=True,defaultaxesvalues=['ra','dec','stokes','freq'])
			os.system('rm -rf '+casaimage_prefix+'_IQUV.fits')
			if keep_wsclean_images==False:
				for i in wsclean_images:
					os.system('rm -rf '+i)
			return imagename
		elif stokes==['XX','YY']:
			for i in wsclean_images:
				if '-XX-' in i:
					data=fits.getdata(i)
					header=fits.getheader(i)
				else:
					data=np.append(data,fits.getdata(i),axis=0)
			header['NAXIS4']=2.
			header['CRVAL4']=-5.
			header['CDELT4']=-1.
			fits.writeto(casaimage_prefix+'_XXYY.fits',data=data,header=header,overwrite=True)
			imagename=casaimage_prefix+'.'+imagetype
			if os.path.isdir(imagename):
				os.system('rm -rf '+imagename)
			importfits(fitsimage=casaimage_prefix+'_XXYY.fits',imagename=imagename,defaultaxes=True,defaultaxesvalues=['ra','dec','stokes','freq'])
			os.system('rm -rf '+casaimage_prefix+'_XXYY.fits')
			if keep_wsclean_images==False:
				for i in wsclean_images:
					os.system('rm -rf '+i)
			return imagename
		elif stokes==['Q','U']:
			for i in wsclean_images:
				if '-Q-' in i:
					data=fits.getdata(i)
					header=fits.getheader(i)
				else:
					data=np.append(data,fits.getdata(i),axis=0)
			header['NAXIS4']=2.
			header['CRVAL4']=2.
			header['CDELT4']=1.
			fits.writeto(casaimage_prefix+'_QU.fits',data=data,header=header,overwrite=True)
			imagename=casaimage_prefix+'.'+imagetype
			if os.path.isdir(imagename):
				os.system('rm -rf '+imagename)
			importfits(fitsimage=casaimage_prefix+'_QU.fits',imagename=imagename,defaultaxes=True,defaultaxesvalues=['ra','dec','stokes','freq'])
			os.system('rm -rf '+casaimage_prefix+'_QU.fits')
			if keep_wsclean_images==False:
				for i in wsclean_images:
					os.system('rm -rf '+i)
			return imagename
		elif stokes==['I','Q']:
			for i in wsclean_images:
				if '-I-' in i:
					data=fits.getdata(i)
					header=fits.getheader(i)
				else:
					data=np.append(data,fits.getdata(i),axis=0)
			header['NAXIS4']=2.
			header['CRVAL4']=1.
			header['CDELT4']=1.
			fits.writeto(casaimage_prefix+'_IQ.fits',data=data,header=header,overwrite=True)
			imagename=casaimage_prefix+'.'+imagetype
			if os.path.isdir(imagename):
				os.system('rm -rf '+imagename)
			importfits(fitsimage=casaimage_prefix+'_IQ.fits',imagename=imagename,defaultaxes=True,defaultaxesvalues=['ra','dec','stokes','freq'])
			os.system('rm -rf '+casaimage_prefix+'_IQ.fits')
			if keep_wsclean_images==False:
				for i in wsclean_images:
					os.system('rm -rf '+i)
			return imagename

	def casaimage_to_wsclean(self,casaimage='',wsclean_imagename='',imagetype='image',wsclean_imsize=1024,keep_casaimage=True,stokes='I'):
		# TODO : include circular basis mores also
		'''
		Function to convert CASA image to WSClean image 
		(Stokes modes : 'I', 'IQUV', 'IQ', 'XXYY')
		Parameters
		----------
		casaimage : str
			Name of the CASA image or model
		wsclean_imagename : str
			WSClean imagename prefix
		imagetype : str
			Image type, 'image', 'model' or 'residual'
		wsclean_imsize : int
			WSClean image size, CASA image size should be larger than that
		keep_casaimage : bool
			Keep CASA image or not
		stokes : str 
			Stokes axis in image ('I', 'IQUV', 'IQ', 'XXYY')
		Returns
		--------
		list
			Output WSClean images
		int
			WSClean image size
		list
			WSClean image stokes axes
		'''
		if casaimage[-1]=='/':
			casaimage=casaimage[:-1]
		if casaimage=='' or wsclean_imagename=='':
			if self.verbose:
				self.log_verbose.error('CASA image name or WSClean imagename is not given.\n')
			else:
				print('CASA image name or WSClean imagename is not given.\n')
			return [],0,[]
		if imagetype!='image' and imagetype!='model' and imagetype!='residual':
			if self.verbose:
				self.log_verbose.error('Image type should be \'image\', \'model\' or \'residual\'\n')
			else:
				print('Image type should be \'image\', \'model\' or \'residual\'\n')
			return [],0,[]
		subimages=[]
		outimages=[]
		outstokes=[]
		x_cen=int(imhead(casaimage,mode='list')['crpix1'])
		y_cen=int(imhead(casaimage,mode='list')['crpix2'])
		casa_imsize=max([2*x_cen,2*y_cen])
		if casa_imsize<wsclean_imsize:
			if self.verbose:
				self.log_verbose.error('Requested WSClean image size : '+str(wsclean_imsize)+' is more than CASA image size : '+str(casa_imsize)+'\n')
			else:
				print('Requested WSClean image size : '+str(wsclean_imsize)+' is more than CASA image size : '+str(casa_imsize)+'\n')
			os.system('rm -rf casa*log *.stokes*')
			return [],0,[]
		else:
			box=str(int(x_cen-wsclean_imsize/2))+','+str(int(y_cen-wsclean_imsize/2))+','+str(int(x_cen+wsclean_imsize/2))+','+str(int(y_cen+wsclean_imsize/2))
		if wsclean_imsize==casa_imsize:
			box=''
		try:
			if stokes=='IQUV':
				subimages=[]
				if self.verbose:
					self.log_verbose.info('imsubimage(imagename=\''+casaimage+'\',outfile=\''+casaimage+'.stokesI\',stokes=\'I\',box=\''+box+'\')\n')
				if os.path.exists(casaimage+'.stokesI'):
					os.system('rm -rf '+casaimage+'.stokesI')
				imsubimage(imagename=casaimage,outfile=casaimage+'.stokesI',stokes='I',box=box)
				subimages.append(casaimage+'.stokesI')
				if self.verbose:
					self.log_verbose.info('imsubimage(imagename=\''+casaimage+'\',outfile=\''+casaimage+'.stokesQ\',stokes=\'Q\',box=\''+box+'\')\n')
				if os.path.exists(casaimage+'.stokesQ'):
					os.system('rm -rf '+casaimage+'.stokesQ')
				imsubimage(imagename=casaimage,outfile=casaimage+'.stokesQ',stokes='Q',box=box)
				subimages.append(casaimage+'.stokesQ')
				if self.verbose:
					self.log_verbose.info('imsubimage(imagename=\''+casaimage+'\',outfile=\''+casaimage+'.stokesU\',stokes=\'U\',box=\''+box+'\')\n')
				if os.path.exists(casaimage+'.stokesU'):
					os.system('rm -rf '+casaimage+'.stokesU')
				imsubimage(imagename=casaimage,outfile=casaimage+'.stokesU',stokes='U',box=box)
				subimages.append(casaimage+'.stokesU')
				if self.verbose:
					self.log_verbose.info('imsubimage(imagename=\''+casaimage+'\',outfile=\''+casaimage+'.stokesV\',stokes=\'V\',box=\''+box+'\')\n')
				if os.path.exists(casaimage+'.stokesV'):
					os.system('rm -rf '+casaimage+'.stokesV')
				imsubimage(imagename=casaimage,outfile=casaimage+'.stokesV',stokes='V',box=box)
				subimages.append(casaimage+'.stokesV')
			elif stokes=='XXYY':
				subimages=[]
				if self.verbose:
					self.log_verbose.info('imsubimage(imagename=\''+casaimage+'\',outfile=\''+casaimage+'.stokesXX\',stokes=\'XX\',box=\''+box+'\')\n')
				if os.path.exists(casaimage+'.stokesXX'):
					os.system('rm -rf '+casaimage+'.stokesXX')
				imsubimage(imagename=casaimage,outfile=casaimage+'.stokesXX',stokes='XX',box=box)
				subimages.append(casaimage+'.stokesXX')
				if self.verbose:
					self.log_verbose.info('imsubimage(imagename=\''+casaimage+'\',outfile=\''+casaimage+'.stokesYY\',stokes=\'YY\',box=\''+box+'\')\n')
				if os.path.exists(casaimage+'.stokesYY'):
					os.system('rm -rf '+casaimage+'.stokesYY')
				imsubimage(imagename=casaimage,outfile=casaimage+'.stokesYY',stokes='YY',box=box)
				subimages.append(casaimage+'.stokesYY')
			elif stokes=='IQ':
				subimages=[]
				if self.verbose:
					self.log_verbose.info('imsubimage(imagename=\''+casaimage+'\',outfile=\''+casaimage+'.stokesI\',stokes=\'I\',box=\''+box+'\')\n')
				if os.path.exists(casaimage+'.stokesI'):
					os.system('rm -rf '+casaimage+'.stokesI')
				imsubimage(imagename=casaimage,outfile=casaimage+'.stokesI',stokes='I',box=box)
				subimages.append(casaimage+'.stokesI')
				if self.verbose:
					self.log_verbose.info('imsubimage(imagename=\''+casaimage+'\',outfile=\''+casaimage+'.stokesQ\',stokes=\'Q\',box=\''+box+'\')\n')
				if os.path.exists(casaimage+'.stokesQ'):
					os.system('rm -rf '+casaimage+'.stokesQ')
				imsubimage(imagename=casaimage,outfile=casaimage+'.stokesQ',stokes='Q',box=box)
				subimages.append(casaimage+'.stokesQ')
			elif stokes=='I':
				subimages=[]
				if self.verbose:
					self.log_verbose.info('imsubimage(imagename=\''+casaimage+'\',outfile=\''+casaimage+'.stokesI\',stokes=\'I\',box=\''+box+'\')\n')
				if os.path.exists(casaimage+'.stokesI'):
					os.system('rm -rf '+casaimage+'.stokesI')
				imsubimage(imagename=casaimage,outfile=casaimage+'.stokesI',stokes='I',box=box)
				subimages.append(casaimage+'.stokesI')
		except Exception as e:	
			if self.verbose:
				self.log_verbose.error('Stokes axes are not '+stokes+'. Error message : '+str(e)+'\n')
			else:
				print('Stokes axes are not '+stokes+'\n')
			os.system('rm -rf casa*log *.stokes*')
			return [],0,[]
		for i in subimages:
			stokes=os.path.basename(i).split('.stokes')[-1]
			outstokes.append(stokes)
			if stokes=='I':
				if len(subimages)==1:
					if os.path.isfile(wsclean_imagename+'-'+imagetype+'.fits')==True:
						os.system('rm -rf '+wsclean_imagename+'-'+imagetype+'.fits')
					exportfits(imagename=i,fitsimage=wsclean_imagename+'-'+imagetype+'.fits',stokeslast=True,history=False)
					outimages.append(wsclean_imagename+'-'+imagetype+'.fits')
				else:
					if os.path.isfile(wsclean_imagename+'-'+str(stokes)+'-'+imagetype+'.fits')==True:
						os.system('rm -rf '+wsclean_imagename+'-'+str(stokes)+'-'+imagetype+'.fits')
					exportfits(imagename=i,fitsimage=wsclean_imagename+'-'+str(stokes)+'-'+imagetype+'.fits',stokeslast=True,history=False)
					outimages.append(wsclean_imagename+'-'+str(stokes)+'-'+imagetype+'.fits')
			else:
				if os.path.exists(wsclean_imagename+'-'+str(stokes)+'-'+imagetype+'.fits')==True:
					os.system('rm -rf '+wsclean_imagename+'-'+str(stokes)+'-'+imagetype+'.fits')
				exportfits(imagename=i,fitsimage=wsclean_imagename+'-'+str(stokes)+'-'+imagetype+'.fits',stokeslast=True,history=False)
				outimages.append(wsclean_imagename+'-'+str(stokes)+'-'+imagetype+'.fits')
		os.system('rm -rf casa*log *.stokes*')
		if keep_casaimage==False:
			os.system('rm -rf '+casaimage)
		return outimages,wsclean_imsize,outstokes

	def subtract_leakage_surface(self,imagename,modelname,sigma=10,do_fluxcal=False):
		'''
		Function to subtract quadratic leakage surface

		Parameters
		----------
		imagename : str 
			Name of the image
		modelname : str 
			Name of the model
		sigma : float 
			N-sigma value above which any emission is considered to be real
		do_fluxcal : bool
			Perform polynomial based flux calibration (See details Kansabanik et al. 2021, submitted to ApJ)
		Returns
		-------
		str
			Leakage surface subtracted image
		str 
			Leakage surface subtracted model
		float 
			Stokes Q fractional change
		float
			Stokes U fractional change
		'''
		if sigma<5:
			sigma=5
		self.log_verbose.info('Correcting Stokes I to Stokes Q, U leakage surface.\n')
		outfile_path=os.path.dirname(os.path.realpath(imagename))
		os.system('rm -rf '+outfile_path+'/I.* '+outfile_path+'/I_*')
		imsubimage(imagename=imagename,outfile=outfile_path+'/I.image',stokes='I',dropdeg=False)
		if do_fluxcal==True:
			fluxcal_image=self.mwa_solar_fluxcal(imagename=outfile_path+'/I.image',outfile=outfile_path+'/I_fluxcal.image')
		else:
			fluxcal_image=outfile_path+'/I.image'
		major=imhead(imagename=fluxcal_image)['restoringbeam']['major']['value'] # In arcsec
		minor=imhead(imagename=fluxcal_image)['restoringbeam']['minor']['value'] # In arcsec
		freq=imhead(imagename=fluxcal_image)['refval'][-1]/10**9 # In GHz
		ref_beam_axes_multi=500000
		beam_axes_multi=major*minor	
		scale_tb_limit=(beam_axes_multi/ref_beam_axes_multi)
		expr='1.222e6*IM0/'+str(freq)+'^2/('+str(major*minor)+')'
		if os.path.exists(outfile_path+'/I_Tb.image'):
			os.system('rm -rf '+outfile_path+'/I_Tb.image')
		immath(imagename=fluxcal_image,outfile=outfile_path+'/I_Tb.image',mode='evalexpr',expr=expr)
		imhead(imagename=outfile_path+'/I_Tb.image', mode='put', hdkey='bunit', hdvalue='K')
		ia=image()
		ia.open(imagename)
		data=ia.getchunk()
		ia.close()
		a=image()
		ia.open(modelname)
		modeldata=ia.getchunk()
		ia.close()
		model_datai=modeldata[:,:,0,0]
		model_dataq=modeldata[:,:,1,0]
		model_datau=modeldata[:,:,2,0]
		rmsq=imstat(imagename=imagename,box=self.rms_box_nearsun,stokes='Q')['rms'][0]
		rmsu=imstat(imagename=imagename,box=self.rms_box_nearsun,stokes='U')['rms'][0]
		rmsi=imstat(imagename=imagename,box=self.rms_box_nearsun,stokes='I')['rms'][0]
		i_flux=imstat(imagename=imagename,stokes='I')['flux'][0]
		dataq=data[:,:,1,0]
		datau=data[:,:,2,0]
		datai=data[:,:,0,0]
		datai_copy=copy.deepcopy(datai)
		dataq_copy=copy.deepcopy(dataq)
		datau_copy=copy.deepcopy(datau)
		dataq1=copy.deepcopy(dataq)
		datau1=copy.deepcopy(datau)
		posi=np.where(datai<(sigma*rmsi))
		posi20=np.where(datai<(sigma*rmsi))
		maxpos=imstat(imagename=imagename,stokes='I')['maxpos']
		box=self.negative_box(maxpos,box_width=5).split(',')
		box_coords=[int(i) for i in box]
		posq=np.where(np.abs(dataq)<(3*rmsq))
		posu=np.where(np.abs(datau)<(3*rmsu))
		ia.open(outfile_path+'/I_Tb.image')
		datatb=ia.getchunk()
		ia.close()
		datai_mask=copy.deepcopy(datai)
		for i in range(datai_mask.shape[0]):
			for j in range(datai_mask.shape[1]):
				if i<box_coords[0] or i>box_coords[2] or j<box_coords[1] or j>box_coords[3]:
					datai_mask[i,j]=np.nan
		tb_nonpos=np.where(datatb[:,:,0,0]<=0)
		datatb[:,:,0,0][tb_nonpos]=np.nan
		datatblog=np.log10(datatb[:,:,0,0])
		datatblog=np.round(datatblog,0)
		pos=np.where(np.isnan(datatblog)==True)
		datatblog[posi20]=np.nan
		mediantb=np.nanmedian(datatblog)
		if mediantb>6 or mediantb<5:
			mediantb=5
		datatblog[pos]=0
		postb=np.where((datatblog>mediantb+2)|(datatblog<mediantb))
		postb2=np.where((datatblog<=mediantb+2) & (datatblog>=mediantb))
		datai[posi]=np.nan
		dataq[posi]=np.nan
		datai[np.isnan(datai_mask)==True]=np.nan
		dataq[np.isnan(datai_mask)==True]=np.nan
		area=np.nansum(np.isnan(datai)==False)
		datai[postb]=np.nan
		dataq[postb]=np.nan
		dataq[posq]=np.nan	
		unmasked_pixels=np.nansum(np.isnan(dataq)==False)
		x=[]
		y=[]
		z=[]
		q_by_i=dataq/datai
		for k in range(datai.shape[0]):
			for l in range(datai.shape[1]):
				if np.isnan(q_by_i[k,l])==False:
					x.append(k)
					y.append(l)
					z.append(q_by_i[k,l])
		q_stack=np.vstack((x,y,z)).T
		del x,y,z
		pre_q_image=np.nanmean(np.abs(q_by_i))
		dataq_backup=copy.deepcopy(dataq_copy)
		if (unmasked_pixels/area)>0.3:
			AQ = np.c_[np.ones(q_stack.shape[0]), q_stack[:,:2], np.prod(q_stack[:,:2], axis=1), q_stack[:,:2]**2]
			CQ,_,_,_ = scipy.linalg.lstsq(AQ, q_stack[:,2])	
			for k in range(datai_copy.shape[0]):
				for l in range(datai_copy.shape[1]):
					if np.isnan(datai_mask[k,l])==False:
						dataq_copy[k,l] -= (CQ[4]*k**2. + CQ[5]*l**2. + CQ[3]*k*l + CQ[1]*k + CQ[2]*l + CQ[0])*datai_copy[k,l]
		else:
			CQ=[0,0,0,0,0,0]
		q_by_i_new=copy.deepcopy(dataq_copy/datai_copy)
		q_by_i_new[np.isnan(datai)==True]=np.nan
		new_q_image=np.nanmean(np.abs(q_by_i_new))
		datai=copy.deepcopy(datai_copy)
		datau=copy.deepcopy(datau_copy)
		datai[posi]=np.nan
		datau[posi]=np.nan
		datai[np.isnan(datai_mask)==True]=np.nan
		datau[np.isnan(datai_mask)==True]=np.nan
		area=np.nansum(np.isnan(datai)==False)
		datau[posu]=np.nan
		datai[postb]=np.nan
		datau[postb]=np.nan
		unmasked_pixels=np.nansum(np.isnan(datau)==False)
		u_by_i=datau/datai
		x=[]
		y=[]
		z=[]
		for k in range(datai.shape[0]):
			for l in range(datai.shape[1]):
				if np.isnan(u_by_i[k,l])==False:
					x.append(k)
					y.append(l)
					z.append(u_by_i[k,l])
		u_stack=np.vstack((x,y,z)).T
		del x,y,z
		pre_u_image=np.nanmean(np.abs(u_by_i))
		if (unmasked_pixels/area)>0.3:
			AU = np.c_[np.ones(u_stack.shape[0]), u_stack[:,:2], np.prod(u_stack[:,:2], axis=1), u_stack[:,:2]**2]
			CU,_,_,_ = scipy.linalg.lstsq(AU, u_stack[:,2])
			for k in range(datai_copy.shape[0]):
				for l in range(datai_copy.shape[1]):
					if np.isnan(datai_mask[k,l])==False:
						datau_copy[k,l] -= (CU[4]*k**2. + CU[5]*l**2. + CU[3]*k*l + CU[1]*k + CU[2]*l + CU[0])*datai_copy[k,l]
		else:
			CU=[0,0,0,0,0,0]
		u_by_i_new=copy.deepcopy(datau_copy/datai_copy)
		u_by_i_new[np.isnan(datai)==True]=np.nan
		new_u_image=np.nanmean(np.abs(u_by_i_new))
		'''
		X,Y = np.meshgrid(np.arange(430, 550, 1), np.arange(430, 550, 1))
		XX = X.flatten()
		YY = Y.flatten()
		fig = plt.figure()
		ax = plt.axes(projection='3d')
		Z=(CU[4]*X**2. + CU[5]*Y**2. + CU[3]*X*Y + CU[1]*X + CU[2]*Y + CU[0])
		ax.plot_surface(X, Y, Z*100, rstride=1, cstride=1, alpha=0.3)
		ax.scatter(u_stack[:,0], u_stack[:,1], u_stack[:,2]*100, c='r', s=10)
		plt.xlabel('RA',fontsize=10)
		plt.ylabel('DEC',fontsize=10)
		ax.set_zlabel('Stokes U (%)',fontsize=10)
		ax.axis('tight')
		plt.show()
		X,Y = np.meshgrid(np.arange(430, 550, 1), np.arange(430, 550, 1))
		XX = X.flatten()
		YY = Y.flatten()
		Z=(CQ[4]*X**2. + CQ[5]*Y**2. + CQ[3]*X*Y + CQ[1]*X + CQ[2]*Y + CQ[0])
		fig = plt.figure()
		ax = plt.axes(projection='3d')
		ax.plot_surface(X, Y, Z*100, rstride=1, cstride=1, alpha=0.3)
		ax.scatter(q_stack[:,0], q_stack[:,1], q_stack[:,2]*100, c='r', s=10)
		plt.xlabel('RA',fontsize=10)
		plt.ylabel('DEC',fontsize=10)
		ax.set_zlabel('Stokes Q (%)',fontsize=10)
		ax.axis('tight')
		plt.show()
		plt.imshow(dataq_copy)
		plt.show()'''
		posq=np.where(np.abs(dataq_copy)<(3*rmsq))
		posu=np.where(np.abs(datau_copy)<(3*rmsu))
		data[:,:,0,0]=datai_copy
		data[:,:,1,0]=dataq_copy
		data[:,:,2,0]=datau_copy
		ia.open(imagename)
		ia.putchunk(data)
		ia.close()
		for k in range(model_datai.shape[0]):
			for l in range(model_datai.shape[1]):
				if np.isnan(datai_mask[k,l])==False:
					model_dataq[k,l] -= (CU[4]*k**2. + CU[5]*l**2. + CU[3]*k*l + CU[1]*k + CU[2]*l + CU[0])*model_datai[k,l]
		for k in range(model_datai.shape[0]):
			for l in range(model_datai.shape[1]):
				if np.isnan(datai_mask[k,l])==False:
					model_datau[k,l] -= (CU[4]*k**2. + CU[5]*l**2. + CU[3]*k*l + CU[1]*k + CU[2]*l + CU[0])*model_datai[k,l]
		modelq=copy.deepcopy(model_dataq)
		modelu=copy.deepcopy(model_datau)
		modeldata[:,:,0,0]=model_datai
		modeldata[:,:,1,0]=modelq
		modeldata[:,:,2,0]=modelu
		dataq_by_i=dataq_copy/datai_copy
		datau_by_i=datau_copy/datai_copy
		dataq_by_i[posi]=np.nan
		dataq_by_i[posq]=np.nan
		dataq_by_i[postb]=np.nan
		datau_by_i[posi]=np.nan
		datau_by_i[posu]=np.nan
		datau_by_i[postb]=np.nan
		model_dataq[posq]=0
		model_datau[posu]=0
		model_dataq[postb2]=0
		model_datau[postb2]=0
		q_change=abs(pre_q_image-new_q_image)
		if q_change>0:
			modeldata[:,:,1,0]=model_dataq
		u_change=abs(pre_u_image-new_u_image)
		if u_change>0:
			modeldata[:,:,2,0]=model_datau
		ia.open(modelname)
		ia.putchunk(modeldata)
		ia.close()
		os.system('rm -rf casa*log I*.image')
		os.system('rm -rf '+outfile_path+'/I.* '+outfile_path+'/I_*')
		del modeldata,data,datai_mask,dataq,datai,datau,model_datai,model_dataq,model_datau,datai_copy,dataq_copy,datau_copy
		return imagename,modelname,round(q_change,2),round(u_change,2)

	def subtract_stokesV_solar_leakage_surface(self,imagename,modelname,sigma=10,do_fluxcal=False):
		'''
		Function to subtract quadratic leakage surface

		Parameters
		----------
		imagename : str 
			Name of the image
		modelname : str 
			Name of the model
		sigma : float 
			N-sigma value above which any emission is considered to be real
		do_fluxcal : bool
			Perform polynomial based flux calibration (See details Kansabanik et al. 2021, submitted to ApJ)
		Returns
		-------
		str
			Leakage surface subtracted image 
		str
			Leakage surface subtracted model
		float
			Stokes V change fraction
		'''
		if sigma<5:
			sigma=5
		self.log_verbose.info('Correcting Stokes I to Stokes V leakage surface.\n')
		outfile_path=os.path.dirname(os.path.realpath(imagename))
		os.system('rm -rf '+outfile_path+'/I.* '+outfile_path+'/I_*')
		imsubimage(imagename=imagename,outfile=outfile_path+'/I.image',stokes='I',dropdeg=False)
		if do_fluxcal==True:
			fluxcal_image=self.mwa_solar_fluxcal(imagename=outfile_path+'/I.image',outfile=outfile_path+'/I_fluxcal.image')
		else:
			fluxcal_image=outfile_path+'/I.image'
		major=imhead(imagename=fluxcal_image)['restoringbeam']['major']['value'] # In arcsec
		minor=imhead(imagename=fluxcal_image)['restoringbeam']['minor']['value'] # In arcsec
		freq=imhead(imagename=fluxcal_image)['refval'][-1]/10**9 # In GHz
		ref_beam_axes_multi=500000
		beam_axes_multi=major*minor	
		scale_tb_limit=(beam_axes_multi/ref_beam_axes_multi)
		expr='1.222e6*IM0/'+str(freq)+'^2/('+str(major*minor)+')'
		if os.path.exists(outfile_path+'/I_Tb.image'):
			os.system('rm -rf '+outfile_path+'/I_Tb.image')
		immath(imagename=fluxcal_image,outfile=outfile_path+'/I_Tb.image',mode='evalexpr',expr=expr)
		imhead(imagename=outfile_path+'/I_Tb.image', mode='put', hdkey='bunit', hdvalue='K')
		ia=image()
		ia.open(imagename)
		data=ia.getchunk()
		ia.close()
		a=image()
		ia.open(modelname)
		modeldata=ia.getchunk()
		ia.close()
		model_datai=modeldata[:,:,0,0]
		model_datav=modeldata[:,:,3,0]
		rmsv=imstat(imagename=imagename,box=self.rms_box_nearsun,stokes='V')['rms'][0]
		rmsi=imstat(imagename=imagename,box=self.rms_box_nearsun,stokes='I')['rms'][0]
		i_flux=imstat(imagename=imagename,stokes='I')['flux'][0]
		datai=data[:,:,0,0]
		datav=data[:,:,3,0]
		datai_copy=copy.deepcopy(datai)
		datav_copy=copy.deepcopy(datav)
		datav1=copy.deepcopy(datav)
		posi=np.where(datai<(sigma*rmsi))
		posi20=np.where(datai<(sigma*rmsi))
		maxpos=imstat(imagename=imagename,stokes='I')['maxpos']
		box=self.negative_box(maxpos,box_width=5).split(',')
		box_coords=[int(i) for i in box]
		posv=np.where(np.abs(datav)<(5*rmsv))
		ia.open(outfile_path+'/I_Tb.image')
		datatb=ia.getchunk()
		ia.close()
		datai_mask=copy.deepcopy(datai)
		for i in range(datai_mask.shape[0]):
			for j in range(datai_mask.shape[1]):
				if i<box_coords[0] or i>box_coords[2] or j<box_coords[1] or j>box_coords[3]:
					datai_mask[i,j]=np.nan
		tb_nonpos=np.where(datatb[:,:,0,0]<=0)
		datatb[:,:,0,0][tb_nonpos]=np.nan
		datatblog=np.log10(datatb[:,:,0,0])
		datatblog=np.round(datatblog,0)
		pos=np.where(np.isnan(datatblog)==True)
		datatblog[posi20]=np.nan
		mediantb=np.nanmedian(datatblog)
		if mediantb>6 or mediantb<5:
			mediantb=5
		datatblog[pos]=0
		postb=np.where((datatblog>mediantb+2)|(datatblog<mediantb))
		postb2=np.where((datatblog<=mediantb+2) & (datatblog>=mediantb))
		datai[posi]=np.nan
		datav[posi]=np.nan	
		datai[np.isnan(datai_mask)==True]=np.nan
		datav[np.isnan(datai_mask)==True]=np.nan
		area=np.nansum(np.isnan(datai)==False)
		datav[posv]=np.nan
		datai[postb]=np.nan
		datav[postb]=np.nan	
		unmasked_pixels=np.nansum(np.isnan(datav)==False)
		v_by_i=datav/datai
		x=[]
		y=[]
		z=[]
		for k in range(datai.shape[0]):
			for l in range(datai.shape[1]):
				if np.isnan(v_by_i[k,l])==False:
					x.append(k)
					y.append(l)
					z.append(v_by_i[k,l])
		v_stack=np.vstack((x,y,z)).T
		del x,y,z
		pre_v_image=np.nanmedian(np.abs(v_by_i))
		if (unmasked_pixels/area)>0.3:
			AV = np.c_[np.ones(v_stack.shape[0]), v_stack[:,:2], np.prod(v_stack[:,:2], axis=1), v_stack[:,:2]**2]
			CV,_,_,_ = scipy.linalg.lstsq(AV, v_stack[:,2])	
			for k in range(datai_copy.shape[0]):
				for l in range(datai_copy.shape[1]):
					if np.isnan(datai_mask[k,l])==False:
						datav_copy[k,l] -= (CV[4]*k**2. + CV[5]*l**2. + CV[3]*k*l + CV[1]*k + CV[2]*l + CV[0])*datai_copy[k,l]
		else:
			CV=[0,0,0,0,0,0]
		v_by_i_new=copy.deepcopy(datav_copy/datai_copy)
		v_by_i_new[np.isnan(datai)==True]=np.nan
		new_v_image=np.nanmedian(np.abs(v_by_i_new))
		'''X,Y = np.meshgrid(np.arange(430, 550, 1), np.arange(430, 550, 1))
		XX = X.flatten()
		YY = Y.flatten()
		fig = plt.figure()
		ax = plt.axes(projection='3d')
		Z=(CV[4]*X**2. + CV[5]*Y**2. + CV[3]*X*Y + CV[1]*X + CV[2]*Y + CV[0])
		ax.plot_surface(X, Y, Z*100, rstride=1, cstride=1, alpha=0.3)
		ax.scatter(v_stack[:,0], v_stack[:,1], v_stack[:,2]*100, c='r', s=10)
		plt.xlabel('RA',fontsize=10)
		plt.ylabel('DEC',fontsize=10)
		ax.set_zlabel('Stokes V (%)',fontsize=10)
		ax.axis('tight')
		plt.show()
		'''
		v_change=abs(pre_v_image-new_v_image)
		posv=np.where(np.abs(datav_copy)<(3*rmsv))
		data[:,:,0,0]=datai_copy
		data[:,:,3,0]=datav_copy
		ia.open(imagename)
		ia.putchunk(data)
		ia.close()
		for k in range(model_datai.shape[0]):
			for l in range(model_datai.shape[1]):
				if np.isnan(datai_mask[k,l])==False:
					model_datav[k,l] -= (CV[4]*k**2. + CV[5]*l**2. + CV[3]*k*l + CV[1]*k + CV[2]*l + CV[0])*model_datai[k,l]
		modelv=copy.deepcopy(model_datav)
		modeldata[:,:,0,0]=model_datai
		modeldata[:,:,3,0]=modelv
		datav_by_i=datav_copy/datai_copy
		datav_by_i[posi]=np.nan
		datav_by_i[posv]=np.nan
		model_datav[posv]=0
		modeldata[:,:,3,0]=model_datav
		ia.open(modelname)
		ia.putchunk(modeldata)
		ia.close()
		os.system('rm -rf casa*log I*.image')
		os.system('rm -rf '+outfile_path+'/I.* '+outfile_path+'/I_*')
		del modeldata,data,datai_mask,datav,datai,model_datai,model_datav,datai_copy,datav_copy
		return imagename,modelname,round(v_change,2)	

	def mwa_solar_fluxcal(self,imagename,outfile,atten=10): #TODO : use reference band shape
		'''
		Function to flux calibrate MWA solar observations using method described in Kansabanik et al. 2021

		Parameters
		----------
		imagename : str 
			Name of the image
		outfile : str
			Output file name
		atten : float 
			Attenuator gain value in dB
		Returns
		-------
		str
			Flux calibrated image
		''' 
		self.log_verbose.info('Perform flux density calibration....\n')
		if imagename[-1]=='/':
			imagename=imagename[:-1]
		if os.path.exists(imagename+'.tempflux'):
			os.system('rm -rf '+imagename+'.tempflux')
		freq_axis=np.where(imhead(imagename=imagename)['axisnames']=='Frequency')[0][0]
		freq=imhead(imagename)['refval'][freq_axis]/10**6 # MHz
		fluxscale_poly=np.poly1d(np.load(datadir+'/flux_scale_polyfit.npy',allow_pickle=True)[0])
		fluxscale=fluxscale_poly(freq)
		# Taking care of beam normalization applied in Kansabanik et al. 2022, robust fux calibration.
		###########################
		metaheader=fits.getheader(self.metafits)
		pointing_ra=metaheader['RA']
		pointing_dec=metaheader['DEC']
		beam=pb.get_pb_radec(pointing_ra,pointing_dec,freq,metafits=self.metafits)
		I_beam=beam[2]
		fluxscale=fluxscale*I_beam
		#############################
		if atten!=10: # Valid for 10 and 14 dB, not tested for any other values and typically other values not used for solar observations
			fluxscale=fluxscale*1.19
		self.log_verbose.info('immath(imagename=\''+imagename+'\',outfile=\''+imagename+'.tempflux\',mode=\'evalexpr\',expr=\'IM0*'+str(fluxscale)+'\')\n')
		immath(imagename=imagename,outfile=imagename+'.tempflux',mode='evalexpr',expr='IM0*'+str(fluxscale))
		if os.path.exists(outfile):
			os.system('rm -rf '+outfile)
		os.system('mv '+imagename+'.tempflux '+outfile)
		return outfile

	def cal_solar_qu_leakage(self,imagename,sigma=10,do_fluxcal=False): 
		'''
		Function to calculate Stokes QU leakage for solar observation (Not vaild for any other astrophysical observation)

		Parameters
		----------
		imagename : str 
			Name of the image
		do_fluxcal : bool 
			Perform flux calibration or not
		outfile_path : str 
			Name of the directory to save leakage data numpy table (default : image directory)
		Returns
		-------
		float
			Stokes Q leakage
		float 
			Stokes U leakage
		'''
		outfile_path=os.path.dirname(os.path.realpath(imagename))
		os.system('rm -rf '+outfile_path+'/I_* '+outfile_path+'/I.*')
		imsubimage(imagename=imagename,outfile=outfile_path+'/I.image',stokes='I',dropdeg=False)
		if do_fluxcal==True:
			atten=float(fits.getheader(self.metafits)['ATTEN_DB'])
			fluxcal_image=self.mwa_solar_fluxcal(imagename=outfile_path+'/I.image',outfile=outfile_path+'/I_fluxcal.image',atten=atten)
		else:
			fluxcal_image=outfile_path+'/I.image'
		major=imhead(imagename=fluxcal_image)['restoringbeam']['major']['value'] # In arcsec
		minor=imhead(imagename=fluxcal_image)['restoringbeam']['minor']['value'] # In arcsec
		freq=imhead(imagename=fluxcal_image)['refval'][-1]/10**9 # In GHz
		ref_beam_axes_multi=500000
		beam_axes_multi=major*minor	
		scale_tb_limit=(beam_axes_multi/ref_beam_axes_multi)
		expr='1.222e6*IM0/'+str(freq)+'^2/('+str(major*minor)+')'
		if os.path.exists(outfile_path+'/I_Tb.image'):
			os.system('rm -rf '+outfile_path+'/I_Tb.image')
		immath(imagename=fluxcal_image,outfile=outfile_path+'/I_Tb.image',mode='evalexpr',expr=expr)
		imhead(imagename=outfile_path+'/I_Tb.image', mode='put', hdkey='bunit', hdvalue='K')
		ia=image()
		ia.open(imagename)
		data=ia.getchunk()
		ia.close()
		imsize=imhead(imagename,mode='list')['shape'][0]
		rmsq=imstat(imagename=imagename,box=self.rms_box,stokes='Q')['rms'][0]
		rmsu=imstat(imagename=imagename,box=self.rms_box,stokes='U')['rms'][0]
		rmsi=imstat(imagename=imagename,box=self.rms_box,stokes='I')['rms'][0]
		i_flux=imstat(imagename=imagename,stokes='I')['flux'][0]
		dataq=copy.deepcopy(data[:,:,1,0])
		datau=copy.deepcopy(data[:,:,2,0])
		datai=copy.deepcopy(data[:,:,0,0])
		posi=np.where(datai<(sigma*rmsi))
		posq=np.where(np.abs(dataq)<(sigma*rmsq))
		posu=np.where(np.abs(datau)<(sigma*rmsu))
		posi20=np.where(datai<(20*rmsi))
		ia.open(outfile_path+'/I_Tb.image')
		datatb=ia.getchunk()
		ia.close()
		tb_nonpos=np.where(datatb[:,:,0,0]<=0)
		datatb[:,:,0,0][tb_nonpos]=np.nan
		datatblog=datatb[:,:,0,0]/10**6
		pos=np.where(np.isnan(datatblog)==True)
		datatblog[posi20]=np.nan
		maxpos=imstat(imagename=imagename,stokes='I')['maxpos']
		box=self.negative_box(maxpos,box_width=3).split(',')
		box_coords=[int(i) for i in box]
		for i in range(datatblog.shape[0]):
			for j in range(datatblog.shape[1]):
				if i<box_coords[0] or i>box_coords[2] or j<box_coords[1] or j>box_coords[3]:
					datatblog[i,j]=np.nan
				else:
					pass
		mediantb=np.nanmedian(datatblog)
		if mediantb>5:
			print ('No quiet Sun region found.\n')
			return 0,0
		elif mediantb<0.5:
			mediantb=0.5
		if scale_tb_limit<1:
			postb=np.where(np.round(datatblog,0)>1)
		else:
			postb=np.where(np.round(datatblog,0)>scale_tb_limit)
		cellsize=np.rad2deg(np.abs(imhead(imagename=imagename,mode='list')['cdelt1'])) # In degree
		circular_mask=self.create_circular_mask(imsize,imsize, center=(int(imsize/2),int(imsize/2)), radius=int(5/cellsize))
		datai[~circular_mask]=np.nan
		dataq[~circular_mask]=np.nan
		# Q leakage calculation
		#######################
		datai[posi]=np.nan
		datai[postb]=np.nan
		unmasked_pixels=np.nansum(np.isnan(datai)==False)
		tot_sun_size=unmasked_pixels*self.cellsize**2
		dataq[posi]=np.nan
		unmasked_pixels=np.nansum(np.isnan(dataq/datai)==False)
		unmasked_source_size=unmasked_pixels*self.cellsize**2
		psf_size=np.pi*self.psf**2
		if (int(unmasked_source_size/psf_size)<1 or unmasked_source_size/tot_sun_size<0.1):
			os.system('rm -rf casa*log I*.image')
			os.system('rm -rf '+outfile_path+'/I_*')
			q_leakage=0
		else:
			q_leakage=np.nanmedian(dataq/datai)
		# U leakage calculation
		#######################
		datau=copy.deepcopy(data[:,:,2,0])
		datai=copy.deepcopy(data[:,:,0,0])
		datai[posi]=np.nan
		datai[postb]=np.nan
		datai[~circular_mask]=np.nan
		datau[~circular_mask]=np.nan
		unmasked_pixels=np.nansum(np.isnan(datai)==False)
		tot_sun_size=unmasked_pixels*self.cellsize**2
		datau[posi]=np.nan
		unmasked_pixels=np.nansum(np.isnan(datau/datai)==False)
		unmasked_source_size=unmasked_pixels*self.cellsize**2
		psf_size=np.pi*self.psf**2
		if (int(unmasked_source_size/psf_size)<1 or unmasked_source_size/tot_sun_size<0.1):
			os.system('rm -rf casa*log I.*image')
			os.system('rm -rf '+outfile_path+'/I.*')
			u_leakage=0
		else:
			u_leakage=np.nanmedian(datau/datai)
		os.system('rm -rf casa*log I.*image')
		os.system('rm -rf '+outfile_path+'/I_*')
		self.log_verbose.info('Calculation of Stokes leakages has been done.\n')
		return q_leakage,u_leakage


	def create_circular_mask(self,h, w, center=None, radius=None):
		'''
		Function to create a circular mask

		Parameters
		----------
		h : int 
			Number of pixels along Y axis
		w : int 
			Number of pixels along X axis
		center : tuple 
			(x_cen,y_cen), center of the mask circle
		radius : float 
			Radius in number of pixels
		Returns
		-------
		numpy.array
			Cirular mask
		'''
		if center is None:
		    center=(int(w/2),int(h/2))
		if radius is None:
		    radius=min(center[0],center[1],w-center[0],h-center[1])
		Y,X=np.ogrid[:h,:w]
		dist_from_center=np.sqrt((X-center[0])**2+(Y-center[1])**2)
		mask=dist_from_center<=radius
		return mask

	def cal_solar_v_leakage(self,imagename,sigma=10,do_fluxcal=False): 
		'''
		Function to calculate Stokes V leakage for solar observation (Not vaild for any other astrophysical observation)

		Parameters
		----------
		imagename : str 
			Name of the image
		sigma : float
			Sigma value for thresholding
		do_fluxcal : bool 
			Perform flux calibration or not
		Returns
		-------
		float
			Stokes V leakage
		'''
		outfile_path=os.path.dirname(os.path.realpath(imagename))
		os.system('rm -rf '+outfile_path+'/I_* '+outfile_path+'/I.*')
		self.log_verbose.info('Correcting Stokes V leakage.\n')
		imsubimage(imagename=imagename,outfile=outfile_path+'/I.image',stokes='I',dropdeg=False)
		if do_fluxcal==True:
			atten=float(fits.getheader(self.metafits)['ATTEN_DB'])
			fluxcal_image=self.mwa_solar_fluxcal(imagename=outfile_path+'/I.image',outfile=outfile_path+'/I_fluxcal.image',atten=atten)
		else:
			fluxcal_image=outfile_path+'/I.image'
		major=imhead(imagename=fluxcal_image)['restoringbeam']['major']['value'] # In arcsec
		minor=imhead(imagename=fluxcal_image)['restoringbeam']['minor']['value'] # In arcsec
		freq=imhead(imagename=fluxcal_image)['refval'][-1]/10**9 # In GHz
		ref_beam_axes_multi=500000
		beam_axes_multi=major*minor	
		scale_tb_limit=(beam_axes_multi/ref_beam_axes_multi)
		expr='1.222e6*IM0/'+str(freq)+'^2/('+str(major*minor)+')'
		if os.path.exists(outfile_path+'/I_Tb.image'):
			os.system('rm -rf '+outfile_path+'/I_Tb.image')
		immath(imagename=fluxcal_image,outfile=outfile_path+'/I_Tb.image',mode='evalexpr',expr=expr)
		imhead(imagename=outfile_path+'/I_Tb.image', mode='put', hdkey='bunit', hdvalue='K')
		ia=image()
		ia.open(imagename)
		data=ia.getchunk()
		ia.close()
		imsize=imhead(imagename,mode='list')['shape'][0]
		cellsize=np.rad2deg(np.abs(imhead(imagename=imagename,mode='list')['cdelt1'])) # In degree
		circular_mask=self.create_circular_mask(imsize,imsize, center=(int(imsize/2),int(imsize/2)), radius=int(5/cellsize))
		rmsv=imstat(imagename=imagename,box=self.rms_box_nearsun,stokes='V')['rms'][0]
		rmsi=imstat(imagename=imagename,box=self.rms_box_nearsun,stokes='I')['rms'][0]
		datav=data[:,:,3,0]
		datai=data[:,:,0,0]
		datav_copy=copy.deepcopy(datav)
		posi=np.where(datai<(sigma*rmsi))
		posv=np.where(np.abs(datav)<(sigma*rmsv))
		posi20=np.where(datai<(20*rmsi))
		datai[~circular_mask]=np.nan
		datav[~circular_mask]=np.nan
		datai[posi]=np.nan
		unmasked_pixels=np.nansum(np.isnan(datai)==False)
		tot_sun_size=unmasked_pixels*self.cellsize**2
		datav[posi]=np.nan
		ia.open(outfile_path+'/I_Tb.image')
		datatb=ia.getchunk()
		ia.close()
		tb_nonpos=np.where(datatb[:,:,0,0]<=0)
		datatb[:,:,0,0][tb_nonpos]=np.nan
		datatblog=datatb[:,:,0,0]/10**6
		pos=np.where(np.isnan(datatblog)==True)
		datatblog[posi20]=np.nan
		maxpos=imstat(imagename=imagename,stokes='I')['maxpos']
		box=self.negative_box(maxpos,box_width=3).split(',')
		box_coords=[int(i) for i in box]
		for i in range(datatblog.shape[0]):
			for j in range(datatblog.shape[1]):
				if i<box_coords[0] or i>box_coords[2] or j<box_coords[1] or j>box_coords[3]:
					datatblog[i,j]=np.nan
				else:
					pass
		mediantb=np.nanmedian(datatblog)
		if mediantb>5:
			print ('No quiet Sun region found.\n')
			return 0
		elif mediantb<0.5:
			mediantb=0.5
		if scale_tb_limit<1:
			postb=np.where(np.round(datatblog,0)<=1)
		else:
			postb=np.where(np.round(datatblog,0)<=scale_tb_limit)
		datai[posi20]=np.nan
		datav[postb]=np.nan
		v_by_i=datav/datai
		unmasked_pixels=np.nansum(np.isnan(datav/datai)==False)
		unmasked_source_size=unmasked_pixels*self.cellsize**2
		psf_size=np.pi*self.psf**2
		if (int(unmasked_source_size/psf_size)<1 or unmasked_source_size/tot_sun_size<0.1):
			os.system('rm -rf casa*log I.*image')
			os.system('rm -rf '+outfile_path+'/I.*')
			return 0
		os.system('rm -rf casa*log I.*image')
		os.system('rm -rf '+outfile_path+'/I_*')
		return np.nanmedian(datav/datai)

	def correct_solar_quv_surface_leakage(self,imagename,modelname,sigma,overwrite=False,stokes='QUV',v_frac_threshold=0.01,do_fluxcal=False,rms_thresh=True):
		'''
		Function to correct for solar Stokes Q, U and V leakage based on second order leakage surface fit 
		(It is based on the fact that we do not expect any linear polarisation from the Quiet Sun emission)

		Parameters
		----------
		imagename : str 
			Name of image
		modelname : str 
			Name of the model
		sigma : float 
			N-sigma threshold to choose Stokes I emission region
		overwrite : bool 
			Overwrite the image and model or not
		stokes : str
			Stokes plane to correct
		do_fluxcal : bool
			Perform flux calibration
		rms_thresh : bool
			Put a cutoff on the final leakage corrected model based on rms 
		Returns
		-------
		str
			Stokes Q ,U, V leakage corrected image name
		str 
			Stokes Q ,U, V leakage corrected model name 
		float
			Stokes Q fractional change
		float 
			Stokes U fractional change
		float
			Stokes V fractional change	 
		'''
		q_change_frac,u_change_frac,v_change_frac=0,0,0
		self.log_verbose.info('Correcting Stokes I to Stokes '+','.join(list(stokes))+'.....\n')
		self.log_verbose.info('Image name : '+imagename+'\n')
		self.log_verbose.info('Model name : '+modelname+'\n')
		if os.path.exists(imagename.split('.image')[0]+'_quvcor_surface.image'):
			os.system('rm -rf '+imagename.split('.image')[0]+'_quvcor_surface.image')
		if os.path.exists(modelname.split('.model')[0]+'_quvcor_surface.model'):
			os.system('rm -rf '+modelname.split('.model')[0]+'_quvcor_surface.model')
		if overwrite==False:
			os.system('cp -r '+imagename+' '+imagename.split('.image')[0]+'_quvcor_surface.image')
			os.system('cp -r '+modelname+' '+modelname.split('.model')[0]+'_quvcor_surface.model')
			imagename=imagename.split('.image')[0]+'_quvcor_surface.image'
			modelname=modelname.split('.model')[0]+'_quvcor_surface.model'
		outfile_path=os.path.dirname(os.path.realpath(imagename))
		imagename,modelname,q_change_frac,u_change_frac=self.subtract_leakage_surface(imagename,modelname,sigma=sigma,do_fluxcal=do_fluxcal)
		os.system('cp -r '+imagename+' '+imagename+'.tempQU')
		os.system('cp -r '+modelname+' '+modelname+'.tempQU')
		if 'V'in stokes:
			imagename,modelname,v_change_frac=self.subtract_stokesV_solar_leakage_surface(imagename,modelname,sigma=sigma,do_fluxcal=do_fluxcal)
		if abs(round(v_change_frac*100,1))<(v_frac_threshold*100):
			self.log_verbose.info('Determined V leakage : '+str(v_change_frac*100)+' % is less than '+str(v_frac_threshold*100)+' %. Do not correcting Stokes V image for leakage.\n')
			v_leakage=0.0
			os.system('mv '+imagename+'.tempQU '+imagename)
			os.system('mv '+modelname+'.tempQU '+modelname)
		else:
			os.system('rm -rf '+imagename+'.tempQU')
			os.system('rm -rf '+modelname+'.tempQU')
		ia=image()
		ia.open(imagename) # Correcting image
		data=ia.getchunk()
		datai=data[:,:,0,:]
		dataq=data[:,:,1,:]
		datau=data[:,:,2,:]
		datav=data[:,:,3,:]
		rmsq=imstat(imagename=imagename,box=self.rms_box_nearsun,stokes='Q')['rms'][0]
		rmsu=imstat(imagename=imagename,box=self.rms_box_nearsun,stokes='U')['rms'][0]
		rmsv=imstat(imagename=imagename,box=self.rms_box_nearsun,stokes='V')['rms'][0]
		rmsi=imstat(imagename=imagename,box=self.rms_box_nearsun,stokes='I')['rms'][0]
		ia.putchunk(data)
		ia.close()
		imsubimage(imagename=imagename,outfile=outfile_path+'/I.image',stokes='I',dropdeg=False)
		if do_fluxcal:
			atten=float(fits.getheader(self.metafits)['ATTEN_DB'])
			fluxcal_image=self.mwa_solar_fluxcal(imagename=outfile_path+'/I.image',outfile=outfile_path+'/I_fluxcal.image',atten=atten)
		else:
			fluxcal_image=outfile_path+'/I.image'
		major=imhead(imagename=fluxcal_image)['restoringbeam']['major']['value'] # In arcsec
		minor=imhead(imagename=fluxcal_image)['restoringbeam']['minor']['value'] # In arcsec
		freq=imhead(imagename=fluxcal_image)['refval'][-1]/10**9 # In GHz
		ref_beam_axes_multi=500000
		beam_axes_multi=major*minor	
		expr='1.222e6*IM0/'+str(freq)+'^2/('+str(major*minor)+')'
		if os.path.exists(outfile_path+'/I_Tb.image'):
			os.system('rm -rf '+outfile_path+'/I_Tb.image')
		immath(imagename=fluxcal_image,outfile=outfile_path+'/I_Tb.image',mode='evalexpr',expr=expr)
		imhead(imagename=outfile_path+'/I_Tb.image', mode='put', hdkey='bunit', hdvalue='K')
		ia.open(outfile_path+'/I_Tb.image')
		datatb=ia.getchunk()
		ia.close()
		posi20=np.where(datai[:,:,0]<(20*rmsi))
		tb_nonpos=np.where(datatb[:,:,0,0]<=0)
		datatb[:,:,0,0][tb_nonpos]=np.nan
		datatblog=np.log10(datatb[:,:,0,0])
		datatblog=np.round(datatblog,0)
		pos=np.where(np.isnan(datatblog)==True)
		datatblog[posi20]=np.nan
		mediantb=np.nanmedian(datatblog)
		if mediantb>6 or mediantb<5:
			mediantb=5
		datatblog[pos]=0
		postb=np.where((datatblog>mediantb+2)|(datatblog<mediantb))
		postb2=np.where((datatblog<=mediantb+2) & (datatblog>=mediantb))
		ia.open(modelname) # Correcting model
		datam=ia.getchunk()
		im=datam[:,:,0,:]
		qm=datam[:,:,1,:]
		um=datam[:,:,2,:]
		vm=datam[:,:,3,:]
		posi=np.where(datai<(sigma*rmsi))
		posq=np.where(np.abs(dataq)<(sigma*rmsq))
		posu=np.where(np.abs(datau)<(sigma*rmsu))
		posv=np.where(np.abs(datav)<(sigma*rmsv))
		if rms_thresh==True:
			posi=np.where(datai<(sigma*rmsi))
			posq=np.where(np.abs(dataq)<(3*rmsq))
			posu=np.where(np.abs(datau)<(3*rmsu))
			posv=np.where(np.abs(datav)<(3*rmsv))
			qm[posi]=0
			qm[posq]=0
			um[posi]=0
			um[posu]=0
			vm[posi]=0
			vm[posv]=0
		qm[postb2]=0
		um[postb2]=0
		if 'Q' in stokes:
			datam[:,:,1,:]=qm
		if 'U' in stokes:
			datam[:,:,2,:]=um
		if 'V' in stokes:
			datam[:,:,3,:]=vm
		ia.putchunk(datam)
		ia.close()
		if overwrite==True:
			if os.path.isdir(imagename.split('.image')[0]+'_quvcor_surface.image')==True:
				os.system('cp -r '+imagename.split('.image')[0]+'_quvcor_surface.image '+imagename)
				os.system('rm -rf '+imagename.split('.image')[0]+'_quvcor_surface.image')
			if os.path.isdir(modelname.split('.model')[0]+'_quvcor_surface.model')==True:
				os.system('cp -r '+modelname.split('.model')[0]+'_quvcor_surface.model '+modelname)
				os.system('rm -rf '+modelname.split('.model')[0]+'_quvcor_surface.model')
		os.system('rm -rf casa*log '+modelname.split('.model')[0]+'*_quvcor_surface.model '+imagename.split('.image')[0]+'*_quvcor_surface.image'+\
				modelname.split('.model')[0]+'*_vcor_surface.model '+imagename.split('.image')[0]+'*_vcor_surface.image')
		os.system('rm -rf '+outfile_path+'/I_* '+outfile_path+'/I.*')
		return imagename,modelname,q_change_frac,u_change_frac,v_change_frac

	def correct_solar_quv_leakage(self,imagename,modelname,sigma,overwrite=False,stokes='QUV',v_frac_threshold=0.01,do_fluxcal=False,rms_thresh=True):
		'''
		Function to correct for solar Stokes Q, U and V leakage based on median leakage 
		(It is based on the fact that we do not expect any linear polarisation from the Quiet Sun emission)

		Parameters
		----------
		imagename : str 
			Name of image
		modelname : str 
			Name of the model
		sigma : float 
			N-sigma threshold to choose Stokes I emission region
		overwrite : bool 
			Overwrite the image and model or not
		stokes : str
			Stokes plane to correct
		v_frac_threshold: float
			Stokes V fraction threshold for quiet Sun, which should be less than 1 (default : 0.01)
		do_fluxcal : bool
			Perform flux calibration
		rms_thresh : bool
			Put a cutoff on the final leakage corrected model based on rms 
		Returns
		-------
		str
			Stokes Q ,U, V leakage corrected image name
		str 
			Stokes Q ,U, V leakage corrected model name 
		float
			Stokes I to Q leakage
		float 
			Stokes I to U leakage
		float
			Stokes I to V leakage 
		'''
		if sigma<5:
			sigma=5
		q_leakage,u_leakage,v_leakage=0,0,0
		self.log_verbose.info('Correcting Stokes I to Stokes '+','.join(list(stokes))+'.....\n')
		self.log_verbose.info('Image name : '+imagename+'\n')
		self.log_verbose.info('Model name : '+modelname+'\n')
		if os.path.exists(imagename.split('.image')[0]+'_quvcor.image'):
			os.system('rm -rf '+imagename.split('.image')[0]+'_quvcor.image')
		if os.path.exists(modelname.split('.model')[0]+'_quvcor.model'):
			os.system('rm -rf '+modelname.split('.model')[0]+'_quvcor.model')
		if overwrite==False:
			os.system('cp -r '+imagename+' '+imagename.split('.image')[0]+'_quvcor.image')
			os.system('cp -r '+modelname+' '+modelname.split('.model')[0]+'_quvcor.model')
			imagename=imagename.split('.image')[0]+'_quvcor.image'
			modelname=modelname.split('.model')[0]+'_quvcor.model'
		outfile_path=os.path.dirname(os.path.realpath(imagename))
		q_leakage,u_leakage=self.cal_solar_qu_leakage(imagename,sigma=sigma,do_fluxcal=do_fluxcal)
		v_leakage=self.cal_solar_v_leakage(imagename,sigma=sigma,do_fluxcal=do_fluxcal)
		if abs(v_leakage)<v_frac_threshold:
			self.log_verbose.info('Determined V leakage : '+str(round(v_leakage*100,3))+' % is less than '+str(round(v_frac_threshold*100,3))+\
				' %. Do not correcting Stokes V image for leakage.\n')
			v_leakage=0.0
		self.log_verbose.info('Stokes I to Q leakage : '+str(round(abs(q_leakage)*100,3))+'%\n')
		self.log_verbose.info('Stokes I to U leakage : '+str(round(abs(u_leakage)*100,3))+'%\n')
		self.log_verbose.info('Stokes I to V leakage : '+str(round(abs(v_leakage)*100,3))+'%\n')
		ia=image()
		ia.open(imagename) # Correcting image
		data=ia.getchunk()
		datai=data[:,:,0,:]
		dataq=data[:,:,1,:]-(q_leakage*datai)
		datau=data[:,:,2,:]-(u_leakage*datai)
		datav=data[:,:,3,:]-(v_leakage*datai)
		data[:,:,1,:]=dataq
		data[:,:,2,:]=datau
		data[:,:,3,:]=datav
		ia.putchunk(data)
		ia.close()
		pol_posv=np.where(np.abs(datav/datai)>1)
		rmsv=imstat(imagename=imagename,box=self.rms_box,stokes='V')['rms'][0]
		rmsi=imstat(imagename=imagename,box=self.rms_box,stokes='I')['rms'][0]
		# Determining Tb map
		####################
		imsubimage(imagename=imagename,outfile=outfile_path+'/I.image',stokes='I',dropdeg=False)
		if os.path.exists(outfile_path+'/QU.image'):
			os.system('rm -rf '+outfile_path+'/QU.image')
		imsubimage(imagename=imagename,outfile=outfile_path+'/QU.image',stokes='QU')
		if os.path.exists(outfile_path+'/polL.image'):
			os.system('rm -rf '+outfile_path+'/polL.image')
		immath(imagename=outfile_path+'/QU.image',outfile=outfile_path+'/polL.image',mode='poli')
		os.system('rm -rf '+outfile_path+'/QU.image')
		rmsl=imstat(imagename=outfile_path+'/polL.image',box=self.rms_box_nearsun)['rms'][0]
		ia.open(outfile_path+'/polL.image')
		datal=ia.getchunk()[:,:,0,:]
		ia.close()
		pol_posl=np.where(np.abs(datal/datai)>1)
		if do_fluxcal:
			atten=float(fits.getheader(self.metafits)['ATTEN_DB'])
			fluxcal_image=self.mwa_solar_fluxcal(imagename=outfile_path+'/I.image',outfile=outfile_path+'/I_fluxcal.image',atten=atten)
		else:
			fluxcal_image=outfile_path+'/I.image'
		major=imhead(imagename=fluxcal_image)['restoringbeam']['major']['value'] # In arcsec
		minor=imhead(imagename=fluxcal_image)['restoringbeam']['minor']['value'] # In arcsec
		freq=imhead(imagename=fluxcal_image)['refval'][-1]/10**9 # In GHz
		ref_beam_axes_multi=500000
		beam_axes_multi=major*minor	
		scale_tb_limit=(beam_axes_multi/ref_beam_axes_multi)
		expr='1.222e6*IM0/'+str(freq)+'^2/('+str(major*minor)+')'
		if os.path.exists(outfile_path+'/I_Tb.image'):
			os.system('rm -rf '+outfile_path+'/I_Tb.image')
		immath(imagename=fluxcal_image,outfile=outfile_path+'/I_Tb.image',mode='evalexpr',expr=expr)
		imhead(imagename=outfile_path+'/I_Tb.image', mode='put', hdkey='bunit', hdvalue='K')
		ia.open(outfile_path+'/I_Tb.image')
		datatb=ia.getchunk()
		ia.close()
		posi20=np.where(datai[:,:,0]<(20*rmsi))	
		tb_nonpos=np.where(datatb[:,:,0,0]<=0)
		datatb[:,:,0,0][tb_nonpos]=np.nan
		datatblog=datatb[:,:,0,0]/10**6
		pos=np.where(np.isnan(datatblog)==True)
		datatblog[posi20]=np.nan
		maxpos=imstat(imagename=imagename,stokes='I')['maxpos']
		box=self.negative_box(maxpos,box_width=3).split(',')
		box_coords=[int(i) for i in box]
		for i in range(datatblog.shape[0]):
			for j in range(datatblog.shape[1]):
				if i<box_coords[0] or i>box_coords[2] or j<box_coords[1] or j>box_coords[3]:
					datatblog[i,j]=np.nan
				else:
					pass	
		mediantb=np.nanmedian(datatblog)
		if mediantb>5:
			print ('No quiet Sun region found.\n')
			q_leakage=u_leakage=v_leakage=0
			if overwrite==True:
				if os.path.isdir(imagename.split('.image')[0]+'_quvcor_surface.image')==True:
					os.system('cp -r '+imagename.split('.image')[0]+'_quvcor_surface.image '+imagename)
					os.system('rm -rf '+imagename.split('.image')[0]+'_quvcor_surface.image')
				if os.path.isdir(modelname.split('.model')[0]+'_quvcor_surface.model')==True:
					os.system('cp -r '+modelname.split('.model')[0]+'_quvcor_surface.model '+modelname)
					os.system('rm -rf '+modelname.split('.model')[0]+'_quvcor_surface.model')
			os.system('rm -rf casa*log '+modelname.split('.model')[0]+'_quvcor_surface.model '+imagename.split('.image')[0]+'_quvcor_surface.image'+\
					modelname.split('.model')[0]+'_vcor_surface.model '+imagename.split('.image')[0]+'_vcor_surface.image')
			os.system('rm -rf '+outfile_path+'/I_* '+outfile_path+'/I.* '+outfile_path+'/polL.image')
			return imagename,modelname,q_leakage,u_leakage,v_leakage
		if scale_tb_limit<1:
			postb=np.where(np.round(datatblog,0)<=1)
		else:
			postb=np.where(np.round(datatblog,0)<=scale_tb_limit)
		#####################
		ia.open(modelname) # Correcting model
		datam=ia.getchunk()
		im=datam[:,:,0,:]
		qm=datam[:,:,1,:]-(q_leakage*im)
		um=datam[:,:,2,:]-(u_leakage*im)
		vm=datam[:,:,3,:]-(v_leakage*im)
		posi=np.where(datai<(sigma*rmsi))
		posq=posu=np.where(np.abs(datal)<(20*rmsl))		
		posv=np.where(np.abs(datav)<(sigma*rmsv))
		qm[posi]=0
		qm[posq]=0
		um[posi]=0
		um[posu]=0
		if rms_thresh==True:
			vm[posi]=0
			vm[posv]=0
		qm[postb]=0
		um[postb]=0
		dataq[posi]=np.nan
		dataq[posq]=np.nan
		dataq[postb]=np.nan
		if 'Q' in stokes:
			datam[:,:,1,:]=qm
		if 'U' in stokes:
			datam[:,:,2,:]=um
		if 'V' in stokes:
			datam[:,:,3,:]=vm
		ia.putchunk(datam)
		ia.close()
		if overwrite==True:
			if os.path.isdir(imagename.split('.image')[0]+'_quvcor_surface.image')==True:
				os.system('cp -r '+imagename.split('.image')[0]+'_quvcor_surface.image '+imagename)
				os.system('rm -rf '+imagename.split('.image')[0]+'_quvcor_surface.image')
			if os.path.isdir(modelname.split('.model')[0]+'_quvcor_surface.model')==True:
				os.system('cp -r '+modelname.split('.model')[0]+'_quvcor_surface.model '+modelname)
				os.system('rm -rf '+modelname.split('.model')[0]+'_quvcor_surface.model')
		os.system('rm -rf casa*log '+modelname.split('.model')[0]+'_quvcor_surface.model '+imagename.split('.image')[0]+'_quvcor_surface.image'+\
				modelname.split('.model')[0]+'_vcor_surface.model '+imagename.split('.image')[0]+'_vcor_surface.image')
		os.system('rm -rf '+outfile_path+'/I_* '+outfile_path+'/I.* '+outfile_path+'/polL.image')
		return imagename,modelname,q_leakage,u_leakage,v_leakage

	def ant_to_baselines(self,ant_list=[],include_baselines=True,total_ant=0,strict=False):
		'''
		Function to generate baseline list originated from list of antennas
		Parameters
		----------
		ant_list : list
			Antenna list
		include_baselines : bool
			Include the baseline in the list or exclude then from all baselines
		total_ant : int
			Total number of antennas (default : 0, calculates from the measurement set)	
		strict : bool
			Use only the baselines between antenna list or not		
		Returns
		-------
		list
			List of baselines
		'''
		if total_ant==0:
			AM=am.AccessMS(self.msname)
			total_ant=AM.get_num_antenna()
		ant_list=[int(i) for i in ant_list]
		if len(ant_list)!=0 or len(ant_list)!=total_ant:
			if include_baselines:
				selected_baselines=[]
				for i in range(total_ant):
					for j in range(total_ant):
						if strict==False:
							if i<j and i in ant_list:
								selected_baselines.append(str(i)+'&'+str(j))
						else:
							if i<j and i in ant_list and j in ant_list:
								selected_baselines.append(str(i)+'&'+str(j))
			else:
				selected_baselines=[]
				for i in range(total_ant):
					for j in range(total_ant):
						if strict==False:
							if i<j and i not in ant_list:
								selected_baselines.append(str(i)+'&'+str(j))
						else:
							if i<j and (i not in ant_list or j not in ant_list):
								selected_baselines.append(str(i)+'&'+str(j))
			return selected_baselines
		else:
			all_baselines=[]
			for i in range(total_ant):
				for j in range(total_ant):
					if i<j:
						all_baselines.append(str(i)+'&'+str(j))
			return all_baselines

	def selfcal_iteration(self,num_iter,rms_thresh,sigma,maskstr,antenna_to_use,startmodel,startmask,ref_ant,minsnr,calmode,maskfile='',want_auto_masking=False,cpus=3,absmem=10,\
		wlayers=1,stokes='I',use_ankflagger=True,do_flag=False,do_bandpass=False,correct_phasecenter=False,ra=0,dec=0,box_width=3,calibrator_caltable=[],maskregion='',\
		weight='briggs',robust=0.5,modify_datacolumn=False,mwapb_cor=False,do_solarqu_cor=False,do_fluxcal=False,solar_imaging=True,restore_flag=True,use_cubical=False,solint='inf',\
		mwapb_cor_restore=False,quvcor_stokes='QUV',applymode='calflag',solmode='R',solnorm=True,normtype='mean'):
		'''
		Function to perform a self-calibration loop, make an image, put the model in the measurement set, and perform the calibration

		Parameters
		----------
		num_iter : int 
			Number of self-calibration iteration
		rms_thresh : list 
			List of rms for threshold
		sigma : float 
			Threshold sigma
		maskstr : str 
			CASA mask string for CLEANing (For CASA only)
		antenna_to_use : list 
			List of antennas for CLEANing
		startmodel : str 
			CASA model to start the CLEANing (For both CASA imager and WSClean, for WSClean imager, CASA models will be converted to WSClean models)
		startmask :  str
			CASA mask to start (For CASA only)
		ref_ant : int
			Reference antenna
		minsnr : float 
			Minimum SNR of gain solutions
		calmode: str 
			'p' for phase-only and 'ap' for amplitude-phase calibration
		maskfile : str 
			Name of the maskfile
		want_auto_masking : bool
			If True use CASA auto-multithresh for auto masking or WSclean auto-masking
		cpus : int
			Number of cpu threads to use
		absmem : float
			Memory in GB to use
		wlayers : int
			Number of w-stacking layers (For wsclean only) 
		stokes : str 
			Stokes plane to image
		use_ankflagger : bool 
			Use aNKflagger for flagging after each selfcal round
		do_flag: bool
			Flag after selfcal round
		do_bandpass : str
			Perform bandpass calibration
		correct_phasecenter : bool 
			Correct the phase center
		ra : float 
			New RA in degree to change phasecenter
		dec : float 
			New DEC in degree to change phasecenter
		box_width : float 
			Negative box width around the maximum pixel in degree (default : 3 degree)
		calibrator_caltable : list 
			List of calilbrator caltables
		maskregion : str 
			Mask region in case of auto-masking (Only for CASA)
		weight : str
			Visibility weighting during imaging , 'uniform', 'natural' or 'briggs'. Default is 'briggs'
		robust : float
			Robust parameter for briggs weighting, default : 0.5
		modify_datacolumn : bool
			Modify data column
		mwapb_cor : bool
			Perform MWA primary beam correction
		do_fluxcal : bool
			Perform flux calibration or not
		do_solarqu_cor : bool
			Correct solar Stokes I to Q,U imaged based leakage correction
		solar_imaging : bool 
			Performing Solar image calibration		
		restore_flag : bool
			Restore previous round flags
		use_cubical : bool
			Use CubiCal for calibration (default : True), otherwise use CASA
		solint : float
			Solution interval in second (default : 'inf')
		mwapb_cor_restore : bool
			Restore MWA PB correction
		quvcor_stokes : str
			Stokes plane to perform image based leakage correction
		applymode : str
			Calibration apply mode ('calonly' or 'calflag')
		solmode : str
			Solver mode (default : 'R')
		solnorm : bool
			Perform normalized gain solutions or not
		normtype : str
			Gain normalisation, mean or median
		Returns
		-------
		float
			rms based dynamic range
		float
			Negative based dynamic range
		int
			Message code
		'''
		os.chdir(self.mspath)
		imagename=self.msname.split('.ms')[0]+'_'+str(num_iter)+'_'+str(self.bandwidth)+'kHz_'+str(self.time_width)+'s' # Imagename prefix
		present_file=glob.glob(imagename+'*')
		if len(present_file)!=0:
			os.system('rm -rf '+imagename+'*')
			self.log_verbose.info('rm -rf '+imagename+'*\n')
		caltable_name=self.msname.split('.ms')[0]+'.cal' # Caltable name
		if os.path.isdir(caltable_name):
			os.system('rm -rf '+caltable_name)
		threshold=[str(rms*sigma)+'Jy' for rms in rms_thresh]
		# Making image
		if do_bandpass:
			selfcalmode='Bandpass'
		else:
			selfcalmode='Intensity'
		self.log_verbose.info('============================\n')
		self.log_verbose.info(selfcalmode+' Selfcal iteration number : '+str(num_iter)+'\n')
		self.log_verbose.info('============================\n')
		if calmode=='p':
			clean_gain=0.05
			smallscalebias=0.8
		else:
			clean_gain=0.1
			smallscalebias=0.6
		if stokes=='I':
			stokes_list=['I']
			pol='i'
		elif stokes=='RR':
			stokes_list=['RR']
			pol='rr'
		elif stokes=='LL':
			stokes_list=['LL']
			pol='ll'
		elif stokes=='XX':
			stokes_list=['XX']
			pol='xx'
		elif stokes=='YY':
			stokes_list=['YY']
			pol='yy'
		elif stokes=='RRLL':
			stokes_list=['RR','LL']
			pol='rrll'
		elif stokes=='XXYY':
			stokes_list=['XX','YY']
			pol='xxyy'
		elif stokes=='IQ':
			stokes_list=['I','Q']
			pol='iq'
		elif stokes=='IQUV':
			stokes_list=['I','Q','U','V']
			pol='iquv'
		else:
			if self.verbose==False:
				print('Specify correct stokes list.\n')
			self.log_verbose.error('Specify correct stokes list.\n')
			return 14
		if self.wsclean==False or self.wsclean_start==None:
			if maskfile=='':
				maskfile=maskstr
			robust=robust*2
			if maskfile!='': # Using the user provided mask file
				self.log_verbose.info('poltclean(vis=\''+self.msname+'\',imagename=\''+imagename+'\',selectdata=True,startmodel=\''+startmodel+'\',startmask=\''+\
						startmask+'\',stokes=\''+stokes+'\',antenna=\''+antenna_to_use+'\',imsize=['+str(self.imsize)+'],cell=\''+str(self.cellsize)+'arcsec'\
						+'\',niter=10000,gain='+str(clean_gain)+',threshold='+str(threshold)+',deconvolver=\'multiscale\',scales='+str(self.multiscale_scales)
						+',uvtaper=\''+self.uvtaper+'\',weighting=\''+weight+'\',robust='+str(robust)+',interactive=False,mask=\''+str(maskfile)\
						+'\',uvrange=\''+str(self.imaging_uvrange_lambda)+'\',smallscalebias='+str(smallscalebias)+')\n')
				poltclean(vis=self.msname,imagename=imagename,selectdata=True,startmodel=startmodel,startmask=startmask,stokes=stokes,antenna=antenna_to_use,imsize=[self.imsize],\
					cell=str(self.cellsize)+'arcsec',niter=10000,gain=clean_gain,threshold=threshold,deconvolver='multiscale',scales=self.multiscale_scales,uvtaper=self.uvtaper,\
					weighting=weight,robust=robust,interactive=False,mask=maskfile,uvrange=self.imaging_uvrange_lambda,smallscalebias=smallscalebias)
			elif want_auto_masking==True and maskfile=='' and maskstr=='': # Use auto-masking
				try_count=0
				if startmodel!='': # Add auto masking safety
					automask_trials=2
				else:
					automask_trials=10
				while True:
					if try_count==0:
						self.log_verbose.info('Normal auto-masking.\n')
						self.log_verbose.info('poltclean(vis=\''+self.msname+'\',imagename=\''+imagename+'\',selectdata=True,startmodel=\''+startmodel+'\',startmask=\''\
							+startmask+'\',stokes=\''+stokes+'\',antenna=\''+antenna_to_use+'\',imsize=['+str(self.imsize)+'],cell=\''+str(self.cellsize)+'arcsec'+\
							'\',niter=10000,gain='+str(clean_gain)+',threshold='+str(threshold)+',deconvolver=\'multiscale\',scales='+str(self.multiscale_scales)+\
							',uvtaper=\''+str(self.uvtaper)+'\',weighting=\''+weight+'\',robust='+str(robust)+',interactive=False,usemask='+\
							'\'auto-multithresh\',mask=\'\',pbmask=0.0,sidelobethreshold=3.0,noisethreshold='+str(float(sigma))+',lownoisethreshold='+str(float(sigma/3.0))+\
							',negativethreshold=0.0,smoothfactor=1.0,minbeamfrac=0.1,growiterations=75,minpercentchange=5.0,uvrange=\''+str(self.imaging_uvrange_lambda)+'\',automask_trials='+\
							str(automask_trials)+',maskregion=\''+maskregion+'\',smallscalebias='+str(smallscalebias)+')\n')
						poltclean(vis=self.msname,imagename=imagename,selectdata=True,startmodel=startmodel,startmask=startmask,stokes=stokes,antenna=antenna_to_use,\
						imsize=[self.imsize],cell=str(self.cellsize)+'arcsec',niter=10000,gain=clean_gain,threshold=threshold,deconvolver='multiscale',scales=self.multiscale_scales,\
						uvtaper=self.uvtaper,weighting=weight,robust=robust,interactive=False,usemask='auto-multithresh',mask='',pbmask=0.0,sidelobethreshold=3.0,\
						noisethreshold=float(sigma),lownoisethreshold=float(sigma/3.0),negativethreshold=0.0,smoothfactor=1.0,minbeamfrac=0.5,growiterations=75,minpercentchange=5.0,\
						uvrange=self.imaging_uvrange_lambda,automask_trials=automask_trials,maskregion=maskregion,smallscalebias=smallscalebias)
					elif try_count==1:
						self.log_verbose.info('Trying with auto-masking with no restriction of minimum beam fraction.\n')
						self.log_verbose.info('poltclean(vis=\''+self.msname+'\',imagename=\''+imagename+'\',selectdata=True,startmodel=\''+startmodel\
							+'\',startmask=\''+startmask+'\',stokes=\''+stokes+',antenna=\''+antenna_to_use+'\',imsize=['+str(self.imsize)+'],cell=\''+str(self.cellsize)+'arcsec'+\
							'\',niter=10000,gain='+str(clean_gain)+',threshold='+str(threshold)+',deconvolver=\'multiscale\',scales='+str(self.multiscale_scales)+\
							',uvtaper=\''+str(self.uvtaper)+'\',weighting=\''+weight+'\',robust='+str(robust)+',interactive=False,usemask=\''+\
							'auto-multithresh\',mask=\'\',pbmask=0.0,sidelobethreshold=3.0,noisethreshold='+str(float(sigma))+\
							',lownoisethreshold='+str(float(sigma/3.0))+',negativethreshold=0.0,smoothfactor=1.0,'+\
							'minbeamfrac=0.1,growiterations=75,minpercentchange=-1.0,uvrange=\''+str(self.imaging_uvrange_lambda)+'\',automask_trials='+str(automask_trials)+\
							',maskregion=\''+maskregion+'\',smallscalebias='+str(smallscalebias)+')\n')
						poltclean(vis=self.msname,imagename=imagename,selectdata=True,startmodel=startmodel,startmask=startmask,stokes=stokes,antenna=antenna_to_use,\
						imsize=[self.imsize],cell=str(self.cellsize)+'arcsec',niter=10000,gain=clean_gain,threshold=threshold,deconvolver='multiscale',scales=self.multiscale_scales,\
						uvtaper=self.uvtaper,weighting=weight,robust=robust,interactive=False,usemask='auto-multithresh',mask='',pbmask=0.0,sidelobethreshold=3.0,\
						noisethreshold=float(sigma),lownoisethreshold=float(sigma/3.0),negativethreshold=0.0,smoothfactor=1.0,minbeamfrac=0.1,growiterations=75,minpercentchange=-1.0,\
						uvrange=self.imaging_uvrange_lambda,automask_trials=automask_trials,maskregion=maskregion,smallscalebias=smallscalebias)
					elif try_count==2:
						self.log_verbose.info('Trying without masking.\n')
						self.log_verbose.info('poltclean(vis=\''+self.msname+'\',imagename=\''+imagename+'\',selectdata=True,startmodel=\''+startmodel+'\',startmask=\''+\
							startmask+'\'stokes=\''+stokes+'\',antenna=\''+antenna_to_use+'\',imsize=['+str(self.imsize)+'],cell=\''+str(self.cellsize)+'arcsec'+\
							'\',niter=10000,gain='+str(clean_gain)+',threshold='+str(threshold)+',deconvolver=\'multiscale\',scales='+str(self.multiscale_scales)+\
							',uvtaper=\''+str(self.uvtaper)+'\',weighting=\''+weight+'\',robust='+\
							str(robust)+',interactive=False,usemask=\'user\',mask=\'\',uvrange=\''+str(self.imaging_uvrange_lambda)+'\',smallscalebias='+str(smallscalebias)+')\n')
						poltclean(vis=self.msname,imagename=imagename,selectdata=True,startmodel=startmodel,startmask=startmask,stokes=stokes,antenna=antenna_to_use,\
						imsize=[self.imsize],cell=str(self.cellsize)+'arcsec',niter=10000,gain=clean_gain,threshold=threshold,deconvolver='multiscale',scales=self.multiscale_scales,\
						uvtaper=self.uvtaper,weighting=weight,robust=robust,interactive=False,usemask='user',mask='',uvrange=self.imaging_uvrange_lambda,smallscalebias=smallscalebias)
					else:
						break
					modelflux=imstat(imagename=imagename+'.model')['sum'][0]
					if modelflux==0.0:
						try_count+=1
					else:
						break
			else: # If no masking
				self.log_verbose.info('poltclean(vis=\''+self.msname+'\',imagename=\''+imagename+'\',selectdata=True,startmodel=\''+startmodel+'\',dtartmask=\''+startmask\
						+'\',stokes=\''+stokes+'\',antenna=\''+antenna_to_use+'\',imsize=['+str(self.imsize)+'],cell=\''+str(self.cellsize)+'arcsec'\
						+'\',niter=10000,gain='+str(clean_gain)+',threshold=\''+str(threshold)+'\',deconvolver=\'multiscale\',scales='+str(self.multiscale_scales)
						+',uvtaper=\''+self.uvtaper+'\',weighting=\''+weight+'\',robust='+str(robust)+',interactive=False,mask='+str([maskfile])+',uvrange=\''\
						+str(self.imaging_uvrange_lambda)+'\',smallscalebias='+str(smallscalebias)+')\n')		
				poltclean(vis=self.msname,imagename=imagename,selectdata=True,startmodel=startmodel,startmask=startmask,stokes=stokes,antenna=antenna_to_use,imsize=[self.imsize],\
					cell=str(self.cellsize)+'arcsec',niter=10000,gain=clean_gain,threshold=threshold,deconvolver='multiscale',scales=self.multiscale_scales,uvtaper=self.uvtaper,\
					weighting=weight,robust=robust,interactive=False,mask='',uvrange=self.imaging_uvrange_lambda,smallscalebias=smallscalebias)
			casa_imagename=imagename+'.image'
			casa_modelname=imagename+'.model'
			casa_residualname=imagename+'.residual'
		else:
			AM=am.AccessMS(self.msname)
			nant=AM.get_num_antenna()
			if antenna_to_use=='':
				ant_list=[]
				for i in range(nant):
					ant_list.append(str(i))
				self.log_verbose.info('Antennas in use : All\n')
			else:
				self.log_verbose.info('Antennas in use : '+antenna_to_use+'\n')
				ant_list=antenna_to_use.split(',')
			baseline_flagged=False
			if len(ant_list)<nant and calmode=='p':
				baselines=self.ant_to_baselines(ant_list=ant_list,include_baselines=False,total_ant=nant)
				baselines_str=';'.join(baselines)
				flagdata(vis=self.msname,mode='manual',antenna=baselines_str,flagbackup=True)
				baseline_flagged=True
			scales=[str(i) for i in self.multiscale_scales]
			if weight=='briggs':
				weight=weight+' '+str(robust)
			taper_tukey=10
			wsclean_args=['-scale '+str(self.cellsize)+'asec','-size '+str(self.imsize)+' '+str(self.imsize),'-no-dirty','-j '+str(cpus),'-abs-mem '+str(absmem),'-weight '+weight,'-taper-edge '+str(self.imaging_maxuv_lambda),\
			'-taper-edge-tukey '+str(taper_tukey),'-name '+imagename,'-nwlayers '+str(wlayers),'-pol '+str(pol),'-minuv-l '+str(self.imaging_minuv_lambda),'-niter 100000','-mgain 0.8','-gain '+str(clean_gain),\
			'-multiscale','-multiscale-scales '+','.join(scales),'-auto-mask '+str(sigma),'-multiscale-gain 0.05']#,'-use-weights-as-taper']
			if calmode=='ap' or do_bandpass==True:
				wsclean_args.append('-auto-threshold '+str(3))
				wsclean_args.append('-multiscale-scale-bias 0.9')
			else:
				wsclean_args.append('-multiscale-scale-bias 1.0')
				wsclean_args.append('-auto-threshold '+str(sigma-0.5))
			if pol=='i' or pol=='xxyy':
				wsclean_args.append('-no-negative')
			if startmodel!='':
				if self.verbose:
					self.log_verbose.info('casaimage_to_wsclean(casaimage=\''+startmodel+'\',wsclean_imagename=\''+imagename+'\',imagetype=\'model\',wsclean-imsize='\
										+str(self.imsize)+',keep_casaimage=True,stokes=\''+stokes+'\')\n')
				wsclean_startmodels,wsclean_imsize,wsclean_stokes=self.casaimage_to_wsclean(casaimage=startmodel,wsclean_imagename=imagename,imagetype='model',\
										wsclean_imsize=self.imsize,keep_casaimage=True,stokes=stokes)
				if len(wsclean_startmodels)==len(stokes_list) and wsclean_imsize==self.imsize and wsclean_stokes==stokes_list:
					wsclean_args.append('-continue')
				elif len(wsclean_startmodels)!=0:
					for i in wsclean_startmodels:
						os.system('rm -rf '+i)
			wsclean_args.append('-quiet')
			if self.verbose:
				self.log_verbose.info('wsclean '+' '.join(wsclean_args)+' '+self.msname+'\n')			
			os.system(self.wsclean_start+' '+' '.join(wsclean_args)+' '+self.msname)	
			if baseline_flagged:
				flaglist=flagmanager(vis=self.msname,mode='list')
				flaglist_keys=list(flaglist.keys())
				flaglist_keys.remove('MS')
				if len(flaglist_keys)>0:
					last_flag_key=len(flaglist_keys)-1
					last_flagversion=flaglist[last_flag_key]['name']
					if 'flagdata' in last_flagversion:  # Restore flagversion if the last flagversion is from applycal
						self.log_verbose.info('flagmanager(vis=\''+self.msname+'\',mode=\'restore\',versionname=\''+str(last_flagversion)+'\',merge=\'replace\')\n')
						self.log_verbose.info('flagmanager(vis=\''+self.msname+'\',mode=\'delete\',versionname=\''+str(last_flagversion)+'\')\n')
						flagmanager(vis=self.msname,mode='restore',versionname=last_flagversion,merge='replace')
						flagmanager(vis=self.msname,mode='delete',versionname=last_flagversion)		
			wsclean_images=glob.glob(imagename+'*image.fits')
			wsclean_models=glob.glob(imagename+'*model.fits')
			wsclean_residuals=glob.glob(imagename+'*residual.fits')
			self.log_verbose.info('wsclean_to_casaimage(wsclean_images='+str(wsclean_images)+',casaimage_prefix=\''+imagename+'\',imagetype=\'image\',keep_wsclean_images=False)\n')
			casa_imagename=self.wsclean_to_casaimage(wsclean_images=wsclean_images,casaimage_prefix=imagename,imagetype='image',keep_wsclean_images=False)
			self.log_verbose.info('wsclean_to_casaimage(wsclean_images='+str(wsclean_models)+',casaimage_prefix=\''+imagename+'\',imagetype=\'model\',keep_wsclean_images=False)\n')
			casa_modelname=self.wsclean_to_casaimage(wsclean_images=wsclean_models,casaimage_prefix=imagename,imagetype='model',keep_wsclean_images=False)
			self.log_verbose.info('wsclean_to_casaimage(wsclean_images='+str(wsclean_residuals)+',casaimage_prefix=\''+imagename+'\',imagetype=\'residual\',keep_wsclean_images=False)\n')
			casa_residualname=self.wsclean_to_casaimage(wsclean_images=wsclean_residuals,casaimage_prefix=imagename,imagetype='residual',keep_wsclean_images=False)	
		out_dict,negative_dyn_range=self.calc_dyn_range(casa_imagename,sigma,box_width=box_width,stokes_list=stokes_list) # Calculating the dynamic range of the image
		out_dict_keys=out_dict.keys()
		if calmode=='p':
			self.log_verbose.info('remove_outside_model(\''+casa_imagename+'\',\''+casa_modelname+'\',box_width=5,overwrite=True)\n')			
			self.remove_outside_model(casa_imagename,casa_modelname,box_width=5,overwrite=True)			
		if 'NAN' in out_dict_keys:
			self.log_verbose.error(B.error_msgs(3))
			os.system('rm -rf casa*log')
			os.chdir(self.cwd)
			os.system('rm -rf casa*log')
			return 3     # If image is not made, no point in continuing
		if os.path.isdir(casa_modelname)==False:
			self.log_verbose.error(B.error_msgs(4)+'\n')
			os.system('rm -rf casa*log')
			os.chdir(self.cwd)
			os.system('rm -rf casa*log')
			return 4	   # If model is not present no point in continuing
		else:
			modelflux=imstat(imagename=casa_modelname)['sum'][0]
			if modelflux==0:
				self.log_verbose.error(B.error_msgs(5)+'\n')
				os.system('rm -rf casa*log')
				os.chdir(self.cwd)
				os.system('rm -rf casa*log')	
				return 5 # If modelflux is 0, no point in continuing
			else:
				if restore_flag:
					flaglist=flagmanager(vis=self.msname,mode='list')
					if len(flaglist)>1:
						flaglist_keys=list(flaglist.keys())
						flaglist_keys.remove('MS')
						last_flag_key=len(flaglist_keys)-1
						last_flagversion=flaglist[last_flag_key]['name']
						# Restore the flag and delete the present flag version
						if 'applycal' in last_flagversion or 'cubical' in last_flagversion or 'aNKflag' in last_flagversion or 'uvsub' in last_flagversion:  
																											# Restore flagversion if the last flagversion is from applycal or uvsub flagging
							self.log_verbose.info('flagmanager(vis=\''+self.msname+'\',mode=\'restore\',versionname=\''+str(last_flagversion)+'\',merge=\'replace\')\n')
							self.log_verbose.info('flagmanager(vis=\''+self.msname+'\',mode=\'delete\',versionname=\''+str(last_flagversion)+'\')\n')
							flagmanager(vis=self.msname,mode='restore',versionname=last_flagversion,merge='replace')
							flagmanager(vis=self.msname,mode='delete',versionname=last_flagversion)
				self.log_verbose.info('clearcal(vis=\''+self.msname+'\')\n')
				clearcal(vis=self.msname)
				self.log_verbose.info('delmod(vis=\''+self.msname+'\',scr=True)\n') 
				delmod(vis=self.msname,scr=True) # Clear the MODEL column
				if do_bandpass==True:
					self.remove_model_negative(casa_imagename,casa_modelname,sigma=3,overwrite=True) # Removing negatives from model and less than 3 sigma regions
				elif do_bandpass==False and calmode=='p':
					self.remove_model_negative(casa_imagename,casa_modelname,sigma=sigma,overwrite=True) # Removing negatives from model and less than 7 regions
				else:
					self.remove_model_negative(casa_imagename,casa_modelname,sigma=5,overwrite=True) # Removing negatives from model and less than 5 sigma regions
				if correct_phasecenter==True:
					if ra==0 or dec==0:
						AM=am.AccessMS(self.msname)
						radec_str,ra,dec=AM.get_phasecenter()
					shifted=self.shift_phasecenter(imagename+'.model',ra,dec)
					if shifted==0:
						self.log_verbose.info('Image phasecenter changed to : RA = '+str(ra)+' deg, DEC = '+str(dec)+' deg\n')
					elif shifted==1:
						self.log_verbose.info('Shifting is not required. Continuing with previous phasecenter.\n')
					else:
						self.log_verbose.info('Error in image phsecenter shifting. Continuing with previous phasecenter.\n')
				if stokes=='IQUV':
					if mwapb_cor==True:
						if self.verbose:
							self.log_verbose.info('mwapb_cor(\''+casa_imagename+'\',\''+'mwapbcor\',\''+self.metafits+'\',verbose='+str(self.verbose)+',nthreads='+str(cpus)+')\n')
						pbcor_casa_imagename=pb.mwapb_cor(casa_imagename,'mwapbcor',self.metafits,verbose=self.verbose,nthreads=cpus)
						if self.verbose:
							self.log_verbose.info('mwapb_cor(\''+casa_modelname+'\',\''+'mwapbcor\',\''+self.metafits+'\',verbose='+str(self.verbose)+',nthreads='+str(cpus)+')\n')
						pbcor_casa_modelname=pb.mwapb_cor(casa_modelname,'mwapbcor',self.metafits,verbose=self.verbose,nthreads=cpus)
					else:
						pbcor_casa_imagename=casa_imagename
						pbcor_casa_modelname=casa_modelname
					if do_solarqu_cor==True and solar_imaging==True:
						if self.verbose:
							self.log_verbose.info('correct_solar_quv_leakage(\''+pbcor_casa_imagename+'\',\''+pbcor_casa_modelname+'\','+str(10)+',overwrite=False,stokes=\''+\
								str(quvcor_stokes)+'\',do_fluxcal='+str(do_fluxcal)+')\n')
						quvcor_imagename,quvcor_modelname,q_leakage,u_leakage,v_leakage=self.correct_solar_quv_leakage(pbcor_casa_imagename,pbcor_casa_modelname,10,overwrite=False,\
																								stokes=quvcor_stokes,do_fluxcal=do_fluxcal)
						self.log_verbose.info('Q leakage : '+str(round(q_leakage*100,3))+', U leakage : '+str(round(u_leakage*100,3))+', V leakage : '+str(round(v_leakage*100,3))+' %\n')
						poldistortion_correction=False
						polmodel_threshold=False
						previous_image=''
						previous_model=''
						if polmodel_threshold==True:
							self.pol_model_threshold(quvcor_imagename,quvcor_modelname,3,1)
					else:
						quvcor_imagename=pbcor_casa_imagename
						quvcor_modelname=pbcor_casa_modelname
					if mwapb_cor_restore==True:
						if self.verbose:
							self.log_verbose.info('mwapb_cor(\''+quvcor_imagename+'\',\''+os.path.basename(casa_imagename).split('.image')[0]+'\',\''+self.metafits+\
										'\',verbose='+str(self.verbose)+',nthreads='+str(cpus)+',restore=True)\n')
						casa_imagename=pb.mwapb_cor(quvcor_imagename,os.path.basename(casa_imagename).split('.image')[0],self.metafits,verbose=self.verbose,nthreads=cpus,restore=True)
						if self.verbose:
							self.log_verbose.info('mwapb_cor(\''+quvcor_modelname+'\',\''+os.path.basename(casa_modelname).split('.model')[0]+'\',\''+self.metafits+\
										'\',verbose='+str(self.verbose)+',nthreads='+str(cpus)+',restore=True)\n')
						casa_modelname=pb.mwapb_cor(quvcor_modelname,os.path.basename(casa_modelname).split('.model')[0],self.metafits,verbose=self.verbose,nthreads=cpus,restore=True)
						if os.path.exists('quvcor.image'):
							os.system('rm -rf quvcor.image')
						os.system('cp -r '+quvcor_imagename+' quvcor.image')
						if os.path.exists('quvcor.model'):
							os.system('rm -rf quvcor.model')
						os.system('cp -r '+quvcor_modelname+' quvcor.model')
					else:
						if quvcor_imagename!=casa_imagename:
							if os.path.exists(casa_imagename):
								os.system('rm -rf '+casa_imagename)
							os.system('cp -r '+quvcor_imagename+' '+casa_imagename)
						if quvcor_modelname!=casa_modelname:
							if os.path.exists(casa_modelname):
								os.system('rm -rf '+casa_modelname)
							os.system('cp -r '+quvcor_modelname+' '+casa_modelname)
						if quvcor_imagename!='' and os.path.exists(quvcor_imagename) and quvcor_imagename!=casa_imagename:
							os.system('rm -rf '+quvcor_imagename)
						if quvcor_modelname!='' and os.path.exists(quvcor_modelname) and quvcor_modelname!=casa_modelname:
							os.system('rm -rf '+quvcor_modelname)
				'''if do_bandpass==True and do_fluxcal==True:
					atten=float(fits.getheader(self.metafits)['ATTEN_DB'])
					casa_modelname=self.mwa_solar_fluxcal(imagename=casa_modelname,outfile=casa_modelname,atten=atten)
					solnorm=False
				else:
					solnorm=True'''
				self.log_verbose.info('ft(vis=\''+self.msname+'\',model=\''+casa_modelname+'\',nterms=1,usescratch=True)\n') 
				ft(vis=self.msname,model=casa_modelname,nterms=1,usescratch=True) # Putting the model into MS
				if use_cubical:
					AM=am.AccessMS(self.msname)
					ref_ant_name=AM.get_antenna_name(antenna_id=int(ref_ant))
					timeres=AM.calc_timeres()
					nchan=AM.get_num_channels()
					if solint!='int' and solint!='inf':
						ntime=int(solint/timeres)
					elif solint=='inf':
						ntime=1
					else:
						ntime=AM.get_num_timestamps()
					if len(calibrator_caltable)>0:
						for i in range(len(calibrator_caltable)):
							if i==0:
								datacolumn='DATA'
							else:
								datacolumn='CORRECTED_DATA'
							cubical_args=['paircars_cubical','--data-ms '+ self.msname,'--out-casa-gaintables 0','--bbc-plot 0','--data-column '+datacolumn,'--sol-jones g',\
							'--flags-save FLAG','--g-xfer-from '+calibrator_caltable[i],'--out-mode ac','--out-overwrite 0','--out-backup 0','--g-clip-low 0','--g-clip-after 0']
							if self.verbose:
								self.log_verbose.info(' '.join(cubical_args)+'\n')
							os.system(' '.join(cubical_args)+' > tmp')
							if os.path.exists('cubical.cc-out'):
								os.system('rm -rf cubical.cc-out* tmp')
							flagdata(vis=self.msname,mode='clip',clipzeros=True,flagbackup=False)
						datacolumn='CORRECTED_DATA'
					else:
						datacolumn='DATA'
					if do_bandpass==True:
						af=agentflagger()
						af.open(self.msname)
						versionlist=af.getflagversionlist()
						if len(versionlist)!=0:
							for version_name in versionlist:
								if 'cubical' in version_name:
									version_num=int(version_name.split(':')[0].split(' ')[0].split('_')[-1])+1
								else:
									version_num=1
						else:
							version_num=1
						now = datetime.now()
						dt_string = now.strftime("%Y-%m-%d %H:%M:%S")
						af.saveflagversion('cubical_'+str(version_num),'Flags autosave on '+dt_string)
						af.done()
						os.system('rm -rf cubical.cc-out* tmp')
						cubical_args=['paircars_cubical','--data-ms '+ self.msname,'--out-casa-gaintables 0','--model-list MODEL_DATA','--bbc-plot 0','--g-type complex-diag','--data-column '+datacolumn,\
						'--sol-min-bl '+str(self.calib_minuv_meter),'--flags-save FLAG ','--g-ref-ant '+str(ref_ant_name),'--out-mode so','--g-update-type diag',' --g-diag-only 1','--out-overwrite 0','--out-backup 0',\
						'--g-save-to '+caltable_name,'--g-freq-int '+str(nchan),'--g-freq-int '+str(ntime)]
						if self.verbose:
							self.log_verbose.info(' '.join(cubical_args)+'\n')
						success=os.system(' '.join(cubical_args)+' > tmp')
						if success!=0:
							self.log_verbose.error('No good solution found. No caltable made.\n')
							os.system('rm -rf casa*log '+caltable_name+'*')
							os.chdir(self.cwd)
							os.system('rm -rf casa*log '+caltable_name+'*')
							if os.path.exists('cubical.cc-out'):
								os.system('rm -rf cubical.cc-out* tmp')
							return 7
						else:
							os.system('cp -r cubical.cc-out/cc.stats.pickle '+caltable_name+'.stat')
							self.log_verbose.info('flag_cubical_caltable(\''+caltable_name+'\',statfile=\''+caltable_name+'.stat\',overwrite=True,verbose='+str(self.verbose)+')\n')
							caltable_name=B.flag_cubical_caltable(caltable_name,statfile=caltable_name+'.stat',sigma=5,overwrite=True,verbose=self.verbose)
							self.log_verbose.info('normalize_cubical_caltable(\''+caltable_name+'\',soltype=\'B\',normtype=\'mean\',overwrite=True,verbose='+str(self.verbose)+')\n')
							caltable_name=B.normalize_cubical_caltable(caltable_name,soltype='B',normtype='mean',overwrite=True,verbose=self.verbose)
							cubical_args=['paircars_cubical','--data-ms '+ self.msname,'--out-casa-gaintables 0','--bbc-plot 0','--data-column '+datacolumn,\
							'--flags-save FLAG','--g-xfer-from '+caltable_name,'--out-mode ac','--out-overwrite 0','--out-backup 0','--sol-jones g']
							if self.verbose:
								self.log_verbose.info(' '.join(cubical_args)+'\n')
							os.system(' '.join(cubical_args)+' > tmp')
							if os.path.exists('cubical.cc-out'):
								os.system('rm -rf cubical.cc-out* tmp')
						flagdata(vis=self.msname,mode='clip',clipzeros=True,flagbackup=False)
					else:
						af=agentflagger()
						af.open(self.msname)
						versionlist=af.getflagversionlist()
						if len(versionlist)!=0:
							for version_name in versionlist:
								if 'cubical' in version_name:
									version_num=int(version_name.split(':')[0].split(' ')[0].split('_')[-1])+1
								else:
									version_num=1
						else:
							version_num=1
						now = datetime.now()
						dt_string = now.strftime("%Y-%m-%d %H:%M:%S")
						af.saveflagversion('cubical_'+str(version_num),'Flags autosave on '+dt_string)
						af.done()
						if calmode=='p':
							g_type='phase-diag'
						else:
							g_type='complex-diag'
						os.system('rm -rf cubical.cc-out* tmp')
						cubical_args=['paircars_cubical','--data-ms '+ self.msname,'--out-casa-gaintables 0','--model-list MODEL_DATA','--bbc-plot 0','--g-type '+g_type,'--data-column '+datacolumn,\
							'--sol-min-bl '+str(self.calib_minuv_meter),'--flags-save FLAG ','--g-ref-ant '+str(ref_ant_name),'--g-save-to '+caltable_name,'--out-mode so','--g-update-type diag',' --g-diag-only 1',\
							'--out-overwrite 0','--out-backup 0','--g-freq-int 1','--g-freq-int '+str(ntime)]
						if self.verbose:
							self.log_verbose.info(' '.join(cubical_args)+'\n')
						success=os.system(' '.join(cubical_args)+' > tmp')
						if success!=0:
							self.log_verbose.error('No good solution found. No caltable made.\n')
							if os.path.exists('cubical.cc-out'):
								os.system('rm -rf cubical.cc-out* tmp')
							os.system('rm -rf casa*log '+caltable_name+'*')
							os.chdir(self.cwd)
							os.system('rm -rf casa*log '+caltable_name+'*')
							return 7
						else:
							if calmode!='p':
								os.system('cp -r cubical.cc-out/cc.stats.pickle '+caltable_name+'.stat')
								if os.path.exists('cubical.cc-out'):
									os.system('rm -rf cubical.cc-out* tmp')
								self.log_verbose.info('flag_cubical_caltable(\''+caltable_name+'\',statfile=\''+caltable_name+'.stat\',overwrite=True,verbose='+str(self.verbose)+')\n')
								caltable_name=B.flag_cubical_caltable(caltable_name,statfile=caltable_name+'.stat',sigma=5,overwrite=True,verbose=self.verbose)
							if calmode!='p' and mwapb_cor==False:
								caltable_name=B.normalize_cubical_caltable(caltable_name,soltype='G',normtype='mean',overwrite=True,verbose=self.verbose)
							os.system('rm -rf cubical.cc-out* tmp')
							cubical_args=['paircars_cubical','--data-ms '+ self.msname,'--out-casa-gaintables 0','--bbc-plot 0','--data-column '+datacolumn,\
							'--flags-save FLAG','--g-xfer-from '+caltable_name,'--out-mode ac','--out-overwrite 0','--out-backup 0','--sol-jones g']
							if self.verbose:
								self.log_verbose.info(' '.join(cubical_args)+'\n')
							os.system(' '.join(cubical_args)+' > tmp')
							if os.path.exists('cubical.cc-out'):
								os.system('rm -rf cubical.cc-out* tmp')
						flagdata(vis=self.msname,mode='clip',clipzeros=True,flagbackup=False)
				else:
					if do_bandpass==True:
						if solint!='int' and solint!='inf':
							solint=str(solint)+'s'
						self.log_verbose.info('Doing bandpass.....\n')
						self.log_verbose.info('bpass_solver(\''+caltable_name+'\',spw=\'\',timerange=\'\',calmode=\''+calmode+'\',uvrange=\''+self.calib_uvrange_lambda+\
						'\',solnorm='+str(solnorm)+',refant=\''+str(ref_ant)+'\',minsnr='+str(minsnr)+',normtype+\''+str(normtype)+'\',gaintable='+str(calibrator_caltable)+\
						',solint=\''+solint+'\')\n')
						self.bpass_solver(caltable_name,spw='',timerange='',calmode=calmode,uvrange=self.calib_uvrange_lambda,solnorm=solnorm,refant=str(ref_ant),minsnr=minsnr,\
										normtype=normtype,gaintable=calibrator_caltable,solint=solint) # Perform bandpass calibration
						if do_fluxcal:
							modify_datacolumn=True
							atten=float(fits.getheader(self.metafits)['ATTEN_DB'])
							self.log_verbose.info('fluxcal_bptable(\''+caltable_name+'\',\''+self.metafits+'\',atten='+str(atten)+')\n')
							caltable_name=B.fluxcal_bptable(caltable_name,self.metafits,atten=atten)
					else:
						if solint!='int' and solint!='inf':
							solint=str(solint)+'s'
						try_count=0
						while try_count<5:
							self.log_verbose.info('Doing gaincal.....\n')
							self.log_verbose.info('gaincal(vis=\''+self.msname+'\',caltable=\''+caltable_name+'\',refant=\''+str(ref_ant)+'\',minsnr='+str(minsnr)
								+',calmode=\''+calmode+'\',solnorm='+str(solnorm)+',uvrange=\''+self.calib_uvrange_lambda+'\',gaintype=\'G\',solmode=\''+solmode+'\',rmsthresh=[10,7,5,3.5],gaintable='+
									str(calibrator_caltable)+',solint=\''+solint+'\',normtype=\''+normtype+'\')\n') 
							gaincal(vis=self.msname,caltable=caltable_name,refant=str(ref_ant),minsnr=minsnr,calmode=calmode,solnorm=solnorm,uvrange=self.calib_uvrange_lambda,\
								gaintype='G',solmode=solmode,rmsthresh=[10,7,5,3.5],gaintable=calibrator_caltable,solint=solint,normtype=normtype) # Performing gain calibration
							calc_flag_frac=B.calc_flag_fraction_caltable(caltable_name)
							if calc_flag_frac>=0.4:
								try_count+=1
								self.log_verbose.info('Large antennas flagged. Retrying to calibrate...\n')
								time.sleep(2)
								continue
							else:
								break
						if calc_flag_frac>=0.4:
							self.log_verbose.info('Performing gaincal without robust calibration due large number of flagged solutions.\n')
							self.log_verbose.info('gaincal(vis=\''+self.msname+'\',caltable=\''+caltable_name+'\',refant=\''+str(ref_ant)+'\',minsnr='+str(3)
								+',calmode=\''+calmode+'\',solnorm='+str(solnorm)+',normtype=\''+normtype+'\',uvrange=\''+self.calib_uvrange_lambda+'\',gaintype=\'G\',gaintable='+str(calibrator_caltable)+')\n') 
							gaincal(vis=self.msname,caltable=caltable_name,refant=str(ref_ant),minsnr=3,calmode=calmode,solnorm=solnorm,uvrange=self.calib_uvrange_lambda,\
								gaintype='G',gaintable=calibrator_caltable,normtype=normtype) # Performing gain calibration
					if os.path.isdir(caltable_name)==False:
						self.log_verbose.error('No good solution found. No caltable made.\n')
						os.system('rm -rf casa*log')
						os.chdir(self.cwd)
						os.system('rm -rf casa*log')
						return 7
					applycal_caltable=copy.deepcopy(calibrator_caltable)
					applycal_caltable.append(caltable_name)
					AM.remove_locks()
					self.log_verbose.info('Applying solutions from :'+str(applycal_caltable)+'\n')
					self.log_verbose.info('applycal(vis=\''+self.msname+'\',gaintable='+str(applycal_caltable)+',applymode=\''+applymode+'\',flagbackup=True,calwt=[True])\n')	
					applycal(vis=self.msname,gaintable=applycal_caltable,applymode=applymode,flagbackup=True,calwt=[True]) # Applying the solution
					if mwapb_cor:
						applycal(vis=self.msname,gaintable=applycal_caltable,applymode=applymode,flagbackup=True,calwt=[True])
				self.log_verbose.info('get_residual_flux_frac(\''+casa_imagename+'\',\''+casa_residualname+'\','+str(sigma)+',stokes_list='+str(stokes_list)+')\n')
				residual_flux,residual_flux_frac=self.get_residual_flux_frac(casa_imagename,casa_residualname,sigma,stokes_list=stokes_list)
				if use_ankflagger==True and do_flag==True and (residual_flux<10 or residual_flux_frac<0.1):
					try:
						self.log_verbose.info('do_uvsub_ankflag(\''+self.msname+'\',nthread=1,verbose='+str(False)+')\n')
						fg.do_uvsub_ankflag(self.msname,nthread=1,verbose=False)
					except Exception as e:
						self.log_verbose.info('Error in aNKflagger : '+str(e)+'\n')
						self.log_verbose.info('Error in running aNKflagger. Using rms threshold flagging.\n')
						self.log_verbose.info('do_uvsub_flagger(\''+self.msname+'\',mode=\'uvsub_flag\',rmsthresh=[10,7,5,3.5])\n')
						fg.do_uvsub_flagger(self.msname,mode='uvsub_flag',rmsthresh=[10,7,5,3.5])
				elif do_flag==True and (residual_flux<10 or residual_flux_frac<0.1):
					self.log_verbose.info('do_uvsub_flagger(\''+self.msname+'\',mode=\'uvsub_flag\',rmsthresh=[10,7,5,3.5])\n')
					fg.do_uvsub_flagger(self.msname,mode='uvsub_flag',rmsthresh=[10,7,5,3.5])
				if modify_datacolumn==True:
					tb=table()
					tb.open(self.msname)
					cor_data=tb.getcol('CORRECTED_DATA')
					data=tb.getcol('DATA')
					tb.close()
					tb.open(self.msname,nomodify=False)
					if use_cubical:
						tb.removecols('BITFLAG')
					tb.putcol('DATA',cor_data)
					tb.flush()
					tb.close()
				self.log_verbose.info('Success.\n')
				os.system('rm -rf casa*log')
				os.chdir(self.cwd)
				os.system('rm -rf casa*log')
				return 0,out_dict,negative_dyn_range
			
#########################################
# Finished Intensity Selfcal Class
#########################################


