import os,glob,subprocess,matplotlib,pandas as pd
import numpy as np,julian,smtplib,imaplib,datetime as dtt,psutil,json,urllib.request,copy,time,requests,matplotlib.pyplot as plt,paircarsdata,warnings
from casatools import msmetadata,table,measures,quanta,agentflagger,image,calibrater,ms,imager
from casatasks import *
from astropy import wcs
from . import access_ms as am
from . import mwapb as pb
from astropy.io import fits
from astropy.coordinates import EarthLocation,SkyCoord,AltAz
from astropy.time import Time
from astropy import units as u
from datetime import datetime
from email.mime.multipart import MIMEMultipart
from email.mime.text import MIMEText
from email.mime.application import MIMEApplication
from numpy.linalg import det
from pathlib import Path
from cubical import param_db
from collections import OrderedDict
matplotlib.use('Agg')

'''
Code is written by Devojyoti Kansabanik, 05 Jan, 2021
'''
############ Basic functions ###############
datadir = os.path.dirname(paircarsdata.__file__)

class ImageBasic:
	'''
	Generic class to calculate different imaging related parameters

	Parameters
	----------
	msname : str 
		Name of the measurement set
	includeflag : bool
			Include flag data or not whilw calculating maximum baseline
	'''
	def __init__(self,msname,includeflag=False):
		self.msname=msname
		self.md=msmetadata()
		self.md.open(msname)
		self.freq=self.md.meanfreq(0)
		self.md.close()
		self.tb=table()
		self.max_baseline=self.calc_max_baseline(includeflag=includeflag)

	############################################
	# Imaging related #
	############################################
	@staticmethod
	def calc_sun_dia(freq):
		'''
		Function to calculate the diameter of the Sun at a given frequency (White 2016)
		Parameters
		----------
		freq : float
			Frequency in MHz
		Returns
		-------
		float
			Diameter of the Sun in arcmin
		'''
		freq=freq/10**3 # Convert in GHz
		dia=32+(2.2*(freq)**(-0.6))
		return dia

	def calc_psf(self,freq=0):
		'''
		Function to calculate PSF size in arcsec

		Parameters
		----------
		freq : float
	 		Frequency in MHz (default : 0, using central frequency of the ms)
		Returns
		-------
		float
			PSF size in arcsec
		'''
		if freq==0:
			wavelength = 299792458.0/self.freq
		else:
			wavelength = 299792458.0/(freq*10**6)
		psf	= (1.2*(wavelength/self.max_baseline))*(180/np.pi*3600.0) # In arcsec
		return psf

	def calc_cellsize(self,num_pixel_in_psf,freq=0):
		'''
		Calculate pixel size in arcsec

		Parameters
		----------
		num_pixel_in_psf : int 
			Number of pixels in one PSF
		freq : float 
			Frequency in MHz (default : 0, using central frequency of the ms)
		Returns
		-------
		float
			Pixel size in arcsec
		'''
		psf	=	self.calc_psf(freq=freq)	
		pixel	=	int(psf/num_pixel_in_psf) 
		return pixel

	def choose_scales(self,num_pixel_in_psf,max_size,freq=0):
		'''
		Function to calculate multiscale scales

		Parameters
		----------
		num_pixel_in_psf : int 
			Number of pixels in one PSF
		max_size : float 
			Maximum source size in arcsec
		freq : float 
			Frequency in MHz (default : 0, using central frequency of the ms)
		Returns
		-------
		list
			List of multiscale lists in number of pixels
		'''	
		if freq==0:
			freq=self.freq/10**6
		psf=self.calc_psf(freq=freq)
		cellsize=self.calc_cellsize(num_pixel_in_psf,freq=freq)
		psf_pix=num_pixel_in_psf
		sun_size=(ImageBasic.calc_sun_dia(freq)/2)*60
		if sun_size<max_size:
			max_size=sun_size
		scale=[0,psf_pix,int(max_size/cellsize)]
		if 3*psf_pix<int(max_size/cellsize):
			for i in range((3*psf_pix)+num_pixel_in_psf,int(max_size/cellsize),psf_pix):
				scale.append(i)
		elif 2*psf_pix<int(max_size/cellsize):
			for i in range((2*psf_pix)+num_pixel_in_psf,int(max_size/cellsize),psf_pix):
				scale.append(i)
		elif psf_pix<int(max_size/cellsize):
			for i in range((psf_pix)+num_pixel_in_psf,int(max_size/cellsize),psf_pix):
				scale.append(i)
		final_scale=[]
		for i in scale:
			if i<15:
				final_scale.append(i)		
		return sorted(final_scale)

	def MWA_field_of_view(self,freq=0,FWHM=True):
		'''
		Calculate optimum field of view in arcsec

		Parameters
		----------
		freq : float 
			Frequency in MHz (default : 0, using central frequency of the ms)
		FWHM : bool
			Upto FWHM, otherwise upto first null
		Returns
		-------
		float
			Field of view in arcsec
		'''
		if FWHM==True:
			if freq==0:
				FOV=np.sqrt(610)*150*10**6/self.freq  # 600 deg^2 is the image FoV at 150MHz for MWA. So extrapolating this to central frequency
			else:
				FOV=np.sqrt(610)*150/freq  # 610 deg^2 is the image FoV at 150MHz for MWA. So extrapolating this to central frequency
		else:
			if freq==0:
				FOV=60*110*10**6/self.freq  # 3600 deg^2 is the image FoV at 110MHz for MWA upto first null. So extrapolating this to central frequency
			else:
				FOV=60*110/freq  # 3600 deg^2 is the image FoV at 110MHz for MWA upto first null. So extrapolating this to central frequency
		return FOV*3600 ### In arcsecs

	def num_pixels(self,num_pixel_in_psf,freq=0,FWHM=True):
		'''
		Number of image pixels
	
		Parameters
		----------
		num_pixel_in_psf : int
			Number of pixels in one PSF
		freq : float 
			Frequency in MHz (default : 0, using central frequency of the ms)
		FWHM : bool
			Upto FWHM, otherwise upto first null
		Returns
		-------
		int
			Number of pixels in the image
		'''
		FOV=self.MWA_field_of_view(freq=freq,FWHM=FWHM)
		cellsize=self.calc_cellsize(num_pixel_in_psf,freq=freq)
		num=FOV/cellsize
		pow2=int(np.log2(num))
		return 2**pow2

	def calc_calib_uvrange(self,max_angular_scale,uvbin=10,includeflag=False):
		'''
		This function calculate the uvrange to be used for calibration. 
	
		Parameters
		----------
		max_angular_scale : float 
			Maximum angular scale to exclude short baselines from calibration (In degree)
		uvbin : float
			Binning in uv-lambda
		includeflag : bool
			Include flag data or not whilw calculating maximum baseline
		Returns
		-------
		str
			uv-range string (uvmin~uvmax lambda)
		float
			uvmin in meter
		float
			uvmax in meter
		float
			uvmin in wavelength unit
		float
			uvmax in wavelength unit
		'''
		wavelength=299792458.0/self.freq
		uvmin_lambda=1.22/np.deg2rad(max_angular_scale)
		uvmin=uvmin_lambda*wavelength
		self.tb.open(self.msname)
		uvw=self.tb.getcol('UVW')
		if includeflag==False:
			flag=self.tb.getcol('FLAG')
			flag=np.prod(np.prod(flag,axis=0,dtype='bool'),axis=0,dtype='bool')
			for i in range(3):
				uvw[i][flag]=np.nan
		self.tb.close()
		u,v,w=[uvw[i, :] for i in range(3)]
		uvdist=np.sqrt(u**2+v**2)
		nanpos=np.where(np.isnan(uvdist)==True)
		uvdist=np.delete(uvdist,nanpos)
		uvlambda=uvdist/wavelength
		uvlambda_hist=np.histogram(uvlambda,bins=int(max(uvlambda)/uvbin))
		max_uvpoints=np.max(uvlambda_hist[1])
		try:
			cutpos1=np.min(np.where(uvlambda_hist[0]<max_uvpoints*0.01))
			uvlambda1=uvlambda_hist[1][cutpos1]
		except:
			uvlambda1=np.nanmax(uvlambda)
		try:
			cutpos2=np.min(np.where(uvlambda_hist[0]==0))
			uvlambda2=uvlambda_hist[1][cutpos2]
		except:
			uvlambda2=uvlambda1
		uvmax=min(uvlambda1,uvlambda2)
		uvmax=uvmax*wavelength
		uvmax_lambda=uvmax/wavelength
		return str(int(uvmin_lambda))+'~'+str(int(uvmax_lambda))+'lambda',int(uvmin),int(uvmax),int(uvmin_lambda),int(uvmax_lambda)

	def calc_uvtaper(self,uvbin=10,includeflag=False):
		'''
		Function return uv-taper value
		Parameters
		----------
		uvbin : float
			Binning in uv-lambda
		includeflag : bool
			Include flag data or not whilw calculating maximum baseline
		Returns
		-------
		str
			UV taper (uvtaper lambda)
		'''
		wavelength=299792458.0/self.freq
		self.tb.open(self.msname)
		uvw=self.tb.getcol('UVW')
		if includeflag==False:
			flag=self.tb.getcol('FLAG')
			flag=np.prod(np.prod(flag,axis=0,dtype='bool'),axis=0,dtype='bool')
			for i in range(3):
				uvw[i][flag]=np.nan
		self.tb.close()
		u,v,w=[uvw[i, :] for i in range(3)]
		uvdist=np.sqrt(u**2+v**2)
		nanpos=np.where(np.isnan(uvdist)==True)
		uvdist=np.delete(uvdist,nanpos)
		uvlambda=uvdist/wavelength
		uvlambda_hist=np.histogram(uvlambda,bins=int(max(uvlambda)/uvbin))
		max_uvpoints=np.max(uvlambda_hist[1])
		try:
			cutpos1=np.min(np.where(uvlambda_hist[0]<max_uvpoints*0.01))
			uvlambda1=uvlambda_hist[1][cutpos1]
		except:
			uvlambda1=np.nanmax(uvlambda)
		try:
			cutpos2=np.min(np.where(uvlambda_hist[0]==0))
			uvlambda2=uvlambda_hist[1][cutpos2]
		except:
			uvlambda2=uvlambda1
		uvmax=min(uvlambda1,uvlambda2)
		uvmax=uvmax*wavelength
		uvmax_lambda=uvmax/wavelength
		return str(uvmax_lambda)+'lambda'

	def calc_suntaper(self,uvbin=10,includeflag=False):
		'''
		Function return uv-taper value to treat Sun as a point source of size 32 arcmin
		Parameters
		----------
		uvbin : float
			Binning in uv-lambda
		includeflag : bool
			Include flag data or not whilw calculating maximum baseline
		Returns
		-------
		str
			UV taper (uvtaper lambda)
		str
			UV taper in meter
		'''
		uvlimstring,umin,umax,uvminlambda,uvmaxlambda=self.calc_calib_uvrange(32/60.0,uvbin=uvbin,includeflag=includeflag)
		return str(uvminlambda)+'lambda',umin

	def calc_max_baseline(self,includeflag=False):
		'''
		Get the maximum baseline in meter
		Parameters
		----------
		includeflag : bool
			Include flag data or not while calculating maximum baseline
		Returns
		-------
		float
			Maximum baseline length in meter
		'''
		self.tb.open(self.msname)
		uvw=self.tb.getcol('UVW')
		if includeflag==False:
			flag=self.tb.getcol('FLAG')
			flag=np.prod(np.prod(flag,axis=0,dtype='bool'),axis=0,dtype='bool')
			for i in range(3):
				uvw[i][flag]=np.nan
		self.tb.close()
		u,v,w=[uvw[i, :] for i in range(3)]
		uvdist=np.sqrt(u**2+v**2)
		return int(np.nanmax(uvdist))
	
	def calc_max_baseline_lambda(self,includeflag=False):
		'''
		Get the maximum baseline in meter
		Parameters
		----------
		includeflag : bool
			Include flag data or not while calculating maximum baseline
		Returns
		-------
		float
			Maximum baseline in lambda
		'''
		self.tb.open(self.msname)
		uvw=self.tb.getcol('UVW')
		if includeflag==False:
			flag=self.tb.getcol('FLAG')
			flag=np.prod(np.prod(flag,axis=0,dtype='bool'),axis=0,dtype='bool')
			for i in range(3):
				uvw[i][flag]=np.nan
		self.tb.close()
		u,v,w=[uvw[i, :] for i in range(3)]
		uvdist=np.sqrt(u**2+v**2)
		wavelength=299792458.0/self.freq
		uvmax_lambda=uvdist/wavelength
		return int(np.nanmax(uvmax_lambda))

	def calc_max_uvw(self):
		'''
		Get the maximum u,v,w

		Returns
		-------
		float
			Maximum U
		float
			Maximum V
		float 
			Maximum W
		'''
		self.tb.open(self.msname)
		uvw=self.tb.getcol('UVW')
		flag=self.tb.getcol('FLAG')
		flag=np.prod(np.prod(flag,axis=0,dtype='bool'),axis=0,dtype='bool')
		self.tb.close()
		for i in range(3):
			uvw[i][flag]=np.nan
		u,v,w=[uvw[i, :] for i in range(3)]
		maxu=np.nanmax(np.abs(u))
		maxv=np.nanmax(np.abs(v))
		maxw=np.nanmax(np.abs(w))
		return maxu,maxv,maxw

	def get_best_chan(self,rmsbox='',chanrange='',timerange='',weight='natural',robust=0.5,refant=1,max_angular_scale=12,docal=True,DRbased=False,caltables=[],verbose=False):
		'''
		Function to calculate rms per channel in a given region of the image.
		Parameters
		----------
		rmsbox : str
			Box in CASA format calculate rms (default : full image)
		chanrange : str
			Channel range to image (startchan~endchan) 
		timerange : str
			Time range 
		weight : str
			Image weighting
		robust : float
			Robust parameter for \'briggs\' weighting
		refant : int
			Reference antenna
		max_angular_scale : float
			Maximum allowed angular scale in degree
		docal : bool
			Perform point source model calibration
		DRbased : bool
			Whether search best time based on DR or rms
		caltables : list
			List of caltable solutions to apply (If docal = False)
		verbose : bool
			Verbose output or not
		Returns
		-------
		int
			Channel number with lowest rms
		list
			List of channels
		list
			List of rms
		'''
		check_thread_availability(4)
		print ('Determining best frequency channel.....\n')
		if chanrange!='' and  '~' in chanrange:	
			start_chan=int(chanrange.split('~')[0])
			end_chan=int(chanrange.split('~')[1])
		else:
			chanrange=''
			start_chan=0
		if os.path.exists(self.msname.split('.ms')[0]+'_chancube.ms'):
			os.system('rm -rf '+self.msname.split('.ms')[0]+'_chancube.ms '+self.msname.split('.ms')[0]+'_chancube.ms.flagversions')
		if os.path.exists(self.msname.split('.ms')[0]+'_chancube.ms.flagversions'):
			os.system('rm -rf '+self.msname.split('.ms')[0]+'_chancube.ms.flagversions')
		if timerange!='' and chanrange=='':	
			if verbose:
				print ('split(vis=\''+self.msname+'\',outputvis=\''+self.msname.split('.ms')[0]+'_chancube.ms\',timerange=\''+timerange+'\',datacolumn=\'data\')\n')
			if os.path.exists(self.msname.split('.ms')[0]+'_chancube.ms'):
				os.system('rm -rf '+self.msname.split('.ms')[0]+'_chancube.ms '+self.msname.split('.ms')[0]+'_chancube.ms.flagversions')
			split(vis=self.msname,outputvis=self.msname.split('.ms')[0]+'_chancube.ms',timerange=timerange,datacolumn='data')		
		elif timerange=='' and chanrange!='':
			if verbose:
				print ('split(vis=\''+self.msname+'\',outputvis=\''+self.msname.split('.ms')[0]+'_chancube.ms\',spw=\'0:'+chanrange+'\',datacolumn=\'data\')\n')
			if os.path.exists(self.msname.split('.ms')[0]+'_chancube.ms'):
				os.system('rm -rf '+self.msname.split('.ms')[0]+'_chancube.ms '+self.msname.split('.ms')[0]+'_chancube.ms.flagversions')
			split(vis=self.msname,outputvis=self.msname.split('.ms')[0]+'_chancube.ms',spw='0:'+chanrange,datacolumn='data')	
		elif timerange!='' and chanrange!='':
			if verbose:
				print ('split(vis=\''+self.msname+'\',outputvis=\''+self.msname.split('.ms')[0]+'_chancube.ms\',timerange=\''+timerange+'\',spw=\'0:'+chanrange+\
						',datacolumn=\'data\')\n')
			if os.path.exists(self.msname.split('.ms')[0]+'_chancube.ms'):
				os.system('rm -rf '+self.msname.split('.ms')[0]+'_chancube.ms '+self.msname.split('.ms')[0]+'_chancube.ms.flagversions')
			split(vis=self.msname,outputvis=self.msname.split('.ms')[0]+'_chancube.ms',timerange=timerange,spw='0:'+chanrange,datacolumn='data')		
		else:
			os.system('cp -r '+self.msname+' '+self.msname.split('.ms')[0]+'_chancube.ms')
		temp_ms=self.msname.split('.ms')[0]+'_chancube.ms'
		if verbose:
			print ('Want point source calibration : '+str(docal)+'\n')
			print ('Supplied caltables : '+str(caltables)+'\n')
		if docal==True and len(caltables)==0:
			if verbose:
				print ('gaincal(vis=\''+temp_ms+'\',caltable=\'point_source_try.cal\',smodel=[1,0,0,0],refant=\''+str(refant)+'\',solnorm=True,solint=\'5s\',calmode=\'p\')\n')
			gaincal(vis=temp_ms,caltable='point_source_try.cal',smodel=[1,0,0,0],refant=str(refant),solnorm=True,solint='5s',calmode='p')
			if verbose:
				print ('applycal(vis=\''+temp_ms+'\',gaintable=[\'point_source_try.cal\'],applymode=\'calonly\',interp=[\'linear\'])\n')
			applycal(vis=temp_ms,gaintable=['point_source_try.cal'],applymode='calonly',interp=['linear'])
			datacolumn_to_split='corrected'
		elif len(caltables)!=0:
			if verbose:
				print ('applycal(vis=\''+temp_ms+'\',gaintable=\''+str(caltables)+'\',applymode=\'calflag\',interp=[\'linear\'])\n')
			applycal(vis=temp_ms,gaintable=caltables,applymode='calflag',interp=['linear'])
			datacolumn_to_split='corrected'
		else:
			datacolumn_to_split='data'	
		os.system('rm -rf point_source_try.cal')
		imsize=self.num_pixels(3,FWHM=False)
		cell=self.calc_cellsize(3)
		uvtaper=self.calc_uvtaper(uvbin=10,includeflag=False)
		if rmsbox=='':
			rmsbox='50,50,'+str(imsize-50)+','+str(int(imsize/4))
		uvrange_list=self.calc_calib_uvrange(max_angular_scale,uvbin=10,includeflag=False)
		uvrange=uvrange_list[0]
		min_uvrange=uvrange_list[3]
		max_uvrange=uvrange_list[4]
		if weight!='natural' and weight!='briggs' and weight!='uniform':
			print ('Wrong weight parameter : '+weight+'. Instead of it using natural weighting.\n')
			weight='natural'
		if weight=='briggs':
			weight=weight+' '+str(robust)
		if os.path.isdir(temp_ms.split('.ms')[0]+'/tempdir')==False:
			os.makedirs(temp_ms.split('.ms')[0]+'/tempdir')
		tmp_dir=temp_ms.split('.ms')[0]+'/tempdir'
		AM=am.AccessMS(temp_ms)
		unflag_chans=AM.get_unflag_chan()
		nchan=AM.get_num_channels()
		wsclean_args=['-scale '+str(cell)+'asec','-size '+str(imsize)+' '+str(imsize),'-weight '+weight,'-taper-edge-tukey '+str(20),'-name '+temp_ms.split('.ms')[0]+'_chancube',\
					'-nwlayers 1','-niter 0','-mgain 0.8','-pol i','-taper-edge '+str(max_uvrange),'-minuv-l '+str(min_uvrange),'-no-reorder','-temp-dir '+tmp_dir,\
					'-channels-out '+str(nchan),'-j 4','-abs-mem 1']
		wsclean_args.append('-quiet')
		wsclean=get_wsclean_path()
		if wsclean!=None:
			if verbose:
				print (wsclean+' '+' '.join(wsclean_args)+' '+temp_ms+'\n')
			os.system(wsclean+' '+' '.join(wsclean_args)+' '+temp_ms)		
			os.system('rm -rf '+temp_ms.split('.ms')[0]+'_chancube*-image.fits '+temp_ms.split('.ms')[0]+'_chancube-MFS*.fits')
			imagelist=sorted(glob.glob(temp_ms.split('.ms')[0]+'_chancube-*-dirty.fits'))
			if len(imagelist)==0:
				print ('No image is made.\n')
				os.system('rm -r '+temp_ms.split('.ms')[0])
				return 
			chanlist=[]
			rmslist=[]
			DRlist=[]
			for image in imagelist:
				rms=imstat(imagename=image,box=rmsbox)['rms'][0]
				max_pix=imstat(imagename=image)['max'][0]
				if rms!=0:
					chan=int(os.path.basename(image).split('chancube-')[-1].split('-dirty')[0])
					if chan in unflag_chans:
						rmslist.append(rms)
						chanlist.append(chan)
						DRlist.append(max_pix/rms)
			os.system('rm -rf '+temp_ms)
		else:
			os.system('rm -rf '+temp_ms.split('.ms')[0]+'_chancube*')
			if verbose:
				print ('tclean(vis=\''+temp_ms+'\',timerange=\''+timerange+'\',imagename=\''+temp_ms.split('.ms')[0]+'_chancube\',imsize='+str(imsize)+\
					',cell=\''+str(cell)+'arcsec\',specmode=\'cube\',niter=0,weighting=\''+weight+'\',robust='+str(robust)+',uvrange=\''+uvrange+'\',uvtaper=\''+uvtaper+'\')\n')
			tclean(vis=temp_ms,timerange=timerange,imagename=temp_ms.split('.ms')[0]+'_chancube',imsize=imsize,\
					cell=str(cell)+'arcsec',specmode='cube',niter=0,weighting=weight,robust=robust,uvrange=uvrange,uvtaper=uvtaper)	
			nchan=imhead(imagename=temp_ms.split('.ms')[0]+'_chancube.image',mode='summary')['shape'][-1]
			chanlist=[]
			rmslist=[]
			DRlist=[]
			for chan in range(nchan):
				rms=imstat(imagename=temp_ms.split('.ms')[0]+'_chancube.image',box=rmsbox,chans=str(chan))['rms'][0]
				max_pix=imstat(imagename=image)['max'][0]
				if rms!=0:
					if chan in unflag_chans:
						chanlist.append(chan)
						rmslist.append(rms)
						DRlist.append(max_pix/rms)
		if DRbased==False:
			if len(rmslist)>0 and len(chanlist)>0:
				minrms=np.argmin(rmslist)
				bestchan=chanlist[minrms]
				bestchan+=start_chan
				chanlist=[i+start_chan for i in chanlist]
				os.system('rm -rf '+temp_ms.split('.ms')[0]+'_chancube* '+tmp_dir+' '+temp_ms+'*')
				os.system('rm -r '+temp_ms.split('.ms')[0])
				return bestchan,chanlist,rmslist
			else:
				os.system('rm -rf '+temp_ms.split('.ms')[0]+'_chancube* '+tmp_dir+' '+temp_ms+'*')
				os.system('rm -r '+temp_ms.split('.ms')[0])
				return 
		else:
			if len(DRlist)>0 and len(chanlist)>0:
				maxDR=np.argmax(DRlist)
				bestchan=chanlist[maxDR]
				os.system('rm -rf '+temp_ms.split('.ms')[0]+'_chancube* '+tmp_dir+' '+temp_ms+'*')
				os.system('rm -r '+temp_ms.split('.ms')[0])
				bestchan+=start_chan
				chanlist=[i+start_chan for i in chanlist]
				return bestchan,chanlist,DRlist
			else:
				os.system('rm -rf '+temp_ms.split('.ms')[0]+'_chancube* '+tmp_dir+' '+temp_ms+'*')
				os.system('rm -r '+temp_ms.split('.ms')[0])
				return 

	def get_best_time(self,rmsbox='',time_list=[],chanrange='',weight='natural',robust=0.5,max_angular_scale=12,refant=1,docal=True,DRbased=True,caltables=[],verbose=False,n_cpu=1):
		'''
		Function to calculate rms per timestamp in a given region of the image.
		Parameters
		----------
		rmsbox : str
			Box in CASA format calculate rms (default : full image)
		timelist : list
			Timestamp list to find the best time (timeformat : 'YYYY/MM/DD/hh:mm:ss') 
		chanrange : str
			Channel range to image (startchan~endchan) 
		weight : str
			Image weighting
		robust : float
			Robust parameter for \'briggs\' weighting
		max_angular_scale : float
			Maximum allowed angular scale in degree
		refant : int
			Reference antenna
		docal : bool
			Perform point source model calibration
		DRbased : bool
			Whether search best time based on DR or rms
		caltables : list
			List of caltable solutions to apply (If docal = False)
		verbose : bool
			Verbose output or not
		n_cpu : int
			Number of CPU threads to use
		Returns
		-------
		int
			Timestamp with lowest rms (timeformat : 'YYYY/MM/DD/hh:mm:ss')
		list
			List of timestamps
		list
			List of rms/DR
		''' 
		check_thread_availability(4)
		print ('Determining best timestamp.....\n')
		chanrange=str(chanrange)
		temp_ms=self.msname.split('.ms')[0]+'_splited.ms'
		if os.path.exists(temp_ms):
			os.system('rm -rf '+temp_ms+' '+temp_ms+'.flagversions')
		if os.path.exists(temp_ms+'.flagversions'):
			os.system('rm -rf '+temp_ms+'.flagversions')
		if len(time_list)!=0:
			timelist_mjd=[timestamp_to_mjdsec(i) for i in time_list]
			timerange=mjdsec_to_timestamp(min(timelist_mjd))+'~'+mjdsec_to_timestamp(max(timelist_mjd))
		else:
			timerange=''
		if chanrange!='':	
			if '~' in chanrange:
				start_chan=int(chanrange.split('~')[0])
				end_chan=int(chanrange.split('~')[1])
				width=end_chan-start_chan
			else:
				width=1
			if timerange!='':
				if verbose:
					print ('split(vis=\''+self.msname+'\',outputvis=\''+temp_ms+'\',spw=\'0:'+chanrange+'\',width='+str(width)+',datacolumn=\'data\',timerange=\''+timerange+'\')\n')
				split(vis=self.msname,outputvis=temp_ms,spw='0:'+chanrange,width=width,datacolumn='data',timerange=timerange)
			else:
				if verbose:
					print ('split(vis=\''+self.msname+'\',outputvis=\''+temp_ms+'\',spw=\'0:'+chanrange+'\',width='+str(width)+',datacolumn=\'data\')\n')
				split(vis=self.msname,outputvis=temp_ms,spw='0:'+chanrange,width=width,datacolumn='data')
		else:
			temp_ms=self.msname
		if verbose:
			print ('Want point source calibration : '+str(docal)+'\n')
			print ('Supplied caltables : '+str(caltables)+'\n')
		if docal==True and len(caltables)==0:
			if verbose:
				print ('gaincal(vis=\''+temp_ms+'\',caltable=\'point_source_try.cal\',smodel=[1,0,0,0],refant=\''+str(refant)+'\',solnorm=True,solint=\'5s\',calmode=\'p\')\n')
			gaincal(vis=temp_ms,caltable='point_source_try.cal',smodel=[1,0,0,0],refant=str(refant),solnorm=True,solint='5s',calmode='p')
			if verbose:
				print ('applycal(vis=\''+temp_ms+'\',gaintable=[\'point_source_try.cal\'],applymode=\'calonly\',interp=[\'nearest\'])\n')
			applycal(vis=temp_ms,gaintable=['point_source_try.cal'],applymode='calonly',interp=['nearest'])
		elif len(caltables)!=0:
			if verbose:
				print ('applycal(vis=\''+temp_ms+'\',gaintable=\''+str(caltables)+'\',applymode=\'calflag\',interp=[\'nearest\'])\n')
			applycal(vis=temp_ms,gaintable=caltables,applymode='calflag',interp=['nearest'])
		os.system('rm -rf point_source_try.cal')
		AM=am.AccessMS(temp_ms)
		timestamps=AM.get_timestamps(format=0)
		unflag_timestamps=AM.get_unflag_timestamps()
		if type(time_list)==np.ndarray:
			time_list=time_list.tolist()
		if len(time_list)==0:
			time_list=unflag_timestamps
		timelist_mjd=np.array([timestamp_to_mjdsec(i) for i in time_list])
		timestamps_mjd=np.array([timestamp_to_mjdsec(i) for i in timestamps])
		uvrange=self.calc_calib_uvrange(max_angular_scale,uvbin=10,includeflag=True)
		min_uvrange=uvrange[3]
		max_uvrange=uvrange[4]
		imsize=self.num_pixels(3,FWHM=False)
		cell=self.calc_cellsize(3)
		if rmsbox=='':
			rmsbox='50,50,'+str(imsize-50)+','+str(int(imsize/4))
		if weight!='natural' and weight!='briggs' and weight!='uniform':
			print ('Wrong weight parameter : '+weight+'. Instead of it using natural weighting.\n')
			weight='natural'
		if weight=='briggs':
			weight=weight+' '+str(robust)
		if os.path.isdir(temp_ms.split('.ms')[0]+'/tempdir')==False:
			os.makedirs(temp_ms.split('.ms')[0]+'/tempdir')
		tmp_dir=temp_ms.split('.ms')[0]+'/tempdir'
		intervals=len(timestamps)
		wsclean_args=['-scale '+str(cell)+'asec','-size '+str(imsize)+' '+str(imsize),'-weight '+weight,'-name '+temp_ms.split('.ms')[0]+'_timecube','-nwlayers 1','-niter 0','-mgain 0.8','-pol i','-minuv-l '+str(min_uvrange),\
					'-no-reorder','-temp-dir '+tmp_dir,'-intervals-out '+str(intervals),'-j '+str(n_cpu)]
		wsclean_args.append('-quiet')
		wsclean=get_wsclean_path()
		if wsclean==None:
			print ('WSClean is not installed.')
			return
		if verbose:
			print (wsclean+' '+' '.join(wsclean_args)+' '+temp_ms+'\n')
		os.system(wsclean+' '+' '.join(wsclean_args)+' '+temp_ms)			
		time_indices=[]
		for i in timelist_mjd:
			try:
				pos=np.argmin(np.abs(timestamps_mjd-i))
				time_indices.append(pos)
			except:
				pass
		os.system('rm -rf '+temp_ms.split('.ms')[0]+'_timecube*-image.fits')
		imagelist=sorted(glob.glob(temp_ms.split('.ms')[0]+'_timecube-t*-dirty.fits'))
		if len(imagelist)==0:
			print ('No image is made.\n')
			return 
		final_timelist=[]
		rmslist=[]
		DRlist=[]
		for image in imagelist:
			rms=imstat(imagename=image,box=rmsbox)['rms'][0]
			max_pix=imstat(imagename=image)['max'][0]
			if rms>0:
				time_index=int(os.path.basename(image).split('-t')[1].split('-dirty')[0])
				if time_index in time_indices:
					rmslist.append(rms)
					DRlist.append(max_pix/rms)
					final_timelist.append(timestamps[time_index])
		if DRbased==False:
			if len(rmslist)>0 and len(final_timelist)>0:
				minrms=np.argmin(rmslist)
				besttime=final_timelist[minrms]
				sorted_pos=np.argsort(np.array(minrms))
				final_timelist=(np.array(final_timelist)[sorted_pos]).tolist()
				rmslist=(np.array(rmslist)[sorted_pos]).tolist()
				os.system('rm -rf '+temp_ms.split('.ms')[0]+'_timecube*-dirty.fits '+tmp_dir+' '+self.msname.split('.ms')[0]+'_splited*')
				return besttime,final_timelist,rmslist
			else:
				os.system('rm -rf '+temp_ms.split('.ms')[0]+'_timecube*-dirty.fits '+tmp_dir+' '+self.msname.split('.ms')[0]+'_splited*')
				return
		else:
			if len(DRlist)>0 and len(final_timelist)>0:
				maxDR=np.argmax(DRlist)
				besttime=final_timelist[maxDR]
				sorted_pos=np.argsort(-1*np.array(DRlist))
				final_timelist=(np.array(final_timelist)[sorted_pos]).tolist()
				DRlist=(np.array(DRlist)[sorted_pos]).tolist()
				os.system('rm -rf '+temp_ms.split('.ms')[0]+'_timecube*-dirty.fits '+tmp_dir+' '+self.msname.split('.ms')[0]+'_splited*')
				return besttime,final_timelist,DRlist
			else:
				os.system('rm -rf '+temp_ms.split('.ms')[0]+'_timecube*-dirty.fits '+tmp_dir+' '+self.msname.split('.ms')[0]+'_splited*')
				return

