Context Navigation

source: trunk/python/scantable.py@ 529

Visit:

Last change on this file since 529 was 524, checked in by kil064, 20 years ago
make 'tint' weigting the default for average_time
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 55.1 KB

Rev	Line
[102]	1	from asap._asap import sdtable
[226]	2	from asap import rcParams
[189]	3	from numarray import ones,zeros
[102]	4	import sys
	5
	6	class scantable(sdtable):
	7	"""
	8	The ASAP container for scans
	9	"""
[113]	10
[484]	11	def __init__(self, filename, average=None, unit=None):
[102]	12	"""
	13	Create a scantable from a saved one or make a reference
	14	Parameters:
[181]	15	filename: the name of an asap table on disk
	16	or
	17	the name of a rpfits/sdfits/ms file
	18	(integrations within scans are auto averaged
	19	and the whole file is read)
	20	or
	21	[advanced] a reference to an existing
[102]	22	scantable
[484]	23	average: average all integrations withinb a scan on read.
	24	The default (True) is taken from .asaprc.
	25	unit: brightness unit; must be consistent with K or Jy.
[340]	26	Over-rides the default selected by the reader
	27	(input rpfits/sdfits/ms) or replaces the value
	28	in existing scantables
[489]	29	"""
	30	if average is None or type(average) is not bool:
	31	average = rcParams['scantable.autoaverage']
	32
[484]	33	varlist = vars()
[226]	34	self._vb = rcParams['verbose']
[113]	35	self._p = None
[484]	36
[181]	37	if isinstance(filename,sdtable):
	38	sdtable.__init__(self, filename)
[340]	39	if unit is not None:
[346]	40	self.set_fluxunit(unit)
[181]	41	else:
[411]	42	import os.path
	43	if not os.path.exists(filename):
	44	print "File '%s' not found." % (filename)
[226]	45	return
[411]	46	filename = os.path.expandvars(filename)
	47	if os.path.isdir(filename):
[181]	48	# crude check if asap table
[411]	49	if os.path.exists(filename+'/table.info'):
[181]	50	sdtable.__init__(self, filename)
[340]	51	if unit is not None:
	52	self.set_fluxunit(unit)
[181]	53	else:
[411]	54	print "The given file '%s'is not a valid asap table." % (filename)
[226]	55	return
[484]	56	else:
[181]	57	from asap._asap import sdreader
[340]	58	ifSel = -1
	59	beamSel = -1
[345]	60	r = sdreader(filename,ifSel,beamSel)
[181]	61	print 'Importing data...'
[411]	62	r._read([-1])
	63	tbl = r._getdata()
[346]	64	if unit is not None:
	65	tbl.set_fluxunit(unit)
[489]	66	if average:
	67	from asap._asap import average as _av
[226]	68	print 'Auto averaging integrations...'
[513]	69	tbl2 = _av((tbl,),(),True,'none')
[226]	70	sdtable.__init__(self,tbl2)
[513]	71	del tbl2
[226]	72	else:
	73	sdtable.__init__(self,tbl)
[513]	74	del r,tbl
[484]	75	self._add_history("scantable", varlist)
[102]	76
[451]	77	def save(self, name=None, format=None, stokes=False, overwrite=False):
[116]	78	"""
[280]	79	Store the scantable on disk. This can be an asap (aips++) Table, SDFITS,
	80	Image FITS or MS2 format.
[116]	81	Parameters:
[196]	82	name: the name of the outputfile. For format="FITS" this
[489]	83	is the directory file name into which all the files
[497]	84	will be written (default is 'asap_FITS'). For format
	85	"ASCII" this is the root file name (data in 'name'.txt
	86	and header in 'name'_header.txt)
[116]	87	format: an optional file format. Default is ASAP.
[280]	88	Allowed are - 'ASAP' (save as ASAP [aips++] Table),
[194]	89	'SDFITS' (save as SDFITS file)
	90	'FITS' (saves each row as a FITS Image)
[200]	91	'ASCII' (saves as ascii text file)
[226]	92	'MS2' (saves as an aips++
	93	MeasurementSet V2)
[489]	94	stokes: Convert to Stokes parameters (only available
	95	currently with FITS and ASCII formats.
	96	Default is False.
[411]	97	overwrite: If the file should be overwritten if it exists.
[256]	98	The default False is to return with warning
[411]	99	without writing the output. USE WITH CARE.
[116]	100	Example:
	101	scan.save('myscan.asap')
	102	scan.save('myscan.sdfits','SDFITS')
	103	"""
[411]	104	from os import path
[226]	105	if format is None: format = rcParams['scantable.save']
[256]	106	suffix = '.'+format.lower()
	107	if name is None or name =="":
	108	name = 'scantable'+suffix
[283]	109	print "No filename given. Using default name %s..." % name
[411]	110	name = path.expandvars(name)
[256]	111	if path.isfile(name) or path.isdir(name):
	112	if not overwrite:
	113	print "File %s already exists." % name
	114	return
[451]	115	format2 = format.upper()
	116	if format2 == 'ASAP':
[116]	117	self._save(name)
	118	else:
	119	from asap._asap import sdwriter as _sw
[451]	120	w = _sw(format2)
[445]	121	w.write(self, name, stokes)
[116]	122	return
	123
[102]	124	def copy(self):
	125	"""
	126	Return a copy of this scantable.
	127	Parameters:
[113]	128	none
[102]	129	Example:
	130	copiedscan = scan.copy()
	131	"""
[113]	132	sd = scantable(sdtable._copy(self))
	133	return sd
	134
[102]	135	def get_scan(self, scanid=None):
	136	"""
	137	Return a specific scan (by scanno) or collection of scans (by
	138	source name) in a new scantable.
	139	Parameters:
[513]	140	scanid: a (list of) scanno or a source name, unix-style
	141	patterns are accepted for source name matching, e.g.
	142	'*_R' gets all 'ref scans
[102]	143	Example:
[513]	144	# get all scans containing the source '323p459'
	145	newscan = scan.get_scan('323p459')
	146	# get all 'off' scans
	147	refscans = scan.get_scan('*_R')
	148	# get a susbset of scans by scanno (as listed in scan.summary())
	149	newscan = scan.get_scan([0,2,7,10])
[102]	150	"""
	151	if scanid is None:
	152	print "Please specify a scan no or name to retrieve from the scantable"
	153	try:
	154	if type(scanid) is str:
[113]	155	s = sdtable._getsource(self,scanid)
[102]	156	return scantable(s)
	157	elif type(scanid) is int:
[381]	158	s = sdtable._getscan(self,[scanid])
	159	return scantable(s)
	160	elif type(scanid) is list:
[113]	161	s = sdtable._getscan(self,scanid)
[102]	162	return scantable(s)
[381]	163	else:
	164	print "Illegal scanid type, use 'int' or 'list' if ints."
[102]	165	except RuntimeError:
	166	print "Couldn't find any match."
	167
	168	def __str__(self):
