source: trunk/python/opacity.py@ 1710

Last change on this file since 1710 was 1689, checked in by Malte Marquarding, 15 years ago

Ticket #177: added skydip function to determine opacities. This resulted in a slight interface change to accept lists of opacities in the scantable.opacity function

File size: 3.0 KB
Line 
1import os
2import math
3from asap import scantable
4from asap import merge
5from asap import fitter
6from asap import selector
7from asap import rcParams
8
9def _import_data(data):
10 if not hasattr(data, "__len__"):
11 if isinstance(data, scantable):
12 return data
13 elif isinstance(data, str):
14 return scantable(data)
15 tables = []
16 for d in data:
17 if isinstance(d, scantable):
18 tables.append(d)
19 elif isinstance(d, str):
20 if os.path.exists(d):
21 tables.append(scantable(d))
22 else:
23 raise IOError("Data file doesn't exists")
24 else:
25 raise TypeError("data is not a scantable or valid file")
26 return merge(tables)
27
28def skydip(data, averagepol=True, tsky=300., plot=False):
29 """Determine the opacity from a set of 'skydip' obervations.
30 This can be any set of observations over a range of elevations,
31 but will ususally be a dedicated (set of) scan(s).
32 Return a list of 'n' opacities for 'n' IFs. In case of averagepol
33 being 'False' a list of 'n*m' elements where 'm' is the number of
34 polarisations, e.g.
35 nIF = 3, nPol = 2 => [if0pol0, if0pol1, if1pol0, if1pol1, if2pol0, if2pol1]
36
37 The opacity is determined by fitting a first order polynomial to:
38
39
40 Tsys(airmass) = p0 + airmass*p1
41
42 where
43
44 airmass = 1/sin(elevation)
45
46 tau = p1/Tsky
47
48 Parameters:
49 data: a list of file names or scantables or a single
50 file name or scantable.
51 averagepol: Return the average of the opacities for the polarisations
52 This might be useful to set to 'False' if one polarisation
53 is corrupted (Mopra). If set to 'False', an opacity value
54 per polarisation is returned.
55 tksy: The sky temperature (default 300.0K). This might
56 be read from the data in the future.
57 plot: Plot each fit (airmass vs. Tsys). Default is 'False'
58 """
59 scan = _import_data(data)
60 f = fitter()
61 f.set_function(poly=1)
62 sel = selector()
63 basesel = scan.get_selection()
64 inos = scan.getifnos()
65 pnos = scan.getpolnos()
66 opacities = []
67 for ino in inos:
68 sel.set_ifs(ino)
69 opacity = []
70 for pno in pnos:
71 sel.set_polarisations(pno)
72 scan.set_selection(basesel+sel)
73 tsys = scan.get_tsys()
74 elev = scan.get_elevation()
75 airmass = [ 1./math.sin(i) for i in elev ]
76 f.set_data(airmass, tsys)
77 f.fit()
78 if plot:
79 f.plot(residual=True, plotparms=True)
80 raw_input("Hit <return> for next fit...")
81 params = f.get_parameters()["params"]
82 opacity.append(params[1]/tsky)
83 if averagepol:
84 opacities.append(sum(opacity)/len(opacity))
85 else:
86 opacities += opacity
87 sel.reset()
88 scan.set_selection(basesel)
89 return opacities
Note: See TracBrowser for help on using the repository browser.