#########################################
# Calibration parameter estimation class
#########################################

class CalcParams:
	'''
	Calculate calibration parameters	

	Parameters
	----------
	msname : str
		Name of the measurement set
	quality_factor : int 
		Imaging quality factor (0,1,2)
	safety_factor : int
		Safety factor for calibration (0,1,2)
	'''
	def __init__(self,msname,quality_factor,safety_factor):
		self.msname=msname
		self.quality_factor=quality_factor
		self.safety_factor=safety_factor
		self.IB=ImageBasic(self.msname)
		self.AM=am.AccessMS(self.msname)
		
	def calc_calibration_params(self):
		'''
		Return calibration parameters automatically based on given quality factor and safety standard
		For details read the documentation of quality_factor and safety_factor

		Returns
		-------
		float
			Starting sigma value for calibration
		float
			Step to reduce sigma value with self-calibration
		float
			Residual flux fraction to stop the calibration
		float
			Minimum allowed sigma for threshold
		float
			Minimum gain SNR
		float
			Increment of DR_rms
		float
			Increment of DR_neg
		float
			Minimum allowed dynamic range
		float
			Maximum dynamic range
		float
			Minimum antenna based SNR for self-calibration
		float
			Time interval for calibration 
		float
			Frequency interval of calibration
		float
			uvrange for calibration
		'''
		if self.quality_factor==0:
			residual_frac	=	0.03
			if self.safety_factor==0:
				start_sigma		=	9.0
				sigma_step		=	1.0
				min_sigma		=	6
				gain_minsnr		=	3.0		
				DR_delta_rms	=	25
				DR_delta_neg	=	20
				min_DR			=	20
				max_DR			=	100
				min_selfcal_snr	=	2.5
				skip_time		=	25
				skip_freq		=	1280
			elif self.safety_factor==1:
				start_sigma		=	9.0
				sigma_step		=	1.0
				min_sigma		=	6.5
				gain_minsnr		=	3.0		
				DR_delta_rms	=	22
				DR_delta_neg	=	18
				min_DR			=	25
				max_DR			=	500
				min_selfcal_snr	=	2.5
				skip_time		=	20
				skip_freq		=	1280
			else:
				start_sigma		=	9.0
				sigma_step		=	1.0
				min_sigma		=	7
				gain_minsnr		=	3.5	
				DR_delta_rms	=	20
				DR_delta_neg	=	15
				min_DR			=	30
				max_DR			=	1000
				min_selfcal_snr	=	2.5
				skip_time		=	15
				skip_freq		=	1280
		elif self.quality_factor==1:
			residual_frac	=	0.015
			if self.safety_factor==0:
				start_sigma		=	10.0
				sigma_step		=	0.5
				min_sigma		=	6
				gain_minsnr		=	3.5		
				DR_delta_rms	=	20
				DR_delta_neg	=	15
				min_DR			=	30
				max_DR			=	1000
				min_selfcal_snr	=	3
				skip_time		=	20
				skip_freq		=	1280
			elif self.safety_factor==1:
				start_sigma		=	10.0
				sigma_step		=	0.5
				min_sigma		=	6.5
				gain_minsnr		=	4		
				DR_delta_rms	=	18
				DR_delta_neg	=	12
				min_DR			=	35
				max_DR			=	5000
				min_selfcal_snr	=	3
				skip_time		=	15
				skip_freq		=	1280
			else:
				start_sigma		=	10.0
				sigma_step		=	0.5
				min_sigma		=	7
				gain_minsnr		=	4		
				DR_delta_rms	=	15
				DR_delta_neg	=	10
				min_DR			=	40
				max_DR			=	10000
				min_selfcal_snr	=	3
				skip_time		=	10
				skip_freq		=	1280
		else:
			residual_frac	=	0.01
			if self.safety_factor==0:
				start_sigma		=	11.0
				sigma_step		=	0.25
				min_sigma		=	6.5
				gain_minsnr		=	4		
				DR_delta_rms	=	18
				DR_delta_neg	=	12
				min_DR			=	40
				max_DR			=	10000
				min_selfcal_snr	=	3
				skip_time		=	15
				skip_freq		=	1280
			elif self.safety_factor==1:
				start_sigma		=	11.0
				sigma_step		=	0.25
				min_sigma		=	7
				gain_minsnr		=	4.5		
				DR_delta_rms	=	15
				DR_delta_neg	=	10
				min_DR			=	45
				max_DR			=	50000
				min_selfcal_snr	=	3
				skip_time		=	10
				skip_freq		=	1280
			else:
				start_sigma		=	11.0
				sigma_step		=	0.25
				min_sigma		=	7
				gain_minsnr		=	4.5		
				DR_delta_rms	=	12
				DR_delta_neg	=	8
				min_DR			=	50
				max_DR			=	100000
				min_selfcal_snr	=	3
				skip_time		=	5
				skip_freq		=	1280
		uvrange_to_cal=''
		if skip_freq<self.AM.calc_freqres():
			skip_freq=self.AM.calc_freqres()
		if skip_time<self.AM.calc_timeres():
			skip_time=self.AM.calc_timeres()
		return start_sigma,sigma_step,residual_frac,min_sigma,gain_minsnr,DR_delta_rms,DR_delta_neg,min_DR,max_DR,min_selfcal_snr,skip_time,skip_freq,uvrange_to_cal