[381]	169	return sdtable._summary(self,True)
[102]	170
[381]	171	def summary(self,filename=None, verbose=None):
[102]	172	"""
	173	Print a summary of the contents of this scantable.
	174	Parameters:
	175	filename: the name of a file to write the putput to
	176	Default - no file output
[381]	177	verbose: print extra info such as the frequency table
	178	The default (False) is taken from .asaprc
[102]	179	"""
[381]	180	info = sdtable._summary(self, verbose)
	181	if verbose is None: verbose = rcParams['scantable.verbosesummary']
[102]	182	if filename is not None:
[256]	183	if filename is "":
	184	filename = 'scantable_summary.txt'
[415]	185	from os.path import expandvars, isdir
[411]	186	filename = expandvars(filename)
[415]	187	if not isdir(filename):
[413]	188	data = open(filename, 'w')
	189	data.write(info)
	190	data.close()
	191	else:
	192	print "Illegal file name '%s'." % (filename)
[102]	193	print info
[484]	194
[256]	195	def set_cursor(self, thebeam=0,theif=0,thepol=0):
[102]	196	"""
	197	Set the spectrum for individual operations.
	198	Parameters:
	199	thebeam,theif,thepol: a number
	200	Example:
[256]	201	scan.set_cursor(0,0,1)
[135]	202	pol1sig = scan.stats(all=False) # returns std dev for beam=0
[181]	203	# if=0, pol=1
[102]	204	"""
[484]	205	varlist = vars()
[102]	206	self.setbeam(thebeam)
	207	self.setpol(thepol)
	208	self.setif(theif)
[484]	209	self._add_history("set_cursor",varlist)
[102]	210	return
	211
[256]	212	def get_cursor(self):
[113]	213	"""
	214	Return/print a the current 'cursor' into the Beam/IF/Pol cube.
	215	Parameters:
	216	none
	217	Returns:
	218	a list of values (currentBeam,currentIF,currentPol)
	219	Example:
	220	none
	221	"""
[102]	222	i = self.getbeam()
	223	j = self.getif()
	224	k = self.getpol()
[113]	225	if self._vb:
[181]	226	print "--------------------------------------------------"
[256]	227	print " Cursor position"
[181]	228	print "--------------------------------------------------"
[226]	229	out = 'Beam=%d IF=%d Pol=%d ' % (i,j,k)
[113]	230	print out
[102]	231	return i,j,k
	232
[436]	233	def stats(self, stat='stddev', mask=None, allaxes=None):
[102]	234	"""
[135]	235	Determine the specified statistic of the current beam/if/pol
[102]	236	Takes a 'mask' as an optional parameter to specify which
	237	channels should be excluded.
	238	Parameters:
[135]	239	stat: 'min', 'max', 'sumsq', 'sum', 'mean'
	240	'var', 'stddev', 'avdev', 'rms', 'median'
	241	mask: an optional mask specifying where the statistic
[102]	242	should be determined.
[436]	243	allaxes: if True apply to all spectra. Otherwise
	244	apply only to the selected (beam/pol/if)spectra only.
	245	The default is taken from .asaprc (True if none)
[102]	246	Example:
[113]	247	scan.set_unit('channel')
[102]	248	msk = scan.create_mask([100,200],[500,600])
[135]	249	scan.stats(stat='mean', mask=m)
[102]	250	"""
[436]	251	if allaxes is None: allaxes = rcParams['scantable.allaxes']
[135]	252	from asap._asap import stats as _stats
[256]	253	from numarray import array,zeros,Float
[102]	254	if mask == None:
	255	mask = ones(self.nchan())
[256]	256	axes = ['Beam','IF','Pol','Time']
	257
[436]	258	beamSel,IFSel,polSel = (self.getbeam(),self.getif(),self.getpol())
	259	if allaxes:
[256]	260	n = self.nbeam()self.nif()self.npol()*self.nrow()
	261	shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
	262	arr = array(zeros(n),shape=shp,type=Float)
	263
	264	for i in range(self.nbeam()):
	265	self.setbeam(i)
	266	for j in range(self.nif()):
	267	self.setif(j)
	268	for k in range(self.npol()):
	269	self.setpol(k)
	270	arr[i,j,k,:] = _stats(self,mask,stat,-1)
	271	retval = {'axes': axes, 'data': arr, 'cursor':None}
	272	tm = [self._gettime(val) for val in range(self.nrow())]
[181]	273	if self._vb:
[256]	274	self._print_values(retval,stat,tm)
[436]	275	self.setbeam(beamSel)
	276	self.setif(IFSel)
	277	self.setpol(polSel)
[256]	278	return retval
[102]	279
	280	else:
[241]	281	statval = _stats(self,mask,stat,-1)
[181]	282	out = ''
	283	for l in range(self.nrow()):
	284	tm = self._gettime(l)
	285	out += 'Time[%s]:\n' % (tm)
[436]	286	if self.nbeam() > 1: out += ' Beam[%d] ' % (beamSel)
	287	if self.nif() > 1: out += ' IF[%d] ' % (IFSel)
	288	if self.npol() > 1: out += ' Pol[%d] ' % (polSel)
[181]	289	out += '= %3.3f\n' % (statval[l])
	290	out += "--------------------------------------------------\n"
[226]	291
[181]	292	if self._vb:
	293	print "--------------------------------------------------"
	294	print " ",stat
	295	print "--------------------------------------------------"
	296	print out
[256]	297	retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
	298	return retval
[102]	299
[436]	300	def stddev(self,mask=None, allaxes=None):
[135]	301	"""
	302	Determine the standard deviation of the current beam/if/pol
	303	Takes a 'mask' as an optional parameter to specify which
	304	channels should be excluded.
	305	Parameters:
	306	mask: an optional mask specifying where the standard
	307	deviation should be determined.
[436]	308	allaxes: optional flag to show all or a cursor selected
[226]	309	spectrum of Beam/IF/Pol. Default is all or taken
	310	from .asaprc
[135]	311
	312	Example:
	313	scan.set_unit('channel')
	314	msk = scan.create_mask([100,200],[500,600])
	315	scan.stddev(mask=m)
	316	"""
[436]	317	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	318	return self.stats(stat='stddev',mask=mask, allaxes=allaxes);
[135]	319
[436]	320	def get_tsys(self, allaxes=None):
[113]	321	"""
	322	Return the System temperatures.
	323	Parameters:
[436]	324	allaxes: if True apply to all spectra. Otherwise
	325	apply only to the selected (beam/pol/if)spectra only.
	326	The default is taken from .asaprc (True if none)
[113]	327	Returns:
	328	a list of Tsys values.
	329	"""
[436]	330	if allaxes is None: allaxes = rcParams['scantable.allaxes']
[256]	331	from numarray import array,zeros,Float
	332	axes = ['Beam','IF','Pol','Time']
	333
