import os
import time,psutil,paircarsdata,logging
from optparse import OptionParser
from casatasks import *
from casatools import msmetadata,table,measures,quanta,agentflagger,image,calibrater,ms,imager
from paircarstools.basic_func import *
from paircarstools.access_ms import *
from paircarsutils.paircars_casatasks.poltclean import *
from paircarstools.flagger import *
from astropy.io import fits
from paircarstools.fullpol_selfcal_LTS import *
from paircarstools.mwapb import mwapb_cor
os.system('rm -rf casa*log')

'''
Code is written by Devojyoti Kansabanik , 13 May, 2021
'''

global imagelog
global wsclean_start
formatter = logging.Formatter('%(asctime)s %(levelname)-8s %(message)s',datefmt='%Y-%m-%d %H:%M:%S')
imagelog = logging.getLogger('final_image_log')
imagelog.setLevel(logging.DEBUG)

def init_logger_console(logdir='',logger_name='',verbose=False):
	'''
	Initial logger
	Parameters
	----------
	working_dir : str
		Name of the working directory
	verbose : bool
		Verbose output or not
	'''
	if verbose==True:
		console=logging.StreamHandler(sys.stdout)
		console.setFormatter(formatter)
		imagelog.addHandler(console)
	if logdir!='' and logger_name!='':
		if os.path.exists(logdir+'/'+logger_name+'.log'):
			os.system('rm -rf '+logdir+'/'+logger_name+'.log')
		if os.path.isdir(logdir)==False:
			os.makedirs(logdir)
		filehandle=logging.FileHandler(logdir+'/'+logger_name+'.log')
		filehandle.setFormatter(formatter)
		imagelog.addHandler(filehandle)
	imagelog.propagate = False
	return

def calc_residual_flux(imagename,residual,nsigma,rms_box,stokes_list=['I']):
	'''
	Function to calculate residual flux.
	Parameters
	----------
	imagename : str 
		Name of the image
	nsigma : float 
		Sigma of present CLEAN cycle
	rms_box : str 
		Box to calculate rms (Only for CASA imaging)
	stokes_list : list 
		Stokes plane list
	Returns
	--------
	float
		Median residual flux over the stokes plane
	'''
	imagename_prefix=imagename.split('.image')[0]
	do_reduce_list=[]
	residual_frac_list=[]
	rms_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 '+imagename_prefix+'.reduce_sigma_*')
	for stokes in stokes_list:
		if stokes=='I' or stokes=='XX' or stokes=='YY':
			if os.path.isdir(imagename_prefix+'.reduce_sigma_I.image')==True:
				os.system('rm -rf '+imagename_prefix+'.reduce_sigma_I.image')
			if os.path.isdir(imagename_prefix+'.reduce_sigma_I.residual')==True:
				os.system('rm -rf '+imagename_prefix+'.reduce_sigma_I.residual')
			if stokes=='pseudoI':
				stokes='I'
			immath(imagename=imagename,mode='evalexpr',stokes=stokes,outfile=imagename_prefix+'.reduce_sigma_I.image')
			immath(imagename=residual,mode='evalexpr',stokes=stokes,outfile=imagename_prefix+'.reduce_sigma_I.residual')
			rms=imstat(imagename=imagename_prefix+'.reduce_sigma_I.image',box=rms_box,stokes=stokes)['rms'][0]
			rms_list.append(rms)
			ia.open(imagename_prefix+'.reduce_sigma_I.image')			
			ia.calcmask('\"'+imagename_prefix+'.reduce_sigma_I.image\">'+str(nsigma*rms),'mymask')
			ia.close()
			makemask(inpimage=imagename_prefix+'.reduce_sigma_I.image',inpmask=imagename_prefix+'.reduce_sigma_I.image:mymask',\
						output=imagename_prefix+'.reduce_sigma_I.residual:mymask',mode='copy')
			try:
				image_pix_sum=imstat(imagename=imagename_prefix+'.reduce_sigma_I.image')['sum'][0]
				residual_pix_sum=imstat(imagename=imagename_prefix+'.reduce_sigma_I.residual')['sum'][0]
			except:
				image_pix_sum=1
				residual_pix_sum=0
		else:
			if stokes=='pseudoI':
				stokes='I'
			immath(imagename=imagename,mode='evalexpr',stokes=stokes,expr='abs(IM0)',outfile=imagename_prefix+'.reduce_sigma_'+stokes+'.image')
			immath(imagename=residual,mode='evalexpr',stokes=stokes,outfile=imagename_prefix+'.reduce_sigma_'+stokes+'.residual')
			rms=imstat(imagename=imagename_prefix+'.reduce_sigma_'+stokes+'.image',box=rms_box,stokes=stokes)['rms'][0]
			rms_list.append(rms)
			ia.open(imagename_prefix+'.reduce_sigma_'+stokes+'.image')			
			ia.calcmask('\"'+imagename_prefix+'.reduce_sigma_'+stokes+'.image\">'+str(nsigma*rms),'mymask')
			ia.close()
			makemask(inpimage=imagename_prefix+'.reduce_sigma_'+stokes+'.image',inpmask=imagename_prefix+'.reduce_sigma_'+stokes+'.image:mymask',\
					output=imagename_prefix+'.reduce_sigma_'+stokes+'.residual:mymask',mode='copy')
			try:
				image_pix_sum=imstat(imagename=imagename_prefix+'.reduce_sigma_'+stokes+'.image')['sum'][0]
				residual_pix_sum=imstat(imagename=imagename_prefix+'.reduce_sigma_'+stokes+'.residual')['sum'][0]
			except:
				image_pix_sum=1
				residual_pix_sum=0
		if image_pix_sum==0:
			residual_frac_list.append(0)
		else:
			residual_frac_list.append(residual_pix_sum/image_pix_sum)
	os.system('rm -rf '+imagename_prefix+'.reduce_sigma_*')
	os.chdir(cwd)
	residual_frac_median=np.median(np.array(residual_frac_list))
	return residual_frac_median,rms_list


def modify_header(imagename,astrometry=False,imaging_cmd=''):
	'''
	Function to include P-AIRCARS speficific metadata in fits files
	Parameters
	----------
	imagename : str
		Name of the fits image
	astrometry : bool
		Precise astrometry or not
	imaging_cmd : str
		Imaging commands used
	'''
	header=fits.getheader(imagename)
	header['PIPELINE']='P-AIRCARS'
	header['WRITTER']='Devojyoti Kansabanik'
	if astrometry==True:
		astrometry=1
	else:
		astrometry=0
	header['ASTRO']=astrometry
	if imaging_cmd!='':
		header['HISTORY']=imaging_cmd
	fits.writeto(imagename,data=fits.getdata(imagename),header=header,overwrite=True)
	return 

def export_images(imagename,OBSID,cell,imsize,imaging_cmd='',savedir='',savemodel=False,saveresidual=False,cutoutbox=[],astrometry=False): 
	'''
	Function to save final images
	Parameters
	----------
	imagename : str 
		Name of the image to export
	OBSID : int 
		OBSID of the observation	
	cell : float
	 	Cell size of the image in arcsecond
	imsize : int 
		Number of pixels in image
	savedir : str 
		Directory to save the final images
	savemodel : str 
		Whether save model images or not
	saveresidual : str 
		Whether save residual images or not
	cutoutbox : list 
		Cutout box for final image [x_width,y_width] in degree
	astrometry : bool 
		Astrometry corrected or not
	Returns
	-------
	list
		List of image,model,residual
	'''
	cwd=os.getcwd()
	if os.path.isdir(savedir+'/All_final_images/'+str(OBSID))==False:	
		os.makedirs(savedir+'/All_final_images/'+str(OBSID))
	imagedir=savedir+'/All_final_images/'+str(OBSID)
	if savemodel==True:
		if os.path.isdir(savedir+'/All_final_models/'+str(OBSID))==False:
			os.makedirs(savedir+'/All_final_models/'+str(OBSID))
	modeldir=savedir+'/All_final_models/'+str(OBSID)
	if saveresidual==True:
		if os.path.isdir(savedir+'/All_final_residuals/'+str(OBSID))==False:
			os.makedirs(savedir+'/All_final_residuals/'+str(OBSID))
	resdir=savedir+'/All_final_residuals/'+str(OBSID)
	output=[] 
	if len(cutoutbox)!=0:
		x_pix=int((float(cutoutbox[0])*3600)/float(cell))
		y_pix=int((float(cutoutbox[1])*3600)/float(cell))
		x_cen=int(imsize/2)
		y_cen=x_cen
		ximsize=imhead(imagename=imagename+'.image',mode='list')['shape'][0]
		yimsize=imhead(imagename=imagename+'.image',mode='list')['shape'][1]
		if x_pix>ximsize:
			x_pix=ximsize-2
		if y_pix>yimsize:
			y_pix=yimsize-2
		box=str(int(x_cen-x_pix/2))+','+str(int(y_cen-y_pix/2))+','+str(int(x_cen+x_pix/2))+','+str(int(y_cen+y_pix/2))
		os.system('rm -rf '+imagename+'.cutout*')
		if os.path.isdir(imagedir):		
			os.system('rm -rf '+imagedir+'/'+os.path.basename(imagename)+'_*.fits '+imagedir+'/'+os.path.basename(imagename)+'_*.junk')
		if os.path.isdir(modeldir):
			os.system('rm -rf '+modeldir+'/'+os.path.basename(imagename)+'_*.fits '+modeldir+'/'+os.path.basename(imagename)+'_*.junk')
		if os.path.isdir(resdir):
			os.system('rm -rf '+resdir+'/'+os.path.basename(imagename)+'_*.fits '+resdir+'/'+os.path.basename(imagename)+'_*.junk')
		imsubimage(imagename=imagename+'.image',outfile=imagename+'.cutout.image',box=box)
		exportfits(imagename=imagename+'.cutout.image',fitsimage=imagedir+'/'+os.path.basename(imagename)+'_image.fits',history=False,overwrite=True)
		output.append(imagedir+'/'+os.path.basename(imagename)+'_image.fits')
		modify_header(imagedir+'/'+os.path.basename(imagename)+'_image.fits',astrometry=astrometry,imaging_cmd=imaging_cmd)
		if savemodel==True:
			imsubimage(imagename=imagename+'.model',outfile=imagename+'.cutout.model',box=box)
			exportfits(imagename=imagename+'.cutout.model',fitsimage=modeldir+'/'+os.path.basename(imagename)+'_model.fits',history=False,overwrite=True)
			output.append(modeldir+'/'+os.path.basename(imagename)+'_model.fits')
			modify_header(modeldir+'/'+os.path.basename(imagename)+'_model.fits',astrometry=astrometry,imaging_cmd=imaging_cmd)
		if saveresidual==True:
			imsubimage(imagename=imagename+'.residual',outfile=imagename+'.cutout.residual',box=box)
			exportfits(imagename=imagename+'.cutout.residual',fitsimage=resdir+'/'+os.path.basename(imagename)+'_res.fits',history=False,overwrite=True)
			output.append(resdir+'/'+os.path.basename(imagename)+'_res.fits')
			modify_header(resdir+'/'+os.path.basename(imagename)+'_res.fits',astrometry=astrometry,imaging_cmd=imaging_cmd)
		os.system('rm -rf '+imagename+'.cutout*')
	else:
		if os.path.isdir(imagedir):		
			os.system('rm -rf '+imagedir+'/'+os.path.basename(imagename)+'_*.fits '+imagedir+'/'+os.path.basename(imagename)+'_*.junk')
		if os.path.isdir(modeldir):
			os.system('rm -rf '+modeldir+'/'+os.path.basename(imagename)+'_*.fits '+modeldir+'/'+os.path.basename(imagename)+'_*.junk')
		if os.path.isdir(resdir):
			os.system('rm -rf '+resdir+'/'+os.path.basename(imagename)+'_*.fits '+resdir+'/'+os.path.basename(imagename)+'_*.junk')
		exportfits(imagename=imagename+'.image',fitsimage=imagedir+'/'+os.path.basename(imagename)+'_image.fits',history=False,overwrite=True)
		output.append(imagedir+'/'+os.path.basename(imagename)+'_image.fits')
		modify_header(imagedir+'/'+os.path.basename(imagename)+'_image.fits',astrometry=astrometry,imaging_cmd=imaging_cmd)
		if savemodel==True:
			exportfits(imagename=imagename+'.model',fitsimage=modeldir+'/'+os.path.basename(imagename)+'_model.fits',history=False,overwrite=True)
			output.append(modeldir+'/'+os.path.basename(imagename)+'_model.fits')
			modify_header(modeldir+'/'+os.path.basename(imagename)+'_model.fits',astrometry=astrometry,imaging_cmd=imaging_cmd)
		if saveresidual==True:
			exportfits(imagename=imagename+'.residual',fitsimage=resdir+'/'+os.path.basename(imagename)+'_res.fits',history=False,overwrite=True)
			output.append(resdir+'/'+os.path.basename(imagename)+'_res.fits')
			modify_header(resdir+'/'+os.path.basename(imagename)+'_res.fits',astrometry=astrometry,imaging_cmd=imaging_cmd)
	os.chdir(imagedir)
	os.system('rm -rf *.pb *.mask *.model *.image *.flux *.sumwt *.residual *psf*')
	os.chdir(cwd)
	return output