###########################################
# Flagging related #
###########################################
def calc_flag_fraction(msname):
	'''
	Function to calculate the fraction of total data flagged

	Parameters
	----------
	msname : str 
		Name of the measurement set
	Returns
	-------
	float
		Fraction of the total data flagged
	'''
	tb=table()
	tb.open(msname)
	flag=tb.getcol('FLAG')
	tb.close()
	flagged_data=np.sum(flag)
	total_data=flag.size
	flagged_fraction=flagged_data/float(total_data)
	return flagged_fraction

def append_cubical_caltable(caltable1,caltable2,outfile,append_axis='time',verbose=False):
	'''
	Function to append multiple CubiCal caltable of same solution type
	Parameters
	----------
	caltable1 : str
		First caltable name
	caltable2 : str
		Seceond caltable name
	outfile : str
		Name of the output table
	append_axis : str
		Append axis , time or frequency
	verbose : bool
		Verbose output
	Returns
	-------
	str
		Merged caltable name
	'''
	if verbose:
		verbose=0
	else:
		verbose=None
	db1=param_db.load(caltable1,verbose=verbose)
	names1=db1.names()
	db2=param_db.load(caltable2,verbose=verbose)
	names2=db2.names()
	if names1!=names2:
		print ('Mismatch in gaintable solution type. Can not merge.\n')
		db1.close()
		db2.close()
		return
	else:
		db=param_db.create(outfile+'.temp',verbose=verbose)
		for i in range(len(names1)):
			g1=db1[names1[i]]
			g2=db2[names1[i]]
			g1_index=g1.grid_index
			g2_index=g2.grid_index
			g1_data=g1.get_cube()
			g2_data=g2.get_cube()
			axis_labels=g1.axis_labels
			g1_shape=g1_data.shape
			g2_shape=g2_data.shape
			db.define_param(names1[i],db1[names1[i]].dtype,db1[names1[i]].axis_labels,interpolation_axes=["time", "freq"])	
			g_index=copy.deepcopy(g1_index)	
			if append_axis=='time':
				for p in list(g2_index['time'].keys()):
					if p not in list(g_index['time'].keys()):
						g_index['time'][p]=0
				g_data=np.ma.append(g1_data,g2_data,axis=1)
			elif append_axis=='freq':
				for p in list(g2_index['freq'].keys()):
					if p not in list(g_index['freq'].keys()):
						g_index['freq'][p]=0
				g_data=np.ma.append(g1_data,g2_data,axis=2)
			else:
				print ('Append axis should be either time or frequency.\n')
				return
			db.add_chunk(names1[i],g_data,grid=g_index)
		db.close()
		if os.path.exists(outfile):
			os.system('rm -rf '+outfile)
		os.system('mv '+outfile+'.temp '+outfile)
		return outfile

def flag_cubical_caltable(caltable,statfile='',sigma=3,overwrite=False,verbose=False):
	'''
	Function to flag bad solutions based on chi-square statistics
	Parameters
	----------
	caltable : str
		Name of the caltable
	statfile : str
		Cubical solution statistics file
	overwrite : bool
		Overwrite the existing caltable
	verbose : bool
		Verbose output
	Returns
	-------
	str
		Output flagged caltable name
	'''
	if verbose:
		verbose=0
	else:
		verbose=None
	if statfile=='' or os.path.exists(statfile)==False:
		print ('CubiCal solutions statistics file does not exist.\n')
		return caltable
	else:
		try:
			stat=np.load(statfile,allow_pickle=True)
			chi2=stat[0]['chi2']
		except Exception as e:
			print ('Error in reading stat file. Error : '+str(e)+'\n')
			return caltable
	zero_chi=np.where(chi2==0)
	chi2[zero_chi]=np.nan
	median_chi=np.nanmedian(chi2)
	std_chi=np.nanstd(chi2)
	pos=np.where(abs(chi2-median_chi)>sigma*std_chi)
	db=param_db.load(caltable,verbose=verbose)
	gb_keys=db.names()
	gain=db.get(gb_keys[0]).get_cube()
	gain_err=db.get(gb_keys[1]).get_cube()
	gain_data=gain.data
	gain_err_data=gain_err.data
	gain_flags=gain.mask
	for t in range(gain_flags.shape[0]):
		for dirc in range(gain_flags.shape[1]):
			gain_flags[t,dirc,:,:,0,0][pos]=True
			gain_flags[t,dirc,:,:,0,1][pos]=True
			gain_flags[t,dirc,:,:,1,0][pos]=True
			gain_flags[t,dirc,:,:,1,1][pos]=True
	flagged_gain=np.ma.array(gain_data,mask=gain_flags,fill_value=(1),dtype=gain_data.dtype)
	if overwrite==False:
		new_caltable=caltable+'.flagged'
	else:
		os.system('rm -rf '+caltable+' '+caltable+'.skel')
		new_caltable=caltable
	# Creating new caltable
	if os.path.exists(new_caltable):
		os.system('rm -rf '+new_caltable+' '+new_caltable+'.skel')
	db_norm=param_db.create(new_caltable,verbose=verbose)
	for gname in gb_keys:
		db_norm.define_param(gname,db[gname].dtype,db[gname].axis_labels,interpolation_axes=["time", "freq"])
		if 'err' not in gname:
			db_norm.add_chunk(gname,flagged_gain,grid=db[gname].grid_index)
		else:
			db_norm.add_chunk(gname,gain_err,grid=db[gname].grid_index)
	db_norm.close()
	return new_caltable

def multiply_cubical_gaintables(caltable1,caltable2,out_gaintable='multiplied_gaintable.cal',diff_threshold=0.001):
	'''
	Function to matrix multiply of two CubiCal gaintables
	Parameters
	----------
	caltable1 : str
		Name of the first caltable
	caltable2: str
		Name of the second caltable
	out_gaintable : str
		Name of the output gaintable
	diff_threshold : float
		If the different of matrix determinant is greater than this value
	Returns
	-------
	str
		Name of the output gaintable
	bool
		Whether the different of matrix determinant is greater than this value
	'''
	db=param_db.load(caltable1,verbose=False)
	g1=db['G:gain']
	g1_err=db['G:gain.err']
	g1_data=g1.get_cube().data
	g1_err=g1_err.get_cube().data
	db.close()
	db=param_db.load(caltable2,verbose=False)
	gb_keys=db.names()
	g2=db['G:gain']
	g2_err=db['G:gain.err']
	g2_data=g2.get_cube().data
	g2_err=g2_err.get_cube().data
	flags1=g1.get_cube().mask
	flags2=g2.get_cube().mask
	flags=flags1+flags2
	g_final=np.empty(g1_data.shape,dtype=g2_data.dtype)
	det_array=[]
	for t in range(g1_data.shape[1]):
		for f in range(g1_data.shape[2]):
			for ant in range(g1_data.shape[3]):
				g_final[0,t,f,ant,...]=np.dot(g1_data[0,t,f,ant,...],g2_data[0,t,f,ant,...])
				if flags[0,t,f,ant,0,0]!=True:
					det_array.append(np.abs(det(g2_data[0,t,f,ant,...])))
	max_det_value=np.abs(1-np.nanmedian(det_array))
	gain_mult=np.ma.array(g_final,mask=flags,fill_value=(1),dtype=g_final.dtype)
	db_mult=param_db.create(out_gaintable,verbose=False)
	for gname in gb_keys:
		db_mult.define_param(gname,g1_data.dtype,db[gname].axis_labels,interpolation_axes=["time", "freq"])
		if 'err' not in gname:
			db_mult.add_chunk(gname,gain_mult,grid=db[gname].grid_index)
		else:
			db_mult.add_chunk(gname,g2_err,grid=db[gname].grid_index)
	db_mult.close()
	db.close()
	print ('Difference value : '+str(max_det_value)+'\n')
	if max_det_value<diff_threshold:
		return out_gaintable,False,max_det_value
	else:
		return out_gaintable,True,max_det_value


#######################
# Calibration table releated
#######################
def normalize_cubical_caltable(caltable,soltype='G',normtype='median',overwrite=False,verbose=False):
	'''
	Function to normalize cubical caltable
	Parameters
	----------
	caltable : str
		Name of the cubical caltable
	soltype : str
		Solution type ('G' for gaincal or 'B' for bandpass)
	normtype : str
		Normalization type, 'mean' or 'median' (default : median)
	overwrite : bool
		Overwrite the exiting caltable
	verbose : bool
		Verbose output
	Returns
	-------
	str
		New caltable name
	'''
	if verbose:
		verbose=0
	else:
		verbose=None
	db=param_db.load(caltable,verbose=verbose)
	gb_keys=db.names()
	gain=db.get(gb_keys[0]).get_cube()
	gain_err=db.get(gb_keys[-1]).get_cube()
	gain_copy=copy.deepcopy(gain)
	gain_data=gain.data
	gain_flags=gain.mask
	gain_copy_data=gain_copy.data
	gain_data[gain_flags]=np.nan
	gainerr_data=gain_err.data
	gainerr_flags=gain_err.mask
	zero_gain=np.where(gain_data==0)
	gain_data[zero_gain]=np.nan
	if soltype=='G':
		gain_data_1d=gain_data.flatten()
		if normtype=='mean':
			norm_factor=np.nanmean(np.abs(gain_data_1d))
		else:	
			norm_factor=np.nanmedian(np.abs(gain_data_1d))
		gain_copy_data=gain_copy_data/norm_factor
		gain_norm=np.ma.array(gain_copy_data,mask=gain_flags,fill_value=(1),dtype=gain_copy_data.dtype)
	elif soltype=='B':
		pols=['XX','YY']
		for t in range(gain_data.shape[0]):
			for ant in range(gain_data.shape[3]):
				for pol in pols:
					if pol=='XX':
						data=gain_data[t,:,:,ant,0,0]
						if normtype=='mean':
							norm_factor=np.nanmean(data.flatten())
						else:
							norm_factor=np.nanmedian(data.flatten())	
						gain_copy_data[t,:,:,ant,0,0]=gain_copy_data[t,:,:,ant,0,0]/norm_factor
					elif pol=='YY':
						data=gain_data[t,:,:,ant,1,1]
						if normtype=='mean':
							norm_factor=np.nanmean(data.flatten())
						else:
							norm_factor=np.nanmedian(data.flatten())
						gain_copy_data[t,:,:,ant,1,1]=gain_copy_data[t,:,:,ant,1,1]/norm_factor
		gain_copy_data[:,:,:,:,0,1]=gain_copy_data[:,:,:,:,0,1]*0
		gain_copy_data[:,:,:,:,1,0]=gain_copy_data[:,:,:,:,1,0]*0
		gain_copy_datanew=copy.deepcopy(gain_copy_data)
		gain_copy_datanew[gain_flags]=np.nan
		gain_norm=np.ma.array(gain_copy_data,mask=gain_flags,fill_value=(1),dtype=gain_copy_data.dtype)
	else:
		print ('Normalization can only be done for solution type \'G\' or \'B\'.\n')
	if overwrite==False:
		new_caltable=caltable+'.norm'
	else:
		os.system('rm -rf '+caltable+'*')
		new_caltable=caltable
	# Creating new caltable
	if os.path.exists(new_caltable):
		os.system('rm -rf '+new_caltable+'*')
	db_norm=param_db.create(new_caltable,verbose=verbose)
	for gname in gb_keys:
		db_norm.define_param(gname,db[gname].dtype,db[gname].axis_labels,interpolation_axes=["time", "freq"])
		if 'err' not in gname:
			db_norm.add_chunk(gname,gain_norm,grid=db[gname].grid_index)
		else:
			db_norm.add_chunk(gname,gain_err,grid=db[gname].grid_index)
	db_norm.close()
	return new_caltable

