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._asap import atmosphere |
---|
9 | |
---|
10 | |
---|
11 | class model(object): |
---|
12 | def _to_pascals(self, val): |
---|
13 | if val > 2000: |
---|
14 | return val |
---|
15 | return val*100 |
---|
16 | |
---|
17 | def __init__(self, temperature=288, pressure=101325., humidity=0.5, |
---|
18 | elevation=700.): |
---|
19 | """ |
---|
20 | This class implements opacity/atmospheric brightness temperature model |
---|
21 | equivalent to the model available in MIRIAD. The actual math is a |
---|
22 | convertion of the Fortran code written by Bob Sault for MIRIAD. |
---|
23 | It implements a simple model of the atmosphere and Liebe's model (1985) |
---|
24 | of the complex refractive index of air. |
---|
25 | |
---|
26 | The model of the atmosphere is one with an exponential fall-off in |
---|
27 | the water vapour content (scale height of 1540 m) and a temperature |
---|
28 | lapse rate of 6.5 mK/m. Otherwise the atmosphere obeys the ideal gas |
---|
29 | equation and hydrostatic equilibrium. |
---|
30 | |
---|
31 | Note, the model includes atmospheric lines up to 800 GHz, but was not |
---|
32 | rigorously tested above 100 GHz and for instruments located at |
---|
33 | a significant elevation. For high-elevation sites it may be necessary to |
---|
34 | adjust scale height and lapse rate. |
---|
35 | |
---|
36 | Parameters: |
---|
37 | temperature: air temperature at the observatory (K) |
---|
38 | pressure: air pressure at the sea level if the observatory |
---|
39 | elevation is set to non-zero value (note, by |
---|
40 | default is set to 700m) or at the observatory |
---|
41 | ground level if the elevation is set to 0. (The |
---|
42 | value is in Pascals or hPa, default 101325 Pa |
---|
43 | humidity: air humidity at the observatory (fractional), |
---|
44 | default is 0.5 |
---|
45 | elevation: observatory elevation about sea level (in meters) |
---|
46 | """ |
---|
47 | self._atm = atmosphere(temp, self._to_pascals(pressure), humidity) |
---|
48 | self.set_observatory_elevation(elevation) |
---|
49 | |
---|
50 | def get_opacities(self, freq, elevation=None): |
---|
51 | """Get the opacity value(s) for the fiven frequency(ies). |
---|
52 | If no elevation is given the opacities for the zenith are returned. |
---|
53 | If an elevation is specified refraction is also taken into account. |
---|
54 | Parameters: |
---|
55 | freq: a frequency value in Hz, or a list of frequency values. |
---|
56 | One opacity value per frequency is returned as a scalar |
---|
57 | or list. |
---|
58 | elevation: the elevation at which to compute the opacity. If `None` |
---|
59 | is given (default) the zenith opacity is returned. |
---|
60 | |
---|
61 | |
---|
62 | """ |
---|
63 | func = None |
---|
64 | if isinstance(freq, (list, tuple)): |
---|
65 | if elevation is None: |
---|
66 | return self._atm.zenith_opacities(freq) |
---|
67 | else: |
---|
68 | elevation *= math.pi/180. |
---|
69 | return self._atm.opacities(freq, elevation) |
---|
70 | else: |
---|
71 | if elevation is None: |
---|
72 | return self._atm.zenith_opacity(freq) |
---|
73 | else: |
---|
74 | elevation *= math.pi/180. |
---|
75 | return self._atm.opacity(freq, elevation) |
---|
76 | |
---|
77 | def set_weather(self, temperature, pressure, humidity): |
---|
78 | """Update the model using the given environmental parameters. |
---|
79 | Parameters: |
---|
80 | temperature: air temperature at the observatory (K) |
---|
81 | pressure: air pressure at the sea level if the observatory |
---|
82 | elevation is set to non-zero value (note, by |
---|
83 | default is set to 700m) or at the observatory |
---|
84 | ground level if the elevation is set to 0. (The |
---|
85 | value is in Pascals or hPa, default 101325 Pa |
---|
86 | humidity: air humidity at the observatory (fractional), |
---|
87 | default is 0.5 |
---|
88 | """ |
---|
89 | pressure = self._to_pascals(pressure) |
---|
90 | self._atm.set_weather(temperature, pressure, humidity) |
---|
91 | |
---|
92 | def set_observatory_elevation(self, elevation): |
---|
93 | """Update the model using the given the observatory elevation |
---|
94 | Parameters: |
---|
95 | elevation: the elevation at which to compute the opacity. If `None` |
---|
96 | is given (default) the zenith opacity is returned. |
---|
97 | """ |
---|
98 | self._atm.set_observatory_elevation(el) |
---|
99 | |
---|
100 | |
---|
101 | def _import_data(data): |
---|
102 | if not isinstance(data, (list,tuple)): |
---|
103 | if isinstance(data, scantable): |
---|
104 | return data |
---|
105 | elif isinstance(data, str): |
---|
106 | return scantable(data) |
---|
107 | tables = [] |
---|
108 | for d in data: |
---|
109 | if isinstance(d, scantable): |
---|
110 | tables.append(d) |
---|
111 | elif isinstance(d, str): |
---|
112 | if os.path.exists(d): |
---|
113 | tables.append(scantable(d)) |
---|
114 | else: |
---|
115 | raise IOError("Data file doesn't exists") |
---|
116 | else: |
---|
117 | raise TypeError("data is not a scantable or valid file") |
---|
118 | return merge(tables) |
---|
119 | |
---|
120 | def skydip(data, averagepol=True, tsky=300., plot=False, |
---|
121 | temperature=288, pressure=101325., humidity=0.5): |
---|
122 | """Determine the opacity from a set of 'skydip' obervations. |
---|
123 | This can be any set of observations over a range of elevations, |
---|
124 | but will ususally be a dedicated (set of) scan(s). |
---|
125 | Return a list of 'n' opacities for 'n' IFs. In case of averagepol |
---|
126 | being 'False' a list of 'n*m' elements where 'm' is the number of |
---|
127 | polarisations, e.g. |
---|
128 | nIF = 3, nPol = 2 => [if0pol0, if0pol1, if1pol0, if1pol1, if2pol0, if2pol1] |
---|
129 | |
---|
130 | The opacity is determined by fitting a first order polynomial to: |
---|
131 | |
---|
132 | |
---|
133 | Tsys(airmass) = p0 + airmass*p1 |
---|
134 | |
---|
135 | where |
---|
136 | |
---|
137 | airmass = 1/sin(elevation) |
---|
138 | |
---|
139 | tau = p1/Tsky |
---|
140 | |
---|
141 | Parameters: |
---|
142 | data: a list of file names or scantables or a single |
---|
143 | file name or scantable. |
---|
144 | averagepol: Return the average of the opacities for the polarisations |
---|
145 | This might be useful to set to 'False' if one polarisation |
---|
146 | is corrupted (Mopra). If set to 'False', an opacity value |
---|
147 | per polarisation is returned. |
---|
148 | tksy: The sky temperature (default 300.0K). This might |
---|
149 | be read from the data in the future. |
---|
150 | plot: Plot each fit (airmass vs. Tsys). Default is 'False' |
---|
151 | """ |
---|
152 | rcsave = rcParams['verbose'] |
---|
153 | rcParams['verbose'] = False |
---|
154 | if plot: |
---|
155 | from matplotlib import pylab |
---|
156 | scan = _import_data(data) |
---|
157 | f = fitter() |
---|
158 | f.set_function(poly=1) |
---|
159 | sel = selector() |
---|
160 | basesel = scan.get_selection() |
---|
161 | inos = scan.getifnos() |
---|
162 | pnos = scan.getpolnos() |
---|
163 | opacities = [] |
---|
164 | om = opacitymodel(temperature, pressure, humidity) |
---|
165 | for ino in inos: |
---|
166 | sel.set_ifs(ino) |
---|
167 | opacity = [] |
---|
168 | fits = [] |
---|
169 | airms = [] |
---|
170 | tsyss = [] |
---|
171 | if plot: |
---|
172 | pylab.cla() |
---|
173 | pylab.ioff() |
---|
174 | pylab.clf() |
---|
175 | pylab.xlabel("Airmass") |
---|
176 | pylab.ylabel(r"$T_{sys}$") |
---|
177 | for pno in pnos: |
---|
178 | sel.set_polarisations(pno) |
---|
179 | scan.set_selection(basesel+sel) |
---|
180 | freq = scan.get_coordinate(0).get_reference_value()/1e9 |
---|
181 | freqstr = "%0.4f GHz" % freq |
---|
182 | tsys = scan.get_tsys() |
---|
183 | elev = scan.get_elevation() |
---|
184 | airmass = [ 1./math.sin(i) for i in elev ] |
---|
185 | airms.append(airmass) |
---|
186 | tsyss.append(tsys) |
---|
187 | f.set_data(airmass, tsys) |
---|
188 | f.fit() |
---|
189 | fits.append(f.get_fit()) |
---|
190 | params = f.get_parameters()["params"] |
---|
191 | opacity.append(params[1]/tsky) |
---|
192 | if averagepol: |
---|
193 | opacities.append(sum(opacity)/len(opacity)) |
---|
194 | else: |
---|
195 | opacities += opacity |
---|
196 | if plot: |
---|
197 | colors = ['b','g','k'] |
---|
198 | n = len(airms) |
---|
199 | for i in range(n): |
---|
200 | pylab.plot(airms[i], tsyss[i], 'o', color=colors[i]) |
---|
201 | pylab.plot(airms[i], fits[i], '-', color=colors[i]) |
---|
202 | pylab.figtext(0.7,0.3-(i/30.0), |
---|
203 | r"$\tau_{fit}=%0.2f$" % opacity[i], |
---|
204 | color=colors[i]) |
---|
205 | if averagepol: |
---|
206 | pylab.figtext(0.7,0.3-(n/30.0), |
---|
207 | r"$\tau_{avg}=%0.2f$" % opacities[-1], |
---|
208 | color='r') |
---|
209 | n +=1 |
---|
210 | pylab.figtext(0.7,0.3-(n/30.0), |
---|
211 | r"$\tau_{model}=%0.2f$" % om.get_opacities(freq*1e9), |
---|
212 | color='grey') |
---|
213 | |
---|
214 | pylab.title("IF%d : %s" % (ino, freqstr)) |
---|
215 | |
---|
216 | pylab.ion() |
---|
217 | pylab.draw() |
---|
218 | raw_input("Hit <return> for next fit...") |
---|
219 | sel.reset() |
---|
220 | |
---|
221 | scan.set_selection(basesel) |
---|
222 | rcParams['verbose'] = rcsave |
---|
223 | if plot: |
---|
224 | pylab.close() |
---|
225 | return opacities |
---|