1 | import os |
---|
2 | import math |
---|
3 | from asap import scantable |
---|
4 | from asap import merge |
---|
5 | from asap import fitter |
---|
6 | from asap import selector |
---|
7 | from asap import rcParams |
---|
8 | from asap import xyplotter |
---|
9 | |
---|
10 | def _import_data(data): |
---|
11 | if not isinstance(data, (list,tuple)): |
---|
12 | if isinstance(data, scantable): |
---|
13 | return data |
---|
14 | elif isinstance(data, str): |
---|
15 | return scantable(data) |
---|
16 | tables = [] |
---|
17 | for d in data: |
---|
18 | if isinstance(d, scantable): |
---|
19 | tables.append(d) |
---|
20 | elif isinstance(d, str): |
---|
21 | if os.path.exists(d): |
---|
22 | tables.append(scantable(d)) |
---|
23 | else: |
---|
24 | raise IOError("Data file doesn't exists") |
---|
25 | else: |
---|
26 | raise TypeError("data is not a scantable or valid file") |
---|
27 | return merge(tables) |
---|
28 | |
---|
29 | def skydip(data, averagepol=True, tsky=300., plot=False): |
---|
30 | """Determine the opacity from a set of 'skydip' obervations. |
---|
31 | This can be any set of observations over a range of elevations, |
---|
32 | but will ususally be a dedicated (set of) scan(s). |
---|
33 | Return a list of 'n' opacities for 'n' IFs. In case of averagepol |
---|
34 | being 'False' a list of 'n*m' elements where 'm' is the number of |
---|
35 | polarisations, e.g. |
---|
36 | nIF = 3, nPol = 2 => [if0pol0, if0pol1, if1pol0, if1pol1, if2pol0, if2pol1] |
---|
37 | |
---|
38 | The opacity is determined by fitting a first order polynomial to: |
---|
39 | |
---|
40 | |
---|
41 | Tsys(airmass) = p0 + airmass*p1 |
---|
42 | |
---|
43 | where |
---|
44 | |
---|
45 | airmass = 1/sin(elevation) |
---|
46 | |
---|
47 | tau = p1/Tsky |
---|
48 | |
---|
49 | Parameters: |
---|
50 | data: a list of file names or scantables or a single |
---|
51 | file name or scantable. |
---|
52 | averagepol: Return the average of the opacities for the polarisations |
---|
53 | This might be useful to set to 'False' if one polarisation |
---|
54 | is corrupted (Mopra). If set to 'False', an opacity value |
---|
55 | per polarisation is returned. |
---|
56 | tksy: The sky temperature (default 300.0K). This might |
---|
57 | be read from the data in the future. |
---|
58 | plot: Plot each fit (airmass vs. Tsys). Default is 'False' |
---|
59 | """ |
---|
60 | rcsave = rcParams['verbose'] |
---|
61 | rcParams['verbose'] = False |
---|
62 | scan = _import_data(data) |
---|
63 | f = fitter() |
---|
64 | f.set_function(poly=1) |
---|
65 | sel = selector() |
---|
66 | basesel = scan.get_selection() |
---|
67 | inos = scan.getifnos() |
---|
68 | pnos = scan.getpolnos() |
---|
69 | opacities = [] |
---|
70 | for ino in inos: |
---|
71 | sel.set_ifs(ino) |
---|
72 | opacity = [] |
---|
73 | fits = [] |
---|
74 | airms = [] |
---|
75 | tsyss = [] |
---|
76 | |
---|
77 | if plot: |
---|
78 | xyplotter.cla() |
---|
79 | xyplotter.ioff() |
---|
80 | xyplotter.clf() |
---|
81 | xyplotter.xlabel("Airmass") |
---|
82 | xyplotter.ylabel(r"$T_{sys}$") |
---|
83 | for pno in pnos: |
---|
84 | sel.set_polarisations(pno) |
---|
85 | scan.set_selection(basesel+sel) |
---|
86 | freq = scan.get_coordinate(0).get_reference_value()/1e9 |
---|
87 | freqstr = "%0.4f GHz" % freq |
---|
88 | tsys = scan.get_tsys() |
---|
89 | elev = scan.get_elevation() |
---|
90 | airmass = [ 1./math.sin(i) for i in elev ] |
---|
91 | airms.append(airmass) |
---|
92 | tsyss.append(tsys) |
---|
93 | f.set_data(airmass, tsys) |
---|
94 | f.fit() |
---|
95 | fits.append(f.get_fit()) |
---|
96 | params = f.get_parameters()["params"] |
---|
97 | opacity.append(params[1]/tsky) |
---|
98 | if averagepol: |
---|
99 | opacities.append(sum(opacity)/len(opacity)) |
---|
100 | else: |
---|
101 | opacities += opacity |
---|
102 | if plot: |
---|
103 | colors = ['b','g','k'] |
---|
104 | for i in range(len(airms)): |
---|
105 | xyplotter.plot(airms[i], tsyss[i], 'o', color=colors[i]) |
---|
106 | xyplotter.plot(airms[i], fits[i], '-', color=colors[i]) |
---|
107 | xyplotter.figtext(0.7,0.3-(i/30.0), |
---|
108 | r"$\tau_{fit}=%0.2f$" % opacity[i], |
---|
109 | color=colors[i]) |
---|
110 | if averagepol: |
---|
111 | xyplotter.figtext(0.7,0.3-(len(airms)/30.0), |
---|
112 | r"$\tau=%0.2f$" % opacities[-1], |
---|
113 | color='r') |
---|
114 | xyplotter.title("IF%d : %s" % (ino, freqstr)) |
---|
115 | |
---|
116 | xyplotter.ion() |
---|
117 | xyplotter.draw() |
---|
118 | raw_input("Hit <return> for next fit...") |
---|
119 | sel.reset() |
---|
120 | |
---|
121 | scan.set_selection(basesel) |
---|
122 | rcParams['verbose'] = rcsave |
---|
123 | if plot: |
---|
124 | xyplotter.close() |
---|
125 | return opacities |
---|