Context Navigation

source: trunk/python/scantable.py@ 713

Visit:

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

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

Download in other formats: