1 //============================================================================== 2 // 3 // This file is part of GPSTk, the GPS Toolkit. 4 // 5 // The GPSTk is free software; you can redistribute it and/or modify 6 // it under the terms of the GNU Lesser General Public License as published 7 // by the Free Software Foundation; either version 3.0 of the License, or 8 // any later version. 9 // 10 // The GPSTk is distributed in the hope that it will be useful, 11 // but WITHOUT ANY WARRANTY; without even the implied warranty of 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 // GNU Lesser General Public License for more details. 14 // 15 // You should have received a copy of the GNU Lesser General Public 16 // License along with GPSTk; if not, write to the Free Software Foundation, 17 // Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA 18 // 19 // This software was developed by Applied Research Laboratories at the 20 // University of Texas at Austin. 21 // Copyright 2004-2020, The Board of Regents of The University of Texas System 22 // 23 //============================================================================== 24 25 //============================================================================== 26 // 27 // This software was developed by Applied Research Laboratories at the 28 // University of Texas at Austin, under contract to an agency or agencies 29 // within the U.S. Department of Defense. The U.S. Government retains all 30 // rights to use, duplicate, distribute, disclose, or release this software. 31 // 32 // Pursuant to DoD Directive 523024 33 // 34 // DISTRIBUTION STATEMENT A: This software has been approved for public 35 // release, distribution is unlimited. 36 // 37 //============================================================================== 38 39 /** 40 * @file TropModel.hpp 41 * Base class for tropospheric models, plus implementations 42 * of several published models 43 */ 44 45 #ifndef TROP_MODEL_HPP 46 #define TROP_MODEL_HPP 47 48 #include "Exception.hpp" 49 #include "ObsEpochMap.hpp" 50 #include "WxObsMap.hpp" 51 #include "Xvt.hpp" 52 #include "Position.hpp" 53 #include "Matrix.hpp" 54 #include "GNSSconstants.hpp" 55 56 #define THROW_IF_INVALID() {if (!valid) {InvalidTropModel e("Invalid model");GPSTK_THROW(e);}} 57 58 // Model of the troposphere, used to compute non-dispersive delay of 59 // satellite signal as function of satellite elevation as seen at the 60 // receiver. Both wet and hydrostatic (dry) components are computed. 61 // 62 // The default model (implemented here) is a simple Black model. 63 // 64 // In this model (and many others), the wet and hydrostatic (dry) components are 65 // independent, the zenith delays depend only on the weather 66 // (temperature, pressure and humidity), and the mapping functions 67 // depend only on the elevation of the satellite as seen at the 68 // receiver. In general, this is not true; other models depend on, 69 // for example, latitude or day of year. 70 // 71 // Other models may be implemented by inheriting this class and 72 // redefining the virtual functions, and (perhaps) adding other 73 // 'set...()' routines as needed. 74 75 namespace gpstk 76 { 77 /** @addtogroup GPSsolutions */ 78 //@{ 79 80 81 /** Thrown when attempting to use a model for which all necessary 82 * parameters have not been specified. 83 * @ingroup exceptiongroup 84 */ 85 NEW_EXCEPTION_CLASS(InvalidTropModel, gpstk::Exception); 86 87 /** Abstract base class for tropospheric models. 88 * The wet and hydrostatic (dry) components of the tropospheric 89 * delay are each the product of a zenith delay and a mapping 90 * function. Usually the zenith delay depends only on the 91 * weather (temperature, pressure and humidity), while the 92 * mapping function depends only on the satellite elevation, 93 * i.e. the geometry of satellite and receiver. This may not be 94 * true in complex models. 95 * 96 * The full tropospheric delay is the sum of the wet and 97 * hydrostatic (dry) components. A TropModel is valid only when 98 * all the necessary information (weather + whatever else the 99 * model requires) is specified; An InvalidTropModel exception 100 * will be thrown when any correction() or zenith_delay() or 101 * mapping_function() routine is called for an invalid 102 * TropModel. 103 */ 104 class TropModel 105 { 106 public: 107 static const double CELSIUS_TO_KELVIN; 108 109 /// Destructor ~TropModel()110 virtual ~TropModel() {} 111 112 /// Return validity of model isValid(void)113 bool isValid(void) 114 { return valid; } 115 116 /// Return the name of the model name(void)117 virtual std::string name(void) 118 { return std::string("Undefined"); } 119 120 /** Compute and return the full tropospheric delay 121 * @param elevation Elevation of satellite as seen at 122 * receiver, in degrees 123 * @throw InvalidTropModel 124 */ 125 virtual double correction(double elevation) const; 126 127 /** 128 * Compute and return the full tropospheric delay, given the 129 * positions of receiver and satellite and the time tag. This 130 * version is most useful within positioning algorithms, 131 * where the receiver position and timetag may vary; it 132 * computes the elevation (and other receiver location 133 * information) and passes them to appropriate set...() 134 * routines and the correction(elevation) routine. 135 * @param RX Receiver position 136 * @param SV Satellite position 137 * @param tt Time tag of the signal 138 * @throw InvalidTropModel 139 */ 140 virtual double correction(const Position& RX, 141 const Position& SV, 142 const CommonTime& tt); 143 144 /** \deprecated 145 * Compute and return the full tropospheric delay, given the 146 * positions of receiver and satellite and the time tag. This 147 * version is most useful within positioning algorithms, 148 * where the receiver position and timetag may vary; it 149 * computes the elevation (and other receiver location 150 * information) and passes them to appropriate set...() 151 * routines and the correction(elevation) routine. 152 * @param RX Receiver position in ECEF cartesian coordinates (meters) 153 * @param SV Satellite position in ECEF cartesian coordinates (meters) 154 * @param tt Time tag of the signal 155 * @throw InvalidTropModel 156 */ correction(const Xvt & RX,const Xvt & SV,const CommonTime & tt)157 virtual double correction(const Xvt& RX, 158 const Xvt& SV, 159 const CommonTime& tt) 160 { Position R(RX),S(SV); return TropModel::correction(R,S,tt); } 161 162 /** Compute and return the zenith delay for hydrostatic (dry) 163 * component of the troposphere 164 * @throw InvalidTropModel 165 */ 166 virtual double dry_zenith_delay(void) const = 0; 167 168 /** Compute and return the zenith delay for wet component of 169 * the troposphere. 170 * @throw InvalidTropModel 171 */ 172 virtual double wet_zenith_delay(void) const = 0; 173 174 /** Compute and return the mapping function for hydrostatic (dry) 175 * component of the troposphere. 176 * @param elevation Elevation of satellite as seen at 177 * receiver, in degrees 178 * @throw InvalidTropModel 179 */ 180 virtual double dry_mapping_function(double elevation) 181 const = 0; 182 183 /** Compute and return the mapping function for wet component of 184 * the troposphere. 185 * @param elevation Elevation of satellite as seen at 186 * receiver, in degrees 187 * @throw InvalidTropModel 188 */ 189 virtual double wet_mapping_function(double elevation) 190 const = 0; 191 192 /** Re-define the tropospheric model with explicit weather data. 193 * Typically called just before correction(). 194 * @param T temperature in degrees Celsius 195 * @param P atmospheric pressure in millibars 196 * @param H relative humidity in percent 197 * @throw InvalidParameter 198 */ 199 virtual void setWeather(const double& T, 200 const double& P, 201 const double& H); 202 203 /** Re-define the tropospheric model with explicit weather data. 204 * Typically called just before correction(). 205 * @param wx the weather to use for this correction 206 * @throw InvalidParameter 207 */ 208 virtual void setWeather(const WxObservation& wx); 209 210 /** Define the receiver height; this required by some models 211 * before calling correction() or any of the zenith_delay or 212 * mapping_function routines. 213 * @param ht Height of the receiver in meters. 214 */ setReceiverHeight(const double & ht)215 virtual void setReceiverHeight(const double& ht) {} 216 217 /** Define the latitude of the receiver; this is required by 218 * some models before calling correction() or any of the 219 * zenith_delay or mapping_function routines. 220 * @param lat Latitude of the receiver in degrees. 221 */ setReceiverLatitude(const double & lat)222 virtual void setReceiverLatitude(const double& lat) {} 223 224 /** Define the receiver longitude; this is required by some models 225 * before calling correction() or any of the zenith_delay routines. 226 * @param lat Longitude of receiver, in degrees East. 227 */ setReceiverLongitude(const double & lon)228 virtual void setReceiverLongitude(const double& lon) {} 229 230 /** Define the day of year; this is required by some models 231 * before calling correction() or any of the zenith_delay or 232 * mapping_function routines. 233 * @param d Day of year. 234 */ setDayOfYear(const int & d)235 virtual void setDayOfYear(const int& d) {} 236 237 /** Saastamoinen hydrostatic zenith delay as modified by 238 * Davis for gravity. 239 * Used by multiple models. 240 * Ref. Leick, 3rd ed, pg 197, Leick, 4th ed, pg 482, and 241 * Saastamoinen 1973 Atmospheric correction for the 242 * troposphere and stratosphere in radio ranging of 243 * satellites. The use of artificial satellites for geodesy, 244 * Geophys. Monogr. Ser. 15, Amer. Geophys. Union, 245 * pp. 274-251, 1972. 246 * Davis, J.L, T.A. Herring, I.I. Shapiro, A.E.E. Rogers, and 247 * G. Elgered, Geodesy by Radio Interferometry: Effects of 248 * Atmospheric Modeling Errors on Estimates of Baseline 249 * Length, Radio Science, Vol. 20, No. 6, pp. 1593-1607, 250 * 1985. 251 * @param press pressure in millibars 252 * @param lat latitude in degrees 253 * @param height ellipsoid height in meters 254 */ SaasDryDelay(const double pr,const double lat,const double ht) const255 double SaasDryDelay(const double pr, const double lat, const double ht) 256 const 257 { 258 return (0.0022768*pr/(1-0.00266*::cos(2*lat*DEG_TO_RAD)-0.00028*ht/1000.)); 259 } 260 261 /** get weather data by a standard atmosphere model 262 * reference to white paper of Bernese 5.0, P243 263 * @param ht height of the receiver in meters. 264 * @param T temperature in degrees Celsius 265 * @param P atmospheric pressure in millibars 266 * @param H relative humidity in percent 267 */ 268 static void weatherByStandardAtmosphereModel( 269 const double& ht, double& T, double& P, double& H); 270 271 protected: 272 bool valid; ///< true only if current model parameters are valid 273 double temp; ///< latest value of temperature (kelvin or celsius) 274 double press; ///< latest value of pressure (millibars) 275 double humid; ///< latest value of relative humidity (percent) 276 277 }; // end class TropModel 278 279 280 //--------------------------------------------------------------------------------- 281 /// The 'zero' trop model, meaning it always returns zero. 282 class ZeroTropModel : public TropModel 283 { 284 public: 285 /// Return the name of the model name(void)286 virtual std::string name(void) 287 { return std::string("Zero"); } 288 289 /** Compute and return the full tropospheric delay 290 * @param elevation Elevation of satellite as seen at receiver, in degrees 291 * @throw InvalidTropModel 292 */ correction(double elevation) const293 virtual double correction(double elevation) const 294 { return 0.0; } 295 296 /** 297 * Compute and return the full tropospheric delay, given the 298 * positions of receiver and satellite and the time tag. This 299 * version is most useful within positioning algorithms, 300 * where the receiver position and timetag may vary; it 301 * computes the elevation (and other receiver location 302 * information) and passes them to appropriate set...() 303 * routines and the correction(elevation) routine. 304 * @param RX Receiver position 305 * @param SV Satellite position 306 * @param tt Time tag of the signal 307 * @throw InvalidTropModel 308 */ correction(const Position & RX,const Position & SV,const CommonTime & tt)309 virtual double correction(const Position& RX, 310 const Position& SV, 311 const CommonTime& tt) 312 { return 0.0; } 313 314 /** \deprecated 315 * Compute and return the full tropospheric delay, given the 316 * positions of receiver and satellite and the time tag. This 317 * version is most useful within positioning algorithms, 318 * where the receiver position and timetag may vary; it 319 * computes the elevation (and other receiver location 320 * information) and passes them to appropriate set...() 321 * routines and the correction(elevation) routine. 322 * @param RX Receiver position in ECEF cartesian coordinates (meters) 323 * @param SV Satellite position in ECEF cartesian coordinates (meters) 324 * @param tt Time tag of the signal 325 * @throw InvalidTropModel 326 */ correction(const Xvt & RX,const Xvt & SV,const CommonTime & tt)327 virtual double correction(const Xvt& RX, 328 const Xvt& SV, 329 const CommonTime& tt) 330 { return 0.0; } 331 332 /** Compute and return the zenith delay for hydrostatic (dry) 333 * component of the troposphere 334 * @throw InvalidTropModel 335 */ dry_zenith_delay(void) const336 virtual double dry_zenith_delay(void) const 337 { return 0.0; } 338 339 /** Compute and return the zenith delay for wet component of 340 * the troposphere 341 * @throw InvalidTropModel 342 */ wet_zenith_delay(void) const343 virtual double wet_zenith_delay(void) const 344 { return 0.0; } 345 346 /** Compute and return the mapping function for hydrostatic (dry) 347 * component of the troposphere. 348 * @param elevation Elevation of satellite as seen at 349 * receiver, in degrees 350 * @throw InvalidTropModel 351 */ dry_mapping_function(double elevation) const352 virtual double dry_mapping_function(double elevation) const 353 { return 0.0; } 354 355 /** Compute and return the mapping function for wet component of 356 * the troposphere. 357 * @param elevation Elevation of satellite as seen at 358 * receiver, in degrees 359 * @throw InvalidTropModel 360 */ wet_mapping_function(double elevation) const361 virtual double wet_mapping_function(double elevation) const 362 { return 0.0; } 363 364 }; // end class ZeroTropModel 365 366 //@} 367 } 368 #endif 369