Context Navigation

source: trunk/python/scantable.py@ 790

Visit:

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

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

Download in other formats: