SaasTropModel.cpp
Go to the documentation of this file.
1 //==============================================================================
2 //
3 // This file is part of GNSSTk, the ARL:UT GNSS Toolkit.
4 //
5 // The GNSSTk 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 GNSSTk 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 GNSSTk; 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-2022, 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 #include "YDSTime.hpp"
41 #include "SaasTropModel.hpp"
42 
43 #define THROW_IF_INVALID_DETAILED() {if (!valid) {\
44  InvalidTropModel e; \
45  if(!validWeather) e.addText("Invalid trop model: Weather"); \
46  if(!validRxHeight) e.addText("Invalid trop model: Rx Height"); \
47  if(!validRxLatitude) e.addText("Invalid trop model: Rx Latitude"); \
48  if(!validDOY) e.addText("Invalid trop model: day of year"); \
49  GNSSTK_THROW(e);}}
50 
51 namespace gnsstk
52 {
53  // constants for wet mapping function
54  static const double SaasWetA[5]=
55  { 0.00058021897, 0.00056794847, 0.00058118019, 0.00059727542, 0.00061641693 };
56  static const double SaasWetB[5]=
57  { 0.0014275268, 0.0015138625, 0.0014572752, 0.0015007428, 0.0017599082 };
58  static const double SaasWetC[5]=
59  { 0.043472961, 0.046729510, 0.043908931, 0.044626982, 0.054736038 };
60 
61  // constants for dry mapping function
62  static const double SaasDryA[5]=
63  { 0.0012769934, 0.0012683230, 0.0012465397, 0.0012196049, 0.0012045996 };
64  static const double SaasDryB[5]=
65  { 0.0029153695, 0.0029152299, 0.0029288445, 0.0029022565, 0.0029024912 };
66  static const double SaasDryC[5]=
67  { 0.062610505, 0.062837393, 0.063721774, 0.063824265, 0.064258455 };
68 
69  static const double SaasDryA1[5]=
70  { 0.0, 0.000012709626, 0.000026523662, 0.000034000452, 0.000041202191 };
71  static const double SaasDryB1[5]=
72  { 0.0, 0.000021414979, 0.000030160779, 0.000072562722, 0.00011723375 };
73  static const double SaasDryC1[5]=
74  { 0.0, 0.000090128400, 0.000043497037, 0.00084795348, 0.0017037206 };
75 
76 
77  SaasTropModel::SaasTropModel()
78  {
79  validWeather = false;
80  validRxLatitude = false;
81  validDOY = false;
82  validRxHeight = false;
83  } // end SaasTropModel::SaasTropModel()
84 
85 
86  SaasTropModel::SaasTropModel(const double& lat,
87  const int& day)
88  {
89  validWeather = false;
90  validRxHeight = false;
91  SaasTropModel::setReceiverLatitude(lat);
92  SaasTropModel::setDayOfYear(day);
93  } // end SaasTropModel::SaasTropModel
94 
95 
96  SaasTropModel::SaasTropModel(const double& lat,
97  const int& day,
98  const WxObservation& wx)
99  {
100  validRxHeight = false;
101  SaasTropModel::setReceiverLatitude(lat);
102  SaasTropModel::setDayOfYear(day);
103  SaasTropModel::setWeather(wx);
104  } // end SaasTropModel::SaasTropModel(weather)
105 
106 
107  SaasTropModel::SaasTropModel(const double& lat,
108  const int& day,
109  const double& T,
110  const double& P,
111  const double& H)
112  {
113  validRxHeight = false;
114  SaasTropModel::setReceiverLatitude(lat);
115  SaasTropModel::setDayOfYear(day);
116  SaasTropModel::setWeather(T,P,H);
117  } // end SaasTropModel::SaasTropModel()
118 
119  // re-define this to get the throws correct
120  double SaasTropModel::correction(double elevation) const
121  {
122  if(!valid) {
123  if(!validWeather) GNSSTK_THROW(
124  InvalidTropModel("Invalid Saastamoinen trop model: weather"));
125  if(!validRxLatitude) GNSSTK_THROW(
126  InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
127  if(!validRxHeight) GNSSTK_THROW(
128  InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
129  if(!validDOY) GNSSTK_THROW(
130  InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
131  GNSSTK_THROW(
132  InvalidTropModel("Valid flag corrupted in Saastamoinen trop model"));
133  }
134 
135  if(elevation < 0.0) return 0.0;
136 
137  double corr=0.0;
138  try {
139  corr = (dry_zenith_delay() * dry_mapping_function(elevation)
140  + wet_zenith_delay() * wet_mapping_function(elevation));
141  }
142  catch(Exception& e) { GNSSTK_RETHROW(e); }
143 
144  return corr;
145 
146  } // end SaasTropModel::correction(elevation)
147 
148 
149  double SaasTropModel::correction(const Position& RX,
150  const Position& SV,
151  const CommonTime& tt)
152  {
153  SaasTropModel::setReceiverHeight(RX.getHeight());
154  SaasTropModel::setReceiverLatitude(RX.getGeodeticLatitude());
155  SaasTropModel::setDayOfYear(int((static_cast<YDSTime>(tt).doy)));
156 
157  if(!valid) {
158  if(!validWeather) GNSSTK_THROW(
159  InvalidTropModel("Invalid Saastamoinen trop model: weather"));
160  if(!validRxLatitude) GNSSTK_THROW(
161  InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
162  if(!validRxHeight) GNSSTK_THROW(
163  InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
164  if(!validDOY) GNSSTK_THROW(
165  InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
166  valid = true;
167  }
168 
169  double corr=0.0;
170  try {
171  corr = SaasTropModel::correction(RX.elevation(SV));
172  }
173  catch(Exception& e) { GNSSTK_RETHROW(e); }
174 
175  return corr;
176 
177  } // end SaasTropModel::correction(RX,SV,TT)
178 
179 
180  double SaasTropModel::correction(const Xvt& RX,
181  const Xvt& SV,
182  const CommonTime& tt)
183  {
184  Position R(RX),S(SV);
185  return SaasTropModel::correction(R,S,tt);
186  }
187 
188 
189  double SaasTropModel::dry_zenith_delay() const
190  {
191  THROW_IF_INVALID_DETAILED();
192 
193  //return (0.0022768*pr/(1-0.00266*::cos(2*lat*DEG_TO_RAD)-0.00028*ht/1000.));
194  return SaasDryDelay(press,latitude,height);
195 
196  } // end SaasTropModel::dry_zenith_delay()
197 
198 
199  double SaasTropModel::wet_zenith_delay() const
200  {
201  THROW_IF_INVALID_DETAILED();
202 
203  double T = temp+CELSIUS_TO_KELVIN;
204 
205  // partial pressure due to water vapor. Leick 4th ed 8.2.4
206  double pwv = 0.01 * humid * ::exp(-37.2465+0.213166*T-0.000256908*T*T);
207  /* IERS2003 Ch 9 pg 99 - very similar to Leick above
208  *double pwv = 0.01*humid
209  * * 0.01*::exp(33.93711047-1.9121316e-2*T+1.2378847e-5*T*T-6.3431645e3/T)
210  * * (1.00062+3.14e-6*press+5.6e-7*temp);
211  */
212 
213  /* Saastamoinen 1973 Atmospheric correction for the troposphere and
214  * stratosphere in radio ranging of satellites. The use of artificial
215  * satellites for geodesy, Geophys. Monogr. Ser. 15, Amer. Geophys. Union,
216  * pp. 274-251, 1972, modified for gravity as in Davis etal 1985
217  */
218  return ( 0.002277*((1255/T)+0.05)*pwv /
219  (1-0.00266*::cos(2*latitude*DEG_TO_RAD)-0.00028*height/1000.) );
220 
221  } // end SaasTropModel::wet_zenith_delay()
222 
223 
224  double SaasTropModel::dry_mapping_function(double elevation) const
225  {
226  THROW_IF_INVALID_DETAILED();
227  if(elevation < 0.0) return 0.0;
228 
229  double lat,t,ct;
230  lat = fabs(latitude); // degrees
231  t = doy - 28.; // mid-winter
232  if(latitude < 0) // southern hemisphere
233  t += 365.25/2.;
234  t *= 360.0/365.25; // convert to degrees
235  ct = ::cos(t*DEG_TO_RAD);
236 
237  double a,b,c;
238  if(lat < 15.) {
239  a = SaasDryA[0];
240  b = SaasDryB[0];
241  c = SaasDryC[0];
242  }
243  else if(lat < 75.) { // coefficients are for 15,30,45,60,75 deg
244  int i=int(lat/15.0)-1;
245  double frac=(lat-15.*(i+1))/15.;
246  a = SaasDryA[i] + frac*(SaasDryA[i+1]-SaasDryA[i]);
247  b = SaasDryB[i] + frac*(SaasDryB[i+1]-SaasDryB[i]);
248  c = SaasDryC[i] + frac*(SaasDryC[i+1]-SaasDryC[i]);
249 
250  a -= ct * (SaasDryA1[i] + frac*(SaasDryA1[i+1]-SaasDryA1[i]));
251  b -= ct * (SaasDryB1[i] + frac*(SaasDryB1[i+1]-SaasDryB1[i]));
252  c -= ct * (SaasDryC1[i] + frac*(SaasDryC1[i+1]-SaasDryC1[i]));
253  }
254  else {
255  a = SaasDryA[4] - ct * SaasDryA1[4];
256  b = SaasDryB[4] - ct * SaasDryB1[4];
257  c = SaasDryC[4] - ct * SaasDryC1[4];
258  }
259 
260  double se = ::sin(elevation*DEG_TO_RAD);
261  double map = (1.+a/(1.+b/(1.+c)))/(se+a/(se+b/(se+c)));
262 
263  a = 0.0000253;
264  b = 0.00549;
265  c = 0.00114;
266  map += (height/1000.0)*(1./se-(1+a/(1.+b/(1.+c)))/(se+a/(se+b/(se+c))));
267 
268  return map;
269 
270  } // end SaasTropModel::dry_mapping_function()
271 
272 
273  double SaasTropModel::wet_mapping_function(double elevation) const
274  {
275  THROW_IF_INVALID_DETAILED();
276  if(elevation < 0.0) return 0.0;
277 
278  double a,b,c,lat;
279  lat = fabs(latitude); // degrees
280  if(lat < 15.) {
281  a = SaasWetA[0];
282  b = SaasWetB[0];
283  c = SaasWetC[0];
284  }
285  else if(lat < 75.) { // coefficients are for 15,30,45,60,75 deg
286  int i=int(lat/15.0)-1;
287  double frac=(lat-15.*(i+1))/15.;
288  a = SaasWetA[i] + frac*(SaasWetA[i+1]-SaasWetA[i]);
289  b = SaasWetB[i] + frac*(SaasWetB[i+1]-SaasWetB[i]);
290  c = SaasWetC[i] + frac*(SaasWetC[i+1]-SaasWetC[i]);
291  }
292  else {
293  a = SaasWetA[4];
294  b = SaasWetB[4];
295  c = SaasWetC[4];
296  }
297 
298  double se = ::sin(elevation*DEG_TO_RAD);
299  double map = (1.+a/(1.+b/(1.+c)))/(se+a/(se+b/(se+c)));
300 
301  return map;
302 
303  }
304 
305 
306  void SaasTropModel::setWeather(const double& T,
307  const double& P,
308  const double& H)
309  {
310  temp = T;
311  press = P;
312  // humid actually stores water vapor partial pressure
313  double exp=7.5*T/(T+237.3);
314  humid = 6.11 * (H/100.) * std::pow(10.0,exp);
315 
316  validWeather = true;
317  valid = (validWeather && validRxHeight && validRxLatitude && validDOY);
318 
319  }
320 
321 
322  void SaasTropModel::setWeather(const WxObservation& wx)
323  {
324  try
325  {
326  SaasTropModel::setWeather(wx.temperature,wx.pressure,wx.humidity);
327  }
328  catch(InvalidParameter& e)
329  {
330  valid = validWeather = false;
331  GNSSTK_RETHROW(e);
332  }
333  }
334 
335 
336  void SaasTropModel::setReceiverHeight(const double& ht)
337  {
338  height = ht;
339  validRxHeight = true;
340  valid = (validWeather && validRxHeight && validRxLatitude && validDOY);
341  }
342 
343 
344  void SaasTropModel::setReceiverLatitude(const double& lat)
345  {
346  latitude = lat;
347  validRxLatitude = true;
348  valid = (validWeather && validRxHeight && validRxLatitude && validDOY);
349  }
350 
351 
352  void SaasTropModel::setDayOfYear(const int& d)
353  {
354  doy = d;
355  if(doy > 0 && doy < 367) validDOY=true; else validDOY = false;
356  valid = (validWeather && validRxHeight && validRxLatitude && validDOY);
357  }
358 }
YDSTime.hpp
se
double se
obliquity cos, T*cos, sin coefficients
Definition: IERS2003NutationData.hpp:48
gnsstk::SaasWetB
static const double SaasWetB[5]
Definition: SaasTropModel.cpp:56
gnsstk::SaasTropModel::validRxLatitude
bool validRxLatitude
flag for valid Rx latitude
Definition: SaasTropModel.hpp:169
example6.day
day
Definition: example6.py:66
gnsstk::SaasDryA
static const double SaasDryA[5]
Definition: SaasTropModel.cpp:62
gnsstk::YDSTime
Definition: YDSTime.hpp:58
gnsstk::SaasTropModel::doy
int doy
day of year
Definition: SaasTropModel.hpp:167
gnsstk::SaasTropModel::SaasTropModel
SaasTropModel()
Empty constructor.
Definition: SaasTropModel.cpp:77
gnsstk::SaasTropModel::height
double height
height (m) of the receiver above the geoid
Definition: SaasTropModel.hpp:165
SaasTropModel.hpp
gnsstk::WxObservation::temperature
float temperature
degrees Centigrade
Definition: WxObsMap.hpp:82
gnsstk::SaasWetA
static const double SaasWetA[5]
Definition: SaasTropModel.cpp:54
gnsstk::SaasTropModel::wet_mapping_function
virtual double wet_mapping_function(double elevation) const
Definition: SaasTropModel.cpp:273
gnsstk::TropModel::SaasDryDelay
double SaasDryDelay(const double pr, const double lat, const double ht) const
Definition: TropModel.hpp:262
gnsstk::SaasTropModel::validDOY
bool validDOY
flag for valid day of year
Definition: SaasTropModel.hpp:171
gnsstk::WxObservation
A Single Weather Observation.
Definition: WxObsMap.hpp:55
THROW_IF_INVALID_DETAILED
#define THROW_IF_INVALID_DETAILED()
Definition: SaasTropModel.cpp:43
gnsstk::SaasTropModel::validWeather
bool validWeather
flag for valid weather
Definition: SaasTropModel.hpp:168
gnsstk
For Sinex::InputHistory.
Definition: BasicFramework.cpp:50
gnsstk::SaasDryC1
static const double SaasDryC1[5]
Definition: SaasTropModel.cpp:73
std::sin
double sin(gnsstk::Angle x)
Definition: Angle.hpp:144
gnsstk::SaasWetC
static const double SaasWetC[5]
Definition: SaasTropModel.cpp:58
gnsstk::SaasDryC
static const double SaasDryC[5]
Definition: SaasTropModel.cpp:66
gnsstk::WxObservation::pressure
float pressure
millibars
Definition: WxObsMap.hpp:83
gnsstk::TropModel::CELSIUS_TO_KELVIN
static const GNSSTK_EXPORT double CELSIUS_TO_KELVIN
Definition: TropModel.hpp:108
gnsstk::Exception
Definition: Exception.hpp:151
gnsstk::TropModel::humid
double humid
latest value of relative humidity (percent)
Definition: TropModel.hpp:282
gnsstk::SaasTropModel::correction
virtual double correction(double elevation) const
Definition: SaasTropModel.cpp:120
gnsstk::SaasDryB1
static const double SaasDryB1[5]
Definition: SaasTropModel.cpp:71
gnsstk::Position::getHeight
double getHeight() const noexcept
return height above ellipsoid (meters)
Definition: Position.hpp:474
gnsstk::SaasDryB
static const double SaasDryB[5]
Definition: SaasTropModel.cpp:64
gnsstk::SaasTropModel::validRxHeight
bool validRxHeight
flag for valid Rx longitude
Definition: SaasTropModel.hpp:170
gnsstk::SaasTropModel::latitude
double latitude
latitude (deg) of receiver
Definition: SaasTropModel.hpp:166
gnsstk::Position::getGeodeticLatitude
double getGeodeticLatitude() const noexcept
return geodetic latitude (deg N)
Definition: Position.hpp:454
gnsstk::CommonTime
Definition: CommonTime.hpp:84
gnsstk::SaasTropModel::setReceiverHeight
void setReceiverHeight(const double &ht)
Definition: SaasTropModel.cpp:336
gnsstk::Xvt
Definition: Xvt.hpp:60
gnsstk::SaasTropModel::setWeather
virtual void setWeather(const WxObservation &wx)
Definition: SaasTropModel.cpp:322
gnsstk::TropModel::valid
bool valid
true only if current model parameters are valid
Definition: TropModel.hpp:279
std::cos
double cos(gnsstk::Angle x)
Definition: Angle.hpp:146
gnsstk::TropModel::press
double press
latest value of pressure (millibars)
Definition: TropModel.hpp:281
GNSSTK_RETHROW
#define GNSSTK_RETHROW(exc)
Definition: Exception.hpp:369
gnsstk::TrackingCode::P
@ P
Legacy GPS precise code.
gnsstk::SaasTropModel::dry_zenith_delay
virtual double dry_zenith_delay() const
Definition: SaasTropModel.cpp:189
gnsstk::Position
Definition: Position.hpp:136
GNSSTK_THROW
#define GNSSTK_THROW(exc)
Definition: Exception.hpp:366
gnsstk::SaasTropModel::setDayOfYear
void setDayOfYear(const int &d)
Definition: SaasTropModel.cpp:352
gnsstk::WxObservation::humidity
float humidity
percent
Definition: WxObsMap.hpp:84
gnsstk::SaasTropModel::dry_mapping_function
virtual double dry_mapping_function(double elevation) const
Definition: SaasTropModel.cpp:224
gnsstk::SaasDryA1
static const double SaasDryA1[5]
Definition: SaasTropModel.cpp:69
gnsstk::TropModel::temp
double temp
latest value of temperature (kelvin or celsius)
Definition: TropModel.hpp:280
gnsstk::SaasTropModel::setReceiverLatitude
void setReceiverLatitude(const double &lat)
Definition: SaasTropModel.cpp:344
gnsstk::Position::elevation
double elevation(const Position &Target) const
Definition: Position.cpp:1308
gnsstk::SaasTropModel::wet_zenith_delay
virtual double wet_zenith_delay() const
Definition: SaasTropModel.cpp:199
gnsstk::DEG_TO_RAD
static const double DEG_TO_RAD
Conversion Factor from degrees to radians (units: degrees^-1)
Definition: GNSSconstants.hpp:76

gnsstk
Author(s):
autogenerated on Wed Oct 25 2023 02:40:41