source: branches/hpc34/python/opacity.py@ 2794

Last change on this file since 2794 was 2447, checked in by Malte Marquarding, 13 years ago

Issue #264: silence skydip.

File size: 9.1 KB
Line 
1__all__ = ["model", "skydip"]
2import os
3import math
4from asap.scantable import scantable
5from asap.asapmath import merge
6from asap.asapfitter import fitter
7from asap.selector import selector
8from asap._asap import atmosphere
9from asap import rcParams
10
11class 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 given 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
102def _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
121def 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 # quiten output
154 verbose = rcParams["verbose"]
155 rcParams["verbose"] = False
156 try:
157 if plot:
158 from matplotlib import pylab
159 scan = _import_data(data)
160 f = fitter()
161 f.set_function(poly=1)
162 sel = selector()
163 basesel = scan.get_selection()
164 inos = scan.getifnos()
165 pnos = scan.getpolnos()
166 opacities = []
167 om = model(temperature, pressure, humidity)
168 for ino in inos:
169 sel.set_ifs(ino)
170 opacity = []
171 fits = []
172 airms = []
173 tsyss = []
174 if plot:
175 pylab.cla()
176 pylab.ioff()
177 pylab.clf()
178 pylab.xlabel("Airmass")
179 pylab.ylabel(r"$T_{sys}$")
180 for pno in pnos:
181 sel.set_polarisations(pno)
182 scan.set_selection(basesel+sel)
183 freq = scan.get_coordinate(0).get_reference_value()/1e9
184 freqstr = "%0.4f GHz" % freq
185 tsys = scan.get_tsys()
186 elev = scan.get_elevation()
187 airmass = [ 1./math.sin(i) for i in elev ]
188 airms.append(airmass)
189 tsyss.append(tsys)
190 f.set_data(airmass, tsys)
191 f.fit()
192 fits.append(f.get_fit())
193 params = f.get_parameters()["params"]
194 opacity.append(params[1]/tsky)
195 if averagepol:
196 opacities.append(sum(opacity)/len(opacity))
197 else:
198 opacities += opacity
199 if plot:
200 colors = ['b','g','k']
201 n = len(airms)
202 for i in range(n):
203 pylab.plot(airms[i], tsyss[i], 'o', color=colors[i])
204 pylab.plot(airms[i], fits[i], '-', color=colors[i])
205 pylab.figtext(0.7,0.3-(i/30.0),
206 r"$\tau_{fit}=%0.2f$" % opacity[i],
207 color=colors[i])
208 if averagepol:
209 pylab.figtext(0.7,0.3-(n/30.0),
210 r"$\tau_{avg}=%0.2f$" % opacities[-1],
211 color='r')
212 n +=1
213 pylab.figtext(0.7,0.3-(n/30.0),
214 r"$\tau_{model}=%0.2f$" % om.get_opacities(freq*1e9),
215 color='grey')
216
217 pylab.title("IF%d : %s" % (ino, freqstr))
218
219 pylab.ion()
220 pylab.draw()
221 raw_input("Hit <return> for next fit...")
222 sel.reset()
223
224 scan.set_selection(basesel)
225 if plot:
226 pylab.close()
227 return opacities
228 finally:
229 rcParams["verbose"] = verbose
Note: See TracBrowser for help on using the repository browser.