def get_stokes(stokes):
	'''
	Get stokes list from stokes string
	Parameters
	----------
	stokes : str
		Stokes string
	Returns
	-------
	list
		Stokes list
	'''
	if stokes=='I':
		return ['I']
	elif stokes=='pseudoI':
		return ['pseudoI']
	elif stokes=='Q':
		return ['Q']
	elif stokes=='U':
		return ['U']
	elif stokes=='V':
		return ['V']
	elif stokes=='IV':
		return ['I','V']
	elif stokes=='QU':
		return ['Q','U']
	elif stokes=='IQ':
		return ['I','Q']
	elif stokes=='UV':
		rerun ['U','V']
	elif stokes=='IQUV':
		return ['I','Q','U','V']
	elif stokes=='RR':
		return ['RR']
	elif stokes=='LL':
		return ['LL']
	elif stokes=='XX':
		return ['XX']
	elif stokes=='YY':
		return ['YY']
	elif stokes=='RRLL':
		return ['RR','LL']
	elif stokes=='XXYY':
		return ['XX','YY']
	else:
		return

def make_image(msname,metafits,basedir,workdir,sigma=5,stokes='I',savedir='',want_automask=False,maskfile='',quality_factor=1,threshold=[5],savemodel=False,\
	saveresidual=False,cutoutbox='',use_ankflagger=False,residual_frac=0.1,use_wsclean=True,cpus=3,absmem=5,clean_beam=[],num_pixel_in_psf=5,imaging_mode='',verbose=False,ref_ant=1,\
	logdir='',logger_name='',input_file='',wsclean_start='',timebin='0.5s',chanwidth=1):
	'''
	Function to make final images
	Parameters
	----------		
	msname : str 
		Name of the measurement set
	metafits : str 
		Name of the metafits file
	basedir : str
		Base directory of P-AIRCARS
	workdir : str 
		Name of the working directory
	sigma : float 
		Sigma value for thresholding
	stokes : str 
		Stokes planes to image
	savedir : str 
		Directory to save final directory
	want_automask : bool 
		Want to use auto masking (Only for CASA)
	maskfile : str 
		Name of any previous mask file or CASA mask string
	quality_factor : int 
		Quality factor for imaging
	threshold : list 
		rms list for all Stokes plane for thresholding
	savemodel : bool 
		Save model images or not
	saveresidual : bool 
		Save residual images or not
	cutoutbox : list 
		Cutout box of the final image [x_width,y_width] in degree (default : [], no cutout, save full image)
	use_ankflagger : bool 
		Whether use aNKflagger for flagging or not
	residual_frac : float 
		Residual flux fraction to stop CLEANing
	use_wsclean : bool
		Use wsclean or not
	cpus : int 
		Number of cpus to use in wsclean
	absmem : float 
		Absolute memory in GB for wsclean
	clean_beam : list
		Clean beam [maj,min,pa]
	num_pixel_in_psf : int
		Number of pixels in a single PSF
	imaging_mode : str
		Imaging mode ('database' or 'final')
	verbose : bool
		Verbose output
	ref_ant : int
		Reference antenna (Only used for differential calibration)
	logdir : str 
		Directory to save logs
	logger_name : str
		Log file name	
	input_file : str
		P-AIRCARS input file	
	wsclean_start : str
		WSClean path
	timebin : str
		Time average width
	freqwidth : int
		Number of channels to average
	Returns
	-------
	int
	 0 or 1 (Success ot Failure)
	'''
	try:
		sys.path.append(os.getcwd())
		import selfcal_inputs as inputs
	except:
		if input_file==None or os.path.exists(input_file)==False:
			print ('Please provide correct input file.\n')
			return 1
		os.system('cp -r '+inputfile+' '+os.getcwd()+'/selfcal_inputs.py')
		sys.path.append(os.getcwd())
		import selfcal_inputs as inputs
	imagelog.info('Making image of ms : '+str(msname)+'.....\n')
	if verbose==False:
		print ('Making image of ms : '+str(msname)+'.....\n')
	if msname[-1]=='/':
		msname=msname[:-1]
	org_msname=msname
	AM=AccessMS(org_msname)
	nchan=AM.get_num_channels()
	total_time=AM.calc_total_time()
	total_bandwidth=AM.calc_bandwidth()/1000
	freq=round(AM.calc_meanfreq()/10**6,2)
	timestamp=AM.get_timestamps()[0]
	tb=table()
	if os.path.exists(inputs.bookkeeping_file):
		try:
			imagelog.info('get_closest_DR(\''+inputs.bookkeeping_file+'\','+str(freq)+',\''+timestamp+'\')\n')
			rms_dr,neg_dr,mean_rms_dr,mean_neg_dr=get_closest_DR(inputs.bookkeeping_file,freq,timestamp)
			if rms_dr!=None:
				rms_dr=min(rms_dr,mean_rms_dr)
			else:
				rms_dr=10**20
			if neg_dr!=None:
				neg_dr=min(neg_dr,mean_neg_dr)	
			else:
				neg_dr=10**20
		except:
			rms_dr=neg_dr=10**20
		try:
			min_rms_dr,max_rms_dr,min_neg_dr,max_neg_dr=get_minmax_DR(inputs.bookkeeping_file)
		except:
			max_rms_dr=max_neg_dr=1000
			min_rms_dr=min_neg_dr=100
	else:
		rms_dr=neg_dr=10**20
		max_rms_dr=max_neg_dr=1000
		min_rms_dr=min_neg_dr=100
	if timebin!='' or chanwidth>1:
		timebin_value=float(timebin[:-1])
		ms_size=0
		num_trials=0	
		ms_file_made=True
		while ms_size==0 or num_trials<5:
			if chanwidth<nchan and timebin_value<total_time:
				if os.path.exists(org_msname.split('.ms')[0]+'_avg.ms'):
					os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
				if os.path.exists(org_msname.split('.ms')[0]+'_avg.ms.flagversions'):
					os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
				imagelog.info('split(vis=\''+org_msname+'\',outputvis=\''+org_msname.split('.ms')[0]+'_avg.ms\',datacolumn=\'data\',width='+str(chanwidth)+\
										',timebin=\''+timebin+'\')\n')
				split(vis=org_msname,outputvis=org_msname.split('.ms')[0]+'_avg.ms',datacolumn='data',width=chanwidth,timebin=timebin)
				ms_size=os.stat(org_msname.split('.ms')[0]+'_avg.ms/table.f21_TSM0').st_size
				ms_file_made=os.path.isdir(org_msname.split('.ms')[0]+'_avg.ms')
				if ms_size==0 or ms_file_made==False:
					num_trials+=1
					os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
					continue
				num_trials+=1
			else:
				break
		if os.path.isdir(org_msname.split('.ms')[0]+'_avg.ms')==False:
			os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
			os.system('cp -r '+org_msname+' '+org_msname.split('.ms')[0]+'_avg.ms')
	else:
		if os.path.exists(org_msname.split('.ms')[0]+'_avg.ms'):
			os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
		os.system('cp -r '+org_msname+' '+org_msname.split('.ms')[0]+'_avg.ms')
	msname=org_msname.split('.ms')[0]+'_avg.ms'
	wsclean_check_count=0
	if use_wsclean==True:
		tempdir=workdir+'/tempdir'
		if os.path.isdir(tempdir)==False:
			os.makedirs(tempdir)
	casalog=False
	cwd=os.getcwd()
	os.chdir(workdir)
	if cutoutbox!='':
		cutoutbox=cutoutbox.split(',')
	else:
		cutoutbox=[]
	cpus=3
	AM=AccessMS(msname)
	try_count=0
	done_imaging=False
	flag_count=0
	count=0
	if use_wsclean==True and wsclean_start=='':
		wsclean_datadir=os.path.abspath(os.path.dirname(paircarsdata.__file__))
		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:
				imagelog.info('WSClean is not installed. Using CASA for imaging.\n')
				use_wsclean=False
			else:
				wsclean_start=wsclean_path+'/wsclean'
				use_wsclean=True
		except:
			imagelog.warning('WSClean is not installed. Using CASA for imaging.\n')
			use_wsclean=False
		os.system('rm -rf wsclean_test')
	org_msname_touch_file=basedir+'/.Imaging_'+os.path.basename(org_msname).split('.ms')[0]
	while done_imaging==False:
		try:
			freqs=AM.calc_meanfreq()/10**6	
			OBSID=str(fits.getheader(metafits)['GPSTIME'])	
			file_str=workdir+'/'+os.path.basename(splited_ms_rename(org_msname,ref_time_chan=False,change_msname=False)).split('.ms')[0]
			already_images=glob.glob(savedir+'/All_final_images/'+str(OBSID)+'/'+os.path.basename(splited_ms_rename(msname,ref_time_chan=False,change_msname=False)).split('.ms')[0]+'*')
			if len(already_images)>0:
				for i in already_images:
					if '.junk' in i:
						imagelog.warning('Imaging tried but failed. Trying again.\n')
						done_imaging=True
						os.system('rm -rf '+i)
					else:
						imagelog.info('Imaging is already done. Image saved at : '+str(i)+'\n')
						os.system('touch '+org_msname_touch_file+'.success')
						done_imaging=True
						return 0
			os.system('rm -rf '+file_str+'.image '+file_str+'.model '+file_str+'.mask '+file_str+'.psf '+file_str+'.pb '+file_str+'.sumwt '+file_str+'.flux '+file_str+'.residual '+\
							file_str+'*.fits')
			if quality_factor==0:
				gain=0.15
				mgain=0.8
				multiscale_gain=0.12
			elif quality_factor==1:
				gain=0.1
				mgain=0.75
				multiscale_gain=0.1
			else:
				gain=0.08	
				mgain=0.7
				multiscale_gain=0.07
			stokes_list=get_stokes(stokes)
			mask=file_str+'.mask'
			PSC=PolSelfcal(org_msname,metafits,16*60,num_pixel_in_psf=num_pixel_in_psf,largest_scale=20,verbose=True,use_wsclean=use_wsclean,savelog=False,includeflag=False) 
			IB=ImageBasic(msname,includeflag=False)
			max_baseline=IB.calc_max_baseline_lambda()
			calc_calib_uvrange=IB.calc_calib_uvrange(20,includeflag=False)
			calib_uvrange_min=calc_calib_uvrange[1]
			calib_uvrange_max=calc_calib_uvrange[2]
			if inputs.calc_image_parameters==True:
				calc_calib_uvrange=IB.calc_calib_uvrange(20,includeflag=False)
				uvrange=calc_calib_uvrange[0]	
				imaging_minuv=calc_calib_uvrange[3]
				imaging_maxuv=calc_calib_uvrange[4]	
				cell=IB.calc_cellsize(num_pixel_in_psf) 
				if quality_factor==2:
					FWHM=False
					gridder='widefield'
					wprojplanes=-1
				else:
					FWHM=True
					gridder='standard'
					wprojplanes=-1
				imsize=IB.num_pixels(num_pixel_in_psf,FWHM=FWHM)
				scales=IB.choose_scales(num_pixel_in_psf,16*60)
				uvtaper=IB.calc_uvtaper()
				if inputs.weight=='' or (inputs.weight!='uniform' and inputs.weight!='natural' and inputs.weight!='briggs'):
					weight='briggs'
				else:
					weight=inputs.weight
				if weight=='briggs':
					if inputs.robust<-1.0 or inputs.robust>1.0:
						robust=1.0
					else:
						robust=inputs.robust
			else:
				weight=inputs.weight
				robust=inputs.robust
				uvrange_to_cal=inputs.uvrange_to_cal
				uvrange=uvrange_to_cal
				if uvrange_to_cal!='':
					if '~' in uvrange_to_cal:
						uvrange=uvrange_to_cal.split('~')
						imaging_minuv=float(uvrange[0])
						imaging_maxuv=float(uvrange[1])
					elif '>' in uvrange_to_cal:
						uvrange=uvrange_to_cal.split('>')
						imaging_minuv=float(uvrange[1])
						imaging_maxuv=AM.get_max_baseline()/AM.calc_meanwavelength()
					elif '<' in uvrange_to_cal:
						uvrange=uvrange_to_cal.split('<')
						imaging_minuv=0
						imaging_maxuv=float(uvrange[1])
				else:
					calc_calib_uvrange=IB.calc_calib_uvrange(20,includeflag=False)
					uvrange=calc_calib_uvrange[0]
					imaging_minuv=calc_calib_uvrange[3]
					imaging_maxuv=calc_calib_uvrange[4]	
				cell=inputs.cellsize
				imsize=inputs.imsize[0]
				scales=inputs.multiscale_scales		
				if inputs.uvtaper!='':
					uvtaper=inputs.uvtaper
				else:
					uvtaper=IB.calc_uvtaper()
		
			rms_box='50,50,'+str(imsize-50)+','+str(int(imsize/4))
			if sigma<5:
				sigma=5
			continue_clean=False
			while True:
				AM=AccessMS(msname)
				flagfrac=AM.calc_flagfrac()
				if flagfrac<=0.95:
					try:
						if use_wsclean==False:
							casa_imagename=file_str+'.image'
							casa_residualname=file_str+'.residual'
							casa_modelname=file_str+'.model'
							if len(threshold)!=len(stokes_list):	
								threshold=[threshold[0]]*len(stokes_list)
							if maskfile!='':
								mask_str=maskfile
							else:
								mask_str=''
							if count==0:
								threshold_list=[str(sigma*rms)+'Jy' for rms in threshold]
							else:
								threshold_list=[str(3*rms)+'Jy' for rms in threshold]
							if len(clean_beam)==3:
								restoring_beam=[str(clean_beam[0])+'arcsec',str(clean_beam[1])+'arcsec',str(clean_beam[2])+'deg']
							else:
								restoring_beam=[]
							check_thread_availability(4)
							if os.path.isdir(mask)==True:
								imaging_cmd='poltclean(vis=\''+msname+'\',selectdata=True,datacolumn="corrected",uvrange=\''+str(uvrange)+'\',imagename=\''+file_str+\
								'\',imsize=['+str(imsize)+'],cell=\''+str(cell)+'\',stokes=\''+str(stokes)+'\',gridder=\''+gridder+'\',pblimit=-1,deconvolver="multiscale",scales='+\
								str(scales)+',nterms=1,weighting=\''+weight+'\',robust='+str(robust)+',uvtaper=\''+str(uvtaper)+'\',casalogger='+\
								str(casalog)+',niter=100000000000,gain='+str(gain)+',threshold='+str(threshold_list)+',interactive=False,usemask="user",startmask='+\
								str(mask)+',savemodel=\'modelcolumn\',restoringbeam='+str(restoring_beam)+'\',wprojplanes='+str(wprojplanes)+')'
								imagelog.info(imaging_cmd+'\n')
								poltclean(vis=msname,selectdata=True,datacolumn="corrected",uvrange=uvrange,imagename=file_str,imsize=[imsize],cell=cell,stokes=stokes,gridder=\
								gridder,pblimit=-1,deconvolver="multiscale",scales=scales,nterms=1,weighting=weight,robust=robust,uvtaper=uvtaper,casalogger=casalog,\
								niter=100000000000,gain=gain,threshold=threshold_list,interactive=False,usemask="user",startmask=mask,savemodel='modelcolumn',\
								restoringbeam=restoring_beam,wprojplanes=wprojplanes)
							elif want_automask==True and maskfile=='':
								imaging_cmd='poltclean(vis=\''+msname+'\',selectdata=True,datacolumn="corrected",imagename=\''+file_str+'\',imsize=['+str(imsize)+'],cell=\''+str(cell)+\
								'\',stokes=\''+stokes+'\',gridder=\''+gridder+'\',pblimit=-1,deconvolver="multiscale",scales='+str(scales)+',nterms=1,weighting=\''+weight+'\',robust='+\
								str(robust)+'uvtaper='+str(uvtaper)+',niter=100000000000,gain='+str(gain)+',casalogger='+str(casalog)+',threshold='+str(threshold_list)+\
								',interactive=False,usemask=\'auto-multithresh\',negativethreshold=3.0,savemodel=\'modelcolumn\',restoringbeam='+str(restoring_beam)+\
								',wprojplanes='+str(wprojplanes)+')'
								imagelog.info(imaging_cmd+'\n')
								poltclean(vis=msname,selectdata=True,datacolumn="corrected",uvrange=uvrange,imagename=file_str,imsize=[imsize],cell=cell,stokes=stokes,\
								gridder=gridder,pblimit=-1,deconvolver="multiscale",scales=scales,nterms=1,weighting=weight,robust=robust,uvtaper=uvtaper,\
								niter=100000000000,gain=gain,casalogger=casalog,threshold=threshold_list,interactive=False,usemask='auto-multithresh',negativethreshold=3.0,\
								savemodel='modelcolumn',restoringbeam=restoring_beam,wprojplanes=wprojplanes)
							else:
								imaging_cmd='poltclean(vis=\''+msname+'\',selectdata=True,datacolumn="corrected",imagename=\''+str(file_str)+'\',imsize=['+str(imsize)+'],cell=\''+\
								str(cell)+'\',stokes=\''+str(stokes)+'\',gridder=\''+gridder+'\',pblimit=-1,deconvolver="multiscale",scales='+str(scales)+\
								',nterms=1,weighting=\''+weight+'\',robust='+str(robust)+',uvtaper=\''+str(uvtaper)+'\',casalogger='+str(casalog)+\
								',niter=100000000000,threshold='+str(threshold_list)+',interactive=False,usemask="user",mask=\''+str(mask_str)+\
								',savemodel=\'modelcolumn\',restoringbeam='+str(restoring_beam)+',wprojplanes='+str(wprojplanes)+')'
								imagelog.info(imaging_cmd+'\n')
								poltclean(vis=msname,selectdata=True,datacolumn="corrected",uvrange=uvrange,imagename=file_str,imsize=[imsize],cell=cell,stokes=stokes,\
								gridder=gridder,pblimit=-1,deconvolver="multiscale",scales=scales,nterms=1,weighting=weight,robust=robust,uvtaper=uvtaper,casalogger=casalog,\
								niter=100000000000,threshold=threshold_list,interactive=False,usemask="user",mask=mask_str,savemodel='modelcolumn',restoringbeam=restoring_beam,\
								wprojplanes=wprojplanes)
						else:
							if stokes=='I' or stokes=='pseudoI': 
								pol='i'
							elif stokes=='RR':
								pol='rr'
							elif stokes=='LL':
								pol='ll'
							elif stokes=='XX':
								pol='xx'
							elif stokes=='YY':
								pol='yy'
							elif stokes=='RRLL':
								pol='rrll'
							elif stokes=='XXYY':
								pol='xxyy'
							elif stokes=='IQUV':
								pol='iquv'
							scales=[str(i) for i in scales]
							if weight=='briggs':
								weight=weight+' '+str(robust)
							taper_tukey=10
							wsclean_args=['-scale '+str(cell)+'asec','-size '+str(imsize)+' '+str(imsize),'-no-dirty','-j '+str(cpus),'-abs-mem '+str(absmem),'-weight '+weight,'-taper-edge-tukey '+str(taper_tukey),\
							'-name '+str(file_str),'-pol '+str(pol),'-taper-edge '+str(imaging_maxuv),'-minuv-l '+str(imaging_minuv),'-niter 1000000','-mgain '+str(mgain),'-multiscale-gain '+str(multiscale_gain),\
							'-multiscale-scale-bias 0.6','-gain '+str(gain),'-multiscale','-multiscale-scales '+','.join(scales),'-temp-dir '+tempdir]
							if continue_clean:
								wsclean_args.append('-continue')
							if count!=0:
								if quality_factor==0:
									a=5
								elif quality_factor==1:
									sigma=4.5
								else:
									sigma=4
							wsclean_args.append('-auto-threshold '+str(sigma))
							wsclean_args.append('-auto-mask '+str(sigma+0.5))
							if gridder=='standard':
								wsclean_args.append('-nwlayers 1')
							wsclean_args.append('-quiet')
							if len(clean_beam)==3:
								wsclean_args.append('-beam-shape '+str(clean_beam[0])+' '+str(clean_beam[1])+' '+str(clean_beam[2]))
							imaging_cmd=wsclean_start+' '+' '.join(wsclean_args)+' '+msname
							imagelog.info(imaging_cmd+'\n')		
							check_thread_availability(4)
							wsclean_try=0
							wsclean_images=wsclean_models=wsclean_residuals=[]
							while wsclean_try<3:
								a=os.system(wsclean_start+' '+' '.join(wsclean_args)+' '+msname)
								wsclean_images=sorted(glob.glob(file_str+'*image.fits'))
								wsclean_models=sorted(glob.glob(file_str+'*model.fits'))
								wsclean_residuals=sorted(glob.glob(file_str+'*residual.fits'))
								if a==0:
									break
									time.sleep(2)
									for i in wsclean_images:
										if 'backup' in i:
											wsclean_images.remove(i)
									for i in wsclean_models:
										if 'backup' in i:
											wsclean_models.remove(i)
									for i in wsclean_residuals:
										if 'backup' in i:
											wsclean_residuals.remove(i)
								elif a!=0 or (len(wsclean_images)==0 or len(wsclean_models)==0 or len(wsclean_residuals)==0):
									wsclean_try+=1
									check_thread_availability(4)
									continue
							if len(wsclean_images)==0 or len(wsclean_models)==0 or len(wsclean_residuals)==0:
								check_thread_availability(4)
								if verbose==False:
									print ('Error occured in final imaging : Image is not made.\n')
								imagelog.error('Error occured in final imaging : Image is not made.\n')
								if try_count>=2:
									if verbose==False:
										print ('Maximum 2 tries failed. Quit.\n')
									imagelog.error('Maximum 2 tries failed. Quit.\n')
									done_imaging=True
									imagedir=savedir+'/All_final_images/'+str(OBSID)
									os.system('touch '+imagedir+'/'+os.path.basename(file_str)+'_image.junk')
									return 2
								else:
									antenna=AM.get_antenna_string()
									flagdata(vis=msname,mode='unflag',antenna=antenna)
									try_count+=1
									imagelog.info('Trying count : '+str(try_count)+'\n')
									continue
						if count==0:
							for i in wsclean_images:
								if '-I-' in i:
									if os.path.exists(i.split('image.fits')[0]+'backup_image.fits'):
										os.system('rm -rf '+i.split('image.fits')[0]+'backup_image.fits')
									imagelog.info('cp -r '+i+' '+i.split('image.fits')[0]+'backup_image.fits\n')
									os.system('cp -r '+i+' '+i.split('image.fits')[0]+'backup_image.fits')
									backup_stokes_I_image=i.split('image.fits')[0]+'backup_image.fits'
							for i in wsclean_models:
								if '-I-' in i:
									if os.path.exists(i.split('model.fits')[0]+'backup_model.fits'):
										os.system('rm -rf '+i.split('model.fits')[0]+'backup_model.fits')
									imagelog.info('cp -r '+i+' '+i.split('model.fits')[0]+'backup_model.fits\n')
									os.system('cp -r '+i+' '+i.split('model.fits')[0]+'backup_model.fits')
									backup_stokes_I_model=i.split('model.fits')[0]+'backup_model.fits'
							for i in wsclean_residuals:
								if '-I-' in i:
									if os.path.exists(i.split('residual.fits')[0]+'backup_residual.fits'):
										os.system('rm -rf '+i.split('residual.fits')[0]+'backup_residual.fits')
									imagelog.info('cp -r '+i+' '+i.split('residual.fits')[0]+'backup_residual.fits\n')
									os.system('cp -r '+i+' '+i.split('residual.fits')[0]+'backup_residual.fits')
									backup_stokes_I_residual=i.split('residual.fits')[0]+'backup_residual.fits'
							imagelog.info('Converting WSClean images : '+','.join([os.path.basename(i) for i in wsclean_images])+' to CASA image : '+file_str+'.image\n')
							casa_imagename=PSC.wsclean_to_casaimage(wsclean_images=wsclean_images,casaimage_prefix=file_str,imagetype='image',keep_wsclean_images=True)
							imagelog.info('Converting WSClean models : '+','.join([os.path.basename(i) for i in wsclean_models])+' to CASA model : '+file_str+'.model\n')
							casa_modelname=PSC.wsclean_to_casaimage(wsclean_images=wsclean_models,casaimage_prefix=file_str,imagetype='model',keep_wsclean_images=True)
							imagelog.info('Converting WSClean residuals : '+','.join([os.path.basename(i) for i in wsclean_residuals])+' to CASA residual : '+file_str+'.residual\n')
							casa_residualname=PSC.wsclean_to_casaimage(wsclean_images=wsclean_residuals,casaimage_prefix=file_str,imagetype='residual',keep_wsclean_images=True)
							# Checking whether differential calibration is required or not
							##############################################################
							DR1,DR2=get_dynamic_range(casa_imagename,stokes='I')
							check_thread_availability(4)
							'''imagelog.info('pb.mwapb_cor(\''+casa_imagename+'\',\''+file_str+'\',\''+metafits+'\',verbose='+str(verbose)+',nthreads=3,restore=False)\n')
							casa_imageanme=pb.mwapb_cor(casa_imagename,file_str,metafits,verbose=verbose,nthreads=3,restore=False)
							imagelog.info('pb.mwapb_cor(\''+casa_modelname+'\',\''+file_str+'\',\''+metafits+'\',verbose='+str(verbose)+',nthreads=3,restore=False)\n')
							casa_modelanme=pb.mwapb_cor(casa_modelname,file_str,metafits,verbose=verbose,nthreads=3,restore=False)'''
							imagelog.info('PSC.pol_model_threshold(\''+casa_imagename+'\',\''+casa_modelname+'\',10,1,stokes=[\'Q\',\'U\'],outside_mask=False)\n')
							PSC.pol_model_threshold(casa_imagename,casa_modelname,10,1,stokes=['Q','U'],outside_mask=False)
							imagelog.info('delmod(vis=\''+org_msname+'\',scr=True,otf=True)\n')
							delmod(vis=org_msname,scr=True,otf=True)
							imagelog.info('ft(vis=\''+org_msname+'\',model=\''+casa_modelname+'\',usescratch=True)\n')
							ft(vis=org_msname,model=casa_modelname,usescratch=True)
							imagelog.info('gaincal(vis=\''+org_msname+'\',caltable=\'round_0.gcal\',solnorm=True,normtype=\'median\',refant=\''+str(ref_ant)+'\',minsnr=3,solmode=\'R\',rmsthresh=[10,7,5,3.5])\n')
							gaincal(vis=org_msname,caltable='round_0.gcal',solnorm=True,normtype='median',refant=str(ref_ant),minsnr=3,solmode='R',rmsthresh=[10,7,5,3.5])
							imagelog.info('applycal(vis=\''+org_msname+'\',gaintable=[\'round_0.gcal\'])\n')
							applycal(vis=org_msname,gaintable=['round_0.gcal'])
							check_thread_availability(4)
							if flag_count==0:
								prev_flag_frac=calc_flag_fraction(org_msname)
								os.system('cp -r '+org_msname+' '+org_msname.split('.ms')[0]+'.aNKflag.tempms')
								if inputs.use_ankflagger:
									imagelog.info('do_uvsub_ankflag(\''+org_msname+'\',model=\''+casa_modelname+'\',verbose='+str(verbose)+',nthread=3,extendpols=False,chantime_minfrac=1.0,casaflag=\'\')\n')
									do_uvsub_ankflag(org_msname,model=casa_modelname,verbose=verbose,nthread=3,extendpols=False,chantime_minfrac=1.0,casaflag='')
								else:
									do_uvsub_flagger(org_msname,model=casa_modelname,rmsthresh=[10,7,5])
								pres_flag_frac=calc_flag_fraction(org_msname)
								imagelog.info('Previous flag fraction : '+str(prev_flag_frac)+'\n')
								imagelog.info('Present flag fraction : '+str(pres_flag_frac)+'\n')
								if round((pres_flag_frac-prev_flag_frac),2)>0.1 or round(pres_flag_frac,2)>0.3:
									imagelog.info('More data are getting flagged. Restoring the flags.\n')
									os.system('rm -rf '+org_msname+' '+org_msname+'.flagversions')
									os.system('mv '+org_msname.split('.ms')[0]+'.aNKflag.tempms '+org_msname)
								else:
									os.system('rm -rf '+org_msname.split('.ms')[0]+'.aNKflag.tempms')
								flag_count+=1

							os.system('rm -rf round_0.gcal')
							datacolumn_to_apply='CORRECTED_DATA'
							split_datacolumn='corrected'
							#count=1
							continue_clean=True
							if timebin!='' or chanwidth>1:
								ms_size=0
								num_trials=0	
								ms_file_made=True
								while ms_size==0 or num_trials<5:
									timebin_value=float(timebin[:-1])
									if chanwidth<nchan and timebin_value<total_time:
										if os.path.exists(org_msname.split('.ms')[0]+'_avg.ms'):
											os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
										if os.path.exists(org_msname.split('.ms')[0]+'_avg.ms.flagversions'):
											os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
										imagelog.info('split(vis=\''+org_msname+'\',outputvis=\''+org_msname.split('.ms')[0]+'_avg.ms\',datacolumn=\''+split_datacolumn+'\',width='+str(chanwidth)+\
											',timebin=\''+str(timebin)+'\')\n')
										split(vis=org_msname,outputvis=org_msname.split('.ms')[0]+'_avg.ms',datacolumn=split_datacolumn,width=chanwidth,timebin=timebin)
									else:
										if os.path.exists(org_msname.split('.ms')[0]+'_avg.ms'):
											os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
										if os.path.exists(org_msname.split('.ms')[0]+'_avg.ms.flagversions'):
											os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
										imagelog.info('split(vis=\''+org_msname+'\',outputvis=\''+org_msname.split('.ms')[0]+'_avg.ms\',datacolumn=\''+split_datacolumn+'\')\n')
										split(vis=org_msname,outputvis=org_msname.split('.ms')[0]+'_avg.ms',datacolumn=split_datacolumn)
									ms_size=os.stat(org_msname.split('.ms')[0]+'_avg.ms/table.f21_TSM0').st_size
									ms_file_made=os.path.isdir(org_msname.split('.ms')[0]+'_avg.ms')
									if num_trials>0:
										imagelog.info('Trying spliting ms for : '+str(num_trials)+' times.\n')
									if ms_file_made==False:
										num_trials+=1
										os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
										continue
									else:
										ms_size=os.stat(org_msname.split('.ms')[0]+'_avg.ms/table.f21_TSM0').st_size
										if ms_size==0:
											num_trials+=1
											os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
											continue
										else:
											break
										num_trials+=1
									num_trials+=1
								if os.path.isdir(org_msname.split('.ms')[0]+'_avg.ms')==False:
									os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
									os.system('cp -r '+org_msname+' '+org_msname.split('.ms')[0]+'_avg.ms')
									tb=table()
									tb.open(org_msname.split('.ms')[0]+'_avg.ms',nomodify=False)
									try:
										cor_data=tb.getcol('CORRECTED_DATA')
										tb.putcol('DATA',cor_data)
										tb.flush()
									except:
										pass
									tb.close()
							else:
								if os.path.exists(org_msname.split('.ms')[0]+'_avg.ms'):
									os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
								os.system('cp -r '+org_msname+' '+org_msname.split('.ms')[0]+'_avg.ms')
								tb=table()
								tb.open(org_msname.split('.ms')[0]+'_avg.ms',nomodify=False)
								try:
									cor_data=tb.getcol('CORRECTED_DATA')
									tb.putcol('DATA',cor_data)
									tb.flush()
								except:
									pass
								tb.close()
							msname=org_msname.split('.ms')[0]+'_avg.ms'
							if continue_clean:
								check_thread_availability(4)
								imagelog.info('delmod(vis=\''+org_msname.split('.ms')[0]+'_avg.ms\',scr=True,otf=True)\n')
								delmod(vis=org_msname.split('.ms')[0]+'_avg.ms',scr=True,otf=True)
								imagelog.info('ft(vis=\''+org_msname.split('.ms')[0]+'_avg.ms\',model=\''+casa_modelname+'\',usescratch=True)\n')
								ft(vis=org_msname.split('.ms')[0]+'_avg.ms',model=casa_modelname,usescratch=True)
							count+=1
							continue
						elif count==1:
							if backup_stokes_I_image in wsclean_images:
								wsclean_images.remove(backup_stokes_I_image) 
							if backup_stokes_I_model in wsclean_models:
								wsclean_models.remove(backup_stokes_I_model) 
							if backup_stokes_I_residual in wsclean_residuals:
								wsclean_residuals.remove(backup_stokes_I_residual) 
							for i in wsclean_images:
								if '-I-' in i and i!=backup_stokes_I_image:
									os.system('rm -rf '+i)
									imagelog.info('cp -r '+backup_stokes_I_image+' '+i+'\n')
									os.system('cp -r '+backup_stokes_I_image+' '+i)
									break
							for i in wsclean_models:
								if '-I-' in i and i!=backup_stokes_I_model:
									os.system('rm -rf '+i)
									imagelog.info('cp -r '+backup_stokes_I_model+' '+i+'\n')
									os.system('cp -r '+backup_stokes_I_model+' '+i)
									break
							for i in wsclean_residuals:
								if '-I-' in i and i!=backup_stokes_I_residual:
									os.system('rm -rf '+i)
									imagelog.info('cp -r '+backup_stokes_I_residual+' '+i+'\n')
									os.system('cp -r '+backup_stokes_I_residual+' '+i)
									break
							imagelog.info('Converting WSClean images : '+','.join([os.path.basename(i) for i in wsclean_images])+' to CASA image : '+file_str+'.image\n')
							casa_imagename=PSC.wsclean_to_casaimage(wsclean_images=wsclean_images,casaimage_prefix=file_str,imagetype='image',keep_wsclean_images=True)
							imagelog.info('Converting WSClean models : '+','.join([os.path.basename(i) for i in wsclean_models])+' to CASA model : '+file_str+'.model\n')
							casa_modelname=PSC.wsclean_to_casaimage(wsclean_images=wsclean_models,casaimage_prefix=file_str,imagetype='model',keep_wsclean_images=True)
							imagelog.info('Converting WSClean residuals : '+','.join([os.path.basename(i) for i in wsclean_residuals])+' to CASA residual : '+file_str+'.residual\n')
							casa_residualname=PSC.wsclean_to_casaimage(wsclean_images=wsclean_residuals,casaimage_prefix=file_str,imagetype='residual',keep_wsclean_images=True)
							median_res_frac,threshold=calc_residual_flux(casa_imagename,casa_residualname,sigma,rms_box,stokes_list=stokes_list)						
							#imagelog.info('PSC.remove_outside_model(\''+casa_imagename+'\',\''+casa_modelname+'\',sigma=5,box_width=5,overwrite=True)\n')
							#PSC.remove_outside_model(casa_imagename,casa_modelname,sigma=5,box_width=5,overwrite=True) 
							if stokes=='IQUV':
								qucor_image,qucor_model,q_leakage,u_leakage,v_leakage=PSC.correct_solar_quv_leakage(casa_imagename,casa_modelname,10,v_frac_threshold=0.02,overwrite=False)
								imagelog.info('Stokes I to leakage in Stokes Q: '+str(q_leakage*100)+', Stokes U: '+str(u_leakage*100)+', Stokes V: '+str(v_leakage*100)+'%\n')
								leakage_threshold=0.01
								if abs(q_leakage)>=leakage_threshold or abs(u_leakage)>=leakage_threshold or abs(v_leakage)>=0.05:
									cal_cause='Stokes leakage decreases.'
									do_cal=True
									continue_clean=False
								else:
									qucor_model=casa_modelname
									qucor_image=casa_imagename
									cal_cause='Differential calibration'
									do_cal=False
									continue_clean=True
								if do_cal:
									imagelog.info('delmod(vis=\''+org_msname+'\',scr=True,otf=True)\n')
									delmod(vis=org_msname,scr=True,otf=True)
									imagelog.info('Performing differential calibration because : '+cal_cause+'\n')
									k=0
									while True:
										if k>3:	
											if verbose==False:
												print ('Model import tried three times. Model import failed.\n')
											imagelog.info('Model import tried three times. Model import failed.\n')
											imagedir=savedir+'/All_final_images/'+str(OBSID)
											os.system('touch '+imagedir+'/'+os.path.basename(file_str)+'_image.junk')
											return 2
										imagelog.info('ft(vis=\''+org_msname+'\',model=\''+qucor_model+'\',usescratch=True)\n')
										ft(vis=org_msname,model=qucor_model,usescratch=True)
										AMmodel=AccessMS(org_msname)
										model_imported=AM.model_imported()
										if model_imported==True:
											break
										else:
											imagelog.info('Try to import model. Trial number : '+str(k)+'\n')
											time.sleep(1.0)
											imagelog.info('delmod(vis=\''+org_msname+'\',scr=True,otf=True)\n')
											delmod(vis=org_msname,scr=True,otf=True)
										k+=1
									AM=AccessMS(org_msname)
									ref_ant_name=AM.get_antenna_name(antenna_id=int(ref_ant))
									tb.open(org_msname)
									cor_data=tb.getcol('CORRECTED_DATA')
									tb.close()
									nodiff=True
									pre_diff_value=1
									caltable_name=file_str+'.cal'
									for i in range(5):
										datacolumn_to_cal='CORRECTED_DATA'
										cubical_args=['paircars_cubical','--data-ms '+org_msname,'--out-casa-gaintables 0','--model-list MODEL_DATA','--bbc-plot 0','--g-type complex-2x2',\
										'--sol-jones g','--data-column CORRECTED_DATA','--sol-min-bl '+str(calib_uvrange_min),'--sol-max-bl '+str(calib_uvrange_max),'--flags-save FLAG ',\
										'--g-ref-ant '+str(ref_ant_name),'--g-robust-npol 4','--g-save-to round_'+str(i)+'.cal','--out-mode sc','--out-overwrite 0','--out-backup 0',\
										'--flags-save-legacy 1','--out-dir '+workdir+'/cubical','--out-overwrite 1']
										check_thread_availability(4)
										if verbose:
											imagelog.info(' '.join(cubical_args)+'\n')
										c=0
										while True:
											success=os.system(' '.join(cubical_args))#+' > tmp')
											os.system('rm -rf '+workdir+'/cubical* tmp')
											if success==0:
												break
											else:
												imagelog.warning('Error in CubiCal. Trying again...\n')
												time.sleep(10)
												c+=1
												if c>5:
													imagelog.error('Error in Cubical. Exiting imaging...\n')
													done_imaging=True
													imagedir=savedir+'/All_final_images/'+str(OBSID)
													os.system('touch '+imagedir+'/'+os.path.basename(file_str)+'_image.junk')
													return 2
									if flag_count==0:
										check_thread_availability(4)
										prev_flag_frac=calc_flag_fraction(org_msname)
										os.system('cp -r '+org_msname+' '+org_msname.split('.ms')[0]+'.aNKflag.tempms')
										if inputs.use_ankflagger:
											imagelog.info('do_uvsub_ankflag(\''+org_msname+'\',model=\''+casa_modelname+'\',verbose='+str(verbose)+',nthread=3,extendpols=False,chantime_minfrac=1.0,casaflag=\'\')\n')
											do_uvsub_ankflag(org_msname,model=casa_modelname,verbose=verbose,nthread=3,extendpols=False,chantime_minfrac=1.0,casaflag='')
										else:
											do_uvsub_flagger(org_msname,model=casa_modelname,rmsthresh=[10,7,5])
										pres_flag_frac=calc_flag_fraction(org_msname)
										imagelog.info('Previous flag fraction : '+str(prev_flag_frac)+'\n')
										imagelog.info('Present flag fraction : '+str(pres_flag_frac)+'\n')
										if round((pres_flag_frac-prev_flag_frac),2)>0.1 or round(pres_flag_frac,2)>0.3:
											imagelog.info('More data are getting flagged. Restoring the flags.\n')
											os.system('rm -rf '+org_msname+' '+org_msname+'.flagversions')
											os.system('mv '+org_msname.split('.ms')[0]+'.aNKflag.tempms '+org_msname)
										else:
											os.system('rm -rf '+org_msname.split('.ms')[0]+'.aNKflag.tempms')
										flag_count+=1
							if do_cal==True or (median_res_frac>=residual_frac and sigma>=inputs.min_sigma):
								imagelog.info('delmod(vis=\''+org_msname+'\',scr=True,otf=True)\n')
								delmod(vis=org_msname,scr=True,otf=True)
								imagelog.info('ft(vis=\''+org_msname+'\',model=\''+casa_modelname+'\',usescratch=True)\n')
								ft(vis=org_msname,model=casa_modelname,usescratch=True)
								if sigma>inputs.min_sigma:
									sigma-=inputs.sigma_step
								if timebin!='' or chanwidth>1:
									ms_size=0
									num_trials=0	
									ms_file_made=True
									while ms_size==0 and num_trials<5:
										timebin_value=float(timebin[:-1])
										if chanwidth<nchan and timebin_value<total_time:
											if os.path.exists(org_msname.split('.ms')[0]+'_avg.ms'):
												os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
											if os.path.exists(org_msname.split('.ms')[0]+'_avg.ms.flagversions'):
												os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
											imagelog.info('split(vis=\''+org_msname+'\',outputvis=\''+org_msname.split('.ms')[0]+'_avg.ms\',datacolumn=\'corrected\',width='+str(chanwidth)+\
												',timebin=\''+timebin+'\')\n')
											split(vis=org_msname,outputvis=org_msname.split('.ms')[0]+'_avg.ms',datacolumn='corrected',width=chanwidth,timebin=timebin)
										else:
											if os.path.exists(org_msname.split('.ms')[0]+'_avg.ms'):
												os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
											if os.path.exists(org_msname.split('.ms')[0]+'_avg.ms.flagversions'):
												os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
											imagelog.info('split(vis=\''+org_msname+'\',outputvis=\''+org_msname.split('.ms')[0]+'_avg.ms\',datacolumn=\'corrected\')\n')
											split(vis=org_msname,outputvis=org_msname.split('.ms')[0]+'_avg.ms',datacolumn='corrected')
										if num_trials>0:
											imagelog.info('Trying spliting ms for : '+str(num_trials)+' times.\n')
										ms_file_made=os.path.isdir(org_msname.split('.ms')[0]+'_avg.ms')
										if ms_file_made==False:
											num_trials+=1
											os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
											continue
										else:
											ms_size=os.stat(org_msname.split('.ms')[0]+'_avg.ms/table.f21_TSM0').st_size
											if ms_size==0:
												num_trials+=1
												os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
												continue
											else:
												break
										num_trials+=1
									if os.path.isdir(org_msname.split('.ms')[0]+'_avg.ms')==False:
										os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
										os.system('cp -r '+org_msname+' '+org_msname.split('.ms')[0]+'_avg.ms')
										tb=table()
										tb.open(org_msname.split('.ms')[0]+'_avg.ms',nomodify=False)
										try:
											cor_data=tb.getcol('CORRECTED_DATA')
											tb.putcol('DATA',cor_data)
											tb.flush()
										except:
											pass
										tb.close()
								else:
									if os.path.exists(org_msname.split('.ms')[0]+'_avg.ms'):
										os.system('rm -rf '+org_msname.split('.ms')[0]+'_avg.ms '+org_msname.split('.ms')[0]+'_avg.ms.flagversions')
									os.system('cp -r '+org_msname+' '+org_msname.split('.ms')[0]+'_avg.ms')
									tb=table()
									tb.open(org_msname.split('.ms')[0]+'_avg.ms',nomodify=False)
									try:
										cor_data=tb.getcol('CORRECTED_DATA')
										tb.putcol('DATA',cor_data)
										tb.flush()
									except:
										pass
									tb.close()
								msname=org_msname.split('.ms')[0]+'_avg.ms'
							else:
								done_imaging=True
								break
							if continue_clean:
								check_thread_availability(4)
								imagelog.info('delmod(vis=\''+org_msname.split('.ms')[0]+'_avg.ms\',scr=True,otf=True)\n')
								delmod(vis=org_msname.split('.ms')[0]+'_avg.ms',scr=True,otf=True)
								imagelog.info('ft(vis=\''+org_msname.split('.ms')[0]+'_avg.ms\',model=\''+casa_modelname+'\',usescratch=True)\n')
								ft(vis=org_msname.split('.ms')[0]+'_avg.ms',model=casa_modelname,usescratch=True)
							count+=1
							continue
						if count>=2:
							done_imaging=True
							break
					except Exception as e:
						if verbose==False:
							print ('Error occured in final imaging : '+str(e)+'\n')
						imagelog.error('Error occured in final imaging : '+str(e)+'\n')
						if try_count>=2:
							if verbose==False:
								print ('Maximum 2 tries failed. Quit.\n')
							imagelog.error('Maximum 2 tries failed. Quit.\n')
							done_imaging=True
							imagedir=savedir+'/All_final_images/'+str(OBSID)
							os.system('touch '+imagedir+'/'+os.path.basename(file_str)+'_image.junk')
							return 2
						else:
							antenna=AM.get_antenna_string()
							flagdata(vis=msname,mode='unflag',antenna=antenna)
							try_count+=1
							imagelog.info('Trying count : '+str(try_count)+'\n')
							continue
			else:
				done_imaging=True
				imagedir=savedir+'/All_final_images/'+str(OBSID)
				os.system('touch '+imagedir+'/'+os.path.basename(file_str)+'_image.junk')
				return 1
		except Exception as e:
			if verbose==False:
				print ('Error occured in final imaging : '+str(e)+'\n')
			imagelog.error('Error occured in final imaging : '+str(e)+'\n')
			if try_count>=2:
				imagelog.info('Maximum 2 tries failed. Quit.\n')
				done_imaging=True
				imagedir=savedir+'/All_final_images/'+str(OBSID)
				os.system('touch '+imagedir+'/'+os.path.basename(file_str)+'_image.junk')
				return 2
			else:
				os.system('rm -rf '+msname+' '+msname+'.flagversions')
				try_count+=1
				imagelog.info('Trying count : '+str(try_count)+'\n')
				continue
			
	# Exporting images
	##################	
	wsclean_images=sorted(glob.glob(file_str+'*image.fits'))
	wsclean_models=sorted(glob.glob(file_str+'*model.fits'))
	wsclean_residuals=sorted(glob.glob(file_str+'*residual.fits'))
	for i in wsclean_images:
		if 'backup' in i:
			wsclean_images.remove(i)
	for i in wsclean_models:
		if 'backup' in i:
			wsclean_models.remove(i)
	for i in wsclean_residuals:
		if 'backup' in i:
			wsclean_residuals.remove(i)
	imagelog.info('Converting WSClean images : '+','.join([os.path.basename(i) for i in wsclean_images])+' to CASA image : '+file_str+'.image\n')
	casa_imagename=PSC.wsclean_to_casaimage(wsclean_images=wsclean_images,casaimage_prefix=file_str,imagetype='image',keep_wsclean_images=False)
	imagelog.info('Converting WSClean models : '+','.join([os.path.basename(i) for i in wsclean_models])+' to CASA model : '+file_str+'.model\n')
	casa_modelname=PSC.wsclean_to_casaimage(wsclean_images=wsclean_models,casaimage_prefix=file_str,imagetype='model',keep_wsclean_images=False)
	imagelog.info('Converting WSClean residuals : '+','.join([os.path.basename(i) for i in wsclean_residuals])+' to CASA residual : '+file_str+'.residual\n')
	casa_residualname=PSC.wsclean_to_casaimage(wsclean_images=wsclean_residuals,casaimage_prefix=file_str,imagetype='residual',keep_wsclean_images=False)	
	if os.path.isdir(casa_imagename):
		DR1,DR2=get_dynamic_range(casa_imagename,stokes='I')
		if (DR1<50 and DR1<(min_rms_dr*0.1)) or (DR2<10 and DR2<(min_neg_dr*0.1)) or np.isnan(DR1)==True or np.isnan(DR2)==True:
			imagelog.info('RMS based dynamic range : '+str(DR1)+', Negative based dynamic range : '+str(DR2)+'\n')
			imagelog.error('Final imaging dynamic range is too low.\n')
			done_imaging=True
			imagedir=savedir+'/All_final_images/'+str(OBSID)
			os.system('touch '+imagedir+'/'+os.path.basename(file_str)+'_image.junk')
			return 2
		paircars_image_bookkeep(casa_imagename,total_bandwidth,total_time,inputs.bookkeeping_file)
		imagelog.info('RMS based dynamic range : '+str(DR1)+', Negative based dynamic range : '+str(DR2)+'\n')
	imaging_cmd=imaging_cmd.split(' ')
	if '-continue' in imaging_cmd:
		imaging_cmd.remove('-continue')
	imaging_cmd='wsclean '+' '.join(imaging_cmd[1:])
	imagelog.info('export_images(\''+file_str+'\','+str(OBSID)+','+str(cell)+','+str(imsize)+',imaging_cmd=\''+str(imaging_cmd)+'\',savedir=\''+\
				str(savedir)+',savemodel='+str(savemodel)+',saveresidual='+str(saveresidual)+',cutoutbox=\''+str(cutoutbox)+'\',astrometry=False)\n')
	output=export_images(file_str,OBSID,cell,imsize,imaging_cmd=imaging_cmd,savedir=savedir,savemodel=savemodel,saveresidual=saveresidual,cutoutbox=cutoutbox,astrometry=False)
	if len(output)==0:
		imagelog.info('No final image is made.\n')
		msg=1
		os.system('touch '+org_msname_touch_file+'.error')
	elif len(output)==1:
		final_image=output[0]
		imagelog.info('Final image is written to : '+final_image+'\n')
		os.system('touch '+org_msname_touch_file+'.success')
		msg=0
	elif len(output)==2:
		final_image=output[0]
		final_model=output[1]
		imagelog.info('Final image is written to : '+final_image+'\n')
		imagelog.info('Final model is written to : '+final_model+'\n')
		os.system('touch '+org_msname_touch_file+'.success')
		msg=0
	else:
		final_image=output[0]
		final_model=output[1]
		final_res=output[2]
		imagelog.info('Final image is written to : '+final_image+'\n')
		imagelog.info('Final model is written to : '+final_model+'\n')
		imagelog.info('Final residual is written to : '+final_res+'\n')
		os.system('touch '+org_msname_touch_file+'.success')
		msg=0
	os.system('rm -rf *.image *.model *pb *psf* *.sumwt *.flux *.fits')
	os.system('cd ../')
	os.chdir(cwd)		
	return msg

