1 | __all__ = ["model", "skydip"] |
---|
2 | import os |
---|
3 | import math |
---|
4 | from asap.scantable import scantable |
---|
5 | from asap.asapmath import merge |
---|
6 | from asap.asapfitter import fitter |
---|
7 | from asap.selector import selector |
---|
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(temperature, self._to_pascals(pressure), |
---|
48 | humidity) |
---|
49 | self.set_observatory_elevation(elevation) |
---|
50 | |
---|
51 | def get_opacities(self, freq, elevation=None): |
---|
52 | """Get the opacity value(s) for the fiven frequency(ies). |
---|
53 | If no elevation is given the opacities for the zenith are returned. |
---|
54 | If an elevation is specified refraction is also taken into account. |
---|
55 | Parameters: |
---|
56 | freq: a frequency value in Hz, or a list of frequency values. |
---|
57 | One opacity value per frequency is returned as a scalar |
---|
58 | or list. |
---|
59 | elevation: the elevation at which to compute the opacity. If `None` |
---|
60 | is given (default) the zenith opacity is returned. |
---|
61 | |
---|
62 | |
---|
63 | """ |
---|
64 | func = None |
---|
65 | if isinstance(freq, (list, tuple)): |
---|
66 | if elevation is None: |
---|
67 | return self._atm.zenith_opacities(freq) |
---|
68 | else: |
---|
69 | elevation *= math.pi/180. |
---|
70 | return self._atm.opacities(freq, elevation) |
---|
71 | else: |
---|
72 | if elevation is None: |
---|
73 | return self._atm.zenith_opacity(freq) |
---|
74 | else: |
---|
75 | elevation *= math.pi/180. |
---|
76 | return self._atm.opacity(freq, elevation) |
---|
77 | |
---|
78 | def set_weather(self, temperature, pressure, humidity): |
---|
79 | """Update the model using the given environmental parameters. |
---|
80 | Parameters: |
---|
81 | temperature: air temperature at the observatory (K) |
---|
82 | pressure: air pressure at the sea level if the observatory |
---|
83 | elevation is set to non-zero value (note, by |
---|
84 | default is set to 700m) or at the observatory |
---|
85 | ground level if the elevation is set to 0. (The |
---|
86 | value is in Pascals or hPa, default 101325 Pa |
---|
87 | humidity: air humidity at the observatory (fractional), |
---|
88 | default is 0.5 |
---|
89 | """ |
---|
90 | pressure = self._to_pascals(pressure) |
---|
91 | self._atm.set_weather(temperature, pressure, humidity) |
---|
92 | |
---|
93 | def set_observatory_elevation(self, elevation): |
---|
94 | """Update the model using the given the observatory elevation |
---|
95 | Parameters: |
---|
96 | elevation: the elevation at which to compute the opacity. If `None` |
---|
97 | is given (default) the zenith opacity is returned. |
---|
98 | """ |
---|
99 | self._atm.set_observatory_elevation(elevation) |
---|
100 | |
---|
101 | |
---|
102 | def _import_data(data): |
---|
103 | if not isinstance(data, (list,tuple)): |
---|
104 | if isinstance(data, scantable): |
---|
105 | return data |
---|
106 | elif isinstance(data, str): |
---|
107 | return scantable(data) |
---|
108 | tables = [] |
---|
109 | for d in data: |
---|
110 | if isinstance(d, scantable): |
---|
111 | tables.append(d) |
---|
112 | elif isinstance(d, str): |
---|
113 | if os.path.exists(d): |
---|
114 | tables.append(scantable(d)) |
---|
115 | else: |
---|
116 | raise IOError("Data file doesn't exists") |
---|
117 | else: |
---|
118 | raise TypeError("data is not a scantable or valid file") |
---|
119 | return merge(tables) |
---|
120 | |
---|
121 | def skydip(data, averagepol=True, tsky=300., plot=False, |
---|
122 | temperature=288, pressure=101325., humidity=0.5): |
---|
123 | """Determine the opacity from a set of 'skydip' obervations. |
---|
124 | This can be any set of observations over a range of elevations, |
---|
125 | but will ususally be a dedicated (set of) scan(s). |
---|
126 | Return a list of 'n' opacities for 'n' IFs. In case of averagepol |
---|
127 | being 'False' a list of 'n*m' elements where 'm' is the number of |
---|
128 | polarisations, e.g. |
---|
129 | nIF = 3, nPol = 2 => [if0pol0, if0pol1, if1pol0, if1pol1, if2pol0, if2pol1] |
---|
130 | |
---|
131 | The opacity is determined by fitting a first order polynomial to: |
---|
132 | |
---|
133 | |
---|
134 | Tsys(airmass) = p0 + airmass*p1 |
---|
135 | |
---|
136 | where |
---|
137 | |
---|
138 | airmass = 1/sin(elevation) |
---|
139 | |
---|
140 | tau = p1/Tsky |
---|
141 | |
---|
142 | Parameters: |
---|
143 | data: a list of file names or scantables or a single |
---|
144 | file name or scantable. |
---|
145 | averagepol: Return the average of the opacities for the polarisations |
---|
146 | This might be useful to set to 'False' if one polarisation |
---|
147 | is corrupted (Mopra). If set to 'False', an opacity value |
---|
148 | per polarisation is returned. |
---|
149 | tsky: The sky temperature (default 300.0K). This might |
---|
150 | be read from the data in the future. |
---|
151 | plot: Plot each fit (airmass vs. Tsys). Default is 'False' |
---|
152 | """ |
---|
153 | if plot: |
---|
154 | from matplotlib import pylab |
---|
155 | scan = _import_data(data) |
---|
156 | f = fitter() |
---|
157 | f.set_function(poly=1) |
---|
158 | sel = selector() |
---|
159 | basesel = scan.get_selection() |
---|
160 | inos = scan.getifnos() |
---|
161 | pnos = scan.getpolnos() |
---|
162 | opacities = [] |
---|
163 | om = model(temperature, pressure, humidity) |
---|
164 | for ino in inos: |
---|
165 | sel.set_ifs(ino) |
---|
166 | opacity = [] |
---|
167 | fits = [] |
---|
168 | airms = [] |
---|
169 | tsyss = [] |
---|
170 | if plot: |
---|
171 | pylab.cla() |
---|
172 | pylab.ioff() |
---|
173 | pylab.clf() |
---|
174 | pylab.xlabel("Airmass") |
---|
175 | pylab.ylabel(r"$T_{sys}$") |
---|
176 | for pno in pnos: |
---|
177 | sel.set_polarisations(pno) |
---|
178 | scan.set_selection(basesel+sel) |
---|
179 | freq = scan.get_coordinate(0).get_reference_value()/1e9 |
---|
180 | freqstr = "%0.4f GHz" % freq |
---|
181 | tsys = scan.get_tsys() |
---|
182 | elev = scan.get_elevation() |
---|
183 | airmass = [ 1./math.sin(i) for i in elev ] |
---|
184 | airms.append(airmass) |
---|
185 | tsyss.append(tsys) |
---|
186 | f.set_data(airmass, tsys) |
---|
187 | f.fit() |
---|
188 | fits.append(f.get_fit()) |
---|
189 | params = f.get_parameters()["params"] |
---|
190 | opacity.append(params[1]/tsky) |
---|
191 | if averagepol: |
---|
192 | opacities.append(sum(opacity)/len(opacity)) |
---|
193 | else: |
---|
194 | opacities += opacity |
---|
195 | if plot: |
---|
196 | colors = ['b','g','k'] |
---|
197 | n = len(airms) |
---|
198 | for i in range(n): |
---|
199 | pylab.plot(airms[i], tsyss[i], 'o', color=colors[i]) |
---|
200 | pylab.plot(airms[i], fits[i], '-', color=colors[i]) |
---|
201 | pylab.figtext(0.7,0.3-(i/30.0), |
---|
202 | r"$\tau_{fit}=%0.2f$" % opacity[i], |
---|
203 | color=colors[i]) |
---|
204 | if averagepol: |
---|
205 | pylab.figtext(0.7,0.3-(n/30.0), |
---|
206 | r"$\tau_{avg}=%0.2f$" % opacities[-1], |
---|
207 | color='r') |
---|
208 | n +=1 |
---|
209 | pylab.figtext(0.7,0.3-(n/30.0), |
---|
210 | r"$\tau_{model}=%0.2f$" % om.get_opacities(freq*1e9), |
---|
211 | color='grey') |
---|
212 | |
---|
213 | pylab.title("IF%d : %s" % (ino, freqstr)) |
---|
214 | |
---|
215 | pylab.ion() |
---|
216 | pylab.draw() |
---|
217 | raw_input("Hit <return> for next fit...") |
---|
218 | sel.reset() |
---|
219 | |
---|
220 | scan.set_selection(basesel) |
---|
221 | if plot: |
---|
222 | pylab.close() |
---|
223 | return opacities |
---|