Context Navigation

source: trunk/python/scantable.py@ 715

Visit:

Last change on this file since 715 was 714, checked in by mar637, 19 years ago
updated logging; added public get_sourcename
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 58.1 KB

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

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

Download in other formats: