import os
import numpy as np,copy,glob,julian,multiprocessing as mp,psutil
from datetime import datetime as dt, timedelta
from casatools import msmetadata,table,measures,quanta,agentflagger,image,calibrater
from casatasks import *
from astropy.time import Time
from astropy.coordinates import get_sun
from multiprocessing import get_context
from . import basic_func as B
import paircarsdata
os.system('rm -rf casa*log')
'''
Code is written by Devojyoti Kansabanik, 06 Jan, 2021
'''
datadir = os.path.dirname(paircarsdata.__file__)
### Functions related to measurement set operations ###
class AccessMS:
	'''
	Generic class to perform measurement set related operation
	
	Parameters
	----------
	msname : str 
		Name of the measurement set
	'''
	def __init__(self,msname):
		if msname[-1]=='/':
			msname=msname[:-1]
		self.msname=msname
		self.md=msmetadata()
		self.tb=table()
		self.me=measures()
		self.af=agentflagger()
		
	def get_nbaseline(self,autocor=True):
		'''		
		Function to get number of baselines

		Parameters
		----------
		autocor : bool 
			Include auto-correlations into account or not
		Returns
		-------
		int
			Number of baselines
		'''
		self.md.open(self.msname)
		nbs=self.md.nbaselines(ac=autocor)
		self.md.close()
		return nbs

	def get_phasecenter(self):
		'''
		Get phasecenter of the measurement set

		Returns
		-------
		str
			Phasecenter of the measurement set in ['RA','DEC'] in hh mm ss dd mm ss format
		float 
			RA in degree
		float 
			DEC in degree 
		'''
		self.md.open(self.msname)
		phasecenter=self.md.phasecenter()
		freq=self.md.meanfreq(0)
		self.md.close()
		radeg=np.rad2deg(phasecenter['m0']['value'])
		if radeg<0:
			radeg+=360
		decdeg=np.rad2deg(phasecenter['m1']['value'])
		radec_str=B.radec_con_deg_to_hhmmss(radeg,decdeg)
		return radec_str,radeg,decdeg

	def get_npol(self):
		'''
		Function to get number of polarisations 
		Returns
		-------
		int
			Number of correlation products
		'''
		self.md.open(self.msname)
		npol=self.md.ncorrforpol()[0]
		self.md.close()
		return npol		

	def convert_mwa_to_iau(self):
		'''
		Convert the data from MWA coordinate (X=E->W and Y=N->S) to IAU coordinate (X=S->N,Y=W->E) system.
		Thus X=>-Y, and Y=>-X. So, XX=>YY,YY=>XX,XY=>YX,YX=>XY

		Returns
		------
		str
 			Confirmation message of coordinate conversion the measurement set
		'''
		code=vishead(vis=self.msname,mode='get',hdkey='fld_code')[0][0]
		code_list=code.split(',')
		if 'IAU' not in code_list:
			self.tb.open(self.msname,nomodify=False)
			data=self.tb.getcol('DATA')
			data_iau=copy.deepcopy(data)
			data_iau[0,:,:]=data[3,:,:]
			data_iau[1,:,:]=data[2,:,:]
			data_iau[2,:,:]=data[1,:,:]
			data_iau[3,:,:]=data[0,:,:]
			self.tb.putcol('DATA',data_iau)
			self.tb.flush()
			self.tb.close()
			del data,data_iau
			code_list.append('IAU')
			code=','.join(code_list)
			vishead(vis=self.msname,mode='put',hdkey='fld_code',hdvalue=np.array([code]))
			msg='Measurement set :'+self.msname+' has been converted in IAU convention'
			return msg
		else:
			msg='Measurement set :'+self.msname+' is already in IAU convention'
			return msg
		
	def calc_ncoarse(self):
		'''
		Calculate the number of MWA coarse channels in the measurement set

		Returns
		-------
		int
			Number of coarse channel in the measurement set
		'''
		import math
		self.md.open(self.msname)
		ncoarse=math.ceil(self.md.bandwidths()[0]/(1.28*10**6))
		self.md.close()
		return ncoarse

	def calc_freqres(self):
		'''
		Return frequency resolution of the data
	
		Returns
		-------
		float
			Frequency resolution in kHz
		'''
		self.md.open(self.msname)
		freqres=self.md.chanres(0)[0]/10**3
		self.md.close()
		return freqres
		
	def calc_bandwidth(self):
		'''
		Function to calculate bandwidth of the data
		
		Returns
		-------
		float
			Bandwidth in Hz
		'''
		self.md.open(self.msname)
		bandwidth=self.md.bandwidths()[0]
		self.md.close()
		return bandwidth

	def calc_total_time(self):
		'''
		Function to calculate total time of the measurement set

		Returns
		-------
		float
			Total time in seceond
		'''
		self.md.open(self.msname)
		times=self.md.timesforfield(0)
		if len(times)==1:
			timewidth=self.md.exposuretime(1)['value']
		else:
			timewidth=(times[-1]-times[0])+self.calc_timeres()
		self.md.close()
		return timewidth

	def calc_timeres(self):
		'''
		Function to calculate time resolution of the measurement set

		Returns
		-------
		float
			Time resolution in second
		'''
		self.md.open(self.msname)
		times=self.md.timesforfield(0)
		if len(times)==1:
			timeres=self.md.exposuretime(1)['value']
		else:
			timeres=times[1]-times[0]
		self.md.close()
		return timeres

	def get_num_timestamps(self):
		'''
		Function to calculate number of timestamps

		Returns
		-------
		int
			Number of timestamps
		'''
		return len(self.get_timestamps_in_mjdsecs()[0])

	def get_num_channels(self):
		'''
		Function to calculate number of channels

		Returns
		-------
		int
			Number of channels
		'''
		self.md.open(self.msname)
		nchan=self.md.nchan(0)
		self.md.close()
		return nchan

	def get_coarse_chan(self):
		'''
		Function to MWA coarse channels of the ms
		Returns
		-------
		list
			Coarse channel list
		'''
		freqs=self.get_freqs()
		coarse_chans=np.unique(np.array([B.freq_to_MWA_coarse(i/10**6) for i in freqs])).tolist()
		return coarse_chans
		
	def get_chans_for_coarsechan(self,coarse_chan=0):
		'''
		Function to get channel corresponding to a coarse channel
		Parameters
		----------
		int
			Coarse channel number
		Returns
		-------
		list		
			Channel numbers 
		'''
		freqs=self.get_freqs()
		all_coarse_chans=np.array([B.freq_to_MWA_coarse(i/10**6) for i in freqs])
		pos=np.where(all_coarse_chans==coarse_chan)[0]
		return pos

	def get_timestamps_in_mjdsecs(self):
		'''
		Function to return list of timestamps
	
		Returns
		-------
		list
			List of timestamps in MJD seconds
		list
			List of timestamps with errors
		'''
		self.md.open(self.msname)
		mjds=self.md.timesforfield(0)
		time_diff=np.ediff1d(mjds)
		pos=np.where(time_diff<self.calc_timeres()/2.0)
		wrong_timestamps=[]
		if len(pos[0])!=0:
			for p in pos[0]:
				try:
					wrong_timestamps.append(B.mjdsec_to_timestamp(mjds[p],includedate=True,format=0)+'~'+B.mjdsec_to_timestamp(mjds[p+1],includedate=True,format=0))
				except:
					wrong_timestamps.append(B.mjdsec_to_timestamp(mjds[p],includedate=True,format=0))			
		mjds=np.delete(mjds,pos)
		mjds=mjds.tolist()
		self.md.close()
		return mjds,wrong_timestamps

	def get_timestamps(self,format=0):
		'''
		Function to return list of timestamps

		Parameters
		----------
		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
		-------
		list
			List of timestamps in UTC at the format 'YYYY/MM/DD/hh:mm:ss'
		'''
		timestamps=[B.mjdsec_to_timestamp(mjdsec,includedate=True,format=format) for mjdsec in self.get_timestamps_in_mjdsecs()[0]]	
		return timestamps	
	
	def get_obs_startend_time(self):
		'''
		Function to get observation start time

		Returns
		-------
		str
			Observation start time in 'YYYY/MM/DD/hh:mm:ss' format
		str 
			Observation end time in 'YYYY/MM/DD/hh:mm:ss' format
		float
			Observation start time in MJD second
		float
			Observation end time in MJD second
		'''
		self.md.open(self.msname)
		start_mjdsec=self.md.timerangeforobs(0)['begin']['m0']['value']*24*3600.0
		end_mjdsec=self.md.timerangeforobs(0)['end']['m0']['value']*24*3600.0
		self.md.close()
		start=B.mjdsec_to_timestamp(start_mjdsec,includedate=True,format=0)
		end=B.mjdsec_to_timestamp(end_mjdsec,includedate=True,format=0)
		return start,end,start_mjdsec,end_mjdsec

	def get_obs_date(self,format=0):
		'''
		Function to get observation date

		Parameters
		----------
		format : str
			Date string format 
				0: 'YYYY/MM/DD'

				1: 'YYYY-MM-DD'
 
				2: 'YYYY_MM_DD'
		Returns
		-------
		str
			Observation start date
		str
			Observation end date
		'''
		start,end,start_mjdsec,end_mjdsec=self.get_obs_startend_time()
		start_mjd=start_mjdsec/(24*3600.0)
		end_mjd=end_mjdsec/(24*3600.0)
		start_year=julian.from_jd(start_mjd+2400000.5,fmt='jd').date().year
		start_month=julian.from_jd(start_mjd+2400000.5,fmt='jd').date().month
		start_day=julian.from_jd(start_mjd+2400000.5,fmt='jd').date().day
		end_year=julian.from_jd(end_mjd+2400000.5,fmt='jd').date().year
		end_month=julian.from_jd(end_mjd+2400000.5,fmt='jd').date().month
		end_day=julian.from_jd(end_mjd+2400000.5,fmt='jd').date().day
		if format==0:
			start_date=str(start_year)+'/'+str(start_month)+'/'+str(start_day)
			end_date=str(end_year)+'/'+str(end_month)+'/'+str(end_day)
		elif format==1:
			start_date=str(start_year)+'-'+str(start_month)+'-'+str(start_day)
			end_date=str(end_year)+'-'+str(end_month)+'-'+str(end_day)
		else:
			start_date=str(start_year)+'_'+str(start_month)+'_'+str(start_day)
			end_date=str(end_year)+'_'+str(end_month)+'_'+str(end_day)
		return start_date,end_date

	def get_scan_startend_time(self):
		'''
		Function to get scan start time

		Returns
		-------
		str
			Scan start time in 'YYYY/MM/DD/hh:mm:ss' format
		str 
			Scan end time in 'YYYY/MM/DD/hh:mm:ss' format
		float
			Scan start time in MJD second
		float
			Scan end time in MJD second
		'''
		self.md.open(self.msname)
		time_list=self.md.timesforfield(0)
		self.md.close()
		start_mjdsec=time_list[0]
		end_mjdsec=time_list[-1]
		start=B.mjdsec_to_timestamp(start_mjdsec,includedate=True,format=0)
		end=B.mjdsec_to_timestamp(end_mjdsec,includedate=True,format=0)
		return start,end,start_mjdsec,end_mjdsec

	def get_num_antenna(self):
		'''
		Function to get total number of antennas in the measurement set

		Returns
		-------
		int
			Number of antennas
		'''
		self.md.open(self.msname)
		nant=int(self.md.nantennas())
		self.md.close()
		return nant

	def get_antenna_string(self):
		'''
		Function to get total antenna string

		Returns
		-------
		str
			Antenna string
		'''
		nant=self.get_num_antenna()
		ant=''
		for i in range(nant):
			ant+=str(i)+','
		ant=ant[:-1]
		return ant

	def get_freqs(self):
		'''
		Function to return list of channel frequencies

		Returns
		-------
		list
			List of frequencies in Hz
		'''
		self.md.open(self.msname)
		freqs=self.md.chanfreqs(0)
		self.md.close()
		return freqs

	def get_unflag_chan(self,flagfrac=1):
		'''
		Function to get the unflagged channels if flag fraction is less than certain value

		Parameters
		----------
		flagfrac : float 
			Flag fraction per channel (default : 1.0)
		Returns
		-------
		list
			List of unflagged channels
		'''
		self.tb.open(self.msname)
		flag=self.tb.getcol('FLAG')
		self.tb.close()
		nchan=self.get_num_channels()
		unflagged_chan=[]
		for chan in range(nchan):
			flagged_data=flag[:,chan,:]
			total_data=float(len(flagged_data.flatten()))
			flagged_data=float(np.sum(flagged_data.flatten()))
			if (flagged_data/total_data)<flagfrac:
				unflagged_chan.append(chan)
		return unflagged_chan
	
	def get_unflag_freqs(self,flagfrac=1):
		'''
		Function to get the unflagger channel frequencies if the flag fraction is less than certain value
		Parameters
		----------
		flagfrac : float 
			Flag fraction per channel (default : 1.0)
		Returns
		-------
		list
			List of unflagged channels
		list
			List of unflagged channel frequencies in Hz
		'''
		freqlist=self.get_freqs()
		unflag_chans=self.get_unflag_chan(flagfrac=flagfrac)
		unflag_freqs=freqlist[unflag_chans].tolist()
		return unflag_chans,unflag_freqs

	def chan_to_freq(self,chan):
		'''
		Function to get frequency of a particular channel
		Parameters
		----------
		chan : int
			Channel number
		Returns
		-------
		float
			Frequency of the channel in MHz
		'''
		freqlist=self.get_freqs()
		return freqlist[chan]/10**6

	def get_pol_types(self):
		'''
		Function to get polarization product types
		Returns
		-------
		list
			Polarization products list
		'''
		self.md.open(self.msname)
		pol_indices=self.md.corrtypesforpol(0)	
		self.md.close()
		pol_types=[]
		pol_dic={5:'RR',6:'RL',7:'LR',8:'LL',9:'XX',10:'XY',11:'YX',12:'YY'}
		for i in pol_indices:
			pol_types.append(pol_dic[i])	
		return pol_types

	def get_field_names(self):
		'''
		Function to return field names
		Returns
		-------
		list
			List of field names
		'''
		self.md.open(self.msname)
		field_name_list=self.md.fieldnames()
		self.md.close()
		return field_name_list

	def get_field_id(self,fieldname):
		'''
		Function to return field id from field name
		Parameters
		----------
		fieldname : str
			Name of the field
		Returns
		-------
		int
			Field if
		'''
		self.md.open(self.msname)
		try:
			field_id=self.md.fieldsforname(fieldname)[0]
			self.md.close()
			return field_id	
		except:
			self.md.close()
			print ('Wrong field name.\n')
			return

	def get_scan_numers(self):
		'''
		Function to get scan numbers
		Returns
		-------
		list
			Scan numbers list
		'''
		self.md.open(self.msname)
		scan_numbers=self.md.scannumbers().tolist()
		self.md.close()
		return scan_numbers

	def get_scan_numers_of_field(self,field_id):
		'''
		Function to get scan numbers of a field
		Parameters
		----------
		int
			Field ID
		Returns
		-------
		list
			Scan numbers list
		'''
		self.md.open(self.msname)
		scan_numbers=self.md.scansforfield(field_id).tolist()
		self.md.close()
		return scan_numbers

	def get_unflag_times_mjd(self,flagfrac=1):
		'''
		Function to get the unflagged timestamps in MJD seconds if flag fraction is less than certain value

		Parameters
		----------
		flagfrac : float 
			Flag fraction per timestamp (default : 1.0)		
		Returns
		-------
		list
			List of unflagged times in MJD second
		'''
		self.tb.open(self.msname)
		flags=self.tb.getcol('FLAG')
		times=self.tb.getcol('TIME')
		self.tb.close()
		nant=self.get_num_antenna()
		n_timestamps=self.get_num_timestamps()
		n_visibility=int(times.size/n_timestamps)
		timestamps=times.reshape(n_timestamps,n_visibility)[:,0]
		flags=flags.reshape(flags.shape[0],flags.shape[1],n_timestamps,n_visibility)
		total_data_per_time=flags.shape[0]*flags.shape[1]*n_visibility
		flag_count=np.nansum(np.nansum(np.nansum(flags,axis=-1),axis=1),axis=0)/total_data_per_time
		pos=np.where(flag_count<flagfrac)
		unflagged_times_mjd=list(timestamps[pos])
		return unflagged_times_mjd

	def get_unflag_timestamps(self,flagfrac=1):
		'''
		Function to get the unflagged timestamps if flag fraction is less than certain value

		Parameters
		----------
		flagfrac : float 
			Flag fraction per timestamp (default : 1.0)
		Returns
		-------
		list
			List of unflagged timestamps in UTC format 'YYYY/MM/DD/hh:mm:ss'
		'''
		unflagged_time=self.get_unflag_times_mjd(flagfrac=flagfrac)
		timestamp_list=[B.mjdsec_to_timestamp(mjd) for mjd in unflagged_time]
		os.system('rm -rf casa*log')
		return timestamp_list

	def get_flag_timestamps(self,flagfrac=1):
		'''
		Function to get the flagged timestamps in MJD seconds if flag fraction is less than certain value

		Parameters
		----------
		flagfrac : float
 			Flag fraction per channel (default : 1.0)		
		Returns
		-------
		list
			List of flagged timestamps
		'''
		unflagged_time=self.get_unflag_timestamps(flagfrac=flagfrac)
		timestamps=self.get_timestamps()
		flagged_times=[]
		for i in timestamps:
			if i not in unflagged_time:
				flagged_times.append(i)
		return flagged_times
	
	def get_flag_times_mjd(self,flagfrac=1):
		'''
		Function to get the flagged timestamps in MJD seconds if flag fraction is less than certain value
	
		Parameters
		----------
		flagfrac : float 
			Flag fraction per channel (default : 1.0)
		Returns
		-------
		list
			List of flagged times in MJD seconds
		'''
		unflagged_time=self.get_unflag_times_mjd(flagfrac=flagfrac)
		timestamps=self.get_timestamps_in_mjdsecs()[0]
		flagged_times=[]
		for i in timestamps:
			if i not in unflagged_time:
				flagged_times.append(float(i))
		return flagged_times

	def get_model_nomodel_chan(self):
		'''
		Function to get the channels wioth no model data

		Returns
		-------
		list
			List of model and nomodel channels
		'''
		self.tb.open(self.msname)
		nchan=self.get_num_channels()
		try:
			model=self.tb.getcol('MODEL_DATA')
			self.tb.close()
		except:
			self.tb.close()
			return [],[i for i in range(nchan)]
		nomodel_chan=[]
		model_chan=[]
		for chan in range(nchan):
			model_data=model[:,chan,:]
			if np.abs(np.sum(model_data))==0.0 or np.abs(np.sum(model_data))==len(model_data.flatten())/2:
				nomodel_chan.append(chan)
			else:
				model_chan.append(chan)
		return model_chan,nomodel_chan

	def get_antenna_id(self,antenna_name=''):
		'''
		Function to get antenna id from antenna name

		Parameters
		----------
		antenna_name : str 
			Name of the antenna
		Returns
		-------
		int
			Antenna id
		'''
		if antenna_name=='':
			return None
		self.md.open(self.msname)
		antenna_names=self.md.antennanames()
		antenna_ids=self.md.antennaids()
		self.md.close()
		if antenna_name not in antenna_names:
			print ('Antenna name not in list.\n')
		else:
			pos=antenna_names.index(antenna_name)
			antenna_id=antenna_ids[pos]
			return antenna_id

	def get_antenna_name(self,antenna_id=0):
		'''
		Function to get antenna name from antenna id

		Parameters
		----------
		antenna_name : int 
			Antenna number
		Returns
		-------
		str
			Antenna name
		'''
		self.md.open(self.msname)
		antenna_names=self.md.antennanames()
		antenna_ids=self.md.antennaids().tolist()
		self.md.close()
		if antenna_id not in antenna_ids:
			print ('Antenna id not in list.\n')
		else:
			pos=antenna_ids.index(antenna_id)
			antenna_name=antenna_names[pos]
			return antenna_name

	def calc_meanfreq(self):
		'''
		Function to return central frequency of the measurement set (Only MS with single SPW)

		Returns
		-------
		float
			Central frequency in Hz
		'''
		self.md.open(self.msname)
		cent_freq=self.md.meanfreq(0)
		self.md.close()
		return cent_freq

	def calc_meanwavelength(self):
		'''
		Function to return central wavelength of the measurement set in meter
	
		Returns
		-------		
		float
			Central wavelength in metre
		'''
		freq=self.calc_meanfreq()
		wavelength=299792458.0/freq
		return wavelength

	def calc_flagfrac(self):
		'''
		Function to calculate flag fraction

		Returns
		-------
		float
			Fraction of the data flagged
		'''
		self.tb.open(self.msname)
		flag=self.tb.getcol('FLAG')
		flag_frac=np.nansum(flag)/len(flag.flatten())
		self.tb.close()
		return flag_frac

	def radec_sun(self):
		'''
		RA DEC of the Sun at the middle of the measurement set

		Returns
		-------
		str
			RA DEC of the Sun in J2000
		float
			RA in degree 
		float
			DEC in degree
		'''
		self.md.open(self.msname)
		start_end_time=self.md.timerangeforobs(0)
		self.md.close()
		start_time_mjd=start_end_time['begin']['m0']['value']
		end_time_mjd=start_end_time['end']['m0']['value']
		mid_time_mjd=(start_time_mjd+end_time_mjd)/2.0
		mid_time_mjd=Time(str(mid_time_mjd),format='mjd')
		sun_coord=get_sun(mid_time_mjd)
		sun_ra_deg=sun_coord.ra.deg
		sun_dec_deg=sun_coord.dec.deg
		if sun_ra_deg<0:
			sun_ra_deg+=360
		sun_radec=B.radec_con_deg_to_hhmmss(sun_ra_deg,sun_dec_deg) 
		sun_radec_string='J2000 '+str(sun_radec[0])+' '+str(sun_radec[1])
		return sun_radec_string,sun_ra_deg,sun_dec_deg
		
	def move_phasecenter_to_sun(self,force=False):
		'''
		Move the phasecenter of the measurement set at the center of the Sun
		Parameters
		----------
		force : bool
			Move phase center to Sun even if header says it is already moved
		Returns
		-------
		str
			Message of the measurement set phase center chaged at the center of the Sun
		'''
		sun_radec_string,sunra_deg,sundec_deg=self.radec_sun()
		radec_str,ra,dec=self.get_phasecenter()
		separation=np.sqrt((ra-sunra_deg)**2+(dec-sundec_deg)**2)*3600
		code=vishead(vis=self.msname,mode='get',hdkey='fld_code')[0][0]
		code_list=code.split(',')
		if 'FIXVIS' not in code_list or force==True:
			IB=B.ImageBasic(self.msname,includeflag=False)
			cell=IB.calc_cellsize(5)
			print ('Phase center is away from the Sun by : '+str(separation)+' arcsec.\n')
			if separation>cell:
				code_list.append('FIXVIS')
				code=','.join(code_list)
				if os.path.islink(self.msname):
					msdir=os.path.realpath(self.msname)
					os.unlink(self.msname)
					linked_path=True
				else:
					msdir=self.msname
					linked_path=False
				if os.path.isdir(msdir+'.phaseshift'):
					os.system('rm -rf '+msdir+'.phaseshift')
				
				phaseshift(vis=msdir,outputvis=msdir+'.phaseshift',phasecenter=sun_radec_string)
				if os.path.isdir(msdir):
					os.system('rm -rf '+msdir)
				os.system('mv '+msdir+'.phaseshift '+msdir)
				if linked_path:
					os.system('ln -s '+msdir+' '+self.msname)
				vishead(vis=self.msname,mode='put',hdkey='fld_code',hdvalue=np.array([code]))
				print ('Phasecenter of the observation is moved to Sun center at :'+sun_radec_string+'\n')
				return 'Phasecenter of the observation is moved to Sun center at :'+sun_radec_string+'\n'
			elif separation<=cell:
				print ('Phasecenter shift is smaller than single pixle size.\n')
				return 'Phasecenter shift is smaller than single pixle size.\n'
		else:
			print ('Phasecenter is already shifted to the Sun\n')
			return 'Phasecenter is already shifted to the Sun\n'

	def move_phasecenter_to_source(self,radec=''):
		'''
		Function to move the phasecenter of measurement set at particular RA-DEC

		Parameters
		----------
		radec : str 
			RA DEC string, format 'J2000 00h00m00s +00d00m00s'
		Returns
		-------
		str
			Message of the phasecenter of the measurement set is moved to the new phasecenter
		'''
		if radec=='':
			return 'No RA-DEC is given'
		else:
			code=vishead(vis=self.msname,mode='get',hdkey='fld_code')[0][0]
			code_list=code.split(',')
			radec_string='_'.join(radec.split(' ')[1:])
			if 'FIXVIS_'+radec_string not in code_list:
				for i in range(len(code_list)):
					if 'FIXVIS' in code_list[i]:
						code_list.remove(code_list[i])
				code=','.join(code_list)
				if os.path.islink(self.msname):
					msdir=os.path.realpath(self.msname)
					os.unlink(self.msname)
					linked_path=True
				else:
					msdir=self.msname
					linked_path=False
				if os.path.isdir(msdir+'.phaseshift'):
					os.system('rm -rf '+msdir+'.phaseshift')
				phaseshift(vis=msdir,outputvis=msdir+'.phaseshift',phasecenter=radec)
				if os.path.isdir(msdir):
					os.system('rm -rf '+msdir)
				os.system('mv '+msdir+'.phaseshift '+msdir)
				if linked_path:
					os.system('ln -s '+msdir+' '+self.msname)
				code_list.append('FIXVIS_'+radec_string)
				code=','.join(code_list)
				vishead(vis=self.msname,mode='put',hdkey='fld_code',hdvalue=np.array([code]))
				return 'Move phasecenter to :'+radec+'.\n'
			else:
				return 'Phasecenter is already at :'+radec+'.\n'
		return
			
	def get_max_baseline(self,includeflag=False):
		'''
		Get the maximum baseline in meter
		Parameters
		----------
		includeflag : bool
			Include flag data or not whilw 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)
		nanpos=np.where(np.isnan(uvdist)==True)
		uvdist=np.delete(uvdist,nanpos)
		return int(np.max(uvdist))

	def get_min_baseline(self,includeflag=False):
		'''
		Get the minimum baseline in meter
		Parameters
		----------
		includeflag : bool
			Include flag data or not whilw calculating minimum baseline
		Returns
		-------
		float
			Minimum baseline length in meter
		float
			Seceond minimum baseline 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)
		nanpos=np.where(np.isnan(uvdist)==True)
		uvdist=np.delete(uvdist,nanpos)
		zeropos=np.where(uvdist==0)
		uvdist=np.delete(uvdist,zeropos)
		return np.sort(uvdist)[0],np.sort(uvdist)[1]

	def get_listobs(self,listobsfile=''):
		'''
		Function to save listobs
	
		Parameters
		----------
		listobsfile : str
			File name to save listobs output
		Returns
		-------
		str
			listobs file name
		'''
		if listobsfile=='':
			listobsfile=self.msname.split('.ms')[0]+'.listobs'
		if os.path.isfile(listobsfile):
			os.system('rm -rf '+listobsfile)
		listobs(vis=self.msname,listfile=listobsfile)
		return listobsfile

	def model_imported(self):
		'''
		Fuction to check whether model imported or model column exists or not
		Returns
		-------
		bool
			True if model exists and False if does not
		''' 
		self.tb.open(self.msname)
		try:
			model=self.tb.getcol('MODEL_DATA')
			self.tb.close()
			if model.shape[0]==2:
				if np.nansum(np.real(model))==model.flatten().shape[0]:
					return False
				else:
					return True
			elif model.shape[0]==4:
				if np.nansum(np.real(model))==model.flatten().shape[0]/2:
					return False
				else:
					return True
			else:
				return False
		except:
			self.tb.close()
			return False	

	
	def make_antenna_list(self,num_bins=5,antenna_list_file=''):
		'''
		Make the antenna addition list splited in number of bins

		Parameters
		----------
		num_bin : int 
			Number of antenna bins
		antenna_list_file : str
 			Name of the file to save antenna list. If an antenna list file given, antenna list will be loaded from that file
		Returns
		-------
		list
			Lists of antennas to be added in steps
		int 
			Number of antenna
		'''
		if antenna_list_file!='':
			try:
				num_ant,antenna_list=np.load(antenna_list_file,allow_pickle=True)
			except:
				num_ant,antenna_list=np.load(antenna_list_file+'.npy',allow_pickle=True)
			return antenna_list,num_ant
		self.tb.open(self.msname+'/ANTENNA')
		coords=self.tb.getcol('POSITION')
		self.tb.close()
		cm_coords=np.repeat(np.median(coords,axis=1)[:,np.newaxis],coords.shape[1],axis=1)
		rel_coords=coords-cm_coords
		dist=np.sqrt(rel_coords[0,:]**2+rel_coords[1,:]**2)	
		pos=np.argsort(dist)
		while num_bins>=1:
			val,bins=np.histogram(dist,num_bins)
			if val[0]<len(dist)*0.3:
				num_bins-=1				
				continue
			else:
				break
		antenna_list=[]
		start=0
		for elem in val:
			antenna_list.append(np.sort(pos[0:start+elem].astype('int')))
			start+=elem
		if antenna_list_file!='':
			antenna_list_file=antenna_list_file.split('.')[0]
			np.save(antenna_list_file,np.array(['Antennafile',len(dist),antenna_list]))
		return antenna_list,len(dist)

	def get_observatory_loc(self):
		'''
		Give the observatory geodetic location

		Returns
		-------
		float
			Latitude in degree
		float 
			Longitude in degree
		float
			Altitude in meter
		'''
		self.md.open(self.msname)
		pos=self.md.observatoryposition()
		radeg=np.rad2deg(pos['m0']['value'])
		if radeg<0:
			radeg+=360
		LON=round(radeg,2)
		LAT=round(np.rad2deg(pos['m1']['value']),2)
		ALT=pos['m2']['value']
		self.md.close()
		return LAT,LON,ALT

	def get_altaz(self,source_field=0,source_scan=1):
		'''
		Calculate alt az of the phasecenter

		Parameters
		----------
		source_field : int 
			FIELD ID of the source
		source_scan : int 
			Scan number 
		Returns
		-------
		float
			Altitude in radian
		float
			Azimuth in radian
		''' 
		LAT,LON,ALT=self.get_observatory_loc()
		radec_str,radeg,decdeg=self.get_phasecenter()
		self.md.open(self.msname)
		scans=self.md.scansforfield(source_field)
		if source_scan not in scans:
			print ('Scan is not in given field ID')
			self.md.close()
			return
		else:
			mjd=np.mean(self.md.timesforscan(source_scan))
			time_string=B.mjdsec_to_timestamp(mjd,includedate=True,format=2)
			alt,az=B.radec_to_altaz(radeg,decdeg,time_string,LAT,LON,ALT)
			alt=np.deg2rad(alt)
			az=np.deg2rad(az)
			self.md.close()
			return alt,az

	def get_phasecenter_parang(self,source_field=0,combine='field'):
		'''
		Calculate parallactic for phasecenter at a given Earth location.
		Note = Parallactic angle is defined as the orientation of the sky in telescope coordinate. All angles are defined positive in IAU defined sky coordiunate (North to East).
		So, the rotation of the sky with respect to telescope is negative in the sky coordinate. To account this effect parallactic angle is given in 360-parang form.

		Parameters
		----------
		source_field : int 
			FIELD id of the source
		combine : str 
			Combine 'field' or 'scan'  for calculating parallactic angle or '' for no combine
		Returns
		-------
		float	
			1. combine = 'field', A single parallactic angle for the entire field in degree
		dict
			2. combine = 'scan', A python dictionary {'scan':parang} format
		list
			3. combine = '', A list of parang for all timestamps
		'''
		LAT,LON,ALT=self.get_observatory_loc()
		radec_str,radeg,decdeg=self.get_phasecenter()
		self.md.open(self.msname)
		if combine=='field':
			mjd=np.mean(self.md.timesforfield(source_field))
			time_string=B.mjdsec_to_timestamp(mjd,includedate=True,format=2)
			alt,az=B.radec_to_altaz(radeg,decdeg,time_string,LAT,LON,ALT)
			alt=np.deg2rad(alt)
			az=np.deg2rad(az)
			lat=np.deg2rad(LAT)
			p=-np.arctan2(np.sin(az)*np.cos(lat),np.cos(alt)*np.sin(lat) - np.sin(alt)*np.cos(lat)*np.cos(az))
			self.md.close()
			return 360-np.rad2deg(p)
		elif combine=='scan':
			scans=self.md.scansforfield(source_field)
			mjds=[np.mean(self.md.timesforscan(source_scan)) for source_scan in scans]
			time_strings=[B.mjdsec_to_timestamp(mjdsec,includedate=True,format=2) for mjdsec in mjds]
			parang={}		
			for i in range(len(scans)):
				obstime=time_strings[i]
				scan=scans[i]
				alt,az=B.radec_to_altaz(radeg,decdeg,obstime,LAT,LON,ALT)
				alt=np.deg2rad(alt)
				az=np.deg2rad(az)
				lat=np.deg2rad(LAT)
				p=-np.arctan2(np.sin(az)*np.cos(lat),np.cos(alt)*np.sin(lat) - np.sin(alt)*np.cos(lat)*np.cos(az))
				parang[scan]=360-np.rad2deg(p)
				self.md.close()
			return parang 
		else:
			mjds=self.md.timesforfield(source_field) 
			time_strings=[B.mjdsec_to_timestamp(mjdsec,includedate=True,format=2) for mjdsec in mjds]
			parang=[]		
			for i in range(len(time_strings)):
				obstime=time_strings[i]
				alt,az=B.radec_to_altaz(radeg,decdeg,obstime,LAT,LON,ALT)
				alt=np.deg2rad(alt)
				az=np.deg2rad(az)
				lat=np.deg2rad(LAT)
				p=-np.arctan2(np.sin(az)*np.cos(lat),np.cos(alt)*np.sin(lat) - np.sin(alt)*np.cos(lat)*np.cos(az))
				parang.append(360-np.rad2deg(p))
				self.md.close()
			return parang 

	def get_parang(self,ra,dec,source_field=0,combine='field'):
		'''
		Calculate parallactic for phasecenter at a given Earth location for a given RA DEC
		Note = Parallactic angle is defined as the orientation of the sky in telescope coordinate. All angles are defined positive in IAU defined sky coordiunate (North to East).
		So, the rotation of the sky with respect to telescope is negative in the sky coordinate. To account this effect parallactic angle is given in 360-parang form.

		Parameters
		----------
		ra : float
			RA in degree
		dec : float
			DEC in degree
		source_field : int 
			FIELD id of the source
		combine : str 
			Combine 'field' or 'scan' for calculating parallactic angle or '' for no combine
		Returns
		-------
		float
			1. combine = 'field', A single parallactic angle for the entire field in degree
		dict
			2. combine = 'scan', A python dictionary {'scan':parang} format
		list
			3. combine = '', A list of parang for all timestamps
		'''
		LAT,LON,ALT=self.get_observatory_loc()
		radeg=ra
		decdeg=dec
		self.md.open(self.msname)
		if combine=='field':
			mjd=np.mean(self.md.timesforfield(source_field))
			time_string=B.mjdsec_to_timestamp(mjd,includedate=True,format=2)
			alt,az=B.radec_to_altaz(radeg,decdeg,time_string,LAT,LON,ALT)
			alt=np.deg2rad(alt)
			az=np.deg2rad(az)
			lat=np.deg2rad(LAT)
			p=-np.arctan2(np.sin(az)*np.cos(lat),np.cos(alt)*np.sin(lat) - np.sin(alt)*np.cos(lat)*np.cos(az))
			self.md.close()
			return 360-np.rad2deg(p)
		elif combine=='scan':
			scans=self.md.scansforfield(source_field)
			mjds=[np.mean(self.md.timesforscan(source_scan)) for source_scan in scans]
			time_strings=[B.mjdsec_to_timestamp(mjdsec,includedate=True,format=2) for mjdsec in mjds]
			parang={}		
			for i in range(len(scans)):
				obstime=time_strings[i]
				scan=scans[i]
				alt,az=B.radec_to_altaz(radeg,decdeg,obstime,LAT,LON,ALT)
				alt=np.deg2rad(alt)
				az=np.deg2rad(az)
				lat=np.deg2rad(LAT)
				p=-np.arctan2(np.sin(az)*np.cos(lat),np.cos(alt)*np.sin(lat) - np.sin(alt)*np.cos(lat)*np.cos(az))
				parang[scan]=360-np.rad2deg(p)
				self.md.close()
			return parang 
		else:
			mjds=self.md.timesforfield(source_field) 
			time_strings=[B.mjdsec_to_timestamp(mjdsec,includedate=True,format=2) for mjdsec in mjds]
			parang=[]		
			for i in range(len(time_strings)):
				obstime=time_strings[i]
				alt,az=B.radec_to_altaz(radeg,decdeg,obstime,LAT,LON,ALT)
				alt=np.deg2rad(alt)
				az=np.deg2rad(az)
				lat=np.deg2rad(LAT)
				p=-np.arctan2(np.sin(az)*np.cos(lat),np.cos(alt)*np.sin(lat) - np.sin(alt)*np.cos(lat)*np.cos(az))				
				parang.append(360-np.rad2deg(p))
				self.md.close()
			return parang 

	def remove_locks(self):
		'''
		Function to remove tablelock of the ms	
		'''
		showtablelock_path=str(np.load(datadir+'/local_library.npy',allow_pickle=True))
		os.path.join(showtablelock_path)
		os.system('showtablelock '+str(self.msname)+' > '+str(self.msname)+'.write_lock')
		fil=open(str(self.msname)+'.write_lock','r')
		line=fil.readline()
		fil.close()
		if 'write-locked in process' in line or 'read-locked in process' in line or 'is opened (but not locked)' in line:
			process_id=int(line.split(' ')[-1])
			try:
				os.kill(process_id,9)
			except:
				pass
		os.system('rm -rf '+str(self.msname)+'.write_lock')
		return
	