[436]	334	if allaxes:
[256]	335	n = self.nbeam()self.nif()self.npol()*self.nrow()
	336	shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
	337	arr = array(zeros(n),shape=shp,type=Float)
	338
	339	for i in range(self.nbeam()):
	340	self.setbeam(i)
	341	for j in range(self.nif()):
	342	self.setif(j)
	343	for k in range(self.npol()):
	344	self.setpol(k)
	345	arr[i,j,k,:] = self._gettsys()
	346	retval = {'axes': axes, 'data': arr, 'cursor':None}
	347	tm = [self._gettime(val) for val in range(self.nrow())]
[181]	348	if self._vb:
[256]	349	self._print_values(retval,'Tsys',tm)
	350	return retval
	351
[102]	352	else:
[256]	353	i,j,k = (self.getbeam(),self.getif(),self.getpol())
	354	statval = self._gettsys()
[181]	355	out = ''
	356	for l in range(self.nrow()):
	357	tm = self._gettime(l)
	358	out += 'Time[%s]:\n' % (tm)
	359	if self.nbeam() > 1: out += ' Beam[%d] ' % (i)
	360	if self.nif() > 1: out += ' IF[%d] ' % (j)
	361	if self.npol() > 1: out += ' Pol[%d] ' % (k)
[256]	362	out += '= %3.3f\n' % (statval[l])
	363	out += "--------------------------------------------------\n"
	364
[181]	365	if self._vb:
	366	print "--------------------------------------------------"
[256]	367	print " TSys"
[181]	368	print "--------------------------------------------------"
	369	print out
[256]	370	retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
	371	return retval
[113]	372
[407]	373	def get_time(self, row=-1):
[113]	374	"""
	375	Get a list of time stamps for the observations.
[181]	376	Return a string for each integration in the scantable.
[113]	377	Parameters:
[407]	378	row: row no of integration. Default -1 return all rows
[113]	379	Example:
	380	none
	381	"""
	382	out = []
[407]	383	if row == -1:
	384	for i in range(self.nrow()):
	385	out.append(self._gettime(i))
	386	return out
	387	else:
	388	if row < self.nrow():
	389	return self._gettime(row)
[102]	390
	391	def set_unit(self, unit='channel'):
	392	"""
	393	Set the unit for all following operations on this scantable
	394	Parameters:
	395	unit: optional unit, default is 'channel'
[113]	396	one of '*Hz','km/s','channel', ''
[102]	397	"""
[484]	398	varlist = vars()
[113]	399	if unit in ['','pixel', 'channel']:
	400	unit = ''
	401	inf = list(self._getcoordinfo())
	402	inf[0] = unit
	403	self._setcoordinfo(inf)
	404	if self._p: self.plot()
[484]	405	self._add_history("set_unit",varlist)
[113]	406
[484]	407	def set_instrument(self, instr):
[358]	408	"""
	409	Set the instrument for subsequent processing
	410	Parameters:
[407]	411	instr: Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
	412	'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
[358]	413	"""
	414	self._setInstrument(instr)
[484]	415	self._add_history("set_instument",vars())
[358]	416
[276]	417	def set_doppler(self, doppler='RADIO'):
	418	"""
	419	Set the doppler for all following operations on this scantable.
	420	Parameters:
	421	doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
	422	"""
[484]	423	varlist = vars()
[276]	424	inf = list(self._getcoordinfo())
	425	inf[2] = doppler
	426	self._setcoordinfo(inf)
	427	if self._p: self.plot()
[484]	428	self._add_history("set_doppler",vars())
[432]	429
[226]	430	def set_freqframe(self, frame=None):
[113]	431	"""
	432	Set the frame type of the Spectral Axis.
	433	Parameters:
	434	frame: an optional frame type, default 'LSRK'.
	435	Examples:
	436	scan.set_freqframe('BARY')
	437	"""
[484]	438	if frame is None: frame = rcParams['scantable.freqframe']
	439	varlist = vars()
[113]	440	valid = ['REST','TOPO','LSRD','LSRK','BARY', \
	441	'GEO','GALACTO','LGROUP','CMB']
	442	if frame in valid:
	443	inf = list(self._getcoordinfo())
	444	inf[1] = frame
	445	self._setcoordinfo(inf)
[484]	446	self._add_history("set_freqframe",varlist)
[102]	447	else:
[113]	448	print "Please specify a valid freq type. Valid types are:\n",valid
	449
	450	def get_unit(self):
	451	"""
	452	Get the default unit set in this scantable
	453	Parameters:
	454	Returns:
	455	A unit string
	456	"""
	457	inf = self._getcoordinfo()
	458	unit = inf[0]
	459	if unit == '': unit = 'channel'
	460	return unit
[102]	461
[158]	462	def get_abcissa(self, rowno=0):
[102]	463	"""
[158]	464	Get the abcissa in the current coordinate setup for the currently
[113]	465	selected Beam/IF/Pol
	466	Parameters:
[226]	467	rowno: an optional row number in the scantable. Default is the
	468	first row, i.e. rowno=0
[113]	469	Returns:
[407]	470	The abcissa values and it's format string (as a dictionary)
[113]	471	"""
[256]	472	abc = self._getabcissa(rowno)
[407]	473	lbl = self._getabcissalabel(rowno)
[158]	474	return abc, lbl
[407]	475	#return {'abcissa':abc,'label':lbl}
[113]	476
	477	def create_mask(self, args, *kwargs):
	478	"""
[102]	479	Compute and return a mask based on [min,max] windows.
[189]	480	The specified windows are to be INCLUDED, when the mask is
[113]	481	applied.
[102]	482	Parameters:
	483	[min,max],[min2,max2],...
	484	Pairs of start/end points specifying the regions
	485	to be masked
[189]	486	invert: optional argument. If specified as True,
	487	return an inverted mask, i.e. the regions
	488	specified are EXCLUDED
[513]	489	row: create the mask using the specified row for
	490	unit conversions, default is row=0
	491	only necessary if frequency varies over rows.
[102]	492	Example:
[113]	493	scan.set_unit('channel')
	494
	495	a)
[102]	496	msk = scan.set_mask([400,500],[800,900])
[189]	497	# masks everything outside 400 and 500
[113]	498	# and 800 and 900 in the unit 'channel'
	499
	500	b)
	501	msk = scan.set_mask([400,500],[800,900], invert=True)
[189]	502	# masks the regions between 400 and 500
[113]	503	# and 800 and 900 in the unit 'channel'
	504
[102]	505	"""
[513]	506	row = 0
	507	if kwargs.has_key("row"):
	508	row = kwargs.get("row")
	509	data = self._getabcissa(row)
[113]	510	u = self._getcoordinfo()[0]
	511	if self._vb:
	512	if u == "": u = "channel"
	513	print "The current mask window unit is", u
[102]	514	n = self.nchan()
[189]	515	msk = zeros(n)
[484]	516	for window in args:
[102]	517	if (len(window) != 2 or window[0] > window[1] ):
	518	print "A window needs to be defined as [min,max]"
	519	return
	520	for i in range(n):
