source: branches/newfiller/python/opacity.py@ 2406

Last change on this file since 2406 was 1770, checked in by Kana Sugimoto, 14 years ago

New Development: Yes

JIRA Issue: No (to merge the alma branch to trunk)

Ready for Test: No

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): single dish package

Description: merged changes in trunk after ASAP3.0.0 release (r1752-1768)


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