def splited_ms_rename(msname,ref_time_chan=True,change_msname=True):
	'''
	Convert the name of a splited measurement set according to its frequency and timestamp

	Parameters
	----------
	msname : str 
		Name of the measurement set
	ref_time_chan : bool 
		Whether the time and frequency slice is refernce
	change_msname : bool
		Change the msname or just return the name
	Returns
	-------
	str
		Modified name of the measurement set
	'''
	md=msmetadata()
	md.open(msname)
	mjdsecond=md.timesforscan(1)#[0]
	meanmjd=np.nanmean(mjdsecond)
	mean_freq=md.meanfreq(0)/10**6 # In MHz
	md.close()
	utcstring=B.mjdsec_to_timestamp(meanmjd,includedate=True,format=0)
	yyyy=utcstring.split('/')[0]
	mm=utcstring.split('/')[1]
	dd=utcstring.split('/')[2]
	timestamp='_'.join(utcstring.split('/')[3].split(':')) 
	ms_path=os.path.dirname(os.path.realpath(msname))
	if ref_time_chan==False:
		filename=ms_path+'/time_'+yyyy+'_'+mm+'_'+dd+'_'+timestamp+'_freq_'+str(mean_freq)+'.ms'
	else:
		filename=ms_path+'/time_'+yyyy+'_'+mm+'_'+dd+'_'+timestamp+'_freq_'+str(mean_freq)+'_ref.ms'
	if change_msname==True and msname!=filename:
		if os.path.exists(filename):
			os.system('rm -rf '+filename)
		os.system('mv '+msname+' ' +filename)
	return filename