[113]	521	if data[i] >= window[0] and data[i] < window[1]:
[189]	522	msk[i] = 1
[113]	523	if kwargs.has_key('invert'):
	524	if kwargs.get('invert'):
	525	from numarray import logical_not
[513]	526	msk = logical_not(msk)
[102]	527	return msk
[113]	528
[256]	529	def get_restfreqs(self):
	530	"""
	531	Get the restfrequency(s) stored in this scantable.
	532	The return value(s) are always of unit 'Hz'
	533	Parameters:
	534	none
	535	Returns:
	536	a list of doubles
	537	"""
	538	return list(self._getrestfreqs())
[102]	539
[402]	540	def lines(self):
[391]	541	"""
[402]	542	Print the list of known spectral lines
	543	"""
	544	sdtable._lines(self)
	545
[484]	546	def set_restfreqs(self, freqs=None, unit='Hz', lines=None, source=None,
	547	theif=None):
[402]	548	"""
[393]	549	Select the restfrequency for the specified source and IF OR
	550	replace for all IFs. If the 'freqs' argument holds a scalar,
	551	then that rest frequency will be applied to the selected
	552	data (and added to the list of available rest frequencies).
	553	In this way, you can set a rest frequency for each
	554	source and IF combination. If the 'freqs' argument holds
	555	a vector, then it MUST be of length the number of IFs
	556	(and the available restfrequencies will be replaced by
	557	this vector). In this case, all data ('source' and
	558	'theif' are ignored) have the restfrequency set per IF according
	559	to the corresponding value you give in the 'freqs' vector.
	560	E.g. 'freqs=[1e9,2e9]' would mean IF 0 gets restfreq 1e9 and
	561	IF 1 gets restfreq 2e9.
[402]	562
	563	You can also specify the frequencies via known line names
	564	in the argument 'lines'. Use 'freqs' or 'lines'. 'freqs'
	565	takes precedence. See the list of known names via function
	566	scantable.lines()
[391]	567	Parameters:
[402]	568	freqs: list of rest frequencies
[393]	569	unit: unit for rest frequency (default 'Hz')
[402]	570	lines: list of known spectral lines names (alternative to freqs).
	571	See possible list via scantable.lines()
[391]	572	source: Source name (blank means all)
	573	theif: IF (-1 means all)
	574	Example:
[393]	575	scan.set_restfreqs(freqs=1.4e9, source='NGC253', theif=2)
	576	scan.set_restfreqs(freqs=[1.4e9,1.67e9])
[391]	577	"""
[484]	578	varlist = vars()
[391]	579	if source is None:
	580	source = ""
	581	if theif is None:
	582	theif = -1
[393]	583	t = type(freqs)
	584	if t is int or t is float:
	585	freqs = [freqs]
[402]	586	if freqs is None:
	587	freqs = []
	588	t = type(lines)
	589	if t is str:
	590	lines = [lines]
	591	if lines is None:
	592	lines = []
	593	sdtable._setrestfreqs(self, freqs, unit, lines, source, theif)
[484]	594	self._add_history("set_restfreqs", varlist)
	595
[391]	596
	597
[102]	598	def flag_spectrum(self, thebeam, theif, thepol):
	599	"""
	600	This flags a selected spectrum in the scan 'for good'.
	601	USE WITH CARE - not reversible.
	602	Use masks for non-permanent exclusion of channels.
	603	Parameters:
	604	thebeam,theif,thepol: all have to be explicitly
	605	specified
	606	Example:
	607	scan.flag_spectrum(0,0,1)
	608	flags the spectrum for Beam=0, IF=0, Pol=1
	609	"""
[407]	610	if (thebeam < self.nbeam() and
	611	theif < self.nif() and
	612	thepol < self.npol()):
	613	sdtable.setbeam(self, thebeam)
	614	sdtable.setif(self, theif)
	615	sdtable.setpol(self, thepol)
	616	sdtable._flag(self)
[484]	617	self._add_history("flag_spectrum", vars())
[102]	618	else:
	619	print "Please specify a valid (Beam/IF/Pol)"
	620	return
[113]	621
	622	def plot(self, what='spectrum',col='Pol', panel=None):
	623	"""
	624	Plot the spectra contained in the scan. Alternatively you can also
	625	Plot Tsys vs Time
	626	Parameters:
	627	what: a choice of 'spectrum' (default) or 'tsys'
[135]	628	col: which out of Beams/IFs/Pols should be colour stacked
[113]	629	panel: set up multiple panels, currently not working.
	630	"""
[283]	631	print "Warning! Not fully functional. Use plotter.plot() instead"
	632
[113]	633	validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
[124]	634
[113]	635	validyax = ['spectrum','tsys']
[189]	636	from asap.asaplot import ASAPlot
[113]	637	if not self._p:
	638	self._p = ASAPlot()
[189]	639	#print "Plotting not enabled"
	640	#return
	641	if self._p.is_dead:
	642	del self._p
	643	self._p = ASAPlot()
[113]	644	npan = 1
	645	x = None
	646	if what == 'tsys':
	647	n = self.nrow()
	648	if n < 2:
	649	print "Only one integration. Can't plot."
	650	return
[124]	651	self._p.hold()
	652	self._p.clear()
[113]	653	if panel == 'Time':
	654	npan = self.nrow()
[124]	655	self._p.set_panels(rows=npan)
[113]	656	xlab,ylab,tlab = None,None,None
[226]	657	self._vb = False
[256]	658	sel = self.get_cursor()
[113]	659	for i in range(npan):
[226]	660	if npan > 1:
	661	self._p.subplot(i)
[113]	662	for j in range(validcol[col]):
	663	x = None
	664	y = None
	665	m = None
	666	tlab = self._getsourcename(i)
[124]	667	import re
	668	tlab = re.sub('_S','',tlab)
[113]	669	if col == 'Beam':
	670	self.setbeam(j)
	671	elif col == 'IF':
	672	self.setif(j)
	673	elif col == 'Pol':
	674	self.setpol(j)
	675	if what == 'tsys':
	676	x = range(self.nrow())
	677	xlab = 'Time [pixel]'
	678	m = list(ones(len(x)))
	679	y = []
	680	ylab = r'$T_{sys}$'
	681	for k in range(len(x)):
	682	y.append(self._gettsys(k))
	683	else:
[158]	684	x,xlab = self.get_abcissa(i)
[256]	685	y = self._getspectrum(i)
[113]	686	ylab = r'Flux'
[256]	687	m = self._getmask(i)
[113]	688	llab = col+' '+str(j)
	689	self._p.set_line(label=llab)
	690	self._p.plot(x,y,m)
[124]	691	self._p.set_axes('xlabel',xlab)
	692	self._p.set_axes('ylabel',ylab)
	693	self._p.set_axes('title',tlab)
