source: trunk/python/asapmath.py@ 1520

Last change on this file since 1520 was 1391, checked in by Malte Marquarding, 17 years ago

merge from alma branch to get alma/GBT support. Commented out fluxUnit changes as they are using a chnaged interface to PKSreader/writer. Also commented out progress meter related code.

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