def main():
	cwd=os.getcwd()
	start_time=time.time()
	usage= ' Perform final imaging\n'
	parser = OptionParser(usage=usage)
	parser.add_option('--msname',dest="msname",default=None,help="Name of measurement set of a single time anf frequency slice",metavar="Measurement Set")
	parser.add_option('--metafits',dest='metafits',default=None,help='Name of the metafits file',metavar='Metafits file')
	parser.add_option('--basedir',dest='basedir',default=None,help='Name of the base directory',metavar='Directory path')
	parser.add_option('--workdir',dest='workdir',default=None,help='Name of the working directory',metavar='Directory path')
	parser.add_option('--savedir',dest='savedir',default=None,help='Directory name to save final images',metavar="Directory path")
	parser.add_option('--savemodel',dest='savemodel',default=False,help='Want to save final models',metavar="Boolean")
	parser.add_option('--saveres',dest="saveresidual",default=False,help="Want to save residual images",metavar="Boolean")
	parser.add_option('--stokes',dest='stokes',default='pseudoI',help='Stokes planes to image',metavar="String")
	parser.add_option('--cutoutbox',dest='cutoutbox',default='',help='Cutout box \'X_width,Y_width\' in degree',metavar="Comma separated string")
	parser.add_option('--threshold',dest='threshold',default=0.1,help='RMS threshold for cleaning for each Stokes plane',metavar="Comma separated string")
	parser.add_option('--sigma',dest='sigma',default=10,help='Sigma value for thresholding',metavar="Float")
	parser.add_option('--want_automask',dest='want_automask',default=False,help='Want auto masking or not',metavar="Boolean")
	parser.add_option('--maskfile',dest='maskfile',default=None,help='Mask for imaging when auto masking is off',metavar="Maskfile or CASA mask string")
	parser.add_option('--quality_factor',dest='quality_factor',default=1,help='Quality factor of imaging',metavar="Integer")
	parser.add_option('--use_ankflag',dest='use_ankflag',default=False,help='Use aNKflag for flagging or not',metavar="Boolean")
	parser.add_option('--residual_frac',dest='resfrac',default=0.1,help='Residual flux fraction',metavar="Float")
	parser.add_option('--caltables',dest='caltables',default='',help='Caltables to apply',metavar="Comma separated string")
	parser.add_option('--wsclean',dest="use_wsclean",default=True,help="Use WSClean for imaging or not",metavar="Boolean")
	parser.add_option('--imaging_mode',dest='mode',default='final',help='Imaging for database or final imaging',metavar="String")
	parser.add_option('--inputfile',dest='inputfile',default=None,help='Path of the P-AIRCARS input file',metavar="File path")
	parser.add_option('--major_axis',dest='maj',default=None,help='Final image restoring beam major axis (FWHM) in arcsec',metavar="Float")
	parser.add_option('--minor_axis',dest='minor',default=None,help='Final image restoring beam minor axis (FWHM) in arcsec',metavar="Float")
	parser.add_option('--pa',dest='pa',default=None,help='Final image restoring beam position angle in degree',metavar="Float")
	parser.add_option('--num_pixel_in_psf',dest='num_pixel_in_psf',default=5,help='Number of pixels in one PSF',metavar="Integer")
	parser.add_option('--verbose',dest='verbose',default=False,help='Verbose output to save',metavar="Boolean")
	parser.add_option('--time_avg',dest='time_avg',default=None,help='Average in time after applying calibration (second)',metavar="Float")
	parser.add_option('--freq_avg',dest='freq_avg',default=None,help='Average in frequency (kHz)',metavar="Float")
	(options, args) = parser.parse_args()


	sys.path.append(cwd)
	try:
		import selfcal_inputs as inputs
	except Exception as e:
		if options.inputfile==None or os.path.isfile(str(options.inputfile))==False:
			print ('P-AIRCARS input file does not exist. Provide correct P-AIRCARS input file.\n')
			return 1
		inputfile=str(options.inputfile)
		if inputfile[-1]=='/':
			inputfile=inputfile[:-1]
		print ('WARNING : '+str(e)+'\n')
		os.system('cp -r '+inputfile+' '+cwd+'/selfcal_inputs.py')
		import selfcal_inputs as inputs
	input_file=cwd+'/selfcal_inputs.py'
	if os.path.exists(options.metafits)==False:
		print ('Metafits file is not given.Exiting....\n')
		count=0
		return 1
	OBSID=str(fits.getheader(options.metafits)['GPSTIME'])
	if inputs.keep_logger:
		logdir=str(options.basedir)+'/verbose_logs/'+str(OBSID)
		logger_name=str(OBSID)+'_'+str(os.path.basename(options.msname).split('.ms')[0])+'_imagelog'	
	else:
		logdir=logger_name=''
	init_logger_console(logdir=logdir,logger_name=logger_name,verbose=eval(str(options.verbose)))

	if eval(str(options.use_wsclean))==True:
		wsclean_datadir=os.path.abspath(os.path.dirname(paircarsdata.__file__))
		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:
				imagelog.info('WSClean is not installed. Using CASA for imaging.\n')
				use_wsclean=False
			else:
				wsclean_start=wsclean_path+'/wsclean'
				use_wsclean=True
		except:
			imagelog.warning('WSClean is not installed. Using CASA for imaging.\n')
			use_wsclean=False
		os.system('rm -rf wsclean_test')
	else:
		use_wsclean=eval(str(options.use_wsclean))
	if str(options.msname)[-1]=='/':
		msname=str(options.msname)[:-1]
	else:
		msname=str(options.msname)
	os.system('rm -rf '+options.basedir+'/.Finished_final_imaging_'+str(options.mode)+'_'+os.path.basename(options.msname)+'*')
	try:
		verbose=eval(str(options.verbose))

		if options.basedir==None:
			basedir=os.path.dirname(os.path.abspath(str(options.msname)))+'/final_imaging_basedir'
		else:
			basedir=str(options.basedir)

		if options.workdir==None:
			workdir=os.path.dirname(os.path.abspath(str(options.msname)))+'/final_imaging_workdir'
		else:
			workdir=str(options.workdir)
			
		if options.savedir==None:
			savedir=basedir
			imagelog.warning('Save directory is not given. Saving final images in base directory : '+basedir+'\n')
			if verbose==False:
				print('Save directory is not given. Saving final images in base directory : '+basedir+'\n')
		else:
			savedir=str(options.savedir)

		if os.path.isdir(basedir)==False:
			os.makedirs(basedir)
		if os.path.isdir(workdir)==False:
			os.makedirs(workdir)
		if os.path.isdir(savedir)==False:
			os.makedirs(savedir)

		if os.path.isdir(str(options.msname))==False or options.msname==None:
			if verbose==False:
				print ('Measurement set is not present.\n')
			imagelog.error('Measurement set is not present.\n')
			os.system('touch '+cwd+'/.Finished_final_imaging_'+str(options.mode)+'_'+os.path.basename(str(msname))+'_noms')
			return 1
		elif os.path.isfile(str(options.metafits))==False or options.metafits==None:
			if verbose==False:
				print ('Metafits file is not present.\n')
			imagelog.error('Metafits file is not present.\n')
			os.system('touch '+cwd+'/.Finished_final_imaging_'+str(options.mode)+'_'+os.path.basename(str(msname))+'_nometa')
			return 1
		else:
			AM=AccessMS(options.msname)
			flagfrac=AM.calc_flagfrac()
			if flagfrac>=0.9:
				imagelog.error('More than 90% data are flagged.\n')
				done_imaging=True
				imagedir=savedir+'/All_final_images/'+str(OBSID)
				file_str=options.workdir+'/'+os.path.basename(splited_ms_rename(options.msname,ref_time_chan=False,change_msname=False)).split('.ms')[0]
				os.system('touch '+imagedir+'/'+os.path.basename(file_str)+'_image.junk')
				if os.path.exists(savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)):
					imagelog.info('rm -rf '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)+' '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)+'.flagversions')
					os.system('rm -rf '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)+' '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)+'.flagversions')
				if os.path.exists(options.msname):
					os.system('rm -rf '+options.msname+' '+options.msname+'.flagversions')
				return 2
			if os.path.isdir(options.workdir+'/imaging_data')==False:
				os.makedirs(options.workdir+'/imaging_data')
			if str(options.caltables)=='':
				caltables=[]
			else:
				caltables=str(options.caltables).split(',')
				if str(options.caltables)=='':
					caltables=[]
				else:
					caltables=str(options.caltables).split(',')
				if verbose==False:
					print ('Applying gain calibration.....\n')
				imagelog.info('Applying gain calibration.....\n')
				imagelog.info('clearcal(vis=\''+options.msname+'\')\n')
				clearcal(vis=options.msname)
				# Applying gaincal and bandpass
				###############################
				AM=AccessMS(options.msname)
				meanfreq=(AM.calc_meanfreq()/10**6)
				coarse_chan=freq_to_MWA_coarse(meanfreq)
				gaintable=[]
				caltables_for_ms=[]
				interp=[]
				polcaltables=[]
				if options.msname[-1]=='/':
					os.system('cp -r '+options.msname+' '+options.msname+'.backup')
				for i in range(len(caltables)):
					if '.pcal' not in caltables[i]:
						gaintable.append(caltables[i])
						if '.gcal' in caltables[i]:
							interp.append('linear')
						else:
							interp.append('linear')
					elif '.pcal' in caltables[i]:
						polcaltables.append(caltables[i])
				if verbose:
					imagelog.info('applycal(vis=\''+options.msname+'\',gaintable='+str(gaintable)+',applymode=\'calflag\',flagbackup=False,interp='+str(interp)+')\n')
				applycal(vis=options.msname,gaintable=gaintable,applymode='calflag',flagbackup=False,interp=interp)
				if len(polcaltables)>0:
					# Applying cross-hand phase correction
					######################################
					PSC=PolSelfcal(options.msname,options.metafits,16*60,num_pixel_in_psf=int(options.num_pixel_in_psf),largest_scale=20,verbose=verbose,use_wsclean=use_wsclean,savelog=False,includeflag=False) 
					mwa_config=get_MWA_phase(options.metafits) # TODO : Include from cross phase cal solutions
					if mwa_config=='MWAPhaseI':
						crossphase=15
					elif mwa_config=='MWAPhaseIILB' or mwa_config=='MWAPhaseIICOMPACT':
						crossphase=135
					imagelog.info('Applying cross hand phase solution. Cross hand phase : '+str(crossphase)+' deg.\n')
					PSC.apply_cross_hand_phase(cross_phase=crossphase,caltable='',polbasis='Linear',modify_datacolumn=False,datacolumn='CORRECTED_DATA')
					del PSC
					# Applying polarization calibration
					###################################
					update_type_list=['full']
					os.system('rm -rf cubical.cc-out* tmp')
					for i in range(len(polcaltables)):
						datacolumn='CORRECTED_DATA'
						if '.pcal' in polcaltables[i]:
							for j in range(len(update_type_list)):
								cubical_args=['paircars_cubical','--data-ms '+options.msname,'--out-casa-gaintables 0','--bbc-plot 0',\
								'--data-column '+datacolumn,'--flags-apply FLAG','--g-update-type '+str(update_type_list[j]),'--g-xfer-from '+polcaltables[i],\
								'--out-mode ac','--out-overwrite 0','--out-backup 0','--flags-save-legacy 1','--out-dir '+options.workdir+'/cubical','--out-overwrite 1']
								if verbose:
									imagelog.info(' '.join(cubical_args))
								c=0
								while True:
									a=os.system(' '.join(cubical_args))#+' > tmp')
									os.system('rm -rf '+options.workdir+'/cubical* tmp')
									if a==0:
										break
									else:
										time.sleep(5)
										imagelog.warning('Error in applying CubiCal solution. Trying again...\n')	
										c+=1
										if c>=5:
											imagelog.info('CubiCal solution could not be applied.\n')
											done_imaging=True
											imagedir=savedir+'/All_final_images/'+str(OBSID)
											if os.path.exists(savedir+'/failed_ms/'+str(OBSID))==False:
												os.makedirs(savedir+'/failed_ms/'+str(OBSID))
											if os.path.exists(savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname))==False:
												imagelog.info('cp -r '+options.msname+' '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname))
												os.system('cp -r '+options.msname+' '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname))
											file_str=options.workdir+'/'+os.path.basename(splited_ms_rename(options.msname,ref_time_chan=False,change_msname=False)).split('.ms')[0]
											os.system('touch '+imagedir+'/'+os.path.basename(file_str)+'_image.junk')
											return 2				
				os.system('rm -rf '+options.workdir+'/imaging_data')
			AM=AccessMS(options.msname)
			flagfrac=AM.calc_flagfrac()
			if flagfrac>=0.9:
				imagelog.error('More than 90% data are flagged.\n')
				done_imaging=True
				imagedir=savedir+'/All_final_images/'+str(OBSID)
				file_str=options.workdir+'/'+os.path.basename(splited_ms_rename(options.msname,ref_time_chan=False,change_msname=False)).split('.ms')[0]
				os.system('touch '+imagedir+'/'+os.path.basename(file_str)+'_image.junk')
				if os.path.exists(savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)):
					imagelog.info('rm -rf '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)+' '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)+'.flagversions')
					os.system('rm -rf '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)+' '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)+'.flagversions')
				if os.path.exists(options.msname):
					os.system('rm -rf '+options.msname+' '+options.msname+'.flagversions')
				return 2
			freqres=AM.calc_freqres()
			timeres=AM.calc_timeres()
			if options.time_avg!=None:
				timebin=str(options.time_avg)+'s'
			else:
				timebin=''
			if options.freq_avg!=None:
				width=int(float(options.freq_avg)/float(freqres))
			else:
				width=1
			if os.path.exists(options.msname.split('.ms')[0]+'_cor.ms'):
				os.system('rm -rf '+options.msname.split('.ms')[0]+'_cor.ms '+options.msname.split('.ms')[0]+'_cor.ms.flagversions')
			if os.path.exists(options.msname.split('.ms')[0]+'_cor.ms.flagversions'):
				os.system('rm -rf '+options.msname.split('.ms')[0]+'_cor.ms.flagversions')
			if freqres<160:
				chanavg_bin=int(160/float(freqres))
			else:
				chanavg_bin=1
			if timeres<2:
				timeavg_bin='2s'
			else:
				timeavg_bin=''
			imagelog.info('split(vis=\''+options.msname+'\',outputvis=\''+options.msname.split('.ms')[0]+'_cor.ms\',datacolumn=\'corrected\',width='+str(chanavg_bin)+',timebin=\''+\
						timeavg_bin+'\')\n')
			split(vis=options.msname,outputvis=options.msname.split('.ms')[0]+'_cor.ms',datacolumn='corrected',width=chanavg_bin,timebin=timeavg_bin)
			msname=options.msname.split('.ms')[0]+'_cor.ms'
			if verbose==False:
				print ('#############################\nStart imaging for ms : '+str(msname)+'\n')
			imagelog.info('#############################\nStart imaging for ms : '+str(msname)+'\n')

			error_files=glob.glob(basedir+'/.Finished_final_imaging_'+str(options.mode)+'_'+os.path.basename(msname)+'*error*')
			for i in error_files:
				os.system('rm -rf '+i)
			touch_file=basedir+'/.Finished_final_imaging_'+str(options.mode)+'_'+os.path.basename(msname)+'_success'
			touch_file_moreflag=basedir+'/.Finished_final_imaging_'+str(options.mode)+'_'+os.path.basename(msname)+'_moreflag'
			if os.path.exists(touch_file) or os.path.exists(touch_file_moreflag):
				if verbose==False:
					print('Imaging is already done or attempted.\n#############################\n')
				imagelog.info('Imaging is already done or attempted.\n#############################\n')
				org_msname_touch_file=basedir+'/.Imaging_'+os.path.basename(options.msname).split('.ms')[0]
				workdir_touch_files=glob.glob(basedir+'/.Imaging_'+os.path.basename(options.msname).split('.ms')[0]+'*')
				if len(workdir_touch_files)==0:
					basedir_touch_files=glob.glob(basedir+'/.Finished_final_imaging_'+str(options.mode)+'_'+os.path.basename(msname)+'*')
					if len(basedir_touch_files)>0:
						basedir_touch_file_msg=basedir_touch_files[0].split('_')[-1]
					else:
						basedir_touch_file_msg='nofound'
					if basedir_touch_file_msg=='success':
						os.system('touch '+org_msname_touch_file+'.success')
					else:
						os.system('touch '+org_msname_touch_file+'.error')
				if verbose:
					if os.path.isdir(savedir+'/verbose_backup_ms')==False:
						os.makedirs(savedir+'/verbose_backup_ms')
					if os.path.exists(str(msname)):
						os.system('cp -r '+str(msname)+' '+savedir+'/verbose_backup_ms/')
						os.system('rm -rf '+str(msname))
					if os.path.exists(str(workdir)):
						os.system('cp -r '+str(workdir)+' '+savedir+'/verbose_backup_ms/')
						os.system('rm -rf '+str(workdir))
					if os.path.exists(str(options.msname)):
						os.system('cp -r '+str(options.msname)+' '+savedir+'/verbose_backup_ms/')
						os.system('rm -rf '+str(options.msname))
					if os.path.exists(str(options.msname+'.flagversions')):
						os.system('cp -r '+str(options.msname+'.flagversions')+' '+savedir+'/verbose_backup_ms/')
						os.system('rm -rf '+str(options.msname+'.flagversions'))
					if os.path.exists(str(msname.split('.ms')[0]+'_avg.ms')):
						os.system('cp -r '+str(msname.split('.ms')[0]+'_avg.ms')+' '+savedir+'/verbose_backup_ms/')
						os.system('rm -rf '+str(msname.split('.ms')[0]+'_avg.ms'))
				os.system('rm -rf '+str(msname)+' '+str(msname)+'* '+workdir+' '+options.msname+' '+options.msname+'.flagversions '+msname.split('.ms')[0]+'_avg.ms '+msname.split('.ms')[0]+'_avg.ms.flagversions')
				imagelog.info('Total run time : '+str(time.time()-start_time)+' s\n#############################\n')
				return 0

			rmsthresh=str(options.threshold).split(',')
			threshold=[float(i) for i in rmsthresh]
			IB=ImageBasic(str(options.msname),includeflag=False)
			AM=AccessMS(str(options.msname))
			cent_freq=AM.calc_meanfreq()/10**6
			coarse_chan_freq=freq_to_MWA_coarse(cent_freq)*1.28

			if options.maskfile==None:
				maskfile=''
			else:
				maskfile=str(options.maskfile)

			if options.maj!=None and options.minor!=None and options.pa!=None:
				clean_beam=[float(options.maj),float(options.minor),float(options.pa)]
			else:
				clean_beam=[]			
			output=make_image(msname,str(options.metafits),str(basedir),str(workdir),sigma=float(options.sigma),stokes=str(options.stokes),savedir=str(savedir),threshold=threshold,\
			want_automask=eval(str(options.want_automask)),maskfile=maskfile,quality_factor=int(options.quality_factor),savemodel=eval(str(options.savemodel)),\
			saveresidual=eval(str(options.saveresidual)),cutoutbox=str(options.cutoutbox),use_ankflagger=eval(str(options.use_ankflag)),residual_frac=float(options.resfrac),\
			use_wsclean=use_wsclean,clean_beam=clean_beam,num_pixel_in_psf=int(options.num_pixel_in_psf),imaging_mode=str(options.mode),\
			verbose=verbose,logdir=logdir,logger_name=logger_name,input_file=input_file,wsclean_start=wsclean_start,timebin=timebin,chanwidth=width)
			if output==0:
				if verbose==False:
					print('Imaging finished.\n#############################\n')
				imagelog.info('Imaging finished.\n#############################\n')
				os.system('touch '+touch_file)
			elif output==1:
				if verbose==False:
					print('More than 95% data are flagged. Image is not made.\n############################\n')
				imagelog.error('More than 95% data are flagged. Image is not made.\n############################\n')
				touch_file=basedir+'/.Finished_final_imaging_'+str(options.mode)+'_'+os.path.basename(msname)+'_moreflag'
				os.system('touch '+touch_file)
			elif output==2:
				if verbose==False:
					print('Error occured during final imaging.\n############################\n')
				imagelog.error('Error occured during final imaging.\n############################\n')
				touch_file=basedir+'/.Finished_final_imaging_'+str(options.mode)+'_'+os.path.basename(msname)+'_error'
				os.system('touch '+touch_file)
			elif output==3:
				if verbose==False:
					print('Image dynamic range is too low.\n############################\n')
				imagelog.error('Image dynamic range is too low.\n############################\n')
				touch_file=basedir+'/.Finished_final_imaging_'+str(options.mode)+'_'+os.path.basename(msname)+'_lowDR'
				os.system('touch '+touch_file)
			if output==0:
				if os.path.exists(savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)):
					imagelog.info('rm -rf '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)+' '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)+'.flagversions')
					os.system('rm -rf '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)+' '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname)+'.flagversions')
				else:
					os.system('rm -rf '+options.msname+' '+options.msname+'.flagversions')
			else:
				if os.path.exists(savedir+'/failed_ms/'+str(OBSID))==False:
					os.makedirs(savedir+'/failed_ms/'+str(OBSID))
				if os.path.exists(savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname))==False:
					imagelog.info('cp -r '+options.msname+' '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname))
					os.system('cp -r '+options.msname+' '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname))
				if options.msname!=savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname):
					os.system('rm -rf '+options.msname+' '+options.msname+'.flagversions')
			if verbose:
				if os.path.isdir(savedir+'/verbose_backup_ms')==False:
					os.makedirs(savedir+'/verbose_backup_ms')
				if os.path.exists(str(msname)):
					os.system('cp -r '+str(msname)+' '+savedir+'/verbose_backup_ms/')
					os.system('rm -rf '+str(msname))
				if os.path.exists(str(workdir+' '+options.msname)):
					os.system('cp -r '+str(workdir+' '+options.msname)+' '+savedir+'/verbose_backup_ms/')
					os.system('rm -rf '+str(workdir+' '+options.msname))
				if os.path.exists(str(options.msname+'.flagversions')):
					os.system('cp -r '+str(options.msname+'.flagversions')+' '+savedir+'/verbose_backup_ms/')
					os.system('rm -rf '+str(options.msname+'.flagversions'))
			os.system('rm -rf '+str(msname)+' '+str(msname)+'* '+workdir+' '+msname.split('.ms')[0]+'_avg.ms '+msname.split('.ms')[0]+'_avg.ms.flagversions')
			block_dir=os.path.dirname(options.msname)
			num_ms_block_dir=glob.glob(block_dir+'/*.ms')
			if len(num_ms_block_dir)==0:
				imagelog.info('All imaging is done for the block. Deleting block directory : '+str(block_dir)+'\n')
				os.system('rm -rf '+block_dir)
			else:
				imagelog.info(str(len(num_ms_block_dir))+' measurement set still exist in block directory.\n')
			imagelog.info('Total run time : '+str(time.time()-start_time)+' s\n#############################\n')
			return 0
	except Exception as e:
		if verbose==False:
			print('Error occured during final imaging. Error : '+str(e)+'\n############################\n')
		imagelog.error('Error occured during final imaging. Error : '+str(e)+'\n############################\n')
		touch_file=options.basedir+'/.Finished_final_imaging_'+str(options.mode)+'_'+os.path.basename(msname)+'_error'
		print (touch_file)
		os.system('touch '+touch_file)
		if os.path.exists(savedir+'/failed_ms/'+str(OBSID))==False:
			os.makedirs(savedir+'/failed_ms/'+str(OBSID))
		if os.path.exists(savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname))==False:
			imagelog.info('cp -r '+options.msname+' '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname))
			os.system('cp -r '+options.msname+' '+savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname))
		if options.msname!=savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname) and msname!=savedir+'/failed_ms/'+str(OBSID)+'/'+os.path.basename(options.msname):
			os.system('rm -rf '+options.msname+' '+options.msname+'.flagversions')
		print (msname,workdir)
		#if os.path.exists(msname):	
		#	os.system('rm -rf '+str(msname)+' '+str(msname)+'.flagversion')
		os.system('rm -rf '+workdir)
		block_dir=os.path.dirname(options.msname)
		num_ms_block_dir=glob.glob(block_dir+'/*.ms')
		if len(num_ms_block_dir)==0:
			imagelog.info('All imaging is done for the block. Deleting block directory : '+str(block_dir)+'\n')
			os.system('rm -rf '+block_dir)
		else:
			imagelog.info(str(len(num_ms_block_dir))+' measurement set still exist in block directory.\n')
		imagelog.info('Total run time : '+str(time.time()-start_time)+' s\n#############################\n')
		return 1

if __name__=='__main__':
	result=main()
	os._exit(result)