[113]	694	self._p.release()
[256]	695	self.set_cursor(sel[0],sel[1],sel[2])
[226]	696	self._vb = rcParams['verbose']
[113]	697	return
[256]	698
	699	print out
	700
	701	def _print_values(self, dat, label='', timestamps=[]):
	702	d = dat['data']
	703	a = dat['axes']
	704	shp = d.getshape()
	705	out = ''
	706	for i in range(shp[3]):
	707	out += '%s [%s]:\n' % (a[3],timestamps[i])
	708	t = d[:,:,:,i]
	709	for j in range(shp[0]):
	710	if shp[0] > 1: out += ' %s[%d] ' % (a[0],j)
	711	for k in range(shp[1]):
	712	if shp[1] > 1: out += ' %s[%d] ' % (a[1],k)
	713	for l in range(shp[2]):
	714	if shp[2] > 1: out += ' %s[%d] ' % (a[2],l)
	715	out += '= %3.3f\n' % (t[j,k,l])
[484]	716	out += "-"*80
	717	out += "\n"
	718	print "-"*80
[256]	719	print " ", label
[484]	720	print "-"*80
[256]	721	print out
[484]	722
	723	def history(self):
	724	hist = list(self._gethistory())
	725	print "-"*80
	726	for h in hist:
[489]	727	if h.startswith("---"):
	728	print h
	729	else:
	730	items = h.split("##")
	731	date = items[0]
	732	func = items[1]
	733	items = items[2:]
	734	print date
	735	print "Function: %s\n Parameters:" % (func)
	736	for i in items:
	737	s = i.split("=")
	738	print " %s = %s" % (s[0],s[1])
	739	print "-"*80
[484]	740	return
	741
[513]	742	#
	743	# Maths business
	744	#
	745
[524]	746	def average_time(self, mask=None, scanav=True, weight='tint'):
[513]	747	"""
	748	Return the (time) average of a scan, or apply it 'insitu'.
	749	Note:
	750	in channels only
	751	The cursor of the output scan is set to 0.
	752	Parameters:
	753	one scan or comma separated scans
	754	mask: an optional mask (only used for 'var' and 'tsys'
	755	weighting)
	756	scanav: True (default) averages each scan separately
	757	False averages all scans together,
[524]	758	weight: Weighting scheme. 'none', 'var' (1/var(spec)
[523]	759	weighted), 'tsys' (1/Tsys**2 weighted), or 'tint'
[524]	760	(integration time weighted). The default is
	761	'tint'
[513]	762	Example:
	763	# time average the scantable without using a mask
	764	newscan = scan.average_time()
	765	"""
	766	varlist = vars()