def fluxcal_cubical_bptable(caltable,atten=10,verbose=False):
	'''
	Function to scale normalized bandpass table
	Parameters
	----------
	caltable : str
		Name of the caltable
	atten : int
		Attenuation in dB
	verbose : bool
		Verbose output
	Returns
	-------
	str
		Flux calibrated caltable
	'''
	if verbose:
		verbose=0
	else:
		verbose=None
	db=param_db.load(caltable,verbose=verbose)
	gb_keys=db.names()
	gain=db[gb_keys[0]]
	grids=gain.grid_index
	gain=db.get(gb_keys[0]).get_cube()
	gain_err=db.get(gb_keys[-1]).get_cube()
	gain_data=gain.data
	gain_flags=gain.mask
	freq_list=np.array(list(grids['freq'].keys()))
	fluxscale_poly=np.poly1d(np.load(datadir+'/flux_scale_polyfit.npy',allow_pickle=True)[0])
	fluxscale=fluxscale_poly(freq_list)
	if atten!=10: # Valid for 10 and 14 dB, not tested for any other values and typically other values not used for solar observations
		a=1/92.0
		b_ref=(np.sqrt(0.04)/np.sqrt(a*fluxscale))
	else:
		a=1/44.0
		b_ref=(np.sqrt(0.1)/np.sqrt(a*fluxscale))
	b_ref=b_ref.flatten()
	for i in range(gain.shape[1]):
		gain[:,i,:]=gain[:,i,:]*b_ref[i]
	gain_fluxcal=np.ma.array(gain_data,mask=gain_flags,fill_value=(1),dtype=gain_data.dtype)
	# Creating new caltable
	if os.path.exists(caltable):
		os.system('rm -rf '+caltable+'*')
	db_fluxcal=param_db.create(caltable,verbose=verbose)
	for gname in gb_keys:
		db_fluxcal.define_param(gname,db[gname].dtype,db[gname].axis_labels,interpolation_axes=["time", "freq"])
		if 'err' not in gname:
			db_fluxcal.add_chunk(gname,gain_fluxcal,grid=db[gname].grid_index)
		else:
			db_fluxcal.add_chunk(gname,gain_err,grid=db[gname].grid_index)
	db_fluxcal.close()
	return caltable

def fluxcal_bptable(caltable,metafits,atten=10):
	'''
	Function to scale normalized bandpass table
	Parameters
	----------
	caltable : str
		Name of the caltable
	atten : int
		Attenuation in dB
	Returns
	-------
	strget_wsclean_path
		Flux calibrated caltable
	'''
	tb=table()
	tb.open(caltable+'/SPECTRAL_WINDOW')
	freqlist=tb.getcol('CHAN_FREQ')/10**6 #In MHz
	tb.close()
	meanfreq=np.nanmean(np.array(freqlist))
	fluxscale_poly=np.poly1d(np.load(datadir+'/flux_scale_polyfit.npy',allow_pickle=True)[0])
	fluxscale=np.array(fluxscale_poly(meanfreq))
	# Taking care of beam normalization applied in Kansabanik et al. 2022, robust fux calibration.
	###########################
	metaheader=fits.getheader(metafits)
	pointing_ra=metaheader['RA']
	pointing_dec=metaheader['DEC']
	beam=pb.get_pb_radec(pointing_ra,pointing_dec,meanfreq,metafits=metafits)
	I_beam=beam[2]
	fluxscale=fluxscale/(2*I_beam)
	fluxscale=np.array([fluxscale]*len(freqlist))
	#############################
	if atten!=10: # Valid for 10 and 14 dB, not tested for any other values and typically other values not used for solar observations
		a=1/92.0
		b_ref=(np.sqrt(0.04)/np.sqrt(a*fluxscale))
	else:
		a=1/44.0
		b_ref=(np.sqrt(0.1)/np.sqrt(a*fluxscale))	
	b_ref=b_ref.flatten()
	tb=table()
	tb.open(caltable,nomodify=False)
	gain=tb.getcol('CPARAM')
	for i in range(gain.shape[1]):
		gain[:,i,:]=gain[:,i,:]*b_ref[i]
	tb.putcol('CPARAM',gain)
	tb.flush()
	tb.close()
	return caltable

def calc_flag_fraction_caltable(caltable):
	'''
	Calculate flaged fraction from caltable
	
	Parameters
	----------
	caltable : str 
		Name of the CASA caltable	
	Returns
	-------
	float
		Fraction of the total solutions are flagged
	'''
	tb=table()
	tb.open(caltable)
	flag=tb.getcol('FLAG')
	tb.close()
	flagged_fraction=np.sum(flag)/float(flag.size)
	return flagged_fraction

def calc_flag_fraction_cubical_caltable(caltable,verbose=False):
	'''
	Calculate flaged fraction from Cubical caltable
	
	Parameters
	----------
	caltable : str 
		Name of the CubiCal caltable	
	verbose : bool
		Verbose output
	Returns
	-------
	float
		Fraction of the total solutions are flagged
	'''
	if verbose:
		verbose=0
	else:
		verbose=None
	if caltable[-1]=='/':
		caltable=caltable[:-1]
	skel_file=os.path.exists(caltable+'skel')
	db=param_db.load(caltable,verbose=verbose)
	gb_keys=db.names()
	gain=db[gb_keys[0]]
	gain=db.get(gb_keys[0]).get_cube()
	gain_flags=gain.mask
	flagged_fraction=np.nansum(gain_flags)/np.size(gain_flags.flatten())
	if skel_file==False and os.path.exists(caltable+'skel'):
		os.system('rm -rf '+caltable+'skel')
	return flagged_fraction

def calc_flag_chans_caltable(caltable,flag_frac=1.0):
	'''
	Calculate flaged channels from caltable

	Parameters
	----------
	caltable : str 
		Name of the CASA caltable
	flag_frac : float 
		Minimum fraction of data flagged for a single channel	
	Returns
	-------
	list
		Flagged channels list
	float
		flag fraction
	'''
	#TODO: As of now only based on CASA caltable, later CALIBRATE caltables also included
	
	tb=table()
	tb.open(caltable)
	flag=tb.getcol('FLAG')
	tb.close()
	y=(flag==True)
	flag_ants=np.nansum(np.nansum(y,axis=-1),axis=0).astype('float')
	flag1=flag+True
	x=(flag1==True)
	x_tot=np.nansum(np.nansum(x,axis=-1),axis=0).astype('float')
	frac=flag_ants/x_tot
	flagged_chans=np.where(frac>=flag_frac)[0].tolist()
	return flagged_chans,frac

def create_mwapb_caltable(msname,metafits,caltable,verbose=False):
	'''
	Function to create CubiCal caltable for MWA PB at the phase center
	Parameters
	----------
	msname : str
		Name of the measurement set
	metafits : str
		Name of the metafits file
	caltable : str
		Name of the output caltable
	verbose : bool
		Verbose output
	Returns
	-------
	str
		Caltable name
	'''
	AM=am.AccessMS(msname)
	radec_str,radeg,decdeg=AM.get_phasecenter()
	freq=AM.calc_meanfreq()/10**6
	nant=AM.get_num_antenna()
	result=pb.get_pb_radec(radeg,decdeg,freq,metafits=metafits)
	if result[0]!=0:
		print ('Error in calculating primary beam.\n')
		return
	jones_array=result[1]
	# Creating new caltable
	if os.path.exists(caltable):
		os.system('rm -rf '+caltable+' '+caltable+'.skel')
	db_pb=param_db.create(caltable,verbose=verbose)
	gb_keys=['G:gain','G:gain.err']
	pb_gain=jones_array[np.newaxis,np.newaxis,np.newaxis,np.newaxis,:]
	pb_gain=np.repeat(pb_gain,nant,3)
	print (pb_gain.shape)
	pb_gain_err=pb_gain/(10**10)
	pb_gain_flag=pb_gain.astype('bool')
	pb_gain=np.ma.array(pb_gain,mask=pb_gain_flag,fill_value=(1),dtype=pb_gain_flag.dtype)
	pb_gain_err=np.ma.array(pb_gain_err,mask=pb_gain_flag,fill_value=(1),dtype=pb_gain_flag.dtype)
	axis_labels=('dir', 'time', 'freq', 'ant', 'corr1', 'corr2')
	for gname in gb_keys:
		db_pb.define_param(gname,pb_gain.dtype,axis_labels,interpolation_axes=["time", "freq"])
		if 'err' not in gname:
			db_pb.add_chunk(gname,pb_gain,grid=db[gname].grid_index)
		else:
			db_pb.add_chunk(gname,gain_err,grid=db[gname].grid_index)
	db_norm.close()


def smooth_bandpass_edges(bptables=[]):
	'''
	Function to smooth the coarse channel bandpass edges
	Parameters
	----------
	bptables : list
		List of bandpass tables
	'''
	warnings.filterwarnings("ignore")
	tb1=table()
	tb2=table()
	coarse_chan_dic={}
	for bp in bptables:	
		coarse_chans=get_MWA_caltable_coarse_chan(bp)
		coarse_chan_dic[coarse_chans[0]]=[bp,coarse_chans[0],coarse_chans[-1]]
	coarse_chan_dic=dict(sorted(coarse_chan_dic.items()))
	final_list=list(coarse_chan_dic.values())
	for i in range(len(final_list)):
		if i<len(final_list)-1:
			caltable1=final_list[i][0]
			coarse_chan_1=final_list[i][-1]
			caltable2=final_list[i+1][0]
			coarse_chan_2=final_list[i+1][1]
			if coarse_chan_2-coarse_chan_1==1:
				print ('Edge smoothing for : '+caltable1+' and '+caltable2+'\n')
				tb1.open(caltable1,nomodify=False)
				gain1=tb1.getcol('CPARAM')
				gain1_copy=copy.deepcopy(gain1)
				flag=tb1.getcol('FLAG')
				tb2.open(caltable2,nomodify=False)
				gain2=tb2.getcol('CPARAM')
				gain2_copy=copy.deepcopy(gain2)
				bp_1=np.mean(np.mean(gain1,axis=0),axis=-1)
				bp_2=np.mean(np.mean(gain2,axis=0),axis=-1)
				chan_flag=flag[0,:,0]
				bp_1[chan_flag]=np.nan
				bp_2[chan_flag]=np.nan
				pos1=np.where(np.isnan(bp_1))
				pos2=np.where(np.isnan(bp_2))
				bp_1=np.delete(bp_1,pos1)
				bp_2=np.delete(bp_2,pos2)
				bp_end_1=bp_1[-1]
				bp_start_2=bp_2[0]
				bp_mult=bp_end_1/bp_start_2
				gain2_copy*=bp_mult
				tb1.putcol('CPARAM',gain1_copy)
				tb1.flush()
				tb1.close()
				tb2.putcol('CPARAM',gain2_copy)
				tb2.flush()
				tb2.close()
	return coarse_chan_dic


##############################################
# General usuage functions #
##############################################

def seconds_to_timestamp(seconds):
	'''
	Function to convert seconds to hh:mm:ss timestamp
	Parameters
	----------
	seconds : float
		Seconds to convert in hh:mm:ss
	Returns
	-------
	str
		Timestamp string
	'''
	a=dtt.timedelta(seconds=seconds)
	return str(a)

def get_wsclean_path():	
	'''
	Function to get WSClean path
	'''
	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:
			print('WSClean is not installed.\n')
			os.system('rm -rf wsclean_test')
			return 
		else:
			wsclean=wsclean_path+'/wsclean'
			os.system('rm -rf wsclean_test')
			return wsclean
	except:
		print('WSClean is not installed.\n')
		os.system('rm -rf wsclean_test')
		return 

def getnearpos(array,value):
	'''
	Function to return index of two elements nearest to a given number
	
	Parameters
	----------
	array : numpy.array 
		Input numpy array
	value : float 
		Value to which nearest element search
	Returns
	-------
	float
		Index of two elements nearest to the value
	'''
	a = abs(array-value)
	b=np.argsort(a)
	del a
	if len(b)>1:
		return b[0],b[1]
	else:
		return b[0],b[0]

def error_msgs(err_code):
	'''
	Function to transform error code to error message
		
	Parameters
	----------
	err_code : int
		Error code
	Returns
	-------
	str
		Error message
	'''
	if err_code==1:
		return "Split problem"
	elif err_code==2:
		return "Dirty image is not produced"
	elif err_code==3:
		return "Failed to make image in a selfcal iteration"
	elif err_code==4:
		return "Model image is not present"
	elif err_code==5:
		return "No flux is CLEANed in the model"
	elif err_code==6:
		return "Start sigma is too low to start"
	elif err_code==7:
		return "No good solution found"
	elif err_code==8:
		return "Dynamic range decreasing"
	elif err_code==9:
		return "Maximum selfcal iteration is reached"
	elif err_code==10:
		return 'SNR is not sufficient for selfcal'
	elif err_code==11:
		return 'No flux above for the present sigma threshold for starting selfcal'
	elif err_code==12:
		return 'No start sigma and threshold information found for reference time channel'
	elif err_code==13:
		return 'Maximum selfcal iteration is reached and final image has DR less than minimum DR'
	elif err_code==14:
		return 'Invalid stokes axes'		
	elif err_code==14:
		return 'Error in applying solutions'
	elif err_code==100:
		return "Reference time and frequency imaging"
	else:
		return "Succeeded."

def radec_con_deg_to_hhmmss(radeg,decdeg):
	'''
	Convert RA-DEC from degree to hh mm ss dd mm ss format

	Parameters
	----------
	radeg : float 
		Value of the RA in degree
	decdeg : float 
		Value of the DEC in degree
	Returns
	-------
	numpy.array
		['RA','DEC'] in hh mm ss dd mm ss format
	'''
	radeg_copy=copy.deepcopy(radeg)
	decdeg_copy=copy.deepcopy(decdeg)
	radeg=abs(radeg)
	decdeg=abs(decdeg)
	ra_sign=int(radeg_copy/radeg)
	dec_sign=int(decdeg_copy/decdeg)
	if ra_sign<0:
		radeg=copy.deepcopy(radeg_copy)+360.0
	rahh=radeg/15.0
	ramm=(rahh-int(rahh))*60.0
	rass=(ramm-int(ramm))*60.0
	ra=str(int(rahh))+'h'+str(int(ramm))+'m'+str('%.2f'%rass)+'s'
	decdd=decdeg
	decmm=(decdd-int(decdd))*60.0
	decss=(decmm-int(decmm))*60.0
	dec=str(int(dec_sign*decdd))+'d'+str(abs(int(decmm)))+'m'+str('%.2f'%abs(decss))+'s'
	return np.array([ra,dec])

def mjdsec_to_timestamp(mjdsec,includedate=True,format=0):
	'''
	Convert CASA MJD seceonds to CASA timestamp

	Parameters
	----------
	mjdsec : float 
		CASA MJD seconds
	includedate : bool 
		Include date in timestamp
	format : int 
		Datetime string format 
			0: 'YYYY/MM/DD/hh:mm:ss'

			1: 'YYYY-MM-DDThh:mm:ss'

			2: 'YYYY-MM-DD hh:mm:ss' 

			3: 'YYYY_MM_DD_hh_mm_ss'
	Returns
	-------
	str
		CASA time stamp in UTC at the format
	'''
	#TODO : Add more formats
	me=measures()
	qa=quanta()
	today = me.epoch('utc','today')
	mjd = np.array(mjdsec)/86400.0
	today['m0']['value'] =  mjd
	hhmmss = qa.time(today['m0'],prec=8)[0]
	date = qa.splitdate(today['m0'])
	qa.done()
	if (includedate):
		if format==0:
			utcstring = "%s/%02d/%02d/%s" % (date['year'],date['month'],date['monthday'],hhmmss)
		elif format==1:
			utcstring = "%s-%02d-%02dT%s" % (date['year'],date['month'],date['monthday'],hhmmss)
		elif format==2:
			utcstring = "%s-%02d-%02d %s" % (date['year'],date['month'],date['monthday'],hhmmss)
		elif format==3:
			utcstring = "%s_%02d_%02d" % (date['year'],date['month'],date['monthday'])+'_'+'_'.join(hhmmss.split(':'))
	else:
		utcstring = hhmmss
	return utcstring

def timestamp_to_mjdsec(timestamp,format=0):
	'''
	Convert timestamp to mjd second

	Parameters
	----------
	timestamp : str 
		Time stamp to convert 
	format : int 
		Datetime string format 
			0: 'YYYY/MM/DD/hh:mm:ss'
 
			1: 'YYYY-MM-DDThh:mm:ss'

			2: 'YYYY-MM-DD hh:mm:ss' 

			3: 'YYYY_MM_DD_hh_mm_ss'	
	Returns
	-------
	float
		Return correspondong MJD second of the day
	'''
	if format==0:
		try:
			timestamp_datetime=datetime.strptime(timestamp,'%Y/%m/%d/%H:%M:%S.%f')
		except:
			timestamp_datetime=datetime.strptime(timestamp,'%Y/%m/%d/%H:%M:%S')
	elif format==1:
		try:
			timestamp_datetime=datetime.strptime(timestamp,'%Y-%m-%dT%H:%M:%S.%f')
		except:
			timestamp_datetime=datetime.strptime(timestamp,'%Y-%m-%dT%H:%M:%S')
	elif format==2:
		try:
			timestamp_datetime=datetime.strptime(timestamp,'%Y-%m-%d %H:%M:%S.%f')
		except:
			timestamp_datetime=datetime.strptime(timestamp,'%Y-%m-%d %H:%M:%S')
	elif format==3:
		try:
			timestamp_datetime=datetime.strptime(timestamp,'%Y_%m_%d_%H_%M_%S.%f')
		except:
			timestamp_datetime=datetime.strptime(timestamp,'%Y_%m_%d_%H_%M_%S')
	mjd=float("{:.2f}".format((julian.to_jd(timestamp_datetime)-2400000.5)*(24.*3600.)))
	return mjd

