source: trunk/python/opacity.py@ 1831

Last change on this file since 1831 was 1826, checked in by Malte Marquarding, 14 years ago

Tidy up of imports (now imported from asap.). Also fixed some whitespace/tab issues

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