[524]	767	if weight is None: weight = 'tint'
[513]	768	if mask is None: mask = ()
	769	from asap._asap import average as _av
	770	s = scantable(_av((self,), mask, scanav, weight))
	771	s._add_history("average_time",varlist)
	772	return s
	773
	774	def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None,
	775	allaxes=None):
	776	"""
	777	Return a scan where all spectra are converted to either
	778	Jansky or Kelvin depending upon the flux units of the scan table.
	779	By default the function tries to look the values up internally.
	780	If it can't find them (or if you want to over-ride), you must
	781	specify EITHER jyperk OR eta (and D which it will try to look up
	782	also if you don't set it). jyperk takes precedence if you set both.
	783	Parameters:
	784	jyperk: the Jy / K conversion factor
	785	eta: the aperture efficiency
	786	d: the geomtric diameter (metres)
	787	insitu: if False a new scantable is returned.
	788	Otherwise, the scaling is done in-situ
	789	The default is taken from .asaprc (False)
	790	allaxes: if True apply to all spectra. Otherwise
	791	apply only to the selected (beam/pol/if)spectra only
	792	The default is taken from .asaprc (True if none)
	793	"""
	794	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	795	if insitu is None: insitu = rcParams['insitu']
	796	varlist = vars()
	797	if jyperk is None: jyperk = -1.0
	798	if d is None: d = -1.0
	799	if eta is None: eta = -1.0
	800	if not insitu:
	801	from asap._asap import convertflux as _convert
	802	s = scantable(_convert(self, d, eta, jyperk, allaxes))
	803	s._add_history("convert_flux", varlist)
	804	return s
	805	else:
	806	from asap._asap import convertflux_insitu as _convert
	807	_convert(self, d, eta, jyperk, allaxes)
	808	self._add_history("convert_flux", varlist)
	809	return
	810
	811	def gain_el(self, poly=None, filename="", method="linear",
	812	insitu=None, allaxes=None):
	813	"""
	814	Return a scan after applying a gain-elevation correction.
	815	The correction can be made via either a polynomial or a
	816	table-based interpolation (and extrapolation if necessary).
	817	You specify polynomial coefficients, an ascii table or neither.
	818	If you specify neither, then a polynomial correction will be made
	819	with built in coefficients known for certain telescopes (an error
	820	will occur if the instrument is not known).
	821	The data and Tsys are divided by the scaling factors.
	822	Parameters:
	823	poly: Polynomial coefficients (default None) to compute a
	824	gain-elevation correction as a function of
	825	elevation (in degrees).
	826	filename: The name of an ascii file holding correction factors.
	827	The first row of the ascii file must give the column
	828	names and these MUST include columns
	829	"ELEVATION" (degrees) and "FACTOR" (multiply data
	830	by this) somewhere.
	831	The second row must give the data type of the
	832	column. Use 'R' for Real and 'I' for Integer.
	833	An example file would be
	834	(actual factors are arbitrary) :
	835
	836	TIME ELEVATION FACTOR
	837	R R R
	838	0.1 0 0.8
	839	0.2 20 0.85
	840	0.3 40 0.9
	841	0.4 60 0.85
	842	0.5 80 0.8
	843	0.6 90 0.75
	844	method: Interpolation method when correcting from a table.
	845	Values are "nearest", "linear" (default), "cubic"
	846	and "spline"
	847	insitu: if False a new scantable is returned.
	848	Otherwise, the scaling is done in-situ
	849	The default is taken from .asaprc (False)
	850	allaxes: If True apply to all spectra. Otherwise
	851	apply only to the selected (beam/pol/if) spectra only
	852	The default is taken from .asaprc (True if none)
	853	"""
	854
	855	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	856	if insitu is None: insitu = rcParams['insitu']
	857	varlist = vars()
	858	if poly is None:
	859	poly = ()
	860	from os.path import expandvars
	861	filename = expandvars(filename)
	862	if not insitu:
	863	from asap._asap import gainel as _gainEl
	864	s = scantable(_gainEl(self, poly, filename, method, allaxes))
	865	s._add_history("gain_el", varlist)
	866	return s
	867	else:
	868	from asap._asap import gainel_insitu as _gainEl
	869	_gainEl(self, poly, filename, method, allaxes)
	870	self._add_history("gain_el", varlist)
	871	return
	872
	873	def freq_align(self, reftime=None, method='cubic', perif=False,
	874	insitu=None):
	875	"""
	876	Return a scan where all rows have been aligned in frequency/velocity.
	877	The alignment frequency frame (e.g. LSRK) is that set by function
	878	set_freqframe.
	879	Parameters:
	880	reftime: reference time to align at. By default, the time of
	881	the first row of data is used.
	882	method: Interpolation method for regridding the spectra.
	883	Choose from "nearest", "linear", "cubic" (default)
	884	and "spline"
	885	perif: Generate aligners per freqID (no doppler tracking) or
	886	per IF (scan-based doppler tracking)
	887	insitu: if False a new scantable is returned.
	888	Otherwise, the scaling is done in-situ
	889	The default is taken from .asaprc (False)
	890	"""
	891	if insitu is None: insitu = rcParams['insitu']
	892	varlist = vars()
	893	if reftime is None: reftime = ''
	894	perfreqid = not perif
	895	if not insitu:
	896	from asap._asap import freq_align as _align
	897	s = scantable(_align(self, reftime, method, perfreqid))
	898	s._add_history("freq_align", varlist)
	899	return s
	900	else:
	901	from asap._asap import freq_align_insitu as _align
	902	_align(self, reftime, method, perfreqid)
	903	self._add_history("freq_align", varlist)
	904	return
	905
	906	def opacity(self, tau, insitu=None, allaxes=None):
	907	"""
	908	Apply an opacity correction. The data
	909	and Tsys are multiplied by the correction factor.
	910	Parameters:
	911	tau: Opacity from which the correction factor is
	912	exp(tau*ZD)
	913	where ZD is the zenith-distance
	914	insitu: if False a new scantable is returned.
	915	Otherwise, the scaling is done in-situ
	916	The default is taken from .asaprc (False)
	917	allaxes: if True apply to all spectra. Otherwise
	918	apply only to the selected (beam/pol/if)spectra only
	919	The default is taken from .asaprc (True if none)
	920	"""
	921	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	922	if insitu is None: insitu = rcParams['insitu']
	923	varlist = vars()
	924	if not insitu:
	925	from asap._asap import opacity as _opacity
	926	s = scantable(_opacity(self, tau, allaxes))
	927	s._add_history("opacity", varlist)
	928	return s
	929	else:
	930	from asap._asap import opacity_insitu as _opacity
	931	_opacity(self, tau, allaxes)
	932	self._add_history("opacity", varlist)
	933	return
	934
	935	def bin(self, width=5, insitu=None):
	936	"""
	937	Return a scan where all spectra have been binned up.
	938	width: The bin width (default=5) in pixels
	939	insitu: if False a new scantable is returned.
	940	Otherwise, the scaling is done in-situ
	941	The default is taken from .asaprc (False)
	942	"""
	943	if insitu is None: insitu = rcParams['insitu']
	944	varlist = vars()
	945	if not insitu:
	946	from asap._asap import bin as _bin
	947	s = scantable(_bin(self, width))
	948	s._add_history("bin",varlist)
	949	return s
	950	else:
	951	from asap._asap import bin_insitu as _bin
	952	_bin(self, width)
	953	self._add_history("bin",varlist)
	954	return
	955
	956
	957	def resample(self, width=5, method='cubic', insitu=None):
	958	"""
	959	Return a scan where all spectra have been binned up
	960	width: The bin width (default=5) in pixels
	961	method: Interpolation method when correcting from a table.
	962	Values are "nearest", "linear", "cubic" (default)
	963	and "spline"
	964	insitu: if False a new scantable is returned.
	965	Otherwise, the scaling is done in-situ
	966	The default is taken from .asaprc (False)
	967	"""
	968	if insitu is None: insitu = rcParams['insitu']
	969	varlist = vars()
	970	if not insitu:
	971	from asap._asap import resample as _resample
	972	s = scantable(_resample(self, method, width))
	973	s._add_history("resample",varlist)
	974	return s
	975	else:
	976	from asap._asap import resample_insitu as _resample
	977	_resample(self, method, width)
	978	self._add_history("resample",varlist)
	979	return
	980
	981	def average_pol(self, mask=None, weight='none', insitu=None):
	982	"""
	983	Average the Polarisations together.
	984	The polarisation cursor of the output scan is set to 0
	985	Parameters:
	986	scan: The scantable
	987	mask: An optional mask defining the region, where the
	988	averaging will be applied. The output will have all
	989	specified points masked.
	990	weight: Weighting scheme. 'none' (default), or
	991	'var' (variance weighted)
	992	insitu: if False a new scantable is returned.
	993	Otherwise, the scaling is done in-situ
	994	The default is taken from .asaprc (False)
	995	Example:
	996	polav = average_pols(myscan)
	997	"""
	998	if insitu is None: insitu = rcParams['insitu']
	999	varlist = vars()
	1000
	1001	if mask is None:
	1002	mask = ()
	1003	if not insitu:
	1004	from asap._asap import averagepol as _avpol
	1005	s = scantable(_avpol(self, mask, weight))
	1006	s._add_history("average_pol",varlist)
	1007	return s
	1008	else:
	1009	from asap._asap import averagepol_insitu as _avpol
	1010	_avpol(self, mask, weight)
	1011	self._add_history("average_pol",varlist)
	1012	return
	1013
	1014	def smooth(self, kernel="hanning", width=5.0, insitu=None, allaxes=None):
	1015	"""
	1016	Smooth the spectrum by the specified kernel (conserving flux).
	1017	Parameters:
	1018	scan: The input scan
	1019	kernel: The type of smoothing kernel. Select from
	1020	'hanning' (default), 'gaussian' and 'boxcar'.
	1021	The first three characters are sufficient.
	1022	width: The width of the kernel in pixels. For hanning this is
	1023	ignored otherwise it defauls to 5 pixels.
	1024	For 'gaussian' it is the Full Width Half
	1025	Maximum. For 'boxcar' it is the full width.
	1026	insitu: if False a new scantable is returned.
	1027	Otherwise, the scaling is done in-situ
	1028	The default is taken from .asaprc (False)
	1029	allaxes: If True (default) apply to all spectra. Otherwise
	1030	apply only to the selected (beam/pol/if)spectra only
	1031	The default is taken from .asaprc (True if none)
	1032	Example:
	1033	none
	1034	"""
	1035	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	1036	if insitu is None: insitu = rcParams['insitu']
	1037	varlist = vars()
	1038	if not insitu:
	1039	from asap._asap import smooth as _smooth
	1040	s = scantable(_smooth(self,kernel,width,allaxes))
	1041	s._add_history("smooth", varlist)
	1042	return s
	1043	else:
	1044	from asap._asap import smooth_insitu as _smooth
	1045	_smooth(self,kernel,width,allaxes)
	1046	self._add_history("smooth", varlist)
	1047	return
	1048
	1049	def poly_baseline(self, mask=None, order=0, insitu=None):
	1050	"""
	1051	Return a scan which has been baselined (all rows) by a polynomial.
	1052	Parameters:
	1053	scan: a scantable
	1054	mask: an optional mask
	1055	order: the order of the polynomial (default is 0)
	1056	insitu: if False a new scantable is returned.
	1057	Otherwise, the scaling is done in-situ
	1058	The default is taken from .asaprc (False)
	1059	Example:
	1060	# return a scan baselined by a third order polynomial,
	1061	# not using a mask
	1062	bscan = scan.poly_baseline(order=3)
	1063	"""
	1064	if insitu is None: insitu = rcParams['insitu']
	1065	varlist = vars()
	1066	if mask is None:
	1067	from numarray import ones
	1068	mask = list(ones(scan.nchan()))
	1069	from asap.asapfitter import fitter
	1070	f = fitter()
	1071	f._verbose(True)
	1072	f.set_scan(self, mask)
	1073	f.set_function(poly=order)
	1074	sf = f.auto_fit(insitu)
	1075	if insitu:
	1076	self._add_history("poly_baseline", varlist)
	1077	return
	1078	else:
	1079	sf._add_history("poly_baseline", varlist)
	1080	return sf
	1081
	1082	def auto_poly_baseline(self, mask=None, edge=(0,0), order=0,
	1083	threshold=3,insitu=None):
	1084	"""
	1085	Return a scan which has been baselined (all rows) by a polynomial.
	1086	Spectral lines are detected first using linefinder and masked out
	1087	to avoid them affecting the baseline solution.
	1088
	1089	Parameters:
	1090	scan: a scantable
	1091	mask: an optional mask retreived from scantable
	1092	edge: an optional number of channel to drop at
	1093	the edge of spectrum. If only one value is
	1094	specified, the same number will be dropped from
	1095	both sides of the spectrum. Default is to keep
	1096	all channels
	1097	order: the order of the polynomial (default is 0)
	1098	threshold: the threshold used by line finder. It is better to
	1099	keep it large as only strong lines affect the
	1100	baseline solution.
	1101	insitu: if False a new scantable is returned.
	1102	Otherwise, the scaling is done in-situ
	1103	The default is taken from .asaprc (False)
	1104
	1105	Example:
[516]	1106	scan2=scan.auto_poly_baseline(order=7)
[513]	1107	"""
	1108	if insitu is None: insitu = rcParams['insitu']
	1109	varlist = vars()
	1110	from asap.asapfitter import fitter
	1111	from asap.asaplinefind import linefinder
	1112	from asap import _is_sequence_or_number as _is_valid
	1113
	1114	if not _is_valid(edge, int):
