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 #include "YDSTime.hpp" 40 #include "NBTropModel.hpp" 41 42 #define THROW_IF_INVALID_DETAILED() {if (!valid) { \ 43 InvalidTropModel e; \ 44 if(!validWeather) e.addText("Invalid trop model: weather"); \ 45 if(!validRxLatitude) e.addText("Invalid trop model: validRxLatitude"); \ 46 if(!validRxHeight) e.addText("Invalid trop model: validRxHeight"); \ 47 if(!validDOY) e.addText("Invalid trop model: day of year"); \ 48 GPSTK_THROW(e);}} 49 50 namespace gpstk 51 { 52 // ------------------------------------------------------------------------ 53 // Tropospheric model developed by University of New Brunswick and described in 54 // "A Tropospheric Delay Model for the User of the Wide Area Augmentation 55 // System," J. Paul Collins and Richard B. Langley, Technical Report No. 187, 56 // Dept. of Geodesy and Geomatics Engineering, University of New Brunswick, 57 // 1997. See particularly Appendix C. 58 // 59 // This model includes a wet and dry component, and was designed for the user 60 // without access to measurements of temperature, pressure and relative humidity 61 // at ground level. It requires input of the latitude, day of year and height 62 // above the ellipsoid, and it interpolates a table of values, using these 63 // inputs, to get the ground level weather parameters plus other parameters and 64 // the mapping function constants. 65 // 66 // NB in this model, units of temp are degrees Celsius, and humid is the water 67 // vapor partial pressure. 68 69 static const double NBRd=287.054; // J/kg/K = m*m*K/s*s 70 static const double NBg=9.80665; // m/s*s 71 static const double NBk1=77.604; // K/mbar 72 static const double NBk3p=382000; // K*K/mbar 73 74 static const double NBLat[5]={ 15.0, 30.0, 45.0, 60.0, 75.0}; 75 76 // zenith delays, averages 77 static const double NBZP0[5]={1013.25,1017.25,1015.75,1011.75,1013.00}; 78 static const double NBZT0[5]={ 299.65, 294.15, 283.15, 272.15, 263.65}; 79 static const double NBZW0[5]={ 26.31, 21.79, 11.66, 6.78, 4.11}; 80 static const double NBZB0[5]={6.30e-3,6.05e-3,5.58e-3,5.39e-3,4.53e-3}; 81 static const double NBZL0[5]={ 2.77, 3.15, 2.57, 1.81, 1.55}; 82 83 // zenith delays, amplitudes 84 static const double NBZPa[5]={ 0.0, -3.75, -2.25, -1.75, -0.50}; 85 static const double NBZTa[5]={ 0.0, 7.0, 11.0, 15.0, 14.5}; 86 static const double NBZWa[5]={ 0.0, 8.85, 7.24, 5.36, 3.39}; 87 static const double NBZBa[5]={ 0.0,0.25e-3,0.32e-3,0.81e-3,0.62e-3}; 88 static const double NBZLa[5]={ 0.0, 0.33, 0.46, 0.74, 0.30}; 89 90 // mapping function, dry, averages 91 static const double NBMad[5]={1.2769934e-3,1.2683230e-3,1.2465397e-3,1.2196049e-3, 92 1.2045996e-3}; 93 static const double NBMbd[5]={2.9153695e-3,2.9152299e-3,2.9288445e-3,2.9022565e-3, 94 2.9024912e-3}; 95 static const double NBMcd[5]={62.610505e-3,62.837393e-3,63.721774e-3,63.824265e-3, 96 64.258455e-3}; 97 98 // mapping function, dry, amplitudes 99 static const double NBMaa[5]={0.0,1.2709626e-5,2.6523662e-5,3.4000452e-5, 100 4.1202191e-5}; 101 static const double NBMba[5]={0.0,2.1414979e-5,3.0160779e-5,7.2562722e-5, 102 11.723375e-5}; 103 static const double NBMca[5]={0.0,9.0128400e-5,4.3497037e-5,84.795348e-5, 104 170.37206e-5}; 105 106 // mapping function, wet, averages (there are no amplitudes for wet) 107 static const double NBMaw[5]={5.8021897e-4,5.6794847e-4,5.8118019e-4,5.9727542e-4, 108 6.1641693e-4}; 109 static const double NBMbw[5]={1.4275268e-3,1.5138625e-3,1.4572752e-3,1.5007428e-3, 110 1.7599082e-3}; 111 static const double NBMcw[5]={4.3472961e-2,4.6729510e-2,4.3908931e-2,4.4626982e-2, 112 5.4736038e-2}; 113 114 // labels for use with the interpolation routine 115 enum TableEntry { ZP=1, ZT, ZW, ZB, ZL, Mad, Mbd, Mcd, Maw, Mbw, Mcw }; 116 117 // the interpolation routine NB_Interpolate(double lat,int doy,TableEntry entry)118 static double NB_Interpolate(double lat, int doy, TableEntry entry) 119 { 120 const double *pave = NULL, *pamp = NULL; 121 double ret, day=double(doy); 122 123 // assign pointer to the right array 124 switch(entry) { 125 case ZP: pave=&NBZP0[0]; pamp=&NBZPa[0]; break; 126 case ZT: pave=&NBZT0[0]; pamp=&NBZTa[0]; break; 127 case ZW: pave=&NBZW0[0]; pamp=&NBZWa[0]; break; 128 case ZB: pave=&NBZB0[0]; pamp=&NBZBa[0]; break; 129 case ZL: pave=&NBZL0[0]; pamp=&NBZLa[0]; break; 130 case Mad: pave=&NBMad[0]; pamp=&NBMaa[0]; break; 131 case Mbd: pave=&NBMbd[0]; pamp=&NBMba[0]; break; 132 case Mcd: pave=&NBMcd[0]; pamp=&NBMca[0]; break; 133 case Maw: pave=&NBMaw[0]; break; 134 case Mbw: pave=&NBMbw[0]; break; 135 case Mcw: pave=&NBMcw[0]; break; 136 } 137 138 // find the index i such that NBLat[i] <= lat < NBLat[i+1] 139 int i = int(ABS(lat)/15.0)-1; 140 141 if(i>=0 && i<=3) { // mid-latitude -> regular interpolation 142 double m=(ABS(lat)-NBLat[i])/(NBLat[i+1]-NBLat[i]); 143 ret = pave[i]+m*(pave[i+1]-pave[i]); 144 if(entry < Maw) 145 ret -= (pamp[i]+m*(pamp[i+1]-pamp[i])) 146 * std::cos(TWO_PI*(day-28.0)/365.25); 147 } 148 else { // < 15 or > 75 -> simpler interpolation 149 if(i<0) i=0; else i=4; 150 ret = pave[i]; 151 if(entry < Maw) 152 ret -= pamp[i]*std::cos(TWO_PI*(day-28.0)/365.25); 153 } 154 155 return ret; 156 157 } // end double NB_Interpolate(lat,doy,entry) 158 159 // Default constructor NBTropModel(void)160 NBTropModel::NBTropModel(void): 161 validWeather(false), validRxLatitude(false), validDOY(false),validRxHeight(false) 162 {} 163 164 // Create a trop model using the minimum information: latitude and doy. 165 // Interpolate the weather unless setWeather (optional) is called. 166 // @param lat Latitude of the receiver in degrees. 167 // @param day Day of year. NBTropModel(const double & lat,const int & day)168 NBTropModel::NBTropModel(const double& lat, 169 const int& day) 170 : validWeather(false), validRxLatitude(false), validDOY(false), 171 validRxHeight(false) 172 { 173 setReceiverLatitude(lat); 174 setDayOfYear(day); 175 setWeather(); 176 } 177 178 // Create a trop model with weather. 179 // @param lat Latitude of the receiver in degrees. 180 // @param day Day of year. 181 // @param wx the weather to use for this correction. NBTropModel(const double & lat,const int & day,const WxObservation & wx)182 NBTropModel::NBTropModel(const double& lat, 183 const int& day, 184 const WxObservation& wx) 185 : validWeather(false), validRxLatitude(false), validDOY(false), 186 validRxHeight(false) 187 { 188 setReceiverLatitude(lat); 189 setDayOfYear(day); 190 setWeather(wx); 191 } 192 193 // Create a tropospheric model from explicit weather data 194 // @param lat Latitude of the receiver in degrees. 195 // @param day Day of year. 196 // @param T temperature in degrees Celsius 197 // @param P atmospheric pressure in millibars 198 // @param H relative humidity in percent NBTropModel(const double & lat,const int & day,const double & T,const double & P,const double & H)199 NBTropModel::NBTropModel(const double& lat, 200 const int& day, 201 const double& T, 202 const double& P, 203 const double& H) 204 : validWeather(false), validRxLatitude(false), validDOY(false), 205 validRxHeight(false) 206 { 207 setReceiverLatitude(lat); 208 setDayOfYear(day); 209 setWeather(T,P,H); 210 } 211 212 // Create a valid model from explicit input (weather will be estimated 213 // internally by this model). 214 // @param ht Height of the receiver in meters. 215 // @param lat Latitude of the receiver in degrees. 216 // @param d Day of year. NBTropModel(const double & ht,const double & lat,const int & day)217 NBTropModel::NBTropModel(const double& ht, 218 const double& lat, 219 const int& day) 220 : validWeather(false), validRxLatitude(false), validDOY(false), 221 validRxHeight(false) 222 { 223 setReceiverHeight(ht); 224 setReceiverLatitude(lat); 225 setDayOfYear(day); 226 setWeather(); 227 } 228 229 // re-define this to get the throws correction(double elevation) const230 double NBTropModel::correction(double elevation) const 231 { 232 THROW_IF_INVALID_DETAILED(); 233 234 if(elevation < 0.0) return 0.0; 235 236 return (dry_zenith_delay() * dry_mapping_function(elevation) 237 + wet_zenith_delay() * wet_mapping_function(elevation)); 238 } 239 240 // Compute and return the full tropospheric delay, given the positions of 241 // receiver and satellite and the time tag. This version is most useful 242 // within positioning algorithms, where the receiver position and timetag 243 // may vary; it computes the elevation (and other receiver location 244 // information) and passes them to appropriate set...() routines 245 // and the correction(elevation) routine. 246 // @param RX Receiver position 247 // @param SV Satellite position 248 // @param tt Time tag of the signal correction(const Position & RX,const Position & SV,const CommonTime & tt)249 double NBTropModel::correction(const Position& RX, 250 const Position& SV, 251 const CommonTime& tt) 252 { 253 THROW_IF_INVALID_DETAILED(); 254 255 // compute height and latitude from RX 256 setReceiverHeight(RX.getHeight()); 257 setReceiverLatitude(RX.getGeodeticLatitude()); 258 259 // compute day of year from tt 260 setDayOfYear(int((static_cast<YDSTime>(tt)).doy)); 261 262 return TropModel::correction(RX.elevation(SV)); 263 } 264 265 // Compute and return the full tropospheric delay, given the positions of 266 // receiver and satellite and the time tag. This version is most useful 267 // within positioning algorithms, where the receiver position and timetag 268 // may vary; it computes the elevation (and other receiver location 269 // information) and passes them to appropriate set...() routines 270 // and the correction(elevation) routine. 271 // @param RX Receiver position in ECEF cartesian coordinates (meters) 272 // @param SV Satellite position in ECEF cartesian coordinates (meters) 273 // @param tt Time tag of the signal 274 // This function is deprecated; use the Position version correction(const Xvt & RX,const Xvt & SV,const CommonTime & tt)275 double NBTropModel::correction(const Xvt& RX, 276 const Xvt& SV, 277 const CommonTime& tt) 278 { 279 Position R(RX),S(SV); 280 return NBTropModel::correction(R,S,tt); 281 } 282 283 // Compute and return the zenith delay for dry component of the troposphere dry_zenith_delay(void) const284 double NBTropModel::dry_zenith_delay(void) const 285 { 286 THROW_IF_INVALID_DETAILED(); 287 288 double beta = NB_Interpolate(latitude,doy,ZB); 289 double gm = 9.784*(1.0-2.66e-3*std::cos(2.0*latitude*DEG_TO_RAD)-2.8e-7*height); 290 291 // scale factors for height above mean sea level 292 // if weather is given, assume it's measured at ht -> kw=kd=1 293 double kd=1, base=std::log(1.0-beta*height/temp); 294 if(interpolateWeather) 295 kd = std::exp(base*NBg/(NBRd*beta)); 296 297 // compute the zenith delay for dry component 298 return ((1.0e-6*NBk1*NBRd/gm) * kd * press); 299 300 } // end NBTropModel::dry_zenith_delay() 301 302 // Compute and return the zenith delay for wet component of the troposphere wet_zenith_delay(void) const303 double NBTropModel::wet_zenith_delay(void) const 304 { 305 THROW_IF_INVALID_DETAILED(); 306 307 double beta = NB_Interpolate(latitude,doy,ZB); 308 double lam = NB_Interpolate(latitude,doy,ZL); 309 double gm = 9.784*(1.0-2.66e-3*std::cos(2.0*latitude*DEG_TO_RAD)-2.8e-7*height); 310 311 // scale factors for height above mean sea level 312 // if weather is given, assume it's measured at ht -> kw=kd=1 313 double kw=1, base=std::log(1.0-beta*height/temp); 314 if(interpolateWeather) 315 kw = std::exp(base*(-1.0+(lam+1)*NBg/(NBRd*beta))); 316 317 // compute the zenith delay for wet component 318 return ((1.0e-6*NBk3p*NBRd/(gm*(lam+1)-beta*NBRd)) * kw * humid/temp); 319 320 } // end NBTropModel::wet_zenith_delay() 321 322 // Compute and return the mapping function for dry component 323 // of the troposphere 324 // @param elevation Elevation of satellite as seen at receiver, 325 // in degrees dry_mapping_function(double elevation) const326 double NBTropModel::dry_mapping_function(double elevation) const 327 { 328 THROW_IF_INVALID_DETAILED(); 329 330 if(elevation < 0.0) return 0.0; 331 332 double a,b,c,se,map; 333 se = std::sin(elevation*DEG_TO_RAD); 334 335 a = NB_Interpolate(latitude,doy,Mad); 336 b = NB_Interpolate(latitude,doy,Mbd); 337 c = NB_Interpolate(latitude,doy,Mcd); 338 map = (1.0+a/(1.0+b/(1.0+c))) / (se+a/(se+b/(se+c))); 339 340 a = 2.53e-5; 341 b = 5.49e-3; 342 c = 1.14e-3; 343 if(ABS(elevation)<=0.001) se=0.001; 344 map += ((1.0/se)-(1.0+a/(1.0+b/(1.0+c)))/(se+a/(se+b/(se+c))))*height/1000.0; 345 346 return map; 347 348 } // end NBTropModel::dry_mapping_function() 349 350 // Compute and return the mapping function for wet component 351 // of the troposphere 352 // @param elevation Elevation of satellite as seen at receiver, 353 // in degrees wet_mapping_function(double elevation) const354 double NBTropModel::wet_mapping_function(double elevation) const 355 { 356 THROW_IF_INVALID_DETAILED(); 357 358 if(elevation < 0.0) return 0.0; 359 360 double a,b,c,se; 361 se = std::sin(elevation*DEG_TO_RAD); 362 a = NB_Interpolate(latitude,doy,Maw); 363 b = NB_Interpolate(latitude,doy,Mbw); 364 c = NB_Interpolate(latitude,doy,Mcw); 365 366 return ( (1.0+a/(1.0+b/(1.0+c))) / (se+a/(se+b/(se+c))) ); 367 368 } // end NBTropModel::wet_mapping_function() 369 370 // Re-define the weather data. 371 // If called, typically called before any calls to correction(). 372 // @param T temperature in degrees Celsius 373 // @param P atmospheric pressure in millibars 374 // @param H relative humidity in percent setWeather(const double & T,const double & P,const double & H)375 void NBTropModel::setWeather(const double& T, 376 const double& P, 377 const double& H) 378 { 379 interpolateWeather=false; 380 TropModel::setWeather(T,P,H); 381 // humid actually stores water vapor partial pressure 382 double th=300./temp; 383 humid = 2.409e9*H*th*th*th*th*std::exp(-22.64*th); 384 validWeather = true; 385 valid = validWeather && validRxHeight && validRxLatitude && validDOY; 386 387 } // end NBTropModel::setWeather() 388 389 // Re-define the tropospheric model with explicit weather data. 390 // Typically called just before correction(). 391 // @param wx the weather to use for this correction setWeather(const WxObservation & wx)392 void NBTropModel::setWeather(const WxObservation& wx) 393 { 394 interpolateWeather = false; 395 try 396 { 397 TropModel::setWeather(wx); 398 // humid actually stores vapor partial pressure 399 double th=300./temp; 400 humid = 2.409e9*humid*th*th*th*th*std::exp(-22.64*th); 401 validWeather = true; 402 valid = validWeather && validRxHeight && validRxLatitude && validDOY; 403 } 404 catch(InvalidParameter& e) 405 { 406 valid = validWeather = false; 407 GPSTK_RETHROW(e); 408 } 409 } 410 411 // configure the model to estimate the weather from the internal model, 412 // using lat and doy setWeather()413 void NBTropModel::setWeather() 414 { 415 interpolateWeather = true; 416 417 if(!validRxLatitude) 418 { 419 valid = validWeather = false; 420 InvalidTropModel e("NBTropModel must have Rx latitude before interpolating weather"); 421 GPSTK_THROW(e); 422 } 423 if(!validDOY) 424 { 425 valid = validWeather = false; 426 InvalidTropModel e("NBTropModel must have day of year before interpolating weather "); 427 GPSTK_THROW(e); 428 } 429 temp = NB_Interpolate(latitude,doy,ZT); 430 press = NB_Interpolate(latitude,doy,ZP); 431 humid = NB_Interpolate(latitude,doy,ZW); 432 validWeather = true; 433 valid = validWeather && validRxHeight && validRxLatitude && validDOY; 434 } 435 436 // Define the receiver height; this required before calling 437 // correction() or any of the zenith_delay or mapping_function routines. setReceiverHeight(const double & ht)438 void NBTropModel::setReceiverHeight(const double& ht) 439 { 440 height = ht; 441 validRxHeight = true; 442 valid = validWeather && validRxHeight && validRxLatitude && validDOY; 443 if(!validWeather && validRxLatitude && validDOY) 444 setWeather(); 445 } 446 447 // Define the latitude of the receiver; this is required before calling 448 // correction() or any of the zenith_delay or mapping_function routines. setReceiverLatitude(const double & lat)449 void NBTropModel::setReceiverLatitude(const double& lat) 450 { 451 latitude = lat; 452 validRxLatitude = true; 453 valid = validWeather && validRxHeight && validRxLatitude && validDOY; 454 if(!validWeather && validRxLatitude && validDOY) 455 setWeather(); 456 } 457 458 // Define the day of year; this is required before calling 459 // correction() or any of the zenith_delay or mapping_function routines. setDayOfYear(const int & d)460 void NBTropModel::setDayOfYear(const int& d) 461 { 462 doy = d; 463 if (doy > 0 && doy < 367) validDOY=true; else validDOY = false; 464 valid = validWeather && validRxHeight && validRxLatitude && validDOY; 465 if(!validWeather && validRxLatitude && validDOY) 466 setWeather(); 467 } 468 469 } 470