def radec_to_altaz(ra,dec,obstime,LAT,LON,ALT):
	'''
	Function to convert radec to altaz for a given Earth location

	Parameters
	----------
	ra : str 
		RA either in degree or 'hh:mm:ss' or '%fh%fm%fs' format
	dec : str 
		DEC either in degree or 'dd:mm:ss' or '%fd%fm%fs'format
	obstime : str 
		Time of the observation in 'yyyy-mm-dd hh:mm:ss' format
	LAT : float
		Latitude of the Earth location in degree
	LON : float 
		Longitude of Earth location in degree 
	ALT : float 
		Altitude of the Earth location in meter
	Returns
	-------
	float
		Altitude in degree
	float
		Azimuth in degree
	'''
	LOCATION=EarthLocation.from_geodetic(lat=LAT*u.deg,lon=LON*u.deg,height=ALT*u.m)
	observing_time=Time(obstime)  
	aa=AltAz(location=LOCATION,obstime=observing_time)
	try:
		ra=float(ra)
		dec=float(dec)
		coord=SkyCoord(ra,dec,frame='icrs',unit='deg')
	except:
		try:
			coord=SkyCoord(ra,dec)
		except:
			coord=SkyCoord(ra,dec,unit=(u.hourangle,u.deg))
	altaz_object=coord.transform_to(aa)
	alt=altaz_object.alt.degree
	az=altaz_object.az.degree
	return alt,az

def altaz_to_radec(alt,az,obstime,LAT,LON,ALT):
	'''
	Function to convert altaz to radec for a given Earth location

	Parameters
	----------
	alt : float 
		Elevation in degree
	az : float 
		Azimuth in degree
	obstime : str 
		Time of the observation in 'yyyy-mm-dd hh:mm:ss' format
	LAT : float 
		Latitude of the Earth location in degree
	LON : float 
		Longitude of Earth location in degree 
	ALT : float 
		Altitude of the Earth location in meter
	Returns
	-------
	numpy.array
		['RA','DEC'] in hh mm ss dd mm ss format
	'''
	LOCATION=EarthLocation.from_geodetic(lat=LAT*u.deg,lon=LON*u.deg,height=ALT*u.m)
	observing_time=Time(obstime) 
	AltAz=SkyCoord(alt=alt*u.deg,az=az*u.deg,obstime=observing_time,frame='altaz',location=LOCATION)
	radec=AltAz.icrs
	radeg=radec.ra.deg
	decdeg=radec.dec.deg
	if radeg<0:
		radeg=360+radeg
	radec_str=radec_con_deg_to_hhmmss(radeg,decdeg)
	return radec_str

def freq_to_MWA_coarse(freq):
	'''
	Frequency to MWA coarse channel conversion

	Parameters
	----------
	freq : float 
		Frequency in MHz
	Returns
	-------
	int
		MWA coarse channel number
	'''
	freq=float(freq)
	coarse_chans=[[(i*1.28)-0.64,(i*1.28)+0.64] for i in range(300)]
	for i in range(len(coarse_chans)):
		ch0=round(coarse_chans[i][0],2)
		ch1=round(coarse_chans[i][1],2)
		if freq>=ch0 and freq<ch1:
			return i 

def update_mwa_obsids(obsid_file='',verbose=False):
	'''
	Function to update MWA OBSIDs 

	Parameters
	----------
	obsid_file : str 
		Name of the file to save MWA OBSIDs
	verbose : bool
		Verbose output
	Returns
	-------
	str
		OBSID file name
	int
		Update success or failure code (0 or 1)
	'''
	if verbose==True:
		print ('Updating local MWA OBSid file......\n')
	if obsid_file=='':
		obsid_file=datadir+'/MWA_OBSids'
	BASEURL='http://ws.mwatelescope.org/'
	temp_array=np.empty(0,dtype='int')
	if os.path.isfile(obsid_file+'.npy')==True:
		try:
			obsids=np.load(obsid_file+'.npy',allow_pickle=True)
			start_obsid=np.max(obsids)
			temp_array=np.append(temp_array,obsids)
		except:
			start_obsid=972654120
	else:
		start_obsid=972654120
	cur_time=dtt.datetime.now()
	cur_time=str(cur_time).split(' ')[0]+'T'+str(cur_time).split(' ')[-1]
	end_obsid_gps=int(Time(cur_time,format='isot',scale='utc').gps) # Till current time
	while True:
		try:
			searchurl=BASEURL+'metadata/find?maxtime='+str(end_obsid_gps)+'&mintime='+str(end_obsid_gps-86400)+'&page=1000'
			end_obsid=json.load(urllib.request.urlopen(searchurl,timeout=1))[-1][0]
			break
		except:
			end_obsid_gps=end_obsid_gps-86400
			if end_obsid_gps<start_obsid:
				break
	try:
		end_obsid=json.load(urllib.request.urlopen(searchurl,timeout=10))[-1][0]
		if verbose==True:
			print ('Last OBSID in MWA metadata server : '+str(end_obsid)+'\n')
		while True:
			searchurl=BASEURL+'metadata/find?mintime='+str(start_obsid)+'&maxtime='+str(start_obsid+432000)
			try:
				OBSid=json.load(urllib.request.urlopen(searchurl,timeout=150))
				OBSid=np.array(OBSid)[:,0].astype('int')
				start_obsid=np.max(OBSid)+235
			except:
				OBSid=np.empty(0,dtype='int')
				start_obsid=start_obsid+3600
			if len(OBSid)!=0:
				temp_array=np.append(temp_array,OBSid)
			if start_obsid>=end_obsid:
				break
		np.save(obsid_file,temp_array)
		if verbose==True:
			print ('Updated successfully.\n')
		os.system('rm -rf casa*log')
		return obsid_file+'.npy',0
	except Exception as e:
		if verbose==True:
			print ('Error in update : '+str(e)+'\n')
			print ('Update not successful.\n')
		os.system('rm -rf casa*log')
		return obsid_file+'.npy',1

def get_OBSID_from_ms(msname):
	'''
	Function to return OBSID of an MWA observation

	Parameters
	----------
	msname : str 
		Name of the measurement set
	Returns
	-------
	int
		MWA Observation ID
	'''
	obsid_file=datadir+'/MWA_OBSids.npy'
	if os.path.isfile(obsid_file)==True:
		obsids=np.load(obsid_file,allow_pickle=True)
	else:
		obsids=np.empty(0)
	if len(obsids)!=0:
		md=msmetadata()
		md.open(msname)
		obs_mjd_ms=md.timerangeforobs(0)['begin']['m0']['value']
		md.close()
		utc_string=mjdsec_to_timestamp(obs_mjd_ms*24*3600,includedate=True,format=1)	
		GPStime=int(Time(utc_string,format='isot',scale='utc').gps)
		OBSid=obsids[np.argmin(abs(GPStime-obsids))]
		os.system('rm -rf casa*log')
		return OBSid
	else:
		c=0
		while c<5:
			try:
				BASEURL='http://ws.mwatelescope.org/'				
				ms_path=os.path.dirname(os.path.realpath(msname))
				searchurl=BASEURL+'metadata/tconv/?utciso='+utc_string
				GPStime=json.load(urllib.request.urlopen(searchurl,timeout=10))
				searchurl=BASEURL+'metadata/find?maxtime='+str(GPStime)+'&mintime='+str(GPStime-500)+'&page=20000000000000'
				OBSid=json.load(urllib.request.urlopen(searchurl,timeout=15))[-1][0]
				os.system('rm -rf casa*log')
				return OBSid
			except Exception as e:
				print ('Error occured : '+str(e)+'\n')
				c+=1
				time.sleep(5.0)
				if c>=5:
					os.system('rm -rf casa*log')
					return 0

def get_OBSID_from_metafits(metafits):
	'''
	Function to return OBSID of an MWA observation

	Parameters
	----------
	metafits : str 
		Name of the metafits file
	Returns
	-------
	int
		MWA Observation ID
	'''
	OBSid=fits.getheader(metafits)['GPSTIME']
	return OBSid 

def download_metafits(msname,outdir):
	'''
	Function to download MWA metafits of a given measurement set (Require internet connection)

	Parameters
	----------
	msname : str 
		Name of the measurement set
	outdir : str
		Name of the outputdir
	Returns
	-------
	str
		Metafits file name
	'''
	BASEURL='http://ws.mwatelescope.org/'
	OBSid=get_OBSID_from_ms(msname)
	if os.path.exists(outdir)==False:
		os.makedirs(outdir)
	if os.path.isfile(outdir+'/'+str(OBSid)+'.metafits')==False:
		try:
			os.system('wget -O '+outdir+'/'+str(OBSid)+'.metafits http://ws.mwatelescope.org/metadata/fits?obs_id='+str(OBSid))
			metafits=outdir+'/'+str(OBSid)+'.metafits'
		except:
			metafits=None
	else:
		metafits=outdir+'/'+str(OBSid)+'.metafits'
	return metafits

def compress_files(filelist,outputfile):
	'''
	Compress a list of numpy files

	Parameters
	----------
	filelist : list 
		List of numpy table files
	outputfile : str 
		Output compressed file name
	Returns
	-------
	str
		Compressed file name (Compressed file have arrays in format [original_filename,data_array])
	'''
	file_array_list=[]
	file_name_list=[]
	for i in filelist:
		a=np.load(i,allow_pickle=True)
		a=np.append(i,a)
		file_array_list.append(a)
	np.savez_compressed(outputfile,a=np.array(file_array_list))
	os.system('mv '+outputfile+'.npz '+outputfile)
	return outputfile

def multifreq_gaincal_interpolate(gaintables=[],output_gaintable='',outputfreq=0):
	'''
	Function to calculate linearly interpolated gain phase from a set of gaintables (at-least two) at multiple frequencies

	Parameters
	----------
	gaintables : list 
		List of gaintables at multiple frequencies
	output_gaintable : str 
		Name of the output gaintable 
	outputfreq : float
		Output frequency in MHz
	Returns
	-------
	str
		Output gaintable name
	'''
	if outputfreq==0:
		print ('Output frequency is not given.\n')
		return
	elif len(gaintables)<2:
		print ('Minimum two gaintables at different frequencies are required.\n')
		return
	else:
		if output_gaintable=='':
			output_gaintable='Gaintable_interp_'+str(outputfreq)+'MHz.gcal'
		if output_gaintable[-1]=='/':
			output_gaintable=output_gaintable[:-1]
		tb=table()
		freqs=[]
		gains=[]
		for g in gaintables:
			tb.open(g+'/SPECTRAL_WINDOW')
			freqs.append(tb.getcol('REF_FREQUENCY')[0])		
			tb.close()
			tb.open(g)
			gains.append(tb.getcol('CPARAM'))
			tb.close()
		gains=np.array(gains)
		if os.path.exists(output_gaintable)==True:
			print ('Removing previous existing caltable :'+output_gaintable+'\n')
			os.system('rm -rf '+output_gaintable)
		os.system('cp -r '+gaintables[0]+' '+output_gaintable)
		tb.open(output_gaintable)
		interpolated_gains=tb.getcol('CPARAM')
		tb.close()
		for i in range(interpolated_gains.shape[0]):
			for j in range(interpolated_gains.shape[-1]):
				x=np.polyfit(freqs,np.angle(gains[:,i,0,j]),deg=1)
				pq=np.poly1d(x,)
				phase=pq(outputfreq*10**6)
				interpolated_gains[i,0,j]=np.cos(phase)+1j*np.sin(phase)
		tb.open(output_gaintable+'/SPECTRAL_WINDOW',nomodify=False)
		chan_freq=tb.getcol('CHAN_FREQ')
		chan_freq=chan_freq+(outputfreq*10**6-chan_freq[0])
		tb.putcol('CHAN_FREQ',chan_freq)
		ref_freq=tb.getcol('REF_FREQUENCY')
		ref_freq[0]=outputfreq*10**6
		tb.putcol('REF_FREQUENCY',ref_freq)
		tb.flush()
		tb.close()
		tb.open(output_gaintable,nomodify=False)
		tb.putcol('CPARAM',interpolated_gains)
		tb.putkeyword('CALTYPE','p')
		tb.flush()
		tb.close()
		return output_gaintable

def get_caltable_metadata(caltable):
	'''
	Function to get caltable metadata

	Parameters
	----------
	caltable : str 
		Name of the caltable
	Returns
	-------
	dict
		A python dictionary with keywords MSNAME, JonesType, Channel 0 frequency (MHz), Central channel frequency (MHz), Channel width (kHz), Bandwidth (MHz), Start time, End time
	'''
	tb=table()
	tb.open(caltable)
	caltype=tb.getkeywords()['VisCal']
	msname=tb.getkeywords()['MSName']
	tb.close()
	tb.open(caltable+'/SPECTRAL_WINDOW')
	ch0=(tb.getcol('REF_FREQUENCY')[0])/10**6 # In MHz
	chanwidth=(tb.getcol('CHAN_WIDTH')[0]/10**3)[0] # In kHz
	freqlist=tb.getcol('CHAN_FREQ')
	chm=(freqlist[int(len(freqlist)/2)]/10**6)[0] # In MHz
	bw=(tb.getcol('TOTAL_BANDWIDTH')[0]/10**6) # In MHz
	tb.close()
	tb.open(caltable)
	timerange= tb.getcol('TIME')
	start_time=mjdsec_to_timestamp(int(np.min(timerange)),includedate=True,format=0)
	end_time=mjdsec_to_timestamp(int(np.max(timerange)),includedate=True,format=0)
	tb.close()
	result={'MSNAME':msname,'JonesType':caltype,'Channel 0 frequency (MHz)':ch0,'Central channel frequency (MHz)':chm,'Channel width (kHz)':chanwidth,'Bandwidth (MHz)':bw,\
			'Start time':start_time,'End time':end_time}
	os.system('rm -rf casa*log')
	return result

def get_MWA_phase(metafits):
	'''
	Function to get MWA phase

	Parameters
	----------
	metafits : str 
		Name of the metafits file
	Returns
	-------
	str
		MWA phase
	'''
	OBSID=get_OBSID_from_metafits(metafits)
	try:
		url='http://ws.mwatelescope.org/metadata/con?obs_id='+str(OBSID)+'&summary'
		config=json.load(urllib.request.urlopen(url,timeout=15))[0]
		if config=='PHASE1':
			config='MWAPhaseI'
		elif config=='LB':
			config='MWAPhaseIILB'
		elif config=='COMPACT':
			config='MWAPhaseIICOMPACT'
	except:
		tilename=fits.getdata(metafits)['TileName']
		LB=['LB' in i for i in tilename]
		Hex=['Hex' in i for i in tilename]
		if np.sum(LB)!=0:
			config='MWAPhaseIILB'
		elif np.sum(Hex)!=0:
			config='MWAPhaseIICOMPACT'
		else:
			config='MWAPhaseI'
	return config

def paircars_instance_runner(cmd,basedir,paircars_basedir,screen_name,finished_touch_file,jobid,prefix_cmds=[]):
	'''
	Function to run a P-AIRCARS instance

	Parameters
	----------
	cmd : str 
		Command to run
	basedir : str
		Base directory of the measurement set
	paircars_basedir : str
		Base directory for a particular P-AIRCARS job
	screen_name : str 
		Name of the screen
	finished_touch_file : str
		Hidden file to generate at the base directory after finishing the instance
	'''
	batch_file=basedir+'/'+screen_name+'.batch'
	cmd_batch=basedir+'/'+screen_name+'_cmd.batch'
	if os.path.isdir(paircars_basedir+'/Logs_and_Errors/Logs')==False:
		os.makedirs(paircars_basedir+'/Logs_and_Errors/Logs')
	if os.path.isdir(paircars_basedir+'/Logs_and_Errors/Errors')==False:
		os.makedirs(paircars_basedir+'/Logs_and_Errors/Errors')
	outputfile=paircars_basedir+'/Logs_and_Errors/Logs/'+screen_name+'.log'
	output_error=paircars_basedir+'/Logs_and_Errors/Errors/'+screen_name+'.error'
	pid_file=paircars_basedir+'/'+str(jobid)+'_pids.log'
	if finished_touch_file!='':
		cmd+=';sleep 2 ;if ! ls '+finished_touch_file+'_* > .tmp_'+screen_name+' 2> '+'.tmp_error_'+screen_name+' < /dev/null &; then  touch '+finished_touch_file+\
		'_error ;  fi; sleep 0.5; rm -rf .tmp_'+screen_name+'; rm -rf .tmp_error_'+screen_name
	else:
		cmd+=';sleep 2'
	cmd='export PYTHONUNBUFFERED=1;echo \"'+cmd+'\" > '+cmd_batch+';sleep 2; chmod a+rwx '+cmd_batch+'; sleep 2; nohup sh '+cmd_batch+' > '+outputfile+' 2>'+output_error+\
		' < /dev/null &; echo $! >> '+pid_file+'; sleep 2'
	cmd=cmd.split(';')
	cmds=[i+'\n' for i in cmd]
	if os.path.exists(cmd_batch):
		os.system('rm -rf '+cmd_batch)
	if os.path.exists(batch_file):
		os.system('rm -rf '+batch_file)
	if os.path.isfile(batch_file):
		fil=open(batch_file,'r+')
	else:
		fil=open(batch_file,'w')
	if len(prefix_cmds)!=0:
		fil.writelines(prefix_cmds)
	fil.writelines(cmds)
	fil.close()
	os.system('chmod a+rwx '+batch_file)
	del cmd,prefix_cmds
	return basedir+'/'+screen_name+'.batch'

