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

Last change on this file since 1612 was 1612, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-729, CAS-1147

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes

Module(s): Module Names change impacts.

Description: Describe your changes here...

I have changed that almost all log messages are output to casapy.log,
not to the terminal window. After this change, asap becomes to depend on casapy
and is not running in standalone, because asap have to import taskinit module
to access casalogger.


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