[1757] | 1 | //#---------------------------------------------------------------------------
|
---|
| 2 | //# STAtmosphere.h: Model of atmospheric opacity
|
---|
| 3 | //#---------------------------------------------------------------------------
|
---|
| 4 | //# Copyright (C) 2004
|
---|
| 5 | //# ATNF
|
---|
| 6 | //#
|
---|
| 7 | //# The code is based on the Fortran code written by Bob Sault for MIRIAD.
|
---|
| 8 | //# Converted to C++ by Max Voronkov. This code uses a simple model of the
|
---|
| 9 | //# atmosphere and Liebe's model (1985) of the complex refractive index of
|
---|
| 10 | //# air.
|
---|
| 11 | //#
|
---|
| 12 | //# The model of the atmosphere is one with an exponential fall-off in
|
---|
| 13 | //# the water vapour content (scale height of 1540 m) and a temperature lapse
|
---|
| 14 | //# rate of 6.5 mK/m. Otherwise the atmosphere obeys the ideal gas equation
|
---|
| 15 | //# and hydrostatic equilibrium.
|
---|
| 16 | //#
|
---|
| 17 | //# This program is free software; you can redistribute it and/or modify it
|
---|
| 18 | //# under the terms of the GNU General Public License as published by the Free
|
---|
| 19 | //# Software Foundation; either version 2 of the License, or (at your option)
|
---|
| 20 | //# any later version.
|
---|
| 21 | //#
|
---|
| 22 | //# This program is distributed in the hope that it will be useful, but
|
---|
| 23 | //# WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 24 | //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
|
---|
| 25 | //# Public License for more details.
|
---|
| 26 | //#
|
---|
| 27 | //# You should have received a copy of the GNU General Public License along
|
---|
| 28 | //# with this program; if not, write to the Free Software Foundation, Inc.,
|
---|
| 29 | //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
|
---|
| 30 | //#
|
---|
| 31 | //# Correspondence concerning this software should be addressed as follows:
|
---|
| 32 | //# Internet email: Malte.Marquarding@csiro.au
|
---|
| 33 | //# Postal address: Malte Marquarding,
|
---|
| 34 | //# Australia Telescope National Facility,
|
---|
| 35 | //# P.O. Box 76,
|
---|
| 36 | //# Epping, NSW, 2121,
|
---|
| 37 | //# AUSTRALIA
|
---|
| 38 | //#
|
---|
| 39 | //# $Id: STAtmosphere.h 1346 2007-04-26 03:24:41Z mar637 $
|
---|
| 40 | //#---------------------------------------------------------------------------
|
---|
| 41 |
|
---|
| 42 | #ifndef STATMOSPHERE_H
|
---|
| 43 | #define STATMOSPHERE_H
|
---|
| 44 |
|
---|
| 45 | // std includes
|
---|
| 46 | #include <vector>
|
---|
| 47 | #include <complex>
|
---|
| 48 |
|
---|
| 49 | namespace asap {
|
---|
| 50 |
|
---|
| 51 | /**
|
---|
| 52 | * This class implements opacity/atmospheric brightness temperature model
|
---|
| 53 | * equivalent to the model available in MIRIAD. The actual math is a
|
---|
| 54 | * convertion of the Fortran code written by Bob Sault for MIRIAD.
|
---|
| 55 | * It implements a simple model of the atmosphere and Liebe's model (1985)
|
---|
| 56 | * of the complex refractive index of air.
|
---|
| 57 | *
|
---|
| 58 | * The model of the atmosphere is one with an exponential fall-off in
|
---|
| 59 | * the water vapour content (scale height of 1540 m) and a temperature lapse
|
---|
| 60 | * rate of 6.5 mK/m. Otherwise the atmosphere obeys the ideal gas equation
|
---|
| 61 | * and hydrostatic equilibrium.
|
---|
| 62 | *
|
---|
| 63 | * Note, the model includes atmospheric lines up to 800 GHz, but was not
|
---|
| 64 | * rigorously tested above 100 GHz and for instruments located at
|
---|
| 65 | * a significant elevation. For high-elevation sites it may be necessary to
|
---|
| 66 | * adjust scale height and lapse rate.
|
---|
| 67 | *
|
---|
| 68 | * @brief The ASAP atmosphere opacity model
|
---|
| 69 | * @author Max Voronkov
|
---|
| 70 | * @date $Date: 2010-03-17 14:55:17 +1000 (Thu, 26 Apr 2007) $
|
---|
| 71 | * @version
|
---|
| 72 | */
|
---|
| 73 | class STAtmosphere {
|
---|
| 74 | public:
|
---|
| 75 | /**
|
---|
| 76 | * Default Constructor (apart from optional parameters).
|
---|
| 77 | * The class set up this way will assume International Standard Atmosphere (ISA) conditions,
|
---|
| 78 | * except for humidity. The latter is assumed to be 50%, which seems more realistic for
|
---|
| 79 | * Australian telescopes than 0%.
|
---|
| 80 | * @param[in] wvScale water vapour scale height (m), default is 1540m to match MIRIAD's model
|
---|
| 81 | * @param[in] maxAlt maximum altitude of the model atmosphere (m), plane parallel layers are spread linearly up to
|
---|
| 82 | * this height, default is 10000m to match MIRIAD.
|
---|
| 83 | * @param[in] nLayers number of plane parallel layers in the model (essentially for a numberical integration),
|
---|
| 84 | * default is 50 to match MIRIAD.
|
---|
| 85 | **/
|
---|
| 86 | explicit STAtmosphere(double wvScale = 1540., double maxAlt = 10000.0, size_t nLayers = 50);
|
---|
| 87 |
|
---|
| 88 | /**
|
---|
| 89 | * Constructor with explicitly given parameters of the atmosphere
|
---|
| 90 | * @param[in] temperature air temperature at the observatory (K)
|
---|
| 91 | * @param[in] pressure air pressure at the sea level if the observatory elevation
|
---|
| 92 | * is set to non-zero value (note, by default is set to 200m) or at the
|
---|
| 93 | * observatory ground level if the elevation is set to 0. (The value is in Pascals)
|
---|
| 94 | * @param[in] pressure air pressure at the observatory (Pascals)
|
---|
| 95 | * @param[in] humidity air humidity at the observatory (fraction)
|
---|
| 96 | * @param[in] lapseRate temperature lapse rate (K/m), default is 0.0065 K/m to match MIRIAD and ISA
|
---|
| 97 | * @param[in] wvScale water vapour scale height (m), default is 1540m to match MIRIAD's model
|
---|
| 98 | * @param[in] maxAlt maximum altitude of the model atmosphere (m), plane parallel layers are spread linearly up to
|
---|
| 99 | * this height, default is 10000m to match MIRIAD.
|
---|
| 100 | * @param[in] nLayers number of plane parallel layers in the model (essentially for a numberical integration),
|
---|
| 101 | * default is 50 to match MIRIAD.
|
---|
| 102 | **/
|
---|
| 103 | STAtmosphere(double temperature, double pressure, double humidity, double lapseRate = 0.0065,
|
---|
| 104 | double wvScale = 1540., double maxAlt = 10000.0, size_t nLayers = 50);
|
---|
| 105 |
|
---|
| 106 | /**
|
---|
| 107 | * Set the new weather station data, recompute the model
|
---|
| 108 | * @param[in] temperature air temperature at the observatory (K)
|
---|
| 109 | * @param[in] pressure air pressure at the sea level if the observatory elevation
|
---|
| 110 | * is set to non-zero value (note, by default is set to 200m) or at the
|
---|
| 111 | * observatory ground level if the elevation is set to 0. (The value is in Pascals)
|
---|
| 112 | * @param[in] humidity air humidity at the observatory (fraction)
|
---|
| 113 | **/
|
---|
| 114 | void setWeather(double temperature, double pressure, double humidity);
|
---|
| 115 |
|
---|
| 116 | /**
|
---|
| 117 | * Set the elevation of the observatory (height above mean sea level)
|
---|
| 118 | *
|
---|
| 119 | * The observatory elevation affects only interpretation of the pressure supplied as part
|
---|
| 120 | * of the weather data, if this value is non-zero, the pressure (e.g. in setWeather or
|
---|
| 121 | * constructor) is that at mean sea level. If the observatory elevation is set to zero,
|
---|
| 122 | * regardless on real elevation, the pressure is that at the observatory ground level.
|
---|
| 123 | *
|
---|
| 124 | * By default, 200m is assumed.
|
---|
| 125 | * @param[in] elev elevation in metres
|
---|
| 126 | **/
|
---|
| 127 | void setObservatoryElevation(double elev);
|
---|
| 128 |
|
---|
| 129 | /**
|
---|
| 130 | * Calculate zenith opacity at the given frequency. This is a simplified version
|
---|
| 131 | * of the routine implemented in MIRIAD, which calculates just zenith opacity and
|
---|
| 132 | * nothing else. Note, that if the opacity is high, 1/sin(el) law is not correct
|
---|
| 133 | * even in the plane parallel case due to refraction.
|
---|
| 134 | * @param[in] freq frequency of interest in Hz
|
---|
| 135 | * @return zenith opacity (nepers, i.e. dimensionless)
|
---|
| 136 | **/
|
---|
| 137 | double zenithOpacity(double freq) const;
|
---|
| 138 |
|
---|
| 139 | /**
|
---|
| 140 | * Calculate zenith opacity for the range of frequencies. Same as zenithOpacity, but
|
---|
| 141 | * for a vector of frequencies.
|
---|
| 142 | * @param[in] freqs vector of frequencies in Hz
|
---|
| 143 | * @return vector of zenith opacities, one value per frequency (nepers, i.e. dimensionless)
|
---|
| 144 | **/
|
---|
| 145 | std::vector<double> zenithOpacities(const std::vector<double> &freqs) const;
|
---|
| 146 |
|
---|
| 147 | /**
|
---|
| 148 | * Calculate opacity at the given frequency and elevation. This is a simplified
|
---|
| 149 | * version of the routine implemented in MIRIAD, which calculates just the opacity and
|
---|
| 150 | * nothing else. In contract to zenithOpacity, this method takes into account refraction
|
---|
| 151 | * and is more accurate than if one assumes 1/sin(el) factor.
|
---|
| 152 | * @param[in] freq frequency of interest in Hz
|
---|
| 153 | * @param[in] el elevation in radians
|
---|
| 154 | * @return zenith opacity (nepers, i.e. dimensionless)
|
---|
| 155 | **/
|
---|
| 156 | double opacity(double freq, double el) const;
|
---|
| 157 |
|
---|
| 158 | /**
|
---|
| 159 | * Calculate opacities for the range of frequencies at the given elevation. Same as
|
---|
| 160 | * opacity, but for a vector of frequencies.
|
---|
| 161 | * @param[in] freqs vector of frequencies in Hz
|
---|
| 162 | * @param[in] el elevation in radians
|
---|
| 163 | * @return vector of opacities, one value per frequency (nepers, i.e. dimensionless)
|
---|
| 164 | **/
|
---|
| 165 | std::vector<double> opacities(const std::vector<double> &freqs, double el) const;
|
---|
| 166 |
|
---|
| 167 | protected:
|
---|
| 168 | /**
|
---|
| 169 | * Build the atmosphere model based on exponential fall-off, ideal gas and hydrostatic
|
---|
| 170 | * equilibrium. The model parameters are taken from the data members of this class.
|
---|
| 171 | **/
|
---|
| 172 | void recomputeAtmosphereModel();
|
---|
| 173 |
|
---|
| 174 | /**
|
---|
| 175 | * Obtain the number of model layers, do consistency check that everything is
|
---|
| 176 | * resized accordingly
|
---|
| 177 | * @retrun number of model layers
|
---|
| 178 | **/
|
---|
| 179 | size_t nLayers() const;
|
---|
| 180 |
|
---|
| 181 | /**
|
---|
| 182 | * Determine the saturation pressure of water vapour for the given temperature.
|
---|
| 183 | *
|
---|
| 184 | * Reference:
|
---|
| 185 | * Waters, Refraction effects in the neutral atmosphere. Methods of
|
---|
| 186 | * Experimental Physics, vol 12B, p 186-200 (1976).
|
---|
| 187 | *
|
---|
| 188 | * @param[in] temperature temperature in K
|
---|
| 189 | * @return vapour saturation pressure (Pascals)
|
---|
| 190 | **/
|
---|
| 191 | static double wvSaturationPressure(double temperature);
|
---|
| 192 |
|
---|
| 193 | /**
|
---|
| 194 | * Compute the complex refractivity of the dry components of the atmosphere
|
---|
| 195 | * (oxygen lines) at the given frequency.
|
---|
| 196 | * @param[in] freq frequency (Hz)
|
---|
| 197 | * @param[in] temperature air temperature (K)
|
---|
| 198 | * @param[in] pDry partial pressure of dry components (Pascals)
|
---|
| 199 | * @param[in] pVapour partial pressure of water vapour (Pascals)
|
---|
| 200 | * @return complex refractivity
|
---|
| 201 | *
|
---|
| 202 | * Reference:
|
---|
| 203 | * Liebe, An updated model for millimeter wave propogation in moist air,
|
---|
| 204 | * Radio Science, 20, 1069-1089 (1985).
|
---|
| 205 | **/
|
---|
| 206 | static std::complex<double> dryRefractivity(double freq, double temperature,
|
---|
| 207 | double pDry, double pVapour);
|
---|
| 208 |
|
---|
| 209 | /**
|
---|
| 210 | * Compute the complex refractivity of the water vapour monomers
|
---|
| 211 | * at the given frequency.
|
---|
| 212 | * @param[in] freq frequency (Hz)
|
---|
| 213 | * @param[in] temperature air temperature (K)
|
---|
| 214 | * @param[in] pDry partial pressure of dry components (Pascals)
|
---|
| 215 | * @param[in] pVapour partial pressure of water vapour (Pascals)
|
---|
| 216 | * @return complex refractivity
|
---|
| 217 | *
|
---|
| 218 | * Reference:
|
---|
| 219 | * Liebe, An updated model for millimeter wave propogation in moist air,
|
---|
| 220 | * Radio Science, 20, 1069-1089 (1985).
|
---|
| 221 | **/
|
---|
| 222 | static std::complex<double> vapourRefractivity(double freq, double temperature,
|
---|
| 223 | double pDry, double pVapour);
|
---|
| 224 |
|
---|
| 225 | private:
|
---|
| 226 |
|
---|
| 227 | // heights of all model layers
|
---|
| 228 | std::vector<double> itsHeights;
|
---|
| 229 |
|
---|
| 230 | // temperatures of all model layers
|
---|
| 231 | std::vector<double> itsTemperatures;
|
---|
| 232 |
|
---|
| 233 | // partial pressures of dry component for all model layers
|
---|
| 234 | std::vector<double> itsDryPressures;
|
---|
| 235 |
|
---|
| 236 | // partial pressure of water vapour for all model layers
|
---|
| 237 | std::vector<double> itsVapourPressures;
|
---|
| 238 |
|
---|
| 239 | /**
|
---|
| 240 | * Atmosphere parameters
|
---|
| 241 | **/
|
---|
| 242 |
|
---|
| 243 | // ground level temperature (K)
|
---|
| 244 | double itsGndTemperature;
|
---|
| 245 |
|
---|
| 246 | // sea level pressure (Pascals)
|
---|
| 247 | double itsPressure;
|
---|
| 248 |
|
---|
| 249 | // ground level humidity (fraction)
|
---|
| 250 | double itsGndHumidity;
|
---|
| 251 |
|
---|
| 252 | // lapse rate (K/m)
|
---|
| 253 | double itsLapseRate;
|
---|
| 254 |
|
---|
| 255 | // water vapour scale height (m)
|
---|
| 256 | double itsWVScale;
|
---|
| 257 |
|
---|
| 258 | // altitude of the highest layer of the model (m)
|
---|
| 259 | double itsMaxAlt;
|
---|
| 260 |
|
---|
| 261 | // observatory elevation (m)
|
---|
| 262 | double itsObsHeight;
|
---|
| 263 | };
|
---|
| 264 |
|
---|
| 265 | } // namespace asap
|
---|
| 266 |
|
---|
| 267 | #endif // #ifndef STATMOSPHERE_H
|
---|
| 268 |
|
---|