def get_available_threads(wait_time=1):
	'''
	Function to get total available threads
	Parameters
	----------
	wait_time : int
		Wait before checking thread availability
	Returns
	-------
	int
		Number of threads
	'''
	available_cpus=int(psutil.cpu_count()*(1-(psutil.cpu_percent(interval=int(wait_time))/100.0)))
	return available_cpus

def check_thread_availability(nthread,wait_time=1):
	'''
	Function to check the thread availability
	Parameters
	----------
	nthread : int
		Number of required free threads
	wait_time : int
		Wait before checking thread availability
	'''
	if nthread>psutil.cpu_count():
		print ('Requested number of threads is more than total available threads in your system. Using half of total threads of your system.\n')
		nthread=int(psutil.cpu_count()/2)
	available_threads=get_available_threads(wait_time=int(wait_time))
	thread_count_availability=0
	while True:	
		if available_threads<nthread:
			if thread_count_availability==0:
				print('Number of available CPU threads is less than '+str(nthread)+'. Waiting for free threads...\n')
			available_threads=get_available_threads(wait_time=int(wait_time))
			thread_count_availability+=1
		else:
			break	

def get_available_paircars_instance(basedir,jobid,total_instances):
	'''
	Function to get available P-AIRCARS instances

	Parameters
	----------
	basedir : str
		Name of the P-AIRCARS base directory
	job_id : int
		P-AIRCARS job ID
	total_instances : int
		Total P-AIRCARS instances
	Returns
	-------
	int
		Available P-AIRCARS instance
	'''
	c=0
	while True:
		if os.path.isfile(basedir+'/'+str(jobid)+'_pids.log')==False:
			running_pids=[]
		else:
			pids=np.loadtxt(basedir+'/'+str(jobid)+'_pids.log',unpack=True).astype('int')
			if pids.shape==():
				pids=np.array([int(pids)])			
			running_pids=[psutil.pid_exists(pid) for pid in pids]
		available_paircars_instance=int(total_instances-np.sum(running_pids))
		if available_paircars_instance<=0:
			time.sleep(5)
			c+=1
			if c>5:
				available_paircars_instance=0
				return available_paircars_instance
			else:
				continue
		else:
			return available_paircars_instance


def run_batch_scripts_cpulist(batch_file,cpu_list=[]):
	'''
	Function to run batch script on a given CPU list
	Parameters
	----------
	batch_file : str
		Name of the batch file
	cpu_list : list
		CPU thread list
	Returns
	-------
	int
		Success message (0 or 1)
	'''
	if len(cpu_list)==0:
		print ('CPU list is empty. Please provide CPU list.\n')
		return 1
	else:
		cpu_str=''
		for i in cpu_list:
			cpu_str+=str(i)+','
		cpu_str=cpu_str[:-1]
		cmd='taskset --cpu-list '+cpu_str+' bash '+batch_file
		print (cmd+'\n')
		os.system(cmd)
		print ('Started batch file : '+batch_file+'\n')
		return 0

def run_batch_scripts(batch_file,num_cpus,cpu_list=[],force_attach=False):
	'''
	Function to run batch script with cpu affinity
	Parameters
	----------
	batch_file : str
		Name of the batch file
	num_cpus : int
		Number of CPU threads to use	
	cpu_list : list
		List of CPU threads to use
	force_attach : bool
		Attach the cpu even if it is not vacant
	'''
	cpu_fracs=np.array(psutil.cpu_percent(percpu=True,interval=5))
	cpus_argsort=np.argsort(cpu_fracs)
	cpu_list_final=[]
	if len(cpu_fracs)<num_cpus:
		num_cpus=len(cpu_fracs)
	print ('Intended CPU threads : '+str(cpu_list)+'\n')
	print ('Total required CPU threads : '+str(num_cpus)+'\n')
	if len(cpu_list)==0:
		print ('Please provide CPU list.\n')
		return 1
	else:
		c=0
		k=0
		while True:
			for cpu in cpus_argsort:	
				frac=cpu_fracs[cpu]
				if cpu in cpu_list and cpu not in cpu_list_final:
					if force_attach or frac<=10:
						cpu_list_final.append(cpu)
						c+=1
				if len(cpu_list_final)>=num_cpus:
					break
			if len(cpu_list_final)>=num_cpus or len(cpu_list_final)>1:
				cpu_str=''
				cpu_list_final=sorted(cpu_list_final)
				for k in cpu_list_final:
					cpu_str+=str(k)+','
				cpu_str=cpu_str[:-1]
				cmd='taskset --cpu-list '+cpu_str+' bash '+batch_file
				print (cmd+'\n')
				os.system(cmd)
				print ('Started batch file : '+batch_file+'\n')
				return 0
			else:
				if k==0:
					print ('Waiting for free CPU threads.....\n')
				cpu_fracs=np.array(psutil.cpu_percent(percpu=True,interval=5))
				cpus_argsort=np.argsort(cpu_fracs)
				k+=1
				time.sleep(0.5)
				continue

def get_free_cpu_list(max_cpu_frac,per_cpu_free_frac=0.1,cpu_list=[]):
	'''
	Function to get free CPU thread list
	Parameters
	----------
	max_cpu_frac : float
		Maximum CPU fraction to use
	cpu_list : list
		List of CPU threads to search for
	Returns
	-------
	list
		CPU thread list
	'''
	total_allowed_cpus=int(psutil.cpu_count()*max_cpu_frac)	
	cpu_list_final=[]
	if len(cpu_list)==0:
		cpu_list=[i for i in range(psutil.cpu_count())]
	c=0
	k=0
	while True:
		cpu_fracs=np.array(psutil.cpu_percent(percpu=True,interval=5))
		cpus_argsort=np.argsort(cpu_fracs)
		for i in cpus_argsort:	
			frac=cpu_fracs[i]
			if frac<=(1-per_cpu_free_frac)*100:
				if i in cpu_list and i not in cpu_list_final:
					cpu_list_final.append(i)
					c+=1
					if len(cpu_list_final)>=total_allowed_cpus:
						break
			if len(cpu_list_final)>=total_allowed_cpus:
				break
		if len(cpu_list_final)>=1:
			break
		else:
			if k==0:
				print ('Waiting for at least 1 cpu.\n')
			k+=1
	return sorted(cpu_list_final)
	
def get_final_psf(msname,weight='briggs',robust=0.5):
	'''
	Function to get a common final psf parameters

	Parameters
	----------
	msname : str
		Name of the measurement set
	weight : str
		Visibility weighting for imaging
	robust : float
		Robust parameter for 'briggs' weighting
	Returns
	-------
	float 
		Major axis (FWHM) in arcsec
	float
		Minor axis (FWHM) in arcsec
	float 
		Position angle in degree
	'''
	imagename=msname.split('.ms')[0]+'_test_image'
	if msname[-1]=='/':
		msname=msname[:-1]
	AM=am.AccessMS(msname)
	timestamps=AM.get_timestamps()
	chans=AM.get_unflag_chan()
	if os.path.exists(msname+'.temp'):
		os.system('rm -rf '+msname+'.temp '+msname+'.temp.flagversions')
	if os.path.exists(msname+'.temp.flagversions'):
		os.system('rm -rf '+msname+'.temp.flagversions')
	split(vis=msname,outputvis=msname+'.temp',timerange=timestamps[0],datacolumn='data',spw='0:'+str(chans[0]))
	flagdata(vis=msname+'.temp',mode='unflag',uvrange='>0lambda',flagbackup=False)
	IB=ImageBasic(msname+'.temp',includeflag=True)
	cell=IB.calc_cellsize(3) # Assuming 3 pixels in one PSF
	imsize=int(IB.num_pixels(3)/4)
	if imsize>100:
		imsize=100
	scales=IB.choose_scales(3,32*60)
	im=imager()
	im.open(msname+'.temp')
	im.defineimage(nx=imsize,ny=imsize,cellx=str(cell)+'arcsec',celly=str(cell)+'arcsec',stokes='I',mode='mfs')
	im.weight(type=weight,robust=robust)
	beam=im.fitpsf(psf='')
	major=beam[1]['value']
	minor=beam[2]['value']
	pa=beam[3]['value']
	os.system('rm -rf *.psf '+msname+'.temp')
	return major,minor,pa

def check_internet(url=''):
	'''
	Function to check internet connection
	Parameters
	----------
	url : str
		An URL to check internet connectivity. (default : https://www.google.com)
	Returns
	-------
	bool
		Whether internet connection is available or not
	'''
	url='https://www.google.com'
	try_count=5
	timeout=5
	c=0
	while c<=try_count:
		try:
			requests.get(url,timeout=timeout)
			return True
		except:
			continue
		c+=1
	os.system('casa*log')
	return False 

def get_MWA_caltable_coarse_chan(caltable):
	'''
	Function to get the MWA coarse channel number for the caltable
	Parameters
	----------
	caltable : str
		Name of the caltable
	Returns
	-------
	int
		Coarse channel list
	'''
	cal_metadata=get_caltable_metadata(caltable)
	ch0_freq=cal_metadata['Channel 0 frequency (MHz)']
	coarse_chan0=freq_to_MWA_coarse(ch0_freq)
	ncoarse=int(cal_metadata['Bandwidth (MHz)']/1.28)
	coarse_chan=[]
	if ncoarse<1:
		ncoarse=1
	for i in range(ncoarse):
		coarse_chan.append(coarse_chan0+i)
	return coarse_chan

def get_MWA_cubical_caltable_coarse_chan(caltable,verbose=False):
	'''
	Function to get the MWA coarse channel number for the CubiCal caltable
	Parameters
	----------
	caltable : str
		Name of the caltable
	verbose : bool
		Verbose output
	Returns
	-------
	int
		Coarse channel list
	'''	
	try:
		c=0
		while os.path.exists(caltable)==False:
			c+=1			
			time.sleep(1)
			if c>5:
				print (caltable+' does not exist. Please try later.\n')
				return []
		db=param_db.load(caltable)
		gb_keys=db.names()
		gain=db[gb_keys[0]]
		grids=gain.grid_index
		gain=db.get(gb_keys[0]).get_cube()
		gain_err=db.get(gb_keys[-1]).get_cube()
		gain_data=gain.data
		gain_flags=gain.mask
		freq_list=np.array(list(grids['freq'].keys()))/10**6
		ch0_freq=freq_list[0]
		coarse_chan0=freq_to_MWA_coarse(ch0_freq)
		ncoarse=int(abs(freq_list[-1]-freq_list[0])/1.28)
		if ncoarse<1:
			ncoarse=1
		coarse_chan=[]
		for i in range(ncoarse):
			coarse_chan.append(coarse_chan0+i)
		return coarse_chan
	except Exception as e:
		print ('Error in calculating coarse channel : '+str(e)+'\n')
		return []

def make_mwa_crossphase_caltable(msname,cross_phase=15,caltable='',polbasis='Linear'): # TODO : Implement circular basis as well
	'''
	Function to make cross hand phase Jones matrices

	Parameters
	----------
	msname : str
		Name of measurement set
	cross_phase : float 
		Cross hand phase different in degree
	caltable : str 
		Name of the caltable to store cross hand Jones matrices
	polbasis : str 
		'Linear' or 'Circular' (Circular basis not implemented)
	Returns
	-------
	str
		Caltable name
	'''
	if caltable=='':
		print ('No caltable name is specified.\n')
		return None
	cross_phase=np.deg2rad(cross_phase)/2.0
	cross_jones=np.array([[np.cos(cross_phase)+1j*np.sin(cross_phase),0],[0,np.cos(cross_phase)-1j*np.sin(cross_phase)]])[np.newaxis,np.newaxis,np.newaxis,np.newaxis,...]
	cross_jones=cross_jones.astype('complex128')
	AM=am.AccessMS(msname)
	nant=AM.get_num_antenna()
	nchan=AM.get_num_channels()
	nint=AM.get_num_timestamps()
	freqs=AM.get_freqs()
	start_freq=freqs[0]
	end_freq=freqs[-1]
	mjdsecs=AM.get_timestamps_in_mjdsecs()[0]
	startmjd=mjdsecs[0]
	endmjd=mjdsecs[-1]
	cross_jones=np.repeat(np.repeat(np.repeat(cross_jones,nint,axis=1),nchan,2),nant,3)
	cross_jones=np.ma.array(cross_jones,mask=False,fill_value=0)
	dummy_err=(cross_jones*0)+(1e-20)
	dummy_err=dummy_err.astype('complex128')
	dummy_err=np.ma.array(dummy_err,mask=False,fill_value=0)
	time_dict=OrderedDict()
	for mjdsec in mjdsecs:
		time_dict[float(mjdsec)]=0
	freq_dict=OrderedDict()
	for freq in freqs:
		freq_dict[float(freq)]=0
	dir_dict=OrderedDict()
	dir_dict[0]=0
	ant_dict=OrderedDict()
	for i in range(nant):
		ant_name=AM.get_antenna_name(antenna_id=i)
		ant_dict[ant_name]=i
	corr_dict=OrderedDict()
	corr_dict['x']=0
	corr_dict['y']=1
	grid_index={'dir':dir_dict,'time':time_dict,'freq':freq_dict,'ant':ant_dict,'corr1':corr_dict,'corr2':corr_dict}
	cwd=os.getcwd()
	cubical_verbose=None
	db_crossphase=param_db.create(caltable,verbose=cubical_verbose)
	gb_keys=['G:gain', 'G:gain.err']
	for gname in gb_keys:
		db_crossphase.define_param(gname,cross_jones.dtype,tuple(['dir','time','freq','ant','corr1','corr2']),interpolation_axes=["time", "freq"])
		if 'err' not in gname:
			db_crossphase.add_chunk(gname,cross_jones,grid=grid_index)
		else:
			db_crossphase.add_chunk(gname,dummy_err,grid=grid_index)
	db_crossphase.close()
	return caltable
	
