source: branches/alma/python/asapmath.py@ 1453

Last change on this file since 1453 was 1446, checked in by TakTsutsumi, 16 years ago

Merged recent updates (since 2007) from nrao-asap

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.8 KB
RevLine 
[1085]1from asap.scantable import scantable
[258]2from asap import rcParams
[720]3from asap import print_log
[1389]4from asap import selector
[101]5
[143]6def average_time(*args, **kwargs):
[101]7 """
[113]8 Return the (time) average of a scan or list of scans. [in channels only]
[305]9 The cursor of the output scan is set to 0
[113]10 Parameters:
[1361]11 one scan or comma separated scans or a list of scans
[143]12 mask: an optional mask (only used for 'var' and 'tsys' weighting)
[558]13 scanav: True averages each scan separately.
14 False (default) averages all scans together,
[1232]15 weight: Weighting scheme.
16 'none' (mean no weight)
17 'var' (1/var(spec) weighted)
18 'tsys' (1/Tsys**2 weighted)
19 'tint' (integration time weighted)
20 'tintsys' (Tint/Tsys**2)
21 'median' ( median averaging)
[930]22 align: align the spectra in velocity before averaging. It takes
23 the time of the first spectrum in the first scantable
24 as reference time.
[113]25 Example:
26 # return a time averaged scan from scana and scanb
27 # without using a mask
[129]28 scanav = average_time(scana,scanb)
[1361]29 # or equivalent
30 # scanav = average_time([scana, scanb])
[113]31 # return the (time) averaged scan, i.e. the average of
32 # all correlator cycles
[558]33 scanav = average_time(scan, scanav=True)
[101]34 """
[930]35 scanav = False
[143]36 if kwargs.has_key('scanav'):
[930]37 scanav = kwargs.get('scanav')
[524]38 weight = 'tint'
[143]39 if kwargs.has_key('weight'):
40 weight = kwargs.get('weight')
41 mask = ()
42 if kwargs.has_key('mask'):
43 mask = kwargs.get('mask')
[930]44 align = False
45 if kwargs.has_key('align'):
46 align = kwargs.get('align')
[1446]47 compel = False
48 if kwargs.has_key('compel'):
49 compel = kwargs.get('compel')
[489]50 varlist = vars()
[665]51 if isinstance(args[0],list):
[981]52 lst = args[0]
[665]53 elif isinstance(args[0],tuple):
[981]54 lst = list(args[0])
[665]55 else:
[981]56 lst = list(args)
[720]57
[489]58 del varlist["kwargs"]
59 varlist["args"] = "%d scantables" % len(lst)
[981]60 # need special formatting here for history...
[720]61
[876]62 from asap._asap import stmath
63 stm = stmath()
[113]64 for s in lst:
[101]65 if not isinstance(s,scantable):
[720]66 msg = "Please give a list of scantables"
67 if rcParams['verbose']:
68 print msg
69 return
70 else:
71 raise TypeError(msg)
[945]72 if scanav: scanav = "SCAN"
73 else: scanav = "NONE"
[981]74 alignedlst = []
75 if align:
76 refepoch = lst[0].get_time(0)
77 for scan in lst:
78 alignedlst.append(scan.freq_align(refepoch,insitu=False))
79 else:
[1059]80 alignedlst = lst
[1232]81 if weight.upper() == 'MEDIAN':
82 # median doesn't support list of scantables - merge first
83 merged = None
84 if len(alignedlst) > 1:
85 merged = merge(alignedlst)
86 else:
87 merged = alignedlst[0]
88 s = scantable(stm._averagechannel(merged, 'MEDIAN', scanav))
89 del merged
90 else:
[1446]91 #s = scantable(stm._average(alignedlst, mask, weight.upper(), scanav))
92 s = scantable(stm._new_average(alignedlst, compel, mask, weight.upper(), scanav))
[489]93 s._add_history("average_time",varlist)
[720]94 print_log()
[489]95 return s
[101]96
[1074]97def quotient(source, reference, preserve=True):
98 """
99 Return the quotient of a 'source' (signal) scan and a 'reference' scan.
100 The reference can have just one scan, even if the signal has many. Otherwise
101 they must have the same number of scans.
102 The cursor of the output scan is set to 0
103 Parameters:
104 source: the 'on' scan
105 reference: the 'off' scan
106 preserve: you can preserve (default) the continuum or
107 remove it. The equations used are
108 preserve: Output = Toff * (on/off) - Toff
109 remove: Output = Toff * (on/off) - Ton
110 """
111 varlist = vars()
112 from asap._asap import stmath
113 stm = stmath()
114 stm._setinsitu(False)
115 s = scantable(stm._quotient(source, reference, preserve))
116 s._add_history("quotient",varlist)
117 print_log()
118 return s
[101]119
[1389]120def dototalpower(calon, caloff, tcalval=0.0):
121 """
122 Do calibration for CAL on,off signals.
123 Adopted from GBTIDL dototalpower
124 Parameters:
125 calon: the 'cal on' subintegration
126 caloff: the 'cal off' subintegration
127 tcalval: user supplied Tcal value
128 """
129 varlist = vars()
130 from asap._asap import stmath
131 stm = stmath()
132 stm._setinsitu(False)
133 s = scantable(stm._dototalpower(calon, caloff, tcalval))
134 s._add_history("dototalpower",varlist)
135 print_log()
136 return s
137
138def dosigref(sig, ref, smooth, tsysval=0.0, tauval=0.0):
139 """
140 Calculate a quotient (sig-ref/ref * Tsys)
141 Adopted from GBTIDL dosigref
142 Parameters:
143 sig: on source data
144 ref: reference data
145 smooth: width of box car smoothing for reference
146 tsysval: user specified Tsys (scalar only)
147 tauval: user specified Tau (required if tsysval is set)
148 """
149 varlist = vars()
150 from asap._asap import stmath
151 stm = stmath()
152 stm._setinsitu(False)
153 s = scantable(stm._dosigref(sig, ref, smooth, tsysval, tauval))
154 s._add_history("dosigref",varlist)
155 print_log()
156 return s
157
158def calps(scantab, scannos, smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
159 """
160 Calibrate GBT position switched data
161 Adopted from GBTIDL getps
162 Currently calps identify the scans as position switched data if they
163 contain '_ps' in the source name. The data must contains 'CAL' signal
164 on/off in each integration. To identify 'CAL' on state, the word, 'calon'
165 need to be present in the source name field.
166 (GBT MS data reading process to scantable automatically append these
167 id names to the source names)
168
169 Parameters:
170 scantab: scantable
171 scannos: list of scan numbers
172 smooth: optional box smoothing order for the reference
173 (default is 1 = no smoothing)
174 tsysval: optional user specified Tsys (default is 0.0,
175 use Tsys in the data)
176 tauval: optional user specified Tau
177 tcalval: optional user specified Tcal (default is 0.0,
178 use Tcal value in the data)
179 """
180 varlist = vars()
181 # check for the appropriate data
182 s = scantab.get_scan('*_ps*')
183 if s is None:
184 msg = "The input data appear to contain no position-switch mode data."
185 if rcParams['verbose']:
186 print msg
187 return
188 else:
189 raise TypeError(msg)
190 ssub = s.get_scan(scannos)
191 if ssub is None:
192 msg = "No data was found with given scan numbers!"
193 if rcParams['verbose']:
194 print msg
195 return
196 else:
197 raise TypeError(msg)
198 ssubon = ssub.get_scan('*calon')
199 ssuboff = ssub.get_scan('*[^calon]')
200 if ssubon.nrow() != ssuboff.nrow():
201 msg = "mismatch in numbers of CAL on/off scans. Cannot calibrate. Check the scan numbers."
202 if rcParams['verbose']:
203 print msg
204 return
205 else:
206 raise TypeError(msg)
207 cals = dototalpower(ssubon, ssuboff, tcalval)
208 sig = cals.get_scan('*ps')
209 ref = cals.get_scan('*psr')
210 if sig.nscan() != ref.nscan():
211 msg = "mismatch in numbers of on/off scans. Cannot calibrate. Check the scan numbers."
212 if rcParams['verbose']:
213 print msg
214 return
215 else:
216 raise TypeError(msg)
217
218 #for user supplied Tsys
219 if tsysval>0.0:
220 if tauval<=0.0:
221 msg = "Need to supply a valid tau to use the supplied Tsys"
222 if rcParams['verbose']:
223 print msg
224 return
225 else:
226 raise TypeError(msg)
227 else:
228 sig.recalc_azel()
229 ref.recalc_azel()
230 #msg = "Use of user specified Tsys is not fully implemented yet."
231 #if rcParams['verbose']:
232 # print msg
233 # return
234 #else:
235 # raise TypeError(msg)
236 # use get_elevation to get elevation and
237 # calculate a scaling factor using the formula
238 # -> tsys use to dosigref
239
240 #ress = dosigref(sig, ref, smooth, tsysval)
241 ress = dosigref(sig, ref, smooth, tsysval, tauval)
242 ress._add_history("calps", varlist)
243 print_log()
244 return ress
245
246def calnod(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
247 """
248 Do full (but a pair of scans at time) processing of GBT Nod data
249 calibration.
250 Adopted from GBTIDL's getnod
251 Parameters:
252 scantab: scantable
253 scannos: a pair of scan numbers, or the first scan number of the pair
254 smooth: box car smoothing order
255 tsysval: optional user specified Tsys value
256 tauval: optional user specified tau value (not implemented yet)
257 tcalval: optional user specified Tcal value
258 """
259 varlist = vars()
260 from asap._asap import stmath
261 stm = stmath()
262 stm._setinsitu(False)
263
264 # check for the appropriate data
265 s = scantab.get_scan('*_nod*')
266 if s is None:
267 msg = "The input data appear to contain no Nod observing mode data."
268 if rcParams['verbose']:
269 print msg
270 return
271 else:
272 raise TypeError(msg)
273
274 # need check correspondance of each beam with sig-ref ...
275 # check for timestamps, scan numbers, subscan id (not available in
276 # ASAP data format...). Assume 1st scan of the pair (beam 0 - sig
277 # and beam 1 - ref...)
278 # First scan number of paired scans or list of pairs of
279 # scan numbers (has to have even number of pairs.)
280
281 #data splitting
282 scan1no = scan2no = 0
283
284 if len(scannos)==1:
285 scan1no = scannos[0]
286 scan2no = scannos[0]+1
287 pairScans = [scan1no, scan2no]
288 else:
289 #if len(scannos)>2:
290 # msg = "calnod can only process a pair of nod scans at time."
291 # if rcParams['verbose']:
292 # print msg
293 # return
294 # else:
295 # raise TypeError(msg)
296 #
297 #if len(scannos)==2:
298 # scan1no = scannos[0]
299 # scan2no = scannos[1]
300 pairScans = list(scannos)
301
302 if tsysval>0.0:
303 if tauval<=0.0:
304 msg = "Need to supply a valid tau to use the supplied Tsys"
305 if rcParams['verbose']:
306 print msg
307 return
308 else:
309 raise TypeError(msg)
310 else:
311 scantab.recalc_azel()
312 resspec = scantable(stm._donod(scantab, pairScans, smooth, tsysval,tauval,tcalval))
313 resspec._add_history("calnod",varlist)
314 print_log()
315 return resspec
316
317def calfs(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
318 """
319 Calibrate GBT frequency switched data.
320 Adopted from GBTIDL getfs.
321 Currently calfs identify the scans as frequency switched data if they
322 contain '_fs' in the source name. The data must contains 'CAL' signal
323 on/off in each integration. To identify 'CAL' on state, the word, 'calon'
324 need to be present in the source name field.
325 (GBT MS data reading via scantable automatically append these
326 id names to the source names)
327
328 Parameters:
329 scantab: scantable
330 scannos: list of scan numbers
331 smooth: optional box smoothing order for the reference
332 (default is 1 = no smoothing)
333 tsysval: optional user specified Tsys (default is 0.0,
334 use Tsys in the data)
335 tauval: optional user specified Tau
336 """
337 varlist = vars()
338 from asap._asap import stmath
339 stm = stmath()
340 stm._setinsitu(False)
341
342# check = scantab.get_scan('*_fs*')
343# if check is None:
344# msg = "The input data appear to contain no Nod observing mode data."
345# if rcParams['verbose']:
346# print msg
347# return
348# else:
349# raise TypeError(msg)
350 s = scantab.get_scan(scannos)
351 del scantab
352
353 resspec = scantable(stm._dofs(s, scannos, smooth, tsysval,tauval,tcalval))
354 resspec._add_history("calfs",varlist)
355 print_log()
356 return resspec
357
[296]358def simple_math(left, right, op='add', tsys=True):
[242]359 """
[720]360 Apply simple mathematical binary operations to two
[242]361 scan tables, returning the result in a new scan table.
362 The operation is applied to both the correlations and the TSys data
[305]363 The cursor of the output scan is set to 0
[242]364 Parameters:
365 left: the 'left' scan
366 right: the 'right' scan
367 op: the operation: 'add' (default), 'sub', 'mul', 'div'
[296]368 tsys: if True (default) then apply the operation to Tsys
369 as well as the data
[242]370 """
[1357]371 print "simple_math is deprecated use +=/* instead."
[918]372
373def merge(*args):
[945]374 """
[1362]375 Merge a list of scanatables, or comma-sperated scantables into one
376 scnatble.
377 Parameters:
378 A list [scan1, scan2] or scan1, scan2.
379 Example:
380 myscans = [scan1, scan2]
381 allscans = merge(myscans)
382 # or equivalent
383 sameallscans = merge(scan1, scan2)
[945]384 """
[918]385 varlist = vars()
386 if isinstance(args[0],list):
387 lst = tuple(args[0])
388 elif isinstance(args[0],tuple):
389 lst = args[0]
390 else:
391 lst = tuple(args)
392 varlist["args"] = "%d scantables" % len(lst)
393 # need special formatting her for history...
394 from asap._asap import stmath
395 stm = stmath()
396 for s in lst:
397 if not isinstance(s,scantable):
398 msg = "Please give a list of scantables"
399 if rcParams['verbose']:
400 print msg
401 return
402 else:
403 raise TypeError(msg)
404 s = scantable(stm._merge(lst))
405 s._add_history("merge", varlist)
406 print_log()
407 return s
[1074]408
Note: See TracBrowser for help on using the repository browser.