[516]	1115	raise RuntimeError, "Parameter 'edge' has to be an integer or a \
	1116	pair of integers specified as a tuple"
[513]	1117
	1118	# setup fitter
	1119	f = fitter()
	1120	f._verbose(True)
	1121	f.set_function(poly=order)
	1122
	1123	# setup line finder
	1124	fl=linefinder()
	1125	fl.set_options(threshold=threshold)
	1126
	1127	if not insitu:
	1128	workscan=self.copy()
	1129	else:
	1130	workscan=self
	1131
	1132	vb=workscan._vb
	1133	# remember the verbose parameter and selection
	1134	workscan._vb=False
	1135	sel=workscan.get_cursor()
	1136	rows=range(workscan.nrow())
	1137	for i in range(workscan.nbeam()):
	1138	workscan.setbeam(i)
	1139	for j in range(workscan.nif()):
	1140	workscan.setif(j)
	1141	for k in range(workscan.npol()):
	1142	workscan.setpol(k)
	1143	if f._vb:
	1144	print "Processing:"
	1145	print 'Beam[%d], IF[%d], Pol[%d]' % (i,j,k)
	1146	for iRow in rows:
	1147	fl.set_scan(workscan,mask,edge)
	1148	fl.find_lines(iRow)
	1149	f.set_scan(workscan, fl.get_mask())
	1150	f.x=workscan._getabcissa(iRow)
	1151	f.y=workscan._getspectrum(iRow)
	1152	f.data=None
	1153	f.fit()
	1154	x=f.get_parameters()
	1155	workscan._setspectrum(f.fitter.getresidual(),iRow)
	1156	workscan.set_cursor(sel[0],sel[1],sel[2])
	1157	workscan._vb = vb
	1158	if not insitu:
	1159	return workscan
	1160
	1161	def rotate_linpolphase(self, angle, allaxes=None):
	1162	"""
	1163	Rotate the phase of the complex polarization O=Q+iU correlation.
	1164	This is always done in situ in the raw data. So if you call this
	1165	function more than once then each call rotates the phase further.
	1166	Parameters:
	1167	angle: The angle (degrees) to rotate (add) by.
	1168	allaxes: If True apply to all spectra. Otherwise
	1169	apply only to the selected (beam/pol/if)spectra only.
	1170	The default is taken from .asaprc (True if none)
	1171	Examples:
	1172	scan.rotate_linpolphase(2.3)
	1173	"""
	1174	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	1175	varlist = vars()
	1176	from asap._asap import _rotate_linpolphase as _rotate
	1177	_rotate(self, angle, allaxes)
	1178	self._add_history("rotate_linpolphase", varlist)
	1179	return
	1180
	1181
	1182	def rotate_xyphase(self, angle, allaxes=None):
	1183	"""
	1184	Rotate the phase of the XY correlation. This is always done in situ
	1185	in the data. So if you call this function more than once
	1186	then each call rotates the phase further.
	1187	Parameters:
	1188	angle: The angle (degrees) to rotate (add) by.
	1189	allaxes: If True apply to all spectra. Otherwise
	1190	apply only to the selected (beam/pol/if)spectra only.
	1191	The default is taken from .asaprc (True if none)
	1192	Examples:
	1193	scan.rotate_xyphase(2.3)
	1194	"""
	1195	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	1196	varlist = vars()
	1197	from asap._asap import rotate_xyphase as _rotate_xyphase
	1198	_rotate_xyphase(self, angle, allaxes)
	1199	self._add_history("rotate_xyphase", varlist)
	1200	return
	1201
	1202
	1203	def add(self, offset, insitu=None, allaxes=None):
	1204	"""
	1205	Return a scan where all spectra have the offset added
	1206	Parameters:
	1207	offset: the offset
	1208	insitu: if False a new scantable is returned.
	1209	Otherwise, the scaling is done in-situ
	1210	The default is taken from .asaprc (False)
	1211	allaxes: if True apply to all spectra. Otherwise
	1212	apply only to the selected (beam/pol/if)spectra only
	1213	The default is taken from .asaprc (True if none)
	1214	"""
	1215	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	1216	if insitu is None: insitu = rcParams['insitu']
	1217	varlist = vars()
	1218	if not insitu:
	1219	from asap._asap import add as _add
	1220	s = scantable(_add(self, offset, allaxes))
	1221	s._add_history("add",varlist)
	1222	return s
	1223	else:
	1224	from asap._asap import add_insitu as _add
	1225	_add(self, offset, allaxes)
	1226	self._add_history("add",varlist)
	1227	return
	1228
	1229	def scale(self, factor, insitu=None, allaxes=None, tsys=True):
	1230	"""
	1231	Return a scan where all spectra are scaled by the give 'factor'
	1232	Parameters:
	1233	factor: the scaling factor
	1234	insitu: if False a new scantable is returned.
	1235	Otherwise, the scaling is done in-situ
	1236	The default is taken from .asaprc (False)
	1237	allaxes: if True apply to all spectra. Otherwise
	1238	apply only to the selected (beam/pol/if)spectra only.
	1239	The default is taken from .asaprc (True if none)
	1240	tsys: if True (default) then apply the operation to Tsys
	1241	as well as the data
	1242	"""
	1243	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	1244	if insitu is None: insitu = rcParams['insitu']
	1245	varlist = vars()
	1246	if not insitu:
	1247	from asap._asap import scale as _scale
	1248	s = scantable(_scale(self, factor, allaxes, tsys))
	1249	s._add_history("scale",varlist)
	1250	return s
	1251	else:
	1252	from asap._asap import scale_insitu as _scale
	1253	_scale(self, factor, allaxes, tsys)
	1254	self._add_history("scale",varlist)
	1255	return
	1256
	1257
	1258	def quotient(self, other, isreference=True, preserve=True):
	1259	"""
	1260	Return the quotient of a 'source' (signal) scan and a 'reference' scan.
	1261	The reference can have just one row, even if the signal has many.
	1262	Otherwise they must have the same number of rows.
	1263	The cursor of the output scan is set to 0
	1264	Parameters:
	1265	source: the 'other' scan
	1266	isreference: if the 'other' scan is the reference (default)
	1267	or source
	1268	preserve: you can preserve (default) the continuum or
	1269	remove it. The equations used are
	1270	preserve - Output = Tref * (sig/ref) - Tref
	1271	remove - Output = Tref * (sig/ref) - Tsig
	1272	Example:
	1273	# src is a scantable for an 'on' scan, ref for an 'off' scantable
	1274	q1 = src.quotient(ref)
	1275	q2 = ref.quotient(src, isreference=False)
	1276	# gives the same result
	1277	"""
	1278	from asap._asap import quotient as _quot
	1279	if isreference:
	1280	return scantable(_quot(self, other, preserve))
	1281	else:
	1282	return scantable(_quot(other, self, preserve))
	1283
	1284	def __add__(self, other):
	1285	varlist = vars()
	1286	s = None
	1287	if isinstance(other, scantable):
	1288	from asap._asap import b_operate as _bop
	1289	s = scantable(_bop(self, other, 'add', True))
	1290	elif isinstance(other, float):
	1291	from asap._asap import add as _add
	1292	s = scantable(_add(self, other, True))
	1293	else:
	1294	print "Other input is not a scantable or float value"
	1295	return
	1296	s._add_history("operator +", varlist)
	1297	return s
	1298
	1299
	1300	def __sub__(self, other):
	1301	"""
	1302	implicit on all axes and on Tsys
	1303	"""
	1304	varlist = vars()
	1305	s = None
	1306	if isinstance(other, scantable):
	1307	from asap._asap import b_operate as _bop
	1308	s = scantable(_bop(self, other, 'sub', True))
	1309	elif isinstance(other, float):
	1310	from asap._asap import add as _add
	1311	s = scantable(_add(self, -other, True))
	1312	else:
	1313	print "Other input is not a scantable or float value"
	1314	return
	1315	s._add_history("operator -", varlist)
	1316	return s
	1317
	1318	def __mul__(self, other):
	1319	"""
	1320	implicit on all axes and on Tsys
	1321	"""
	1322	varlist = vars()
	1323	s = None
	1324	if isinstance(other, scantable):
	1325	from asap._asap import b_operate as _bop
	1326	s = scantable(_bop(self, other, 'mul', True))
	1327	elif isinstance(other, float):
	1328	if other == 0.0:
	1329	print "Multiplying by zero is not recommended."
	1330	return
	1331	from asap._asap import scale as _sca
	1332	s = scantable(_sca(self, other, True))
	1333	else:
	1334	print "Other input is not a scantable or float value"
	1335	return
	1336	s._add_history("operator *", varlist)
	1337	return s
	1338
	1339
	1340	def __div__(self, other):
	1341	"""
	1342	implicit on all axes and on Tsys
	1343	"""
	1344	varlist = vars()
	1345	s = None
	1346	if isinstance(other, scantable):
	1347	from asap._asap import b_operate as _bop
	1348	s = scantable(_bop(self, other, 'div', True))
	1349	elif isinstance(other, float):
	1350	if other == 0.0:
	1351	print "Dividing by zero is not recommended"
	1352	return
	1353	from asap._asap import scale as _sca
	1354	s = scantable(_sca(self, 1.0/other, True))
	1355	else:
	1356	print "Other input is not a scantable or float value"
	1357	return
	1358	s._add_history("operator /", varlist)
	1359	return s
	1360
[484]	1361	def _add_history(self, funcname, parameters):
	1362	# create date
	1363	sep = "##"
	1364	from datetime import datetime
	1365	dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
	1366	hist = dstr+sep
	1367	hist += funcname+sep#cdate+sep
	1368	if parameters.has_key('self'): del parameters['self']
	1369	for k,v in parameters.iteritems():
	1370	if type(v) is dict:
	1371	for k2,v2 in v.iteritems():
	1372	hist += k2
	1373	hist += "="
	1374	if isinstance(v2,scantable):
	1375	hist += 'scantable'
	1376	elif k2 == 'mask':
[513]	1377	if isinstance(v2,list) or isinstance(v2,tuple):
	1378	hist += str(self._zip_mask(v2))
	1379	else:
	1380	hist += str(v2)
[484]	1381	else:
[513]	1382	hist += str(v2)
[484]	1383	else:
	1384	hist += k
	1385	hist += "="
	1386	if isinstance(v,scantable):
	1387	hist += 'scantable'
	1388	elif k == 'mask':
[513]	1389	if isinstance(v,list) or isinstance(v,tuple):
	1390	hist += str(self._zip_mask(v))
	1391	else:
	1392	hist += str(v)
[484]	1393	else:
	1394	hist += str(v)
	1395	hist += sep
	1396	hist = hist[:-2] # remove trailing '##'
	1397	self._addhistory(hist)
	1398
	1399
	1400	def _zip_mask(self, mask):
	1401	mask = list(mask)
	1402	i = 0
	1403	segments = []
	1404	while mask[i:].count(1):
	1405	i += mask[i:].index(1)
	1406	if mask[i:].count(0):
	1407	j = i + mask[i:].index(0)
	1408	else:
	1409	j = len(mask)
	1410	segments.append([i,j])
[513]	1411	i = j
[484]	1412	return segments

Note: See TracBrowser for help on using the repository browser.

Download in other formats: