[1689] | 1 | import os |
---|
| 2 | import math |
---|
| 3 | from asap import scantable |
---|
| 4 | from asap import merge |
---|
| 5 | from asap import fitter |
---|
| 6 | from asap import selector |
---|
| 7 | from asap import rcParams |
---|
[1725] | 8 | from asap._asap import atmosphere |
---|
[1689] | 9 | |
---|
[1725] | 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 |
---|
[1726] | 21 | equivalent to the model available in MIRIAD. The actual math is a |
---|
[1725] | 22 | convertion of the Fortran code written by Bob Sault for MIRIAD. |
---|
[1726] | 23 | It implements a simple model of the atmosphere and Liebe's model (1985) |
---|
[1725] | 24 | of the complex refractive index of air. |
---|
| 25 | |
---|
| 26 | The model of the atmosphere is one with an exponential fall-off in |
---|
[1726] | 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 |
---|
[1725] | 29 | equation and hydrostatic equilibrium. |
---|
| 30 | |
---|
[1726] | 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 |
---|
[1725] | 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) |
---|
[1726] | 38 | pressure: air pressure at the sea level if the observatory |
---|
| 39 | elevation is set to non-zero value (note, by |
---|
[1725] | 40 | default is set to 700m) or at the observatory |
---|
[1726] | 41 | ground level if the elevation is set to 0. (The |
---|
[1725] | 42 | value is in Pascals or hPa, default 101325 Pa |
---|
[1726] | 43 | humidity: air humidity at the observatory (fractional), |
---|
[1725] | 44 | default is 0.5 |
---|
| 45 | elevation: observatory elevation about sea level (in meters) |
---|
| 46 | """ |
---|
[1754] | 47 | self._atm = atmosphere(temperature, self._to_pascals(pressure), |
---|
| 48 | humidity) |
---|
[1726] | 49 | self.set_observatory_elevation(elevation) |
---|
[1725] | 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) |
---|
[1726] | 82 | pressure: air pressure at the sea level if the observatory |
---|
| 83 | elevation is set to non-zero value (note, by |
---|
[1725] | 84 | default is set to 700m) or at the observatory |
---|
[1726] | 85 | ground level if the elevation is set to 0. (The |
---|
[1725] | 86 | value is in Pascals or hPa, default 101325 Pa |
---|
[1726] | 87 | humidity: air humidity at the observatory (fractional), |
---|
[1725] | 88 | default is 0.5 |
---|
| 89 | """ |
---|
| 90 | pressure = self._to_pascals(pressure) |
---|
| 91 | self._atm.set_weather(temperature, pressure, humidity) |
---|
| 92 | |
---|
[1726] | 93 | def set_observatory_elevation(self, elevation): |
---|
[1725] | 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 | """ |
---|
[1754] | 99 | self._atm.set_observatory_elevation(elevation) |
---|
[1725] | 100 | |
---|
| 101 | |
---|
[1689] | 102 | def _import_data(data): |
---|
[1722] | 103 | if not isinstance(data, (list,tuple)): |
---|
[1689] | 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 | |
---|
[1725] | 121 | def skydip(data, averagepol=True, tsky=300., plot=False, |
---|
| 122 | temperature=288, pressure=101325., humidity=0.5): |
---|
[1689] | 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 | """ |
---|
[1722] | 153 | rcsave = rcParams['verbose'] |
---|
| 154 | rcParams['verbose'] = False |
---|
[1725] | 155 | if plot: |
---|
| 156 | from matplotlib import pylab |
---|
[1689] | 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 = [] |
---|
[1754] | 165 | om = model(temperature, pressure, humidity) |
---|
[1689] | 166 | for ino in inos: |
---|
| 167 | sel.set_ifs(ino) |
---|
| 168 | opacity = [] |
---|
[1722] | 169 | fits = [] |
---|
| 170 | airms = [] |
---|
| 171 | tsyss = [] |
---|
| 172 | if plot: |
---|
[1725] | 173 | pylab.cla() |
---|
| 174 | pylab.ioff() |
---|
| 175 | pylab.clf() |
---|
| 176 | pylab.xlabel("Airmass") |
---|
| 177 | pylab.ylabel(r"$T_{sys}$") |
---|
[1689] | 178 | for pno in pnos: |
---|
| 179 | sel.set_polarisations(pno) |
---|
| 180 | scan.set_selection(basesel+sel) |
---|
[1722] | 181 | freq = scan.get_coordinate(0).get_reference_value()/1e9 |
---|
| 182 | freqstr = "%0.4f GHz" % freq |
---|
[1689] | 183 | tsys = scan.get_tsys() |
---|
| 184 | elev = scan.get_elevation() |
---|
| 185 | airmass = [ 1./math.sin(i) for i in elev ] |
---|
[1722] | 186 | airms.append(airmass) |
---|
| 187 | tsyss.append(tsys) |
---|
[1689] | 188 | f.set_data(airmass, tsys) |
---|
| 189 | f.fit() |
---|
[1722] | 190 | fits.append(f.get_fit()) |
---|
[1689] | 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 |
---|
[1722] | 197 | if plot: |
---|
| 198 | colors = ['b','g','k'] |
---|
[1725] | 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), |
---|
[1722] | 204 | r"$\tau_{fit}=%0.2f$" % opacity[i], |
---|
| 205 | color=colors[i]) |
---|
| 206 | if averagepol: |
---|
[1725] | 207 | pylab.figtext(0.7,0.3-(n/30.0), |
---|
| 208 | r"$\tau_{avg}=%0.2f$" % opacities[-1], |
---|
[1722] | 209 | color='r') |
---|
[1725] | 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') |
---|
[1726] | 214 | |
---|
[1725] | 215 | pylab.title("IF%d : %s" % (ino, freqstr)) |
---|
[1722] | 216 | |
---|
[1725] | 217 | pylab.ion() |
---|
| 218 | pylab.draw() |
---|
[1722] | 219 | raw_input("Hit <return> for next fit...") |
---|
[1689] | 220 | sel.reset() |
---|
[1722] | 221 | |
---|
[1689] | 222 | scan.set_selection(basesel) |
---|
[1722] | 223 | rcParams['verbose'] = rcsave |
---|
| 224 | if plot: |
---|
[1725] | 225 | pylab.close() |
---|
[1689] | 226 | return opacities |
---|