def apply_cross_hand_phase(msname,cross_phase=15,caltable='',polbasis='Linear',datacolumn='DATA',modify_datacolumn=False,verbose=False):
	'''
	Function to applycal cross hand phase solution

	Parameters
	----------
	msname : str
		Name of the measurement set
	cross_phase : float 
		Cross hand phase different in degree
	caltable : str 
		Name of the caltable to store cross hand Jones matrices
	polbasis : str 
		'Linear' or 'Circular' (Circular basis not implemented)
	datacolumn : str
		Datacolumn to apply the solution ('DATA', or 'CORRECTED_DATA')
	modify_datacolumn : bool 
		Modify the DATA column or not
	Returns
	-------
	str
		Caltable name
	'''	
	if caltable=='':
		caltable=os.getcwd()+'/'+os.path.basename(msname).split('.ms')[0]+'_cross_phase.temp'
	if verbose:
		print('make_mwa_crossphase_caltable(cross_phase='+str(cross_phase)+',caltable=\''+caltable+'\',polbasis=\''+polbasis+'\')\n')
	if os.path.exists('cubical.cc-out'):
		os.system('rm -rf cubical.cc-out* tmp')
	caltable_name=make_mwa_crossphase_caltable(msname,cross_phase=cross_phase,caltable=caltable,polbasis=polbasis)
	cubical_args=['paircars_cubical','--data-ms '+msname,'--out-casa-gaintables 0','--bbc-plot 0','--data-column '+datacolumn,'--g-load-from '+caltable_name,'--out-mode ac',\
					'--out-overwrite 0','--out-backup 0','--flags-save-legacy 1']
	if verbose:
		print(' '.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')
	if modify_datacolumn==True:
		tb=table()
		tb.open(msname,nomodify=False)
		cor_data=tb.getcol('CORRECTED_DATA')
		tb.putcol('DATA',cor_data)
		tb.flush()
		tb.close()
	os.system('rm -rf '+caltable+'*')
	return caltable_name

def generate_tb_map(imagename,outfile):
	'''
	Function to generate brightness temperature map
	Parameters
	----------
	imagename : str
		Name of the flux calibrated image
	outfile : str
		Basename of the output brightess temperature image
	Returns
	-------
	str
		Output image name
	'''
	try:
		if os.path.exists(imagename+'.image'):
			os.system('rm -rf '+imagename+'.image')
		importfits(fitsimage=imagename,imagename=imagename+'.image')
		imagename=imagename+'.image'
		image_header=imhead(imagename=imagename)
		imagetype='fits'
	except:
		try:
			image_header=imhead(imagename=imagename)
			imagetype='CASA'
		except:
			print ('Image is either not in CASA or FITS format.\n')
			return 
	major=image_header['restoringbeam']['major']['value'] # In arcsec
	minor=image_header['restoringbeam']['minor']['value'] # In arcsec
	freq_axis=np.where(image_header['axisunits']=='Hz')[0][0]
	freq=image_header['refval'][freq_axis]/10**9 # In GHz
	beam_axes_multi=major*minor	
	if imagetype=='CASA':
		outfile=outfile+'.'+imagename.split('.')[-1]
	else:
		outfile=outfile+'.fits'
	expr='1.222e6*IM0/'+str(freq)+'^2/('+str(major*minor)+')'
	if os.path.exists(outfile):
		os.system('rm -rf '+outfile)
	immath(imagename=imagename,outfile=outfile,mode='evalexpr',expr=expr)
	if imagetype=='fits':
		os.system('rm -rf '+imagename)
	imhead(imagename=outfile, mode='put', hdkey='bunit', hdvalue='K')
	if imagetype=='fits':
		if os.path.exists(outfile+'.temp'):
			os.system('rm -rf '+outfile+'.temp')
		exportfits(imagename=outfile,fitsimage=outfile+'.temp',stokeslast=False)
		os.system('rm -rf '+outfile)
		os.system('mv '+outfile+'.temp '+outfile)
	os.system('rm -rf casa*log')
	return outfile

def negative_box(max_pix,imsize,cellsize,box_width=3):
	'''
	Create a box about the maximum pixel of image to search negative
	Parameters
	----------
	max_pix : list 
		Maximum pixel [xxmax,yymax]
	imsize : int
		Image size (Number of pixels)
	cellsize : float
		Pixelsize in arcsec
	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)/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>=imsize:
		xrtc=imsize-1
	if yrtc>=imsize:
		yrtc=imsize-1
	os.system('rm -rf casa*log')
	return str(int(xblc))+','+str(int(yblc))+','+str(int(xrtc))+','+str(int(yrtc))

def get_image_rms(imagename,stokes='I'):
	'''
	Function to calculate rms of the image
	Parameters	
	----------
	imagename : str
		Name of the image
	stokes : str
		Stokes plane
	Returns
	-------
	float
		Mean of the background sky
	float
		RMS of the image
	'''
	if imagename[-1]=='/':
		imagename=imagename[:-1]
	org_imagename=imagename
	os.system('cp -r '+imagename+' '+imagename+'.temp')
	imagename=imagename+'.temp'
	imageheader=imhead(imagename=imagename,mode='summary')
	imsize=imageheader['shape'][0]
	cellsize=np.rad2deg(imageheader['incr'][1])*3600
	rms_box='50,50,'+str(imsize-50)+','+str(int(imsize/4)) 
	try:
		rms=imstat(imagename=imagename,box=rms_box,stokes=stokes)['rms'][0]
	except:
		os.system('rm -rf '+imagename)
		imcollapse(imagename=org_imagename,outfile=imagename,overwrite=True,axes=2,function='mean')
		stokes=''
		rms=imstat(imagename=imagename,box=rms_box,stokes=stokes)['rms'][0]
	image_mean=imstat(imagename=imagename,box=rms_box,stokes=stokes)['mean'][0]
	os.system('rm -rf '+imagename+' casa*log')
	return image_mean,rms

def determine_quiet_sun_condition(msname,suntaper=200,sigma=3,spw=''):
	'''
	This function determine the presence of non-quiet component from the visibility data.
	This function assumes if the Sun is completely Quiet over the timerange, then beyond the 200 lambda (at metre-wavelength or any user given value at the observing frequency),
	there will be no strong signal in the visibility. It calculates the S/N of the 200-300 lambda range, and provides the Quiet condition based on the signal strength.
	Parameters
	----------
	msname : str
		Name of the measurement set
	suntaper : floar
		Tapering of the quiet Sun at uv space (in lambda)
	sigma : float
		SNR value for classify the quiet Condition
	spw : str
		CASA format spectral window
	Returns
	------
	dict
		A python dictionary with timestamp and corresponding quiet Sun condition {'timestamp1': [SNR,True/False],....}
	'''
	if os.path.islink(msname):
		msname=os.path.realpath(msname)
	if os.path.isdir(msname)==False:
		print ('Measurment set : '+msname+' does not exist.\n')
		return {}
	AM=am.AccessMS(msname)
	timestamps=AM.get_timestamps()
	result_dict={}
	uvrange=str(suntaper)+'~'+str(suntaper+100)+'lambda'
	datacolumn='corrected'
	for t in timestamps:
		try:
			print ('visstat(vis=\''+msname+'\',datacolumn=\''+datacolumn+'\',axis=\'amp\',uvrange=\''+str(uvrange)+'\',timerange=\''+t+'\')\n')
			stat_vis=visstat(vis=msname,datacolumn=datacolumn,axis='amp',uvrange=uvrange,timerange=t)
		except:
			if datacolumn=='corrected':
				datacolumn='data'
			elif datacolumn=='data':
				datacolumn='corrected'
			print ('visstat(vis=\''+msname+'\',datacolumn=\''+datacolumn+'\',axis=\'amp\',uvrange=\''+str(uvrange)+'\',timerange=\''+t+'\')\n')
			stat_vis=visstat(vis=msname,datacolumn=datacolumn,axis='amp',uvrange=uvrange,timerange=t)
		stat_vis=stat_vis['DATA_DESC_ID=0']
		mean_amp=stat_vis['mean']
		std_amp=stat_vis['stddev']
		snr=(mean_amp/std_amp)
		if round(mean_amp,0)>10:
			result_dict[t]=[mean_amp,snr,False]
		else:
			result_dict[t]=[mean_amp,snr,True]
	return result_dict			

def get_dynamic_range(imagename,stokes='I'):
	'''
	Function to calculate dynamic range of the image
	Parameters	
	----------
	imagename : str
		Name of the image
	stokes : str
		Stokes plane
	Returns
	-------
	float
		RMS based dynamic range
	float
		Negative based dynamic range
	'''
	if imagename[-1]=='/':
		imagename=imagename[:-1]
	org_imagename=imagename
	os.system('cp -r '+imagename+' '+imagename+'.temp')
	imagename=imagename+'.temp'
	imageheader=imhead(imagename=imagename,mode='summary')
	imsize=imageheader['shape'][0]
	cellsize=np.rad2deg(imageheader['incr'][1])*3600
	rms_box='50,50,'+str(imsize-50)+','+str(int(imsize/4)) 
	try:
		max_flux=imstat(imagename=imagename,stokes=stokes)['max'][0]
	except:
		os.system('rm -rf '+imagename)
		imcollapse(imagename=org_imagename,outfile=imagename,overwrite=True,axes=2,function='mean')
		stokes=''
		max_flux=imstat(imagename=imagename,stokes=stokes)['max'][0]
	rms=imstat(imagename=imagename,box=rms_box,stokes=stokes)['rms'][0]
	max_pix=imstat(imagename=imagename,stokes=stokes)['maxpos']
	neg_box=negative_box(max_pix,imsize,cellsize,box_width=3)
	min_flux=imstat(imagename=imagename,stokes=stokes,box=neg_box)['min'][0]
	rms_dr=max_flux/rms
	if min_flux!=0:
		neg_dr=max_flux/abs(min_flux)
	else:
		neg_dr=rms_dr
	os.system('rm -rf '+imagename+' casa*log')
	return rms_dr,neg_dr

def bookkeeping_DR(imagename,timestamp,freq,freq_width,time_width,bookkeeping_file):
	'''
	Function for bookkeeping DR of the images 
	Note : All DR is scaled to 0.5s and 40 kHz frequency resolution 
	Parameters
	----------
	imagename : str
		Name of the image
	timestamp : str
		Timestamp
	freq : float
		Frequency in MHz
	freq_width : float
		Frequency width in kHz
	time_width : float
		Temporal width in seconds
	bookkeeping_file : str
		Bookkeeping file name
	Returns
	-------
	str	
		Bookkeeping file name
	'''
	if os.path.exists(bookkeeping_file)==False:
		os.system('touch '+bookkeeping_file)
	rms_dr,neg_dr=get_dynamic_range(imagename,stokes='I')
	scaling=np.sqrt((0.5/time_width)*(40.0/freq_width))
	rms_dr=rms_dr*scaling
	neg_dr=neg_dr*scaling
	if os.path.exists(bookkeeping_file):
		try:
			data=pd.read_pickle(bookkeeping_file)
			data[(timestamp,freq)]=np.array([rms_dr,neg_dr])
		except:
			data=pd.DataFrame(np.array([rms_dr,neg_dr]),columns=[(timestamp,freq)],index=['RMS_DR','NEG_DR'])
	else:
		data=pd.DataFrame(np.array([rms_dr,neg_dr]),columns=[(timestamp,freq)],index=['RMS_DR','NEG_DR'])
	data.to_pickle(bookkeeping_file)
	return bookkeeping_file
		
def paircars_image_bookkeep(imagename,freq_width,time_width,bookkeeping_file):
	'''
	Function to bookkeep P-AIRCARS generated image (Image name should be of P-AIRCARS format)
	Parameters
	----------
	imagename : str
		Name of the image
	freq_width : float
		Frequency width in kHz
	time_width : float
		Temporal width in seconds
	bookkeeping_file : str
		Name of the bookkeeping file
	Returns
	-------
	str
		Name of the bookkeeping file
	'''
	if os.path.exists(bookkeeping_file)==False:
		os.system('touch '+bookkeeping_file)
	timestamp_list=imagename.split('time_')[-1].split('_freq')[0].split('_')
	timestamp='/'.join(timestamp_list[:3])+'/'+':'.join(timestamp_list[3:])
	freq=imagename.split('_freq_')[-1].split('_')[0].split('.image')[0]
	bookkeeping_file=bookkeeping_DR(imagename,timestamp,freq,freq_width,time_width,bookkeeping_file)
	return bookkeeping_file

def get_closest_DR(bookkeeping_file,freq,timestamp):
	'''
	Function to get closest DR from the bookkeeping file
	Parameters
	----------
	bookkeeping_file: str
		Name of the bookkeeping file
	freq : float
		Frequency in MHz
	timestamp : str
		Timestamp (YYYY/MM/DD/hh:mm:ss.ff)
	Returns
	-------
	float
		RMS based DR
	float
		Negative based DR
	float
		Mean rms DR
	float
		Mean negative DR
	'''
	if os.path.exists(bookkeeping_file):
		bookkept_data=pd.read_pickle(bookkeeping_file)
		numpy_data=bookkept_data.to_numpy()
		median_rms_dr=np.nanmedian(numpy_data[0])
		median_neg_dr=np.nanmedian(numpy_data[1])
		columns=list(bookkept_data.columns)
		org_column=copy.deepcopy(columns)
		for i in range(len(columns)):
			columns[i]=list(columns[i])
		for i in range(len(columns)):
			columns[i][0]=timestamp_to_mjdsec(columns[i][0])
		columns=np.array(columns).astype('float')
		dist=lambda a, b: (a[0]-b[0])**2 + (a[1]-b[1])**2
		closest_coord=min(columns, key=lambda co: dist(co,(timestamp_to_mjdsec(timestamp),freq)))
		pos=None
		for i in range(len(columns)):
			if list(columns[i])==list(closest_coord):
				pos=i
		if pos!=None:
			rms_dr=bookkept_data[org_column[pos]].RMS_DR	
			neg_dr=bookkept_data[org_column[pos]].NEG_DR
		else:
			rms_dr=neg_dr=None
		return rms_dr,neg_dr,median_rms_dr,median_neg_dr
	else:
		return None,None,None,None

def get_minmax_DR(bookkeeping_file):
	'''
	Function to get minimum and maximum DR from the bookkeeping file
	Parameters
	----------
	bookkeeping_file: str
		Name of the bookkeeping file
	Returns
	-------
	float
		Minimum RMS based DR
	float
		Maximum RMS based DR
	float
		Minimum negative based DR
	float
		Maximum negative based DR
	'''
	if os.path.exists(bookkeeping_file):
		bookkept_data=pd.read_pickle(bookkeeping_file)
		numpy_data=bookkept_data.to_numpy()
		rms_dr_array=numpy_data[0]
		neg_dr_array=numpy_data[1]
		min_rms_dr=np.nanmin(rms_dr_array)
		max_rms_dr=np.nanmax(rms_dr_array)
		min_neg_dr=np.nanmin(neg_dr_array)
		max_neg_dr=np.nanmax(neg_dr_array)
		return min_rms_dr,max_rms_dr,min_neg_dr,max_neg_dr
	else:
		return None,None,None,None

def get_quicklook_image(imagename,outfile,freq,timestamp,DR_rms,DR_neg,field_of_view=2,imagetype='casa'):
	'''
	Function to get a quick look image
	Parameters
	----------
	imagename : str 
		Name of the CASA image
	outfile : str 
		Output file name
	freq : float 
		Frequency in MHz
	timestamp : str 
		Timestamp string
	DR_rms : float 
		Dynamic range of the image based on rms
	DR_neg : float 
		Dynamic range of the image based on negative flux near the source
	field_of_view : float 
		Field of view to cut the image in degree (default : 2)
	imagetype : str
		Type of the image, 'casa' or 'fits'
	Returns
	-------
	str
		Outfile name
	'''
	if imagename[-1]=='/':
		imagename=imagename[:-1]
	if imagetype=='fits':
		try:
			print ('importfits(fitsimage=\''+imagename+'\',imagename=\''+imagename+'.image\')\n')
			importfits(fitsimage=imagename,imagename=imagename+'.image')
			temp_image=imagename+'.image'
		except Exception as e:
			print ('Can not convert fits to CASA image. Error occured : '+str(e)+'\n')
			return 
	else:
		temp_image=imagename+'_temp.fits'
		os.system('cp -r '+imagename+' '+temp_image)
	header=imhead(imagename=temp_image,mode='list')
	xcent=int(header['shape'][0]/2)
	ycent=int(header['shape'][1]/2)
	cell=np.rad2deg(abs(header['cdelt2'])) # In degree
	freq="{:.2f}".format(float(freq))
	xwidth=ywidth=int((field_of_view)/cell)
	box=str(xcent-int(xwidth/2))+','+str(ycent-int(ywidth/2))+','+str(xcent+int(xwidth/2))+','+str(ycent+int(ywidth/2))
	os.system('rm -rf temp*')
	try:
		if os.path.exists(temp_image+'_temp.image'):
			os.system('rm -rf '+temp_image+'_temp.image')
		imsubimage(imagename=temp_image,outfile=temp_image+'_temp.image',box=box,stokes='I')
	except Exception as e:
		print ('WARNING : '+str(e)+'. Trying to collapse Stokes axes.\n')
		try:
			if os.path.exists(temp_image+'_temp.image'):
				os.system('rm -rf '+temp_image+'_temp.image')
			if os.path.exists(temp_image+'_temp1.image'):
				os.system('rm -rf '+temp_image+'_temp1.image')
			imcollapse(imagename=temp_image,outfile=temp_image+'_temp1.image',axes=2,function='mean')
			imsubimage(imagename=temp_image+'_temp1.image',outfile=temp_image+'_temp.image',box=box)
		except Exception as e:
			print ('WARNING : '+str(e)+'. Continuing with the full image.\n')
			if os.path.exists(temp_image+'_temp.image'):
				os.system('rm -rf '+temp_image+'_temp.image')
			os.system('cp -r '+temp_image+' '+temp_image+'_temp.image')	
	try:
		if os.path.exists(temp_image+'_temp.fits'):
			os.system('rm -rf '+temp_image+'_temp.fits')
		exportfits(imagename=temp_image+'_temp.image',fitsimage=temp_image+'_temp.fits',dropdeg=True)
		data=fits.getdata(temp_image+'_temp.fits')
		wlist=fits.getheader(temp_image+'_temp.fits')
		w = wcs.WCS(wlist)
		fig = plt.figure(figsize=(8,8))
		ax = fig.add_subplot(111, projection = w)
		cax = fig.add_axes([0.91, 0.11, 0.02, 0.77])
		im=ax.imshow(data,cmap='hot',origin='lower')
		fig.colorbar(im,cax=cax,orientation='vertical')
		ax.set_xlabel('Right Ascension',fontsize=12)
		ax.set_ylabel('Declination',fontsize=12)
		title='Frequency : '+str(freq)+' MHz, Timestamp : '+str(timestamp)+' UTC\n'
		if DR_rms!=0 and DR_neg!=0: 
			title+='Dynamic range (rms) : '+str(int(DR_rms))+', Dynamic range (negative) : '+str(int(DR_neg))
		ax.set_title(title)
		ax.xaxis.set_tick_params(labelsize=10)
		ax.yaxis.set_tick_params(labelsize=10)
		cwd=os.getcwd()
		outfile_dir=os.path.dirname(outfile)
		if outfile_dir=='':
			outfile=cwd+'/'+outfile
		plt.savefig(outfile)
		os.system('rm -rf '+temp_image+'* temp* casa*log')
		return outfile
	except Exception as e:	
		print ('Error occured in producing quicklook image. Error : '+str(e)+'\n')
		os.system('rm -rf '+temp_image+'* temp* casa*log')
		return

def get_quicklook_polimage(imagename,outfile,freq,timestamp,DR_rms,DR_neg,field_of_view=2,imagetype='casa'): 
	'''
	Function to get a quick look image
	Parameters
	----------
	imagename : str 
		Name of the image	
	outfile : str
		Output file name
	freq : float
		Frequency in MHz
	timestamp : str 
		Timestamp string
	DR_rms : float 
		Dynamic range of the image based on rms
	DR_neg : float 
		Dynamic range of the image based on negative flux near the source
	field_of_view : float 
		Field of view to cut the image in degree (default : 2)
	imagetype : str
		Image type, 'casa' or 'fits' 
	Returns
	-------
	str
		Outfile name
	'''
	if imagename[-1]=='/':
		imagename=imagename[:-1]
	if imagetype=='fits':
		try:
			if os.path.exists(imagename+'.image'):
				os.system('rm -rf '+imagename+'.image')
			c=0
			while True:
				if os.path.exists(imagename) and os.path.getsize(imagename)>0:
					print ('importfits(fitsimage=\''+imagename+'\',imagename=\''+imagename+'.image\')\n')
					try:
						importfits(fitsimage=imagename,imagename=imagename+'.image')
						temp_image=imagename+'.image'
					except Exception as e:
						print ('Error in importfits '+str(e)+'\n')
						c+=1
						time.sleep(5)
						pass
					break
				elif c>5:
					os.system('rm -rf '+temp_image+'* temp_*.image')
					print ('Errror in producing quicklook images 0 : '+str(e)+'\n')
					return	
				else:
					c+=1
					time.sleep(5)			
		except Exception as e:
			print ('Can not convert fits to CASA image. Error occured : '+str(e)+'\n')
			return 
	else:
		temp_image=imagename+'_temp.image'
		os.system('cp -r '+imagename+' '+temp_image)
	if os.path.exists(temp_image)==False:
		print ('Image file could not be located.\n')
		os.system('rm -rf '+temp_image+'* temp_*.image casa*log')
		return
	try:
		header=imhead(imagename=temp_image,mode='list')
		xcent=int(header['shape'][0]/2)
		ycent=int(header['shape'][1]/2)
		cell=np.rad2deg(abs(header['cdelt2'])) # In degree
		freq="{:.2f}".format(float(freq))
		xwidth=ywidth=int((field_of_view)/cell)
		box=str(xcent-int(xwidth/2))+','+str(ycent-int(ywidth/2))+','+str(xcent+int(xwidth/2))+','+str(ycent+int(ywidth/2))
		stokes_list=['I','Q','U','V']
		fig = plt.figure(figsize=(8,8))
		plt.subplots_adjust(wspace=0.45, hspace=0.1)
	except Exception as e:
		os.system('rm -rf '+temp_image+'* temp_*.image')
		print ('Errror in producing quicklook images 0: '+str(e)+'\n')
		return
	for i in range(len(stokes_list)):
		stokes=stokes_list[i]
		try:
			if os.path.exists(temp_image+'_temp_'+stokes+'.image'):
				os.system('rm -rf '+temp_image+'_temp_'+stokes+'.image')
			imsubimage(imagename=temp_image,outfile=temp_image+'_temp_'+stokes+'.image',box=box,stokes=stokes)
		except Exception as e:
			os.system('rm -rf '+temp_image+'* temp_*.image')
			print ('Errror in producing quicklook images 1: '+str(e)+'\n')
			return 
		try:
			if os.path.exists(temp_image+'_temp_'+stokes+'.fits'):
				os.system('rm -rf '+temp_image+'_temp_'+stokes+'.fits')	
			print ('exportfits(imagename=\''+temp_image+'_temp_'+stokes+'.image\',fitsimage=\''+temp_image+'_temp_'+stokes+'.fits\',dropdeg=True,dropstokes=True)\n')		
			exportfits(imagename=temp_image+'_temp_'+stokes+'.image',fitsimage=temp_image+'_temp_'+stokes+'.fits',dropdeg=True,dropstokes=True)
			if os.path.exists(temp_image+'_temp_'+stokes+'.fits'):
				try:
					data=fits.getdata(temp_image+'_temp_'+stokes+'.fits')
					wlist=fits.getheader(temp_image+'_temp_'+stokes+'.fits')
					w = wcs.WCS(wlist)
					if i==0:
						ax1 = fig.add_subplot(221, projection = w)
						im=ax1.imshow(data,cmap='seismic',origin='lower')
						fig.colorbar(im, ax=ax1)
						ax1.set_title('Stokes : '+stokes)
						ax1.set_xlabel('RA')
						ax1.set_ylabel('DEC')
					elif i==1:		
						ax2 = fig.add_subplot(222, projection = w)
						im=ax2.imshow(data,cmap='seismic',origin='lower')
						fig.colorbar(im, ax=ax2)
						ax2.set_title('Stokes : '+stokes)
						ax2.set_xlabel('RA')
						ax2.set_ylabel('DEC')
					elif i==2:
						ax3 = fig.add_subplot(223, projection = w)
						im=ax3.imshow(data,cmap='seismic',origin='lower')
						fig.colorbar(im, ax=ax3)
						ax3.set_title('Stokes : '+stokes)
						ax3.set_xlabel('RA')
						ax3.set_ylabel('DEC')
					elif i==3:
						ax4 = fig.add_subplot(224, projection = w)
						im=ax4.imshow(data,cmap='seismic',origin='lower')
						fig.colorbar(im, ax=ax4)
						ax4.set_title('Stokes : '+stokes)
						ax4.set_xlabel('RA')
						ax4.set_ylabel('DEC')
				except Exception as e:
					os.system('rm -rf '+temp_image+'* temp_*.image')
					print ('Errror in producing quicklook images 2: '+str(e)+'\n')
					return		
			else:
				os.system('rm -rf '+temp_image+'* temp_*.image casa*log')
				return
		except Exception as e:
			os.system('rm -rf '+temp_image+'* temp_*.image')
			print ('Errror in producing quicklook images 3: '+str(e)+'\n')
			return
	title='Frequency : '+str(freq)+' MHz, Timestamp : '+str(timestamp)+' UTC\n'
	if DR_rms!=0 and DR_neg!=0: 
		title+='Dynamic range I (rms) : '+str(int(DR_rms))+', Dynamic range I (negative) : '+str(int(DR_neg))
	plt.suptitle(title,fontsize=12)	
	cwd=os.getcwd()
	outfile_dir=os.path.dirname(outfile)
	if outfile_dir=='':
		outfile=cwd+'/'+outfile
	plt.savefig(outfile)
	os.system('rm -rf '+temp_image+'* temp_*.image casa*log')
	return outfile

def get_OBSID_for_timerange(timerange='',caldata=True,project_ID='G0002'):
	'''
	Function to get MWA observation IDs for a given time range
	Parameters
	----------
	timerange : str
		Time range to search, format 'yy0/mm0/dd0/hh0:mm0:ss0.ff0~yy1/mm1/dd1/hh1:mm1:ss1.ff1,yy2/mm2/dd2/hh2:mm2:ss2.ff2~....'
	caldata : bool
		Whether search for calibrator observations also or not
	project_ID : str
		Project ID to search data
	Returns
	-------
	list
		OBSID list
	list
		Other calibration OBSID list
	list
		Solar calibration OBSID list
	'''	
	BASEURL='http://ws.mwatelescope.org/'
	obsid_list=[]
	other_cal_obsid_list=[]
	solar_cal_obsid_list=[]
	if check_internet(url=BASEURL)==False:
		print ('Can not connect to the MWA database.\n')
		os.system('rm -rf casa*log')
		return [],[],[]
	else:
		timelist=timerange.split(',')
		for i in timelist:
			starttime=mjdsec_to_timestamp(timestamp_to_mjdsec(i.split('~')[0],format=0),includedate=True,format=1)		
			endtime=mjdsec_to_timestamp(timestamp_to_mjdsec(i.split('~')[-1],format=0),includedate=True,format=1)
			cal_starttime=mjdsec_to_timestamp(int(timestamp_to_mjdsec(i.split('~')[0],format=0))-(24*3600),includedate=True,format=1)		
			cal_endtime=mjdsec_to_timestamp(int(timestamp_to_mjdsec(i.split('~')[-1],format=0))+(24*3600),includedate=True,format=1)	
			searchurl=BASEURL+'metadata/tconv/?utciso='+starttime
			startGPStime=json.load(urllib.request.urlopen(searchurl,timeout=10))
			searchurl=BASEURL+'metadata/tconv/?utciso='+endtime
			endGPStime=json.load(urllib.request.urlopen(searchurl,timeout=10))
			print ('Searching data for timerange : '+starttime+' to '+endtime+'\n')
			print ('Searching calibrator observation for timerange : '+cal_starttime+' to '+cal_endtime+'\n')
			searchurl=BASEURL+'metadata/find?maxtime='+str(endGPStime)+'&mintime='+str(startGPStime)+'&dataquality=1&notdeleted=on&projectid='+str(project_ID)+'&page=20000000000000'
			try:
				OBSids=json.load(urllib.request.urlopen(searchurl,timeout=15))
				for i in range(len(OBSids)):
					obsid_list.append(int(OBSids[i][0]))
			except:
				pass
			searchurl=BASEURL+'metadata/tconv/?utciso='+cal_starttime
			startGPStime=json.load(urllib.request.urlopen(searchurl,timeout=10))
			searchurl=BASEURL+'metadata/tconv/?utciso='+cal_endtime
			endGPStime=json.load(urllib.request.urlopen(searchurl,timeout=10))
			searchurl=BASEURL+'metadata/find?maxtime='+str(endGPStime)+'&mintime='+str(startGPStime)+'&dataquality=1&calibration=on&notdeleted=on&page=20000000000000'
			try:
				OBSids=json.load(urllib.request.urlopen(searchurl,timeout=15))
				for i in range(len(OBSids)):
					if int(OBSids[i][0]) not in cal_obsid_list:
						other_cal_obsid_list.append(int(OBSids[i][0]))
			except:
				pass
			searchurl=BASEURL+'metadata/find?maxtime='+str(endGPStime)+'&mintime='+str(startGPStime)+'&dataquality=1&calibration=on&notdeleted=on&projectid='\
					+str(project_ID)+'&page=20000000000000'
			try:
				OBSids=json.load(urllib.request.urlopen(searchurl,timeout=15))
				for i in range(len(OBSids)):
					if int(OBSids[i][0]) not in solar_cal_obsid_list:
						solar_cal_obsid_list.append(int(OBSids[i][0]))
			except:
				pass
			if len(solar_cal_obsid_list)>0:
				for i in solar_cal_obsid_list:
					if i in other_cal_obsid_list:
						other_cal_obsid_list.remove(i)
		os.system('rm -rf casa*log')
		return obsid_list,other_cal_obsid_list,solar_cal_obsid_list

def plotms(msname='',savedir='',plotflag=False,savelog=True,xaxis='TIME',yaxis='amp',xdatacolumn='DATA',ydatacolumn='DATA',\
			iteration='',ants='',baseline='',spw='',chan='',field=0,scan=1,corr=''):
	'''
	Function to MS data using shadems
	Note : error checks for wrong input parameters has not been implemented, so please give correct inputs.
	Parameters
	----------
	msname : str
		Measurement set name
	savedir : str
		Directory to save final png files (default : current directory)
	plotflag : bool
		Plot flagged data or not (default : False)
	savelog : bool
		Save log file or not (default : False)
	xaxis : str
		X axis (TIME, CHAN, FREQ, CORR, U, V, W, UV, 'amp', 'phase', 'real', 'imag')
	yaxis : str
		Y axis (TIME, CHAN, FREQ, CORR, U, V, W, UV, 'amp', 'phase', 'real', 'imag')
	xdatacolumn : str
		X data column ('DATA','CORRECTED_DATA','MODEL_DATA','DATA-MODEL_DATA','CORRECTED_DATA-MODEL_DATA')	
	iteration : str
		Iteration axis ('field','scan','spw','corr','ant','baseline')
	ants : str
		Antenna numbers; comma separated
	baseline : str
		Baselines to plot, 'ant1-ant2', comma separated
	spw : str
		Spectral windows (comma separated)
	chan : str
		Channels start:stop:step format
	field : int
		Field ID 
	scan : int
		Scan number
	corr : str
		Correction or Stokes, 'XX,XY,YX,YY,RR,RL,LR,LL,I,Q,U,V'	
	Returns
	-------
	list
		Output image name list
	int
		Success message (0 or 1)
	'''	
	cwd=os.getcwd()
	if msname=='':
		print ('No measurement set is given.\n')
		return [],1
	else:
		os.chdir(msname) 
	if xaxis=='amp' or xaxis=='phase' or  xaxis=='real' or xaxis=='imag':
		shadems_args=['shadems','--xaxis='+xdatacolumn+':'+xaxis]
	else:
		shadems_args=['shadems','--xaxis='+xaxis]
	if yaxis=='amp' or yaxis=='phase' or  yaxis=='real' or yaxis=='imag':
		shadems_args.append('--yaxis='+ydatacolumn+':'+yaxis)
	else:
		shadems_args.append('--yaxis='+yaxis)
	if iteration!='' and iteration!=None:
		if iteration!='None':
			shadems_args.append('--iter-'+iteration)
	if ants!='':
		shadems_args.append('--ant-num='+ants)
	if baseline!='':
		shadems_args.append('--baseline='+baseline)
	if spw!='':
		shadems_args.append('--spw='+spw)
	if chan!='':
		shadems_args.append('--chan='+chan)
	if field!=0:
		shadems_args.append('--field='+str(field))
	if scan!=1:
		shadems_args.append('--scan='+str(scan))
	if corr!='':
		shadems_args.append('--corr='+corr)
	md=msmetadata()
	md.open(msname)
	nrows=int(md.nrows())
	nchan=md.nchan(0)
	chunk_size=int(nrows/3)
	shadems_args.append('-j0')
	shadems_args.append('--cmap=kb')
	shadems_args.append('--no-lim-save')
	shadems_args.append('-z'+str(chunk_size))
	if xaxis=='CHAN' or xaxis=='FREQ' or yaxis=='CHAN' or yaxis=='FREQ':
		if nchan<=1 and xaxis=='CHAN':
			shadems_args.append('--xmin=-50')
			shadems_args.append('--xmax=50')
		elif nchan<=1 and xaxis=='CHAN':
			shadems_args.append('--ymin=-50')
			shadems_args.append('--ymax=50')
		elif nchan<=1 and xaxis=='FREQ':
			freq=md.chanfreqs(0)[0]
			freqres=md.chanres(0)[0]
			shadems_args.append('--xmin='+str(freq-5*freqres))
			shadems_args.append('--xmax='+str(freq+5*freqres))
		elif nchan<=1 and yaxis=='FREQ':
			freq=md.chanfreqs(0)[0]
			freqres=md.chanres(0)[0]
			shadems_args.append('--ymin='+str(freq-5*freqres))
			shadems_args.append('--ymax='+str(freq+5*freqres))
		if xaxis=='CHAN' and nchan>1:
			shadems_args.append('--colour-by=CHAN')
			shadems_args.append('--cnum='+str(nchan))
		elif xaxis=='FREQ' and nchan>1:
			minfreq=np.min(md.chanfreqs(0))
			maxfreq=np.max(md.chanfreqs(0))
			freqres=md.chanres(0)[0]
			shadems_args.append('--xmin='+str(minfreq-5*freqres))
			shadems_args.append('--xmax='+str(maxfreq+5*freqres))
	md.close()
	if savedir!='':
		if os.path.isdir(savedir)==False:
			try:
				os.makedirs(savedir)
				savedir=savedir
			except:
				savedir=cwd
				pass
		else:
			savedir=savedir
	else:
		savedir=cwd
	shadems_args.append('--dir='+savedir)
	if plotflag==True:
		shadems_args.append('--noflags')
	shadems_args.append('--fontsize=30')
	shadems_args.append(msname)
	shadems_cmd=' '.join(shadems_args)
	if os.path.exists('log-shadems.txt'):
		os.system('rm -rf log-shadems.txt')
	a=os.system(shadems_cmd)
	if a!=0:
		return [],1
	outputname=[]
	with open('log-shadems.txt','r') as fil:
		for line in fil:
			if 'wrote ' in line:
				outputname.append(line.split('wrote ')[-1].split('\n')[0])
	if savelog==False:
		os.system('rm -rf '+savedir+'/log-shadems.txt '+savedir+'/casa*log')
	os.chdir(cwd)
	return outputname,0










