source: tags/asap-3.1.0/python/opacity.py@ 2803

Last change on this file since 2803 was 1927, checked in by Takeshi Nakazato, 14 years ago

New Development: No

JIRA Issue: No

Ready for Test: No

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Fixed typo in help text.


File size: 8.7 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
9
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 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
Note: See TracBrowser for help on using the repository browser.