| 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
|
|---|
| 8 | from asap._asap import atmosphere
|
|---|
| 9 |
|
|---|
| 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
|
|---|
| 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 |
|
|---|
| 102 | def _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 |
|
|---|
| 121 | def 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
|
|---|
