TropModel.cpp

Go to the documentation of this file.
00001 #pragma ident "$Id: TropModel.cpp 3140 2012-06-18 15:03:02Z susancummins $"
00002 
00003 //============================================================================
00004 //
00005 //  This file is part of GPSTk, the GPS Toolkit.
00006 //
00007 //  The GPSTk is free software; you can redistribute it and/or modify
00008 //  it under the terms of the GNU Lesser General Public License as published
00009 //  by the Free Software Foundation; either version 2.1 of the License, or
00010 //  any later version.
00011 //
00012 //  The GPSTk is distributed in the hope that it will be useful,
00013 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 //  GNU Lesser General Public License for more details.
00016 //
00017 //  You should have received a copy of the GNU Lesser General Public
00018 //  License along with GPSTk; if not, write to the Free Software Foundation,
00019 //  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA
00020 //
00021 //  Copyright 2004, The University of Texas at Austin
00022 //
00023 //============================================================================
00024 
00025 //============================================================================
00026 //
00027 //This software developed by Applied Research Laboratories at the University of
00028 //Texas at Austin, under contract to an agency or agencies within the U.S.
00029 //Department of Defense. The U.S. Government retains all rights to use,
00030 //duplicate, distribute, disclose, or release this software.
00031 //
00032 //Pursuant to DoD Directive 523024
00033 //
00034 // DISTRIBUTION STATEMENT A: This software has been approved for public
00035 //                           release, distribution is unlimited.
00036 //
00037 //=============================================================================
00038 
00045 #include "TropModel.hpp"
00046 #include "EphemerisRange.hpp"             // for Elevation()
00047 #include "MathBase.hpp"                   // SQRT
00048 #include "geometry.hpp"                   // DEG_TO_RAD
00049 #include "GPSEllipsoid.hpp"               // ell.a() = R earth
00050 #include "GNSSconstants.hpp"          // TWO_PI
00051 #include "YDSTime.hpp"
00052 
00053 namespace gpstk
00054 {
00055       // for temperature conversion from Celcius to Kelvin
00056    static const double CELSIUS_TO_KELVIN = 273.15;
00057 
00058       // Compute and return the full tropospheric delay. Typically call
00059       // setWeather(T,P,H) before making this call.
00060       // @param elevation Elevation of satellite as seen at receiver, in degrees
00061    double TropModel::correction(double elevation) const
00062       throw(TropModel::InvalidTropModel)
00063    {
00064       if(!valid)
00065          GPSTK_THROW(InvalidTropModel("Invalid model"));
00066 
00067       if(elevation < 0.0)
00068          return 0.0;
00069 
00070       return (dry_zenith_delay() * dry_mapping_function(elevation)
00071             + wet_zenith_delay() * wet_mapping_function(elevation));
00072 
00073    }  // end TropModel::correction(elevation)
00074 
00075       // Compute and return the full tropospheric delay, given the positions of
00076       // receiver and satellite and the time tag. This version is most useful
00077       // within positioning algorithms, where the receiver position and timetag may
00078       // vary; it computes the elevation (and other receiver location information)
00079       // and passes them to appropriate set...() routines and the
00080       // correction(elevation) routine.
00081       // @param RX  Receiver position in ECEF cartesian coordinates (meters)
00082       // @param SV  Satellite position in ECEF cartesian coordinates (meters)
00083       // @param tt  Time tag of the signal
00084    double TropModel::correction(const Position& RX,
00085                                 const Position& SV,
00086                                 const CommonTime& tt)
00087       throw(TropModel::InvalidTropModel)
00088    {
00089       if(!valid)
00090          GPSTK_THROW(InvalidTropModel("Invalid model"));
00091 
00092       double c;
00093       try
00094       {
00095          c = correction(RX.elevation(SV));
00096       }
00097       catch(InvalidTropModel& e)
00098       {
00099          GPSTK_RETHROW(e);
00100       }
00101       return c;
00102    }  // end TropModel::correction(RX,SV,TT)
00103 
00104       // Re-define the tropospheric model with explicit weather data.
00105       // Typically called just before correction().
00106       // @param T temperature in degrees Celsius
00107       // @param P atmospheric pressure in millibars
00108       // @param H relative humidity in percent
00109    void TropModel::setWeather(const double& T,
00110                               const double& P,
00111                               const double& H)
00112       throw(InvalidParameter)
00113    {
00114       temp = T + CELSIUS_TO_KELVIN;
00115       press = P;
00116       humid = H;
00117       if (temp < 0.0)
00118       {
00119          valid = false;
00120          GPSTK_THROW(InvalidParameter("Invalid temperature parameter."));
00121       }
00122       if (press < 0.0)
00123       {
00124          valid = false;
00125          GPSTK_THROW(InvalidParameter("Invalid pressure parameter."));
00126       }
00127       if (humid < 0.0 || humid > 100.0)
00128       {
00129          valid = false;
00130          GPSTK_THROW(InvalidParameter("Invalid humidity parameter."));
00131       }
00132    }  // end TropModel::setWeather(T,P,H)
00133 
00134       // Re-define the tropospheric model with explicit weather data.
00135       // Typically called just before correction().
00136       // @param wx the weather to use for this correction
00137    void TropModel::setWeather(const WxObservation& wx)
00138       throw(InvalidParameter)
00139    {
00140       if (wx.isAllValid())
00141       {
00142          try
00143          {
00144             setWeather(wx.temperature, wx.pressure, wx.humidity);
00145             valid = true;
00146          }
00147          catch(InvalidParameter& e)
00148          {
00149             valid = false;
00150             GPSTK_RETHROW(e);
00151          }
00152       }
00153       else
00154       {
00155          valid = false;
00156          GPSTK_THROW(InvalidParameter("Invalid weather data"));
00157       }
00158    }
00165    void TropModel::weatherByStandardAtmosphereModel(const double& ht, double& T, double& P, double& H)
00166    {
00167 
00168          // reference height and it's relate weather(T P H)
00169       const double h0  = 0.0;                      // meter
00170       const double Tr  = +18.0;       // Celsius
00171       const double pr  = 1013.25;               // millibarc
00172       const double Hr  = 50;                       // humidity
00173 
00174       T = Tr-0.0065*(ht-h0);
00175       P = pr * std::pow((1 - 0.0000226 * (ht - h0)), 5.225);
00176       H = Hr * std::exp(-0.0006396 * (ht - h0));
00177 
00178    }
00179 
00180    // -----------------------------------------------------------------------
00181    // Simple Black model. This has been used as the 'default' for many years.
00182 
00183       // Default constructor
00184    SimpleTropModel::SimpleTropModel(void)
00185    {
00186       setWeather(20.0, 980.0, 50.0);
00187       Cwetdelay = 0.122382715318184;
00188       Cdrydelay = 2.235486646978727;
00189       Cwetmap = 1.000282213715744;
00190       Cdrymap = 1.001012704615527;
00191       valid = true;
00192    }
00193 
00194       // Creates a trop model from a weather observation
00195       // @param wx the weather to use for this correction.
00196    SimpleTropModel::SimpleTropModel(const WxObservation& wx)
00197       throw(InvalidParameter)
00198    {
00199       setWeather(wx);
00200       valid = true;
00201    }
00202 
00203       // Create a tropospheric model from explicit weather data
00204       // @param T temperature in degrees Celsius
00205       // @param P atmospheric pressure in millibars
00206       // @param H relative humidity in percent
00207    SimpleTropModel::SimpleTropModel(const double& T,
00208                                     const double& P,
00209                                     const double& H)
00210       throw(InvalidParameter)
00211    {
00212       setWeather(T,P,H);
00213       valid = true;
00214    }
00215 
00216       // Re-define the tropospheric model with explicit weather data.
00217       // Typically called just before correction().
00218       // @param T temperature in degrees Celsius
00219       // @param P atmospheric pressure in millibars
00220       // @param H relative humidity in percent
00221    void SimpleTropModel::setWeather(const double& T,
00222                                     const double& P,
00223                                     const double& H)
00224       throw(InvalidParameter)
00225    {
00226       TropModel::setWeather(T,P,H);
00227       GPSEllipsoid ell;
00228       Cdrydelay = 2.343*(press/1013.25)*(temp-3.96)/temp;
00229       double tks = temp * temp;
00230       Cwetdelay = 8.952/tks*humid*std::exp(-37.2465+0.213166*temp-(0.256908e-3)*tks);
00231       Cdrymap =1.0+(0.15)*148.98*(temp-3.96)/ell.a();
00232       Cwetmap =1.0+(0.15)*12000.0/ell.a();
00233       valid = true;
00234    }  // end SimpleTropModel::setWeather(T,P,H)
00235 
00236       // Re-define the tropospheric model with explicit weather data.
00237       // Typically called just before correction().
00238       // @param wx the weather to use for this correction
00239    void SimpleTropModel::setWeather(const WxObservation& wx)
00240       throw(InvalidParameter)
00241    {
00242       TropModel::setWeather(wx);
00243    }
00244 
00245       // Compute and return the zenith delay for dry component of the troposphere
00246    double SimpleTropModel::dry_zenith_delay(void) const
00247       throw(TropModel::InvalidTropModel)
00248    {
00249       if(!valid)
00250          GPSTK_THROW(InvalidTropModel("Invalid model"));
00251 
00252       return Cdrydelay;
00253    }
00254 
00255       // Compute and return the zenith delay for wet component of the troposphere
00256    double SimpleTropModel::wet_zenith_delay(void) const
00257       throw(TropModel::InvalidTropModel)
00258    {
00259       if(!valid)
00260          GPSTK_THROW(InvalidTropModel("Invalid model"));
00261 
00262       return Cwetdelay;
00263    }
00264 
00265       // Compute and return the mapping function for dry component
00266       // of the troposphere
00267       // @param elevation is the Elevation of satellite as seen at receiver,
00268       //                  in degrees
00269    double SimpleTropModel::dry_mapping_function(double elevation) const
00270       throw(TropModel::InvalidTropModel)
00271    {
00272       if(!valid)
00273          GPSTK_THROW(InvalidTropModel("Invalid model"));
00274 
00275       if(elevation < 0.0) return 0.0;
00276 
00277       double d = std::cos(elevation*DEG_TO_RAD);
00278       d /= Cdrymap;
00279       return (1.0/SQRT(1.0-d*d));
00280    }
00281 
00282       // Compute and return the mapping function for wet component
00283       // of the troposphere
00284       // @param elevation is the Elevation of satellite as seen at receiver,
00285       //                  in degrees
00286    double SimpleTropModel::wet_mapping_function(double elevation) const
00287       throw(TropModel::InvalidTropModel)
00288    {
00289       if(!valid)
00290          GPSTK_THROW(InvalidTropModel("Invalid model"));
00291 
00292       if(elevation < 0.0) return 0.0;
00293 
00294       double d = std::cos(elevation*DEG_TO_RAD);
00295       d /= Cwetmap;
00296       return (1.0/SQRT(1.0-d*d));
00297    }
00298 
00299    // ------------------------------------------------------------------------
00300    // Tropospheric model based on Goad and Goodman(1974),
00301    // "A Modified Hopfield Tropospheric Refraction Correction Model," Paper
00302    // presented at the Fall Annual Meeting of the American Geophysical Union,
00303    // San Francisco, December 1974.
00304    // See Leick, "GPS Satellite Surveying," Wiley, NY, 1990, Chapter 9,
00305    // particularly Table 9.1.
00306    // ------------------------------------------------------------------------
00307 
00308    static const double GGdryscale = 8594.777388436570600;
00309    static const double GGwetscale = 2540.042008403690900;
00310 
00311       // Default constructor
00312    GGTropModel::GGTropModel(void)
00313    {
00314       TropModel::setWeather(20.0, 980.0, 50.0);
00315       Cdrydelay = 2.59629761092150147e-4;    // zenith delay, dry
00316       Cwetdelay = 4.9982784999977412e-5;     // zenith delay, wet
00317       Cdrymap = 42973.886942182834900;       // height for mapping, dry
00318       Cwetmap = 12700.210042018454260;       // height for mapping, wet
00319       valid = true;
00320    }  // end GGTropModel::GGTropModel()
00321 
00322       // Creates a trop model from a weather observation
00323       // @param wx the weather to use for this correction.
00324    GGTropModel::GGTropModel(const WxObservation& wx)
00325       throw(InvalidParameter)
00326    {
00327       setWeather(wx);
00328       valid = true;
00329    }  // end GGTropModel::GGTropModel(weather)
00330 
00331       // Create a tropospheric model from explicit weather data
00332       // @param T temperature in degrees Celsius
00333       // @param P atmospheric pressure in millibars
00334       // @param H relative humidity in percent
00335    GGTropModel::GGTropModel(const double& T,
00336                             const double& P,
00337                             const double& H)
00338       throw(InvalidParameter)
00339    {
00340       setWeather(T,P,H);
00341       valid = true;
00342    } // end GGTropModel::GGTropModel()
00343 
00344    double GGTropModel::dry_zenith_delay(void) const
00345       throw(TropModel::InvalidTropModel)
00346    {
00347       if(!valid)
00348          GPSTK_THROW(InvalidTropModel("Invalid model"));
00349 
00350       return (Cdrydelay * GGdryscale);
00351    }  // end GGTropModel::dry_zenith_delay()
00352 
00353    double GGTropModel::wet_zenith_delay(void) const
00354       throw(TropModel::InvalidTropModel)
00355    {
00356       if(!valid)
00357          GPSTK_THROW(InvalidTropModel("Invalid model"));
00358 
00359       return (Cwetdelay * GGwetscale);
00360    }  // end GGTropModel::wet_zenith_delay()
00361 
00362    double GGTropModel::dry_mapping_function(double elevation) const
00363       throw(TropModel::InvalidTropModel)
00364    {
00365       if(!valid)
00366          GPSTK_THROW(InvalidTropModel("Invalid model"));
00367 
00368       if(elevation < 0.0) return 0.0;
00369 
00370       GPSEllipsoid ell;
00371       double ce=std::cos(elevation*DEG_TO_RAD), se=std::sin(elevation*DEG_TO_RAD);
00372       double ad = -se/Cdrymap;
00373       double bd = -ce*ce/(2.0*ell.a()*Cdrymap);
00374       double Rd = SQRT((ell.a()+Cdrymap)*(ell.a()+Cdrymap)
00375                 - ell.a()*ell.a()*ce*ce) - ell.a()*se;
00376 
00377       double Ad[9], ad2=ad*ad, bd2=bd*bd;
00378       Ad[0] = 1.0;
00379       Ad[1] = 4.0*ad;
00380       Ad[2] = 6.0*ad2 + 4.0*bd;
00381       Ad[3] = 4.0*ad*(ad2+3.0*bd);
00382       Ad[4] = ad2*ad2 + 12.0*ad2*bd + 6.0*bd2;
00383       Ad[5] = 4.0*ad*bd*(ad2+3.0*bd);
00384       Ad[6] = bd2*(6.0*ad2+4.0*bd);
00385       Ad[7] = 4.0*ad*bd*bd2;
00386       Ad[8] = bd2*bd2;
00387 
00388          // compute dry component of the mapping function
00389       double sumd=0.0;
00390       for(int j=9; j>=1; j--) {
00391          sumd += Ad[j-1]/double(j);
00392          sumd *= Rd;
00393       }
00394       return sumd/GGdryscale;
00395 
00396    }  // end GGTropModel::dry_mapping_function(elev)
00397 
00398       // compute wet component of the mapping function
00399    double GGTropModel::wet_mapping_function(double elevation) const
00400       throw(TropModel::InvalidTropModel)
00401    {
00402       if(!valid)
00403          GPSTK_THROW(InvalidTropModel("Invalid model"));
00404 
00405       if(elevation < 0.0) return 0.0;
00406 
00407       GPSEllipsoid ell;
00408       double ce = std::cos(elevation*DEG_TO_RAD), se = std::sin(elevation*DEG_TO_RAD);
00409       double aw = -se/Cwetmap;
00410       double bw = -ce*ce/(2.0*ell.a()*Cwetmap);
00411       double Rw = SQRT((ell.a()+Cwetmap)*(ell.a()+Cwetmap)
00412                 - ell.a()*ell.a()*ce*ce) - ell.a()*se;
00413 
00414       double Aw[9], aw2=aw*aw, bw2=bw*bw;
00415       Aw[0] = 1.0;
00416       Aw[1] = 4.0*aw;
00417       Aw[2] = 6.0*aw2 + 4.0*bw;
00418       Aw[3] = 4.0*aw*(aw2+3.0*bw);
00419       Aw[4] = aw2*aw2 + 12.0*aw2*bw + 6.0*bw2;
00420       Aw[5] = 4.0*aw*bw*(aw2+3.0*bw);
00421       Aw[6] = bw2*(6.0*aw2+4.0*bw);
00422       Aw[7] = 4.0*aw*bw*bw2;
00423       Aw[8] = bw2*bw2;
00424 
00425       double sumw=0.0;
00426       for(int j=9; j>=1; j--) {
00427          sumw += Aw[j-1]/double(j);
00428          sumw *= Rw;
00429       }
00430       return sumw/GGwetscale;
00431 
00432    }  // end GGTropModel::wet_mapping_function(elev)
00433 
00434    void GGTropModel::setWeather(const double& T,
00435                                 const double& P,
00436                                 const double& H)
00437       throw(InvalidParameter)
00438    {
00439       TropModel::setWeather(T,P,H);
00440       double th=300./temp;
00441          // water vapor partial pressure (mb)
00442          // this comes from Leick and is not good.
00443          // double wvpp=6.108*(RHum*0.01)*exp((17.15*Tk-4684.0)/(Tk-38.45));
00444       double wvpp=2.409e9*humid*th*th*th*th*std::exp(-22.64*th);
00445       Cdrydelay = 7.7624e-5*press/temp;
00446       Cwetdelay = 1.0e-6*(-12.92+3.719e+05/temp)*(wvpp/temp);
00447       Cdrymap = (5.0*0.002277*press)/Cdrydelay;
00448       Cwetmap = (5.0*0.002277/Cwetdelay)*(1255.0/temp+0.5)*wvpp;
00449       valid = true;
00450    }  // end GGTropModel::setWeather(T,P,H)
00451 
00452       // Re-define the tropospheric model with explicit weather data.
00453       // Typically called just before correction().
00454       // @param wx the weather to use for this correction
00455    void GGTropModel::setWeather(const WxObservation& wx)
00456       throw(InvalidParameter)
00457    {
00458       TropModel::setWeather(wx);
00459    }
00460 
00461    // ------------------------------------------------------------------------
00462    // Tropospheric model with heights based on Goad and Goodman(1974),
00463    // "A Modified Hopfield Tropospheric Refraction Correction Model," Paper
00464    // presented at the Fall Annual Meeting of the American Geophysical Union,
00465    // San Francisco, December 1974.
00466    // (Not the same as GGTropModel because this has height dependence,
00467    // and the computation of this model does not break cleanly into
00468    // wet and dry components.)
00469 
00470       // Default constructor
00471    GGHeightTropModel::GGHeightTropModel(void)
00472    {
00473       validWeather = false; //setWeather(20.0,980.0,50.0);
00474       validHeights = false; //setHeights(0.0,0.0,0.0);
00475       validRxHeight = false;
00476    }
00477 
00478       // Creates a trop model from a weather observation
00479       // @param wx the weather to use for this correction.
00480    GGHeightTropModel::GGHeightTropModel(const WxObservation& wx)
00481       throw(InvalidParameter)
00482    {
00483       valid = validRxHeight = validHeights = false;
00484       setWeather(wx);
00485    }
00486 
00487       // Create a tropospheric model from explicit weather data
00488       // @param T temperature in degrees Celsius
00489       // @param P atmospheric pressure in millibars
00490       // @param H relative humidity in percent
00491    GGHeightTropModel::GGHeightTropModel(const double& T,
00492                                         const double& P,
00493                                         const double& H)
00494       throw(InvalidParameter)
00495    {
00496       validRxHeight = validHeights = false;
00497       setWeather(T,P,H);
00498    }
00499 
00500       // Create a valid model from explicit input.
00501       // @param T temperature in degrees Celsius
00502       // @param P atmospheric pressure in millibars
00503       // @param H relative humidity in percent
00504       // @param hT height at which temperature applies in meters.
00505       // @param hP height at which atmospheric pressure applies in meters.
00506       // @param hH height at which relative humidity applies in meters.
00507    GGHeightTropModel::GGHeightTropModel(const double& T,
00508                                         const double& P,
00509                                         const double& H,
00510                                         const double hT,
00511                                         const double hP,
00512                                         const double hH)
00513       throw(InvalidParameter)
00514    {
00515       validRxHeight = false;
00516       setWeather(T,P,H);
00517       setHeights(hT,hP,hH);
00518    }
00519 
00520       // re-define this to get the throws
00521    double GGHeightTropModel::correction(double elevation) const
00522       throw(TropModel::InvalidTropModel)
00523    {
00524       if(!valid)
00525       {
00526          if(!validWeather)
00527             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
00528          if(!validHeights)
00529             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
00530          if(!validRxHeight)
00531             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
00532       }
00533       if(elevation < 0.0) return 0.0;
00534 
00535       return (dry_zenith_delay() * dry_mapping_function(elevation)
00536             + wet_zenith_delay() * wet_mapping_function(elevation));
00537 
00538    }  // end GGHeightTropModel::correction(elevation)
00539 
00540       // Compute and return the full tropospheric delay, given the positions of
00541       // receiver and satellite and the time tag. This version is most useful
00542       // within positioning algorithms, where the receiver position and timetag
00543       // may vary; it computes the elevation (and other receiver location
00544       // information) and passes them to appropriate set...() routines and
00545       // the correction(elevation) routine.
00546       // @param RX  Receiver position
00547       // @param SV  Satellite position
00548       // @param tt  Time tag of the signal
00549    double GGHeightTropModel::correction(const Position& RX,
00550                                         const Position& SV,
00551                                         const CommonTime& tt)
00552       throw(TropModel::InvalidTropModel)
00553    {
00554       if(!valid)
00555       {
00556          if(!validWeather)
00557             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
00558          if(!validHeights)
00559             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
00560          if(!validRxHeight)
00561             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
00562       }
00563 
00564       // compute height from RX
00565       setReceiverHeight(RX.getHeight());
00566 
00567       return TropModel::correction(RX.elevation(SV));
00568 
00569    }  // end GGHeightTropModel::correction(RX,SV,TT)
00570 
00571       // Compute and return the full tropospheric delay, given the positions of
00572       // receiver and satellite and the time tag. This version is most useful
00573       // within positioning algorithms, where the receiver position and timetag
00574       // may vary; it computes the elevation (and other receiver location
00575       // information) and passes them to appropriate set...() routines and
00576       // the correction(elevation) routine.
00577       // @param RX  Receiver position in ECEF cartesian coordinates (meters)
00578       // @param SV  Satellite position in ECEF cartesian coordinates (meters)
00579       // @param tt  Time tag of the signal
00580       // This function is deprecated; use the Position version
00581    double GGHeightTropModel::correction(const Xvt& RX,
00582                                         const Xvt& SV,
00583                                         const CommonTime& tt)
00584       throw(TropModel::InvalidTropModel)
00585    {
00586       Position R(RX),S(SV);
00587       return GGHeightTropModel::correction(R,S,tt);
00588    }
00589 
00590       // Compute and return the zenith delay for dry component of the troposphere
00591    double GGHeightTropModel::dry_zenith_delay(void) const
00592       throw(TropModel::InvalidTropModel)
00593    {
00594       if(!valid) {
00595          if(!validWeather)
00596             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
00597          if(!validHeights)
00598             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
00599          if(!validRxHeight)
00600             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
00601       }
00602       double hrate=6.5e-3;
00603       double Ts=temp+hrate*height;
00604       double em=978.77/(2.8704e4*hrate);
00605       double Tp=Ts-hrate*hpress;
00606       double ps=press*std::pow(Ts/Tp,em)/1000.0;
00607       double rs=77.624e-3/Ts;
00608       double ho=11.385/rs;
00609       rs *= ps;
00610       double zen=(ho-height)/ho;
00611       zen = rs*zen*zen*zen*zen;
00612          // normalize
00613       zen *= (ho-height)/5;
00614       return zen;
00615 
00616    }  // end GGHeightTropModel::dry_zenith_delay()
00617 
00618       // Compute and return the zenith delay for wet component of the troposphere
00619    double GGHeightTropModel::wet_zenith_delay(void) const
00620       throw(TropModel::InvalidTropModel)
00621    {
00622       if(!valid) {
00623          if(!validWeather)
00624             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
00625          if(!validHeights)
00626             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
00627          if(!validRxHeight)
00628             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
00629       }
00630       double hrate=6.5e-3; //   deg K / m
00631       double Th=temp-273.15-hrate*(hhumid-htemp);
00632       double Ta=7.5*Th/(237.3+Th);
00633          // water vapor partial pressure
00634       double e0=6.11e-5*humid*std::pow(10.0,Ta);
00635       double Ts=temp+hrate*htemp;
00636       double em=978.77/(2.8704e4*hrate);
00637       double Tk=Ts-hrate*hhumid;
00638       double es=e0*std::pow(Ts/Tk,4.0*em);
00639       double rs=(371900.0e-3/Ts-12.92e-3)/Ts;
00640       double ho=11.385*(1255/Ts+0.05)/rs;
00641       double zen=(ho-height)/ho;
00642       zen = rs*es*zen*zen*zen*zen;
00643          //normalize
00644       zen *= (ho-height)/5;
00645       return zen;
00646 
00647    }  // end GGHeightTropModel::wet_zenith_delay()
00648 
00649       // Compute and return the mapping function for dry component
00650       // of the troposphere
00651       // @param elevation Elevation of satellite as seen at receiver,
00652       //                  in degrees
00653    double GGHeightTropModel::dry_mapping_function(double elevation) const
00654       throw(TropModel::InvalidTropModel)
00655    {
00656       if(!valid) {
00657          if(!validWeather)
00658             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
00659          if(!validHeights)
00660             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
00661          if(!validRxHeight)
00662             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
00663       }
00664       if(elevation < 0.0) return 0.0;
00665 
00666       double hrate=6.5e-3;
00667       double Ts=temp+hrate*htemp;
00668       double ho=(11.385/77.624e-3)*Ts;
00669       double se=std::sin(elevation*DEG_TO_RAD);
00670       if(se < 0.0) se=0.0;
00671 
00672       GPSEllipsoid ell;
00673       double rt,a,b,rn[8],al[8],er=ell.a();
00674       rt = (er+ho)/(er+height);
00675       rt = rt*rt - (1.0-se*se);
00676       if(rt < 0) rt=0.0;
00677       rt = (er+height)*(SQRT(rt)-se);
00678       a = -se/(ho-height);
00679       b = -(1.0-se*se)/(2.0*er*(ho-height));
00680       rn[0] = rt*rt;
00681       for(int j=1; j<8; j++) rn[j]=rn[j-1]*rt;
00682       al[0] = 2*a;
00683       al[1] = 2*a*a+4*b/3;
00684       al[2] = a*(a*a+3*b);
00685       al[3] = a*a*a*a/5+2.4*a*a*b+1.2*b*b;
00686       al[4] = 2*a*b*(a*a+3*b)/3;
00687       al[5] = b*b*(6*a*a+4*b)*0.1428571;
00688       if(b*b > 1.0e-35) {
00689          al[6] = a*b*b*b/2;
00690          al[7] = b*b*b*b/9;
00691       } else {
00692          al[6] = 0.0;
00693          al[7] = 0.0;
00694       }
00695       double map=rt;
00696       for(int k=0; k<8; k++) map += al[k]*rn[k];
00697          // normalize
00698       double norm=(ho-height)/5;
00699       return map/norm;
00700 
00701    }  // end GGHeightTropModel::dry_mapping_function(elevation)
00702 
00703       // Compute and return the mapping function for wet component
00704       // of the troposphere
00705       // @param elevation Elevation of satellite as seen at receiver,
00706       //                  in degrees
00707    double GGHeightTropModel::wet_mapping_function(double elevation) const
00708       throw(TropModel::InvalidTropModel)
00709    {
00710       if(!valid) {
00711          if(!validWeather)
00712             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
00713          if(!validHeights)
00714             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
00715          if(!validRxHeight)
00716             GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
00717       }
00718       if(elevation < 0.0) return 0.0;
00719 
00720       double hrate=6.5e-3;
00721       double Ts=temp+hrate*htemp;
00722       double rs=(371900.0e-3/Ts-12.92e-3)/Ts;
00723       double ho=11.385*(1255/Ts+0.05)/rs;
00724       double se=std::sin(elevation*DEG_TO_RAD);
00725       if(se < 0.0) se=0.0;
00726 
00727       GPSEllipsoid ell;
00728       double rt,a,b,rn[8],al[8],er=ell.a();
00729       rt = (er+ho)/(er+height);
00730       rt = rt*rt - (1.0-se*se);
00731       if(rt < 0) rt=0.0;
00732       rt = (er+height)*(SQRT(rt)-se);
00733       a = -se/(ho-height);
00734       b = -(1.0-se*se)/(2.0*er*(ho-height));
00735       rn[0] = rt*rt;
00736       for(int i=1; i<8; i++) rn[i]=rn[i-1]*rt;
00737       al[0] = 2*a;
00738       al[1] = 2*a*a+4*b/3;
00739       al[2] = a*(a*a+3*b);
00740       al[3] = a*a*a*a/5+2.4*a*a*b+1.2*b*b;
00741       al[4] = 2*a*b*(a*a+3*b)/3;
00742       al[5] = b*b*(6*a*a+4*b)*0.1428571;
00743       if(b*b > 1.0e-35) {
00744          al[6] = a*b*b*b/2;
00745          al[7] = b*b*b*b/9;
00746       } else {
00747          al[6] = 0.0;
00748          al[7] = 0.0;
00749       }
00750       double map=rt;
00751       for(int j=0; j<8; j++) map += al[j]*rn[j];
00752          // normalize map function
00753       double norm=(ho-height)/5;
00754       return map/norm;
00755 
00756    }  // end GGHeightTropModel::wet_mapping_function(elevation)
00757 
00758       // Re-define the weather data.
00759       // Typically called just before correction().
00760       // @param T temperature in degrees Celsius
00761       // @param P atmospheric pressure in millibars
00762       // @param H relative humidity in percent
00763    void GGHeightTropModel::setWeather(const double& T,
00764                                       const double& P,
00765                                       const double& H)
00766       throw(InvalidParameter)
00767    {
00768       try
00769       {
00770          TropModel::setWeather(T,P,H);
00771          validWeather = true;
00772          valid = validWeather && validHeights && validRxHeight;
00773       }
00774       catch(InvalidParameter& e)
00775       {
00776          valid = validWeather = false;
00777          GPSTK_RETHROW(e);
00778       }
00779    }  // end GGHeightTropModel::setWeather(T,P,H)
00780 
00781       // Re-define the tropospheric model with explicit weather data.
00782       // Typically called just before correction().
00783       // @param wx the weather to use for this correction
00784    void GGHeightTropModel::setWeather(const WxObservation& wx)
00785       throw(InvalidParameter)
00786    {
00787       try
00788       {
00789          TropModel::setWeather(wx);
00790          validWeather = true;
00791          valid = validWeather && validHeights && validRxHeight;
00792       }
00793       catch(InvalidParameter& e)
00794       {
00795          valid = validWeather = false;
00796          GPSTK_RETHROW(e);
00797       }
00798    }
00799 
00800 
00801       // Re-define the heights at which the weather parameters apply.
00802       // Typically called just before correction().
00803       // @param hT height (m) at which temperature applies
00804       // @param hP height (m) at which atmospheric pressure applies
00805       // @param hH height (m) at which relative humidity applies
00806    void GGHeightTropModel::setHeights(const double& hT,
00807                                       const double& hP,
00808                                       const double& hH)
00809    {
00810       htemp = hT;                 // height (m) at which temp applies
00811       hpress = hP;                // height (m) at which press applies
00812       hhumid = hH;                // height (m) at which humid applies
00813       validHeights = true;
00814       valid = validWeather && validHeights && validRxHeight;
00815    }  // end GGHeightTropModel::setHeights(hT,hP,hH)
00816 
00817       // Define the receiver height; this required before calling
00818       // correction() or any of the zenith_delay or mapping_function routines.
00819    void GGHeightTropModel::setReceiverHeight(const double& ht)
00820    {
00821       height = ht;
00822       validRxHeight = true;
00823       if(!validHeights) {
00824          htemp = hpress = hhumid = ht;
00825          validHeights = true;
00826       }
00827       valid = validWeather && validHeights && validRxHeight;
00828    }  // end GGHeightTropModel::setReceiverHeight(const double& ht)
00829 
00830    // ------------------------------------------------------------------------
00831    // Tropospheric model developed by University of New Brunswick and described in
00832    // "A Tropospheric Delay Model for the User of the Wide Area Augmentation
00833    // System," J. Paul Collins and Richard B. Langley, Technical Report No. 187,
00834    // Dept. of Geodesy and Geomatics Engineering, University of New Brunswick,
00835    // 1997. See particularly Appendix C.
00836    //
00837    // This model includes a wet and dry component, and was designed for the user
00838    // without access to measurements of temperature, pressure and relative humidity
00839    // at ground level. It requires input of the latitude, day of year and height
00840    // above the ellipsoid, and it interpolates a table of values, using these
00841    // inputs, to get the ground level weather parameters plus other parameters and
00842    // the mapping function constants.
00843    //
00844    // NB in this model, units of temp are degrees Celsius, and humid is the water
00845    // vapor partial pressure.
00846 
00847    static const double NBRd=287.054;   // J/kg/K = m*m*K/s*s
00848    static const double NBg=9.80665;    // m/s*s
00849    static const double NBk1=77.604;    // K/mbar
00850    static const double NBk3p=382000;   // K*K/mbar
00851 
00852    static const double NBLat[5]={   15.0,   30.0,   45.0,   60.0,   75.0};
00853 
00854    // zenith delays, averages
00855    static const double NBZP0[5]={1013.25,1017.25,1015.75,1011.75,1013.00};
00856    static const double NBZT0[5]={ 299.65, 294.15, 283.15, 272.15, 263.65};
00857    static const double NBZW0[5]={  26.31,  21.79,  11.66,   6.78,   4.11};
00858    static const double NBZB0[5]={6.30e-3,6.05e-3,5.58e-3,5.39e-3,4.53e-3};
00859    static const double NBZL0[5]={   2.77,   3.15,   2.57,   1.81,   1.55};
00860 
00861    // zenith delays, amplitudes
00862    static const double NBZPa[5]={    0.0,  -3.75,  -2.25,  -1.75,  -0.50};
00863    static const double NBZTa[5]={    0.0,    7.0,   11.0,   15.0,   14.5};
00864    static const double NBZWa[5]={    0.0,   8.85,   7.24,   5.36,   3.39};
00865    static const double NBZBa[5]={    0.0,0.25e-3,0.32e-3,0.81e-3,0.62e-3};
00866    static const double NBZLa[5]={    0.0,   0.33,   0.46,   0.74,   0.30};
00867 
00868    // mapping function, dry, averages
00869    static const double NBMad[5]={1.2769934e-3,1.2683230e-3,1.2465397e-3,1.2196049e-3,
00870                                  1.2045996e-3};
00871    static const double NBMbd[5]={2.9153695e-3,2.9152299e-3,2.9288445e-3,2.9022565e-3,
00872                                  2.9024912e-3};
00873    static const double NBMcd[5]={62.610505e-3,62.837393e-3,63.721774e-3,63.824265e-3,
00874                                  64.258455e-3};
00875 
00876    // mapping function, dry, amplitudes
00877    static const double NBMaa[5]={0.0,1.2709626e-5,2.6523662e-5,3.4000452e-5,
00878                                  4.1202191e-5};
00879    static const double NBMba[5]={0.0,2.1414979e-5,3.0160779e-5,7.2562722e-5,
00880                                  11.723375e-5};
00881    static const double NBMca[5]={0.0,9.0128400e-5,4.3497037e-5,84.795348e-5,
00882                                  170.37206e-5};
00883 
00884    // mapping function, wet, averages (there are no amplitudes for wet)
00885    static const double NBMaw[5]={5.8021897e-4,5.6794847e-4,5.8118019e-4,5.9727542e-4,
00886                            6.1641693e-4};
00887    static const double NBMbw[5]={1.4275268e-3,1.5138625e-3,1.4572752e-3,1.5007428e-3,
00888                            1.7599082e-3};
00889    static const double NBMcw[5]={4.3472961e-2,4.6729510e-2,4.3908931e-2,4.4626982e-2,
00890                            5.4736038e-2};
00891 
00892    // labels for use with the interpolation routine
00893    enum TableEntry { ZP=1, ZT, ZW, ZB, ZL, Mad, Mbd, Mcd, Maw, Mbw, Mcw };
00894 
00895    // the interpolation routine
00896    static double NB_Interpolate(double lat, int doy, TableEntry entry)
00897    {
00898       const double *pave, *pamp;
00899       double ret, day=double(doy);
00900 
00901          // assign pointer to the right array
00902       switch(entry) {
00903          case ZP:  pave=&NBZP0[0]; pamp=&NBZPa[0]; break;
00904          case ZT:  pave=&NBZT0[0]; pamp=&NBZTa[0]; break;
00905          case ZW:  pave=&NBZW0[0]; pamp=&NBZWa[0]; break;
00906          case ZB:  pave=&NBZB0[0]; pamp=&NBZBa[0]; break;
00907          case ZL:  pave=&NBZL0[0]; pamp=&NBZLa[0]; break;
00908          case Mad: pave=&NBMad[0]; pamp=&NBMaa[0]; break;
00909          case Mbd: pave=&NBMbd[0]; pamp=&NBMba[0]; break;
00910          case Mcd: pave=&NBMcd[0]; pamp=&NBMca[0]; break;
00911          case Maw: pave=&NBMaw[0];                 break;
00912          case Mbw: pave=&NBMbw[0];                 break;
00913          case Mcw: pave=&NBMcw[0];                 break;
00914       }
00915 
00916          // find the index i such that NBLat[i] <= lat < NBLat[i+1]
00917       int i = int(ABS(lat)/15.0)-1;
00918 
00919       if(i>=0 && i<=3) {               // mid-latitude -> regular interpolation
00920          double m=(ABS(lat)-NBLat[i])/(NBLat[i+1]-NBLat[i]);
00921          ret = pave[i]+m*(pave[i+1]-pave[i]);
00922          if(entry < Maw)
00923             ret -= (pamp[i]+m*(pamp[i+1]-pamp[i]))
00924                * std::cos(TWO_PI*(day-28.0)/365.25);
00925       }
00926       else {                           // < 15 or > 75 -> simpler interpolation
00927          if(i<0) i=0; else i=4;
00928          ret = pave[i];
00929          if(entry < Maw)
00930             ret -= pamp[i]*std::cos(TWO_PI*(day-28.0)/365.25);
00931       }
00932 
00933       return ret;
00934 
00935    }  // end double NB_Interpolate(lat,doy,entry)
00936 
00937       // Default constructor
00938    NBTropModel::NBTropModel(void)
00939    {
00940       validWeather = false;
00941       validRxLatitude = false;
00942       validDOY = false;
00943       validRxHeight = false;
00944    } // end NBTropModel::NBTropModel()
00945 
00946       // Create a trop model using the minimum information: latitude and doy.
00947       // Interpolate the weather unless setWeather (optional) is called.
00948       // @param lat Latitude of the receiver in degrees.
00949       // @param day Day of year.
00950    NBTropModel::NBTropModel(const double& lat,
00951                             const int& day)
00952    {
00953       validRxHeight = false;
00954       setReceiverLatitude(lat);
00955       setDayOfYear(day);
00956       setWeather();
00957    }  // end NBTropModel::NBTropModel(lat, day);
00958 
00959       // Create a trop model with weather.
00960       // @param lat Latitude of the receiver in degrees.
00961       // @param day Day of year.
00962       // @param wx the weather to use for this correction.
00963    NBTropModel::NBTropModel(const double& lat,
00964                             const int& day,
00965                             const WxObservation& wx)
00966       throw(InvalidParameter)
00967    {
00968       validRxHeight = false;
00969       setReceiverLatitude(lat);
00970       setDayOfYear(day);
00971       setWeather(wx);
00972    }  // end NBTropModel::NBTropModel(weather)
00973 
00974       // Create a tropospheric model from explicit weather data
00975       // @param lat Latitude of the receiver in degrees.
00976       // @param day Day of year.
00977       // @param T temperature in degrees Celsius
00978       // @param P atmospheric pressure in millibars
00979       // @param H relative humidity in percent
00980    NBTropModel::NBTropModel(const double& lat,
00981                             const int& day,
00982                             const double& T,
00983                             const double& P,
00984                             const double& H)
00985       throw(InvalidParameter)
00986    {
00987       validRxHeight = false;
00988       setReceiverLatitude(lat);
00989       setDayOfYear(day);
00990       setWeather(T,P,H);
00991    } // end NBTropModel::NBTropModel()
00992 
00993       // Create a valid model from explicit input (weather will be estimated
00994       // internally by this model).
00995       // @param ht Height of the receiver in meters.
00996       // @param lat Latitude of the receiver in degrees.
00997       // @param d Day of year.
00998    NBTropModel::NBTropModel(const double& ht,
00999                             const double& lat,
01000                             const int& day)
01001    {
01002       setReceiverHeight(ht);
01003       setReceiverLatitude(lat);
01004       setDayOfYear(day);
01005       setWeather();
01006    }
01007 
01008       // re-define this to get the throws
01009    double NBTropModel::correction(double elevation) const
01010       throw(TropModel::InvalidTropModel)
01011    {
01012       if(!valid) {
01013          if(!validWeather)
01014             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: weather"));
01015          if(!validRxLatitude)
01016             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
01017          if(!validRxHeight)
01018             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
01019          if(!validDOY)
01020             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
01021       }
01022 
01023       if(elevation < 0.0) return 0.0;
01024 
01025       return (dry_zenith_delay() * dry_mapping_function(elevation)
01026             + wet_zenith_delay() * wet_mapping_function(elevation));
01027    }  // end NBTropModel::correction(elevation)
01028 
01029       // Compute and return the full tropospheric delay, given the positions of
01030       // receiver and satellite and the time tag. This version is most useful
01031       // within positioning algorithms, where the receiver position and timetag
01032       // may vary; it computes the elevation (and other receiver location
01033       // information) and passes them to appropriate set...() routines
01034       // and the correction(elevation) routine.
01035       // @param RX  Receiver position
01036       // @param SV  Satellite position
01037       // @param tt  Time tag of the signal
01038    double NBTropModel::correction(const Position& RX,
01039                                   const Position& SV,
01040                                   const CommonTime& tt)
01041       throw(TropModel::InvalidTropModel)
01042    {
01043       if(!valid) {
01044          if(!validWeather)
01045 
01046          if(!validRxLatitude)
01047             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
01048          if(!validRxHeight)
01049             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
01050          if(!validDOY)
01051             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
01052       }
01053 
01054          // compute height and latitude from RX
01055       setReceiverHeight(RX.getHeight());
01056       setReceiverLatitude(RX.getGeodeticLatitude());
01057 
01058          // compute day of year from tt
01059       setDayOfYear(int((static_cast<YDSTime>(tt)).doy));
01060 
01061       return TropModel::correction(RX.elevation(SV));
01062 
01063    }  // end NBTropModel::correction(RX,SV,TT)
01064 
01065       // Compute and return the full tropospheric delay, given the positions of
01066       // receiver and satellite and the time tag. This version is most useful
01067       // within positioning algorithms, where the receiver position and timetag
01068       // may vary; it computes the elevation (and other receiver location
01069       // information) and passes them to appropriate set...() routines
01070       // and the correction(elevation) routine.
01071       // @param RX  Receiver position in ECEF cartesian coordinates (meters)
01072       // @param SV  Satellite position in ECEF cartesian coordinates (meters)
01073       // @param tt  Time tag of the signal
01074       // This function is deprecated; use the Position version
01075    double NBTropModel::correction(const Xvt& RX,
01076                                   const Xvt& SV,
01077                                   const CommonTime& tt)
01078       throw(TropModel::InvalidTropModel)
01079    {
01080       Position R(RX),S(SV);
01081       return NBTropModel::correction(R,S,tt);
01082    }
01083 
01084       // Compute and return the zenith delay for dry component of the troposphere
01085    double NBTropModel::dry_zenith_delay(void) const
01086       throw(TropModel::InvalidTropModel)
01087    {
01088       if(!valid) {
01089          if(!validWeather)
01090             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: weather"));
01091          if(!validRxLatitude)
01092             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
01093          if(!validRxHeight)
01094             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
01095          if(!validDOY)
01096             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
01097       }
01098       double beta = NB_Interpolate(latitude,doy,ZB);
01099       double gm = 9.784*(1.0-2.66e-3*std::cos(2.0*latitude*DEG_TO_RAD)-2.8e-7*height);
01100 
01101          // scale factors for height above mean sea level
01102          // if weather is given, assume it's measured at ht -> kw=kd=1
01103       double kd=1, base=std::log(1.0-beta*height/temp);
01104       if(interpolateWeather)
01105          kd = std::exp(base*NBg/(NBRd*beta));
01106 
01107          // compute the zenith delay for dry component
01108       return ((1.0e-6*NBk1*NBRd/gm) * kd * press);
01109 
01110    }  // end NBTropModel::dry_zenith_delay()
01111 
01112       // Compute and return the zenith delay for wet component of the troposphere
01113    double NBTropModel::wet_zenith_delay(void) const
01114       throw(TropModel::InvalidTropModel)
01115    {
01116       if(!valid) {
01117          if(!validWeather)
01118             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: weather"));
01119          if(!validRxLatitude)
01120             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
01121          if(!validRxHeight)
01122             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
01123          if(!validDOY)
01124             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
01125       }
01126       double beta = NB_Interpolate(latitude,doy,ZB);
01127       double lam = NB_Interpolate(latitude,doy,ZL);
01128       double gm = 9.784*(1.0-2.66e-3*std::cos(2.0*latitude*DEG_TO_RAD)-2.8e-7*height);
01129 
01130          // scale factors for height above mean sea level
01131          // if weather is given, assume it's measured at ht -> kw=kd=1
01132       double kw=1, base=std::log(1.0-beta*height/temp);
01133       if(interpolateWeather)
01134          kw = std::exp(base*(-1.0+(lam+1)*NBg/(NBRd*beta)));
01135 
01136          // compute the zenith delay for wet component
01137       return ((1.0e-6*NBk3p*NBRd/(gm*(lam+1)-beta*NBRd)) * kw * humid/temp);
01138 
01139    }  // end NBTropModel::wet_zenith_delay()
01140 
01141       // Compute and return the mapping function for dry component
01142       // of the troposphere
01143       // @param elevation Elevation of satellite as seen at receiver,
01144       //                  in degrees
01145    double NBTropModel::dry_mapping_function(double elevation) const
01146       throw(TropModel::InvalidTropModel)
01147    {
01148       if(!valid) {
01149          if(!validWeather)
01150             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: weather"));
01151          if(!validRxLatitude)
01152             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
01153          if(!validRxHeight)
01154             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
01155          if(!validDOY)
01156             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
01157       }
01158       if(elevation < 0.0) return 0.0;
01159 
01160       double a,b,c,se,map;
01161       se = std::sin(elevation*DEG_TO_RAD);
01162 
01163       a = NB_Interpolate(latitude,doy,Mad);
01164       b = NB_Interpolate(latitude,doy,Mbd);
01165       c = NB_Interpolate(latitude,doy,Mcd);
01166       map = (1.0+a/(1.0+b/(1.0+c))) / (se+a/(se+b/(se+c)));
01167 
01168       a = 2.53e-5;
01169       b = 5.49e-3;
01170       c = 1.14e-3;
01171       if(ABS(elevation)<=0.001) se=0.001;
01172       map += ((1.0/se)-(1.0+a/(1.0+b/(1.0+c)))/(se+a/(se+b/(se+c))))*height/1000.0;
01173 
01174       return map;
01175 
01176    }  // end NBTropModel::dry_mapping_function()
01177 
01178       // Compute and return the mapping function for wet component
01179       // of the troposphere
01180       // @param elevation Elevation of satellite as seen at receiver,
01181       //                  in degrees
01182    double NBTropModel::wet_mapping_function(double elevation) const
01183       throw(TropModel::InvalidTropModel)
01184    {
01185       if(!valid) {
01186          if(!validWeather)
01187             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: weather"));
01188          if(!validRxLatitude)
01189             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
01190          if(!validRxHeight)
01191             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
01192          if(!validDOY)
01193             GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
01194       }
01195       if(elevation < 0.0) return 0.0;
01196 
01197       double a,b,c,se;
01198       se = std::sin(elevation*DEG_TO_RAD);
01199       a = NB_Interpolate(latitude,doy,Maw);
01200       b = NB_Interpolate(latitude,doy,Mbw);
01201       c = NB_Interpolate(latitude,doy,Mcw);
01202 
01203       return ( (1.0+a/(1.0+b/(1.0+c))) / (se+a/(se+b/(se+c))) );
01204 
01205    }  // end NBTropModel::wet_mapping_function()
01206 
01207       // Re-define the weather data.
01208       // If called, typically called before any calls to correction().
01209       // @param T temperature in degrees Celsius
01210       // @param P atmospheric pressure in millibars
01211       // @param H relative humidity in percent
01212    void NBTropModel::setWeather(const double& T,
01213                                 const double& P,
01214                                 const double& H)
01215       throw(InvalidParameter)
01216    {
01217       interpolateWeather=false;
01218       TropModel::setWeather(T,P,H);
01219             // humid actually stores water vapor partial pressure
01220       double th=300./temp;
01221       humid = 2.409e9*H*th*th*th*th*std::exp(-22.64*th);
01222       validWeather = true;
01223       valid = validWeather && validRxHeight && validRxLatitude && validDOY;
01224 
01225    }  // end NBTropModel::setWeather()
01226 
01227       // Re-define the tropospheric model with explicit weather data.
01228       // Typically called just before correction().
01229       // @param wx the weather to use for this correction
01230    void NBTropModel::setWeather(const WxObservation& wx)
01231       throw(InvalidParameter)
01232    {
01233       interpolateWeather = false;
01234       try
01235       {
01236          TropModel::setWeather(wx);
01237             // humid actually stores vapor partial pressure
01238          double th=300./temp;
01239          humid = 2.409e9*humid*th*th*th*th*std::exp(-22.64*th);
01240          validWeather = true;
01241          valid = validWeather && validRxHeight && validRxLatitude && validDOY;
01242       }
01243       catch(InvalidParameter& e)
01244       {
01245          valid = validWeather = false;
01246          GPSTK_RETHROW(e);
01247       }
01248    }
01249 
01250       // configure the model to estimate the weather from the internal model,
01251       // using lat and doy
01252    void NBTropModel::setWeather()
01253       throw(TropModel::InvalidTropModel)
01254    {
01255       interpolateWeather = true;
01256       if(!validRxLatitude)
01257       {
01258          valid = validWeather = false;
01259          GPSTK_THROW(InvalidTropModel(
01260             "NBTropModel must have Rx latitude before interpolating weather"));
01261       }
01262       if(!validDOY)
01263       {
01264          valid = validWeather = false;
01265          GPSTK_THROW(InvalidTropModel(
01266             "NBTropModel must have day of year before interpolating weather"));
01267       }
01268       temp = NB_Interpolate(latitude,doy,ZT);
01269       press = NB_Interpolate(latitude,doy,ZP);
01270       humid = NB_Interpolate(latitude,doy,ZW);
01271       validWeather = true;
01272       valid = validWeather && validRxHeight && validRxLatitude && validDOY;
01273    }
01274 
01275       // Define the receiver height; this required before calling
01276       // correction() or any of the zenith_delay or mapping_function routines.
01277    void NBTropModel::setReceiverHeight(const double& ht)
01278    {
01279       height = ht;
01280       validRxHeight = true;
01281       valid = validWeather && validRxHeight && validRxLatitude && validDOY;
01282       if(!validWeather && validRxLatitude && validDOY)
01283          setWeather();
01284    }  // end NBTropModel::setReceiverHeight()
01285 
01286       // Define the latitude of the receiver; this is required before calling
01287       // correction() or any of the zenith_delay or mapping_function routines.
01288    void NBTropModel::setReceiverLatitude(const double& lat)
01289    {
01290       latitude = lat;
01291       validRxLatitude = true;
01292       valid = validWeather && validRxHeight && validRxLatitude && validDOY;
01293       if(!validWeather && validRxLatitude && validDOY)
01294          setWeather();
01295    }  // end NBTropModel::setReceiverLatitude(lat)
01296 
01297       // Define the day of year; this is required before calling
01298       // correction() or any of the zenith_delay or mapping_function routines.
01299    void NBTropModel::setDayOfYear(const int& d)
01300    {
01301       doy = d;
01302       if(doy > 0 && doy < 367) validDOY=true; else validDOY = false;
01303       valid = validWeather && validRxHeight && validRxLatitude && validDOY;
01304       if(!validWeather && validRxLatitude && validDOY)
01305          setWeather();
01306    }  // end NBTropModel::setDayOfYear(doy)
01307 
01308    // ------------------------------------------------------------------------
01309    // Saastamoinen tropospheric model.
01310    // This model needs work; it is not the Saastamoinen model, but appears to be
01311    // a combination of the Neill mapping functions and an unknown delay model.
01312    // Based on Saastamoinen, J., 'Atmospheric
01313    // Correction for the Troposphere and Stratosphere in Radio Ranging of
01314    // Satellites,' Geophysical Monograph 15, American Geophysical Union, 1972,
01315    // and Ch. 9 of McCarthy, D. and Petit, G., IERS Conventions (2003), IERS
01316    // Technical Note 32, IERS, 2004. The mapping functions are from
01317    // Neill, A.E., 1996, 'Global Mapping Functions for the Atmosphere Delay of
01318    // Radio Wavelengths,' J. Geophys. Res., 101, pp. 3227-3246 (also see IERS TN 32).
01319    //
01320    // This model includes a wet and dry component, and requires input of the
01321    // geodetic latitude, day of year and height above the ellipsoid of the receiver.
01322    //
01323    // Usually, the caller will set the latitude and day of year at the same
01324    // time the weather is set
01325    //   SaasTropModel stm;
01326    //   stm.setReceiverLatitude(lat);
01327    //   stm.setDayOfYear(doy);
01328    //   stm.setWeather(T,P,H);
01329    // Then, when the correction (and/or delay and map) is computed, receiver height
01330    // should be set before the call to correction(elevation):
01331    //   stm.setReceiverHeight(height);
01332    //   trop_corr = stm.correction(elevation);
01333    //
01334    // NB in this model, units of 'temp' are degrees Celcius and
01335    // humid actually stores water vapor partial pressure in mbars
01336    //
01337 
01338    // constants for wet mapping function
01339    static const double SaasWetA[5]=
01340      { 0.00058021897, 0.00056794847, 0.00058118019, 0.00059727542, 0.00061641693 };
01341    static const double SaasWetB[5]=
01342      { 0.0014275268, 0.0015138625, 0.0014572752, 0.0015007428, 0.0017599082 };
01343    static const double SaasWetC[5]=
01344      { 0.043472961, 0.046729510, 0.043908931, 0.044626982, 0.054736038 };
01345 
01346    // constants for dry mapping function
01347    static const double SaasDryA[5]=
01348      { 0.0012769934, 0.0012683230, 0.0012465397, 0.0012196049, 0.0012045996 };
01349    static const double SaasDryB[5]=
01350      { 0.0029153695, 0.0029152299, 0.0029288445, 0.0029022565, 0.0029024912 };
01351    static const double SaasDryC[5]=
01352      { 0.062610505, 0.062837393, 0.063721774, 0.063824265, 0.064258455 };
01353 
01354    static const double SaasDryA1[5]=
01355      { 0.0, 0.000012709626, 0.000026523662, 0.000034000452, 0.000041202191 };
01356    static const double SaasDryB1[5]=
01357      { 0.0, 0.000021414979, 0.000030160779, 0.000072562722, 0.00011723375 };
01358    static const double SaasDryC1[5]=
01359      { 0.0, 0.000090128400, 0.000043497037, 0.00084795348, 0.0017037206 };
01360 
01361       // Default constructor
01362    SaasTropModel::SaasTropModel(void)
01363    {
01364       validWeather = false;
01365       validRxLatitude = false;
01366       validDOY = false;
01367       validRxHeight = false;
01368    } // end SaasTropModel::SaasTropModel()
01369 
01370       // Create a trop model using the minimum information: latitude and doy.
01371       // Interpolate the weather unless setWeather (optional) is called.
01372       // @param lat Latitude of the receiver in degrees.
01373       // @param day Day of year.
01374    SaasTropModel::SaasTropModel(const double& lat,
01375                                 const int& day)
01376    {
01377       validWeather = false;
01378       validRxHeight = false;
01379       SaasTropModel::setReceiverLatitude(lat);
01380       SaasTropModel::setDayOfYear(day);
01381    } // end SaasTropModel::SaasTropModel
01382 
01383       // Create a trop model with weather.
01384       // @param lat Latitude of the receiver in degrees.
01385       // @param day Day of year.
01386       // @param wx the weather to use for this correction.
01387    SaasTropModel::SaasTropModel(const double& lat,
01388                                 const int& day,
01389                                 const WxObservation& wx)
01390       throw(InvalidParameter)
01391    {
01392       validRxHeight = false;
01393       SaasTropModel::setReceiverLatitude(lat);
01394       SaasTropModel::setDayOfYear(day);
01395       SaasTropModel::setWeather(wx);
01396    }  // end SaasTropModel::SaasTropModel(weather)
01397 
01398       // Create a tropospheric model from explicit weather data
01399       // @param lat Latitude of the receiver in degrees.
01400       // @param day Day of year.
01401       // @param T temperature in degrees Celsius
01402       // @param P atmospheric pressure in millibars
01403       // @param H relative humidity in percent
01404    SaasTropModel::SaasTropModel(const double& lat,
01405                                 const int& day,
01406                                 const double& T,
01407                                 const double& P,
01408                                 const double& H)
01409       throw(InvalidParameter)
01410    {
01411       validRxHeight = false;
01412       SaasTropModel::setReceiverLatitude(lat);
01413       SaasTropModel::setDayOfYear(day);
01414       SaasTropModel::setWeather(T,P,H);
01415    } // end SaasTropModel::SaasTropModel()
01416 
01417       // re-define this to get the throws correct
01418    double SaasTropModel::correction(double elevation) const
01419       throw(TropModel::InvalidTropModel)
01420    {
01421       if(!valid) {
01422          if(!validWeather) GPSTK_THROW(
01423             InvalidTropModel("Invalid Saastamoinen trop model: weather"));
01424          if(!validRxLatitude) GPSTK_THROW(
01425             InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
01426          if(!validRxHeight) GPSTK_THROW(
01427             InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
01428          if(!validDOY) GPSTK_THROW(
01429             InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
01430          GPSTK_THROW(
01431             InvalidTropModel("Valid flag corrupted in Saastamoinen trop model"));
01432       }
01433 
01434       if(elevation < 0.0) return 0.0;
01435 
01436       double corr=0.0;
01437       try {
01438          corr = (dry_zenith_delay() * dry_mapping_function(elevation)
01439             + wet_zenith_delay() * wet_mapping_function(elevation));
01440       }
01441       catch(Exception& e) { GPSTK_RETHROW(e); }
01442 
01443       return corr;
01444 
01445    }  // end SaasTropModel::correction(elevation)
01446 
01447       // Compute and return the full tropospheric delay, given the positions of
01448       // receiver and satellite and the time tag. This version is most useful
01449       // within positioning algorithms, where the receiver position and timetag
01450       // may vary; it computes the elevation (and other receiver location
01451       // information) and passes them to appropriate set...() routines
01452       // and the correction(elevation) routine.
01453       // @param RX  Receiver position
01454       // @param SV  Satellite position
01455       // @param tt  Time tag of the signal
01456    double SaasTropModel::correction(const Position& RX,
01457                                     const Position& SV,
01458                                     const CommonTime& tt)
01459       throw(TropModel::InvalidTropModel)
01460    {
01461       SaasTropModel::setReceiverHeight(RX.getHeight());
01462       SaasTropModel::setReceiverLatitude(RX.getGeodeticLatitude());
01463       SaasTropModel::setDayOfYear(int((static_cast<YDSTime>(tt).doy)));
01464 
01465       if(!valid) {
01466          if(!validWeather) GPSTK_THROW(
01467             InvalidTropModel("Invalid Saastamoinen trop model: weather"));
01468          if(!validRxLatitude) GPSTK_THROW(
01469             InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
01470          if(!validRxHeight) GPSTK_THROW(
01471             InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
01472          if(!validDOY) GPSTK_THROW(
01473             InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
01474          valid = true;
01475       }
01476 
01477       double corr=0.0;
01478       try {
01479          corr = SaasTropModel::correction(RX.elevation(SV));
01480       }
01481       catch(Exception& e) { GPSTK_RETHROW(e); }
01482 
01483       return corr;
01484 
01485    }  // end SaasTropModel::correction(RX,SV,TT)
01486 
01487    double SaasTropModel::correction(const Xvt& RX,
01488                                     const Xvt& SV,
01489                                     const CommonTime& tt)
01490       throw(TropModel::InvalidTropModel)
01491    {
01492       Position R(RX),S(SV);
01493       return SaasTropModel::correction(R,S,tt);
01494    }
01495 
01496       // Compute and return the zenith delay for dry component of the troposphere
01497    double SaasTropModel::dry_zenith_delay(void) const
01498       throw(TropModel::InvalidTropModel)
01499    {
01500       if(!valid) {
01501          if(!validWeather) GPSTK_THROW(
01502             InvalidTropModel("Invalid Saastamoinen trop model: weather"));
01503          if(!validRxLatitude) GPSTK_THROW(
01504             InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
01505          if(!validRxHeight) GPSTK_THROW(
01506             InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
01507          if(!validDOY) GPSTK_THROW(
01508             InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
01509          GPSTK_THROW(
01510             InvalidTropModel("Valid flag corrupted in Saastamoinen trop model"));
01511       }
01512 
01513       // correct pressure for height
01514       //double press_at_h =
01515       //   press * std::pow((temp+273.16-4.5*height/1000.0)/(temp+273.16),34.1/4.5);
01516       // humid is zero for the dry component
01517       double delay = 0.0022768 * press //_at_h
01518             / (1 - 0.00266 * ::cos(2*latitude*DEG_TO_RAD) - 0.00028 * height/1000.);
01519 
01520       return delay;
01521 
01522    }  // end SaasTropModel::dry_zenith_delay()
01523 
01524       // Compute and return the zenith delay for wet component of the troposphere
01525    double SaasTropModel::wet_zenith_delay(void) const
01526       throw(TropModel::InvalidTropModel)
01527    {
01528       if(!valid) {
01529          if(!validWeather) GPSTK_THROW(
01530             InvalidTropModel("Invalid Saastamoinen trop model: weather"));
01531          if(!validRxLatitude) GPSTK_THROW(
01532             InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
01533          if(!validRxHeight) GPSTK_THROW(
01534             InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
01535          if(!validDOY) GPSTK_THROW(
01536             InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
01537          GPSTK_THROW(
01538             InvalidTropModel("Valid flag corrupted in Saastamoinen trop model"));
01539       }
01540 
01541       // press is zero for the wet component
01542       double delay = 0.0022768 * humid * 1255/((temp+CELSIUS_TO_KELVIN) + 0.05)
01543             / (1 - 0.00266 * ::cos(2*latitude*DEG_TO_RAD) - 0.00028 * height/1000.);
01544 
01545       return delay;
01546 
01547    }  // end SaasTropModel::wet_zenith_delay()
01548 
01549       // Compute and return the mapping function for dry component of the troposphere
01550       // @param elevation Elevation of satellite as seen at receiver, in degrees
01551    double SaasTropModel::dry_mapping_function(double elevation) const
01552       throw(TropModel::InvalidTropModel)
01553    {
01554       if(!valid) {
01555          if(!validWeather) GPSTK_THROW(
01556             InvalidTropModel("Invalid Saastamoinen trop model: weather"));
01557          if(!validRxLatitude) GPSTK_THROW(
01558             InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
01559          if(!validRxHeight) GPSTK_THROW(
01560             InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
01561          if(!validDOY) GPSTK_THROW(
01562             InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
01563          GPSTK_THROW(
01564             InvalidTropModel("Valid flag corrupted in Saastamoinen trop model"));
01565       }
01566       if(elevation < 0.0) return 0.0;
01567 
01568       double lat,t,ct;
01569       lat = fabs(latitude);         // degrees
01570       t = doy - 28.;                // mid-winter
01571       if(latitude < 0)              // southern hemisphere
01572          t += 365.25/2.;
01573       t *= 360.0/365.25;            // convert to degrees
01574       ct = ::cos(t*DEG_TO_RAD);
01575 
01576       double a,b,c;
01577       if(lat < 15.) {
01578          a = SaasDryA[0];
01579          b = SaasDryB[0];
01580          c = SaasDryC[0];
01581       }
01582       else if(lat < 75.) {          // coefficients are for 15,30,45,60,75 deg
01583          int i=int(lat/15.0)-1;
01584          double frac=(lat-15.*(i+1))/15.;
01585          a = SaasDryA[i] + frac*(SaasDryA[i+1]-SaasDryA[i]);
01586          b = SaasDryB[i] + frac*(SaasDryB[i+1]-SaasDryB[i]);
01587          c = SaasDryC[i] + frac*(SaasDryC[i+1]-SaasDryC[i]);
01588 
01589          a -= ct * (SaasDryA1[i] + frac*(SaasDryA1[i+1]-SaasDryA1[i]));
01590          b -= ct * (SaasDryB1[i] + frac*(SaasDryB1[i+1]-SaasDryB1[i]));
01591          c -= ct * (SaasDryC1[i] + frac*(SaasDryC1[i+1]-SaasDryC1[i]));
01592       }
01593       else {
01594          a = SaasDryA[4] - ct * SaasDryA1[4];
01595          b = SaasDryB[4] - ct * SaasDryB1[4];
01596          c = SaasDryC[4] - ct * SaasDryC1[4];
01597       }
01598 
01599       double se = ::sin(elevation*DEG_TO_RAD);
01600       double map = (1.+a/(1.+b/(1.+c)))/(se+a/(se+b/(se+c)));
01601 
01602       a = 0.0000253;
01603       b = 0.00549;
01604       c = 0.00114;
01605       map += (height/1000.0)*(1./se-(1+a/(1.+b/(1.+c)))/(se+a/(se+b/(se+c))));
01606 
01607       return map;
01608 
01609    }  // end SaasTropModel::dry_mapping_function()
01610 
01611       // Compute and return the mapping function for wet component of the troposphere
01612       // @param elevation Elevation of satellite as seen at receiver, in degrees.
01613    double SaasTropModel::wet_mapping_function(double elevation) const
01614       throw(TropModel::InvalidTropModel)
01615    {
01616       if(!valid) {
01617          if(!validWeather) GPSTK_THROW(
01618             InvalidTropModel("Invalid Saastamoinen trop model: weather"));
01619          if(!validRxLatitude) GPSTK_THROW(
01620             InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
01621          if(!validRxHeight) GPSTK_THROW(
01622             InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
01623          if(!validDOY) GPSTK_THROW(
01624             InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
01625          GPSTK_THROW(
01626             InvalidTropModel("Valid flag corrupted in Saastamoinen trop model"));
01627       }
01628       if(elevation < 0.0) return 0.0;
01629 
01630       double a,b,c,lat;
01631       lat = fabs(latitude);         // degrees
01632       if(lat < 15.) {
01633          a = SaasWetA[0];
01634          b = SaasWetB[0];
01635          c = SaasWetC[0];
01636       }
01637       else if(lat < 75.) {          // coefficients are for 15,30,45,60,75 deg
01638          int i=int(lat/15.0)-1;
01639          double frac=(lat-15.*(i+1))/15.;
01640          a = SaasWetA[i] + frac*(SaasWetA[i+1]-SaasWetA[i]);
01641          b = SaasWetB[i] + frac*(SaasWetB[i+1]-SaasWetB[i]);
01642          c = SaasWetC[i] + frac*(SaasWetC[i+1]-SaasWetC[i]);
01643       }
01644       else {
01645          a = SaasWetA[4];
01646          b = SaasWetB[4];
01647          c = SaasWetC[4];
01648       }
01649 
01650       double se = ::sin(elevation*DEG_TO_RAD);
01651       double map = (1.+a/(1.+b/(1.+c)))/(se+a/(se+b/(se+c)));
01652 
01653       return map;
01654 
01655    }  // end SaasTropModel::wet_mapping_function()
01656 
01657       // Re-define the weather data.
01658       // If called, typically called before any calls to correction().
01659       // @param T temperature in degrees Celsius
01660       // @param P atmospheric pressure in millibars
01661       // @param H relative humidity in percent
01662    void SaasTropModel::setWeather(const double& T,
01663                                   const double& P,
01664                                   const double& H)
01665       throw(InvalidParameter)
01666    {
01667       temp = T;
01668       press = P;
01669          // humid actually stores water vapor partial pressure
01670       double exp=7.5*T/(T+237.3);
01671       humid = 6.11 * (H/100.) * std::pow(10.0,exp);
01672 
01673       validWeather = true;
01674       valid = (validWeather && validRxHeight && validRxLatitude && validDOY);
01675 
01676    }  // end SaasTropModel::setWeather()
01677 
01678       // Re-define the tropospheric model with explicit weather data.
01679       // Typically called just before correction().
01680       // @param wx the weather to use for this correction
01681    void SaasTropModel::setWeather(const WxObservation& wx)
01682       throw(InvalidParameter)
01683    {
01684       try
01685       {
01686          SaasTropModel::setWeather(wx.temperature,wx.pressure,wx.humidity);
01687       }
01688       catch(InvalidParameter& e)
01689       {
01690          valid = validWeather = false;
01691          GPSTK_RETHROW(e);
01692       }
01693    }
01694 
01695       // Define the receiver height; this required before calling
01696       // correction() or any of the zenith_delay or mapping_function routines.
01697    void SaasTropModel::setReceiverHeight(const double& ht)
01698    {
01699       height = ht;
01700       validRxHeight = true;
01701       valid = (validWeather && validRxHeight && validRxLatitude && validDOY);
01702    }  // end SaasTropModel::setReceiverHeight()
01703 
01704       // Define the latitude of the receiver; this is required before calling
01705       // correction() or any of the zenith_delay or mapping_function routines.
01706    void SaasTropModel::setReceiverLatitude(const double& lat)
01707    {
01708       latitude = lat;
01709       validRxLatitude = true;
01710       valid = (validWeather && validRxHeight && validRxLatitude && validDOY);
01711    }  // end SaasTropModel::setReceiverLatitude(lat)
01712 
01713       // Define the day of year; this is required before calling
01714       // correction() or any of the zenith_delay or mapping_function routines.
01715    void SaasTropModel::setDayOfYear(const int& d)
01716    {
01717       doy = d;
01718       if(doy > 0 && doy < 367) validDOY=true; else validDOY = false;
01719       valid = (validWeather && validRxHeight && validRxLatitude && validDOY);
01720    }  // end SaasTropModel::setDayOfYear(doy)
01721 
01722 
01723    //-----------------------------------------------------------------------
01724       /* Tropospheric model implemented in "GPS Code Analysis Tool" (GCAT)
01725        * software.
01726        *
01727        * This model is described in the book "GPS Data processing: code and
01728        * phase Algorithms, Techniques and Recipes" by Hernandez-Pajares, M.,
01729        * J.M. Juan-Zornoza and Sanz-Subirana, J. See Chapter 5.
01730        *
01731        * This book and associated software are freely available at:
01732        *
01733        * http://gage152.upc.es/~manuel/tdgps/tdgps.html
01734        *
01735        * This is a simple but efective model composed of the wet and dry
01736        * vertical tropospheric delays as defined in Gipsy/Oasis-II GPS
01737        * analysis software, and the mapping function as defined by Black and
01738        * Eisner (H. D. Black, A. Eisner. Correcting Satellite Doppler
01739        * Data for Tropospheric Effects. Journal of  Geophysical Research.
01740        * Vol 89. 1984.) and used in MOPS (RTCA/DO-229C) standards.
01741        *
01742        * Usually, the caller will set the receiver height using
01743        * setReceiverHeight() method, and then call the correction() method
01744        * with the satellite elevation as parameter.
01745        *
01746        * @code
01747        *   GCATTropModel gcatTM();
01748        *   ...
01749        *   gcatTM.setReceiverHeight(150.0);
01750        *   trop = gcatTM.correction(elevation);
01751        * @endcode
01752        *
01753        * Another posibility is to set the receiver height when calling
01754        * the constructor.
01755        *
01756        * @code
01757        *   GCATTropModel gcatTM(150.0);    // Receiver height is 150.0 meters
01758        *   ...
01759        *   trop = gcatTM.correction(elevation);
01760        * @endcode
01761        */
01762 
01763       // Constructor to create a GCAT trop model providing  the height of the
01764       // receiver above mean sea level (as defined by ellipsoid model).
01765       //
01766       // @param ht Height of the receiver above mean sea level, in meters.
01767    GCATTropModel::GCATTropModel(const double& ht)
01768    {
01769       setReceiverHeight(ht);
01770       valid = true;
01771    }
01772 
01773 
01774       /* Compute and return the full tropospheric delay. The receiver height
01775        * must has been provided before, whether using the appropriate
01776        * constructor or with the setReceiverHeight() method.
01777        *
01778        * @param elevation  Elevation of satellite as seen at receiver, in
01779        *                   degrees
01780        */
01781    double GCATTropModel::correction(double elevation) const
01782       throw(TropModel::InvalidTropModel)
01783    {
01784       if(!valid) throw InvalidTropModel("Invalid model");
01785 
01786       if(elevation < 5.0) return 0.0;
01787 
01788       return ( (dry_zenith_delay() + wet_zenith_delay()) *
01789                mapping_function(elevation));
01790    }  // end GCATTropModel::correction(elevation)
01791 
01792 
01793       /* Compute and return the full tropospheric delay, given the positions of
01794        * receiver and satellite. This version is most useful within positioning
01795        * algorithms, where the receiver position may vary; it computes the
01796        * elevation and the receiver height and passes them to appropriate
01797        * set...() routines and the correction(elevation) routine.
01798        *
01799        * @param RX  Receiver position in ECEF cartesian coordinates (meters)
01800        * @param SV  Satellite position in ECEF cartesian coordinates (meters)
01801        */
01802    double GCATTropModel::correction( const Position& RX,
01803                                      const Position& SV )
01804       throw(TropModel::InvalidTropModel)
01805    {
01806 
01807       try
01808       {
01809          setReceiverHeight( RX.getAltitude() );
01810       }
01811       catch(GeometryException& e)
01812       {
01813          valid = false;
01814       }
01815 
01816       if(!valid) throw InvalidTropModel("Invalid model");
01817 
01818       double c;
01819       try
01820       {
01821          c = correction(RX.elevationGeodetic(SV));
01822       }
01823       catch(InvalidTropModel& e)
01824       {
01825          GPSTK_RETHROW(e);
01826       }
01827 
01828       return c;
01829 
01830    }  // end GCATTropModel::correction(RX,SV,TT)
01831 
01832 
01833       /* Compute and return the full tropospheric delay, given the positions of
01834        * receiver and satellite and the time tag. This version is most useful
01835        * within positioning algorithms, where the receiver position and timetag
01836        * may vary; it computes the elevation (and other receiver location
01837        * information) and passes them to appropriate set...() routines and the
01838        * correction(elevation) routine.
01839        *
01840        * @param RX  Receiver position in ECEF cartesian coordinates (meters)
01841        * @param SV  Satellite position in ECEF cartesian coordinates (meters)
01842        * @param tt  Time. In this model, tt is a dummy parameter kept only for
01843        *            consistency
01844        *
01845        * This function is deprecated; use the Position version
01846        */
01847    double GCATTropModel::correction( const Xvt& RX,
01848                                      const Xvt& SV,
01849                                      const CommonTime& tt )
01850       throw(TropModel::InvalidTropModel)
01851    {
01852 
01853       Position R(RX),S(SV);
01854 
01855       return GCATTropModel::correction(R,S);
01856    }  // end GCATTropModel::correction(RX,SV,tt)
01857 
01858 
01859       /* Compute and return the zenith delay for the dry component of the
01860        * troposphere.
01861        */
01862    double GCATTropModel::dry_zenith_delay(void) const
01863       throw(TropModel::InvalidTropModel)
01864    {
01865       if(!valid) throw InvalidTropModel("Invalid model");
01866 
01867       double ddry(2.29951*std::exp((-0.000116 * gcatHeight) ));
01868 
01869       return ddry;
01870    }  // end GCATTropModel::dry_zenith_delay()
01871 
01872 
01873       /* Compute and return the mapping function of the troposphere
01874        * @param elevation  Elevation of satellite as seen at receiver,
01875        *                   in degrees
01876        */
01877    double GCATTropModel::mapping_function(double elevation) const
01878       throw(TropModel::InvalidTropModel)
01879    {
01880       if(!valid) throw InvalidTropModel("Invalid model");
01881 
01882       if(elevation < 5.0) return 0.0;
01883 
01884       double d = std::sin(elevation*DEG_TO_RAD);
01885       d = SQRT(0.002001+(d*d));
01886 
01887       return (1.001/d);
01888    }  // end GCATTropModel::mapping_function(elevation)
01889 
01890 
01891       /* Define the receiver height; this is required before calling
01892        * correction() or any of the zenith_delay or mapping_function routines.
01893        * @param ht Height of the receiver above mean sea level, in meters.
01894        */
01895    void GCATTropModel::setReceiverHeight(const double& ht)
01896    {
01897       gcatHeight = ht;
01898       valid = true;
01899    }
01900 
01901 
01902    //---------------------------------------------------------------
01903       /* Tropospheric model implemented in the RTCA "Minimum Operational
01904        * Performance Standards" (MOPS), version C.
01905        *
01906        * This model is described in the RTCA "Minimum Operational Performance
01907        * Standards" (MOPS), version C (RTCA/DO-229C), in Appendix A.4.2.4.
01908        * Although originally developed for SBAS systems (EGNOS, WAAS), it may
01909        * be suitable for other uses as well.
01910        *
01911        * This model needs the day of year, satellite elevation (degrees),
01912        * receiver height over mean sea level (meters) and receiver latitude in
01913        * order to start computing.
01914        *
01915        * On the other hand, the outputs are the tropospheric correction (in
01916        * meters) and the sigma-squared of tropospheric delay residual error
01917        * (meters^2).
01918        *
01919        * A typical way to use this model follows:
01920        *
01921        * @code
01922        *   MOPSTropModel mopsTM;
01923        *   mopsTM.setReceiverLatitude(lat);
01924        *   mopsTM.setReceiverHeight(height);
01925        *   mopsTM.setDayOfYear(doy);
01926        * @endcode
01927        *
01928        * Once all the basic model parameters are set (latitude, height and day
01929        * of year), then we are able to compute the tropospheric correction as
01930        * a function of elevation:
01931        *
01932        * @code
01933        *   trop = mopsTM.correction(elevation);
01934        * @endcode
01935        *
01936        */
01937 
01938       // Some specific constants
01939    static const double MOPSg=9.80665;
01940    static const double MOPSgm=9.784;
01941    static const double MOPSk1=77.604;
01942    static const double MOPSk2=382000.0;
01943    static const double MOPSRd=287.054;
01944 
01945 
01946       // Empty constructor
01947    MOPSTropModel::MOPSTropModel(void)
01948    {
01949       validHeight = false;
01950       validLat    = false;
01951       validTime   = false;
01952       valid       = false;
01953    }
01954 
01955 
01956       /* Constructor to create a MOPS trop model providing the height of
01957        *  the receiver above mean sea level (as defined by ellipsoid model),
01958        *  its latitude and the day of year.
01959        *
01960        * @param ht   Height of the receiver above mean sea level, in meters.
01961        * @param lat  Latitude of receiver, in degrees.
01962        * @param doy  Day of year.
01963        */
01964    MOPSTropModel::MOPSTropModel( const double& ht,
01965                                  const double& lat,
01966                                  const int& doy )
01967    {
01968       setReceiverHeight(ht);
01969       setReceiverLatitude(lat);
01970       setDayOfYear(doy);
01971    }
01972 
01973 
01974       /* Constructor to create a MOPS trop model providing the position of
01975        *  the receiver and current time.
01976        *
01977        * @param RX   Receiver position.
01978        * @param time Time.
01979        */
01980    MOPSTropModel::MOPSTropModel(const Position& RX, const CommonTime& time)
01981    {
01982       setReceiverHeight(RX.getAltitude());
01983       setReceiverLatitude(RX.getGeodeticLatitude());
01984       setDayOfYear(time);
01985    }
01986 
01987 
01988       // Compute and return the full tropospheric delay. The receiver height,
01989       // latitude and Day oy Year must has been set before using the
01990       // appropriate constructor or the provided methods.
01991       // @param elevation Elevation of satellite as seen at receiver, in
01992       // degrees
01993    double MOPSTropModel::correction(double elevation) const
01994       throw(TropModel::InvalidTropModel)
01995    {
01996       if(!valid)
01997       {
01998          if(!validLat)
01999             throw InvalidTropModel("Invalid MOPS trop model: Rx Latitude");
02000          if(!validHeight)
02001             throw InvalidTropModel("Invalid MOPS trop model: Rx Height");
02002          if(!validTime)
02003             throw InvalidTropModel("Invalid MOPS trop model: day of year");
02004       }
02005 
02006       if(elevation < 5.0) return 0.0;
02007 
02008       double map = MOPSTropModel::mapping_function(elevation);
02009 
02010          // Compute tropospheric delay
02011       double tropDelay = ( MOPSTropModel::dry_zenith_delay() +
02012                            MOPSTropModel::wet_zenith_delay() ) * map;
02013 
02014       return tropDelay;
02015 
02016    }  // end MOPSTropModel::correction(elevation)
02017 
02018 
02019       // Compute and return the full tropospheric delay, given the positions of
02020       // receiver and satellite. This version is most useful within
02021       // positioning algorithms, where the receiver position may vary; it
02022       // computes the elevation (and other receiver location information as
02023       // height and latitude) and passes them to appropriate methods. You must
02024       // set time using method setDayOfYear() before calling this method.
02025       // @param RX  Receiver position
02026       // @param SV  Satellite position
02027    double MOPSTropModel::correction( const Position& RX,
02028                                      const Position& SV )
02029       throw(TropModel::InvalidTropModel)
02030    {
02031       try
02032       {
02033          setReceiverHeight( RX.getAltitude() );
02034          setReceiverLatitude(RX.getGeodeticLatitude());
02035          setWeather();
02036       }
02037       catch(GeometryException& e)
02038       {
02039          valid = false;
02040       }
02041 
02042       if(!valid) throw InvalidTropModel("Invalid model");
02043 
02044       double c;
02045       try
02046       {
02047          c = MOPSTropModel::correction(RX.elevationGeodetic(SV));
02048       }
02049       catch(InvalidTropModel& e)
02050       {
02051          GPSTK_RETHROW(e);
02052       }
02053 
02054       return c;
02055 
02056    }  // end MOPSTropModel::correction(RX,SV)
02057 
02058 
02059       // Compute and return the full tropospheric delay, given the positions of
02060       // receiver and satellite and the time tag. This version is most useful
02061       // within positioning algorithms, where the receiver position may vary;
02062       // it computes the elevation (and other receiver location information as
02063       // height  and latitude) and passes them to appropriate methods.
02064       // @param RX  Receiver position in ECEF cartesian coordinates (meters)
02065       // @param SV  Satellite position in ECEF cartesian coordinates (meters)
02066       // @param tt  Time (CommonTime object).
02067    double MOPSTropModel::correction( const Position& RX,
02068                                      const Position& SV,
02069                                      const CommonTime& tt )
02070       throw(TropModel::InvalidTropModel)
02071    {
02072       setDayOfYear(tt);
02073 
02074       return MOPSTropModel::correction(RX,SV);
02075    }  // end MOPSTropModel::correction(RX,SV,TT)
02076 
02077 
02078       // Compute and return the full tropospheric delay, given the positions of
02079       // receiver and satellite and the day of the year. This version is most
02080       // useful within positioning algorithms, where the receiver position may
02081       // vary; it computes the elevation (and other receiver location
02082       // information as height and latitude) and passes them to appropriate
02083       // methods.
02084       // @param RX  Receiver position in ECEF cartesian coordinates (meters)
02085       // @param SV  Satellite position in ECEF cartesian coordinates (meters)
02086       // @param doy Day of year.
02087    double MOPSTropModel::correction( const Position& RX,
02088                                      const Position& SV,
02089                                      const int& doy )
02090       throw(TropModel::InvalidTropModel)
02091    {
02092       setDayOfYear(doy);
02093 
02094       return MOPSTropModel::correction(RX,SV);
02095    }  // end MOPSTropModel::correction(RX,SV,doy)
02096 
02097 
02098 
02099       // deprecated
02100       // Compute and return the full tropospheric delay, given the positions of
02101       // receiver and satellite. You must set time using method setDayOfYear()
02102       // before calling this method.
02103       // @param RX  Receiver position in ECEF cartesian coordinates (meters)
02104       // @param SV  Satellite position in ECEF cartesian coordinates (meters)
02105       // This function is deprecated; use the Position version
02106    double MOPSTropModel::correction( const Xvt& RX,
02107                                      const Xvt& SV )
02108       throw(TropModel::InvalidTropModel)
02109    {
02110       Position R(RX),S(SV);
02111 
02112       return MOPSTropModel::correction(R,S);
02113    }  // end MOPSTropModel::correction(RX,SV)
02114 
02115 
02116       // deprecated
02117       // Compute and return the full tropospheric delay, given the positions of
02118       // receiver and satellite and the time tag. This version is most useful
02119       // within positioning algorithms, where the receiver position may vary;
02120       // it computes the elevation (and other receiver location information as
02121       // height and latitude) and passes them to appropriate methods.
02122       // @param RX  Receiver position in ECEF cartesian coordinates (meters)
02123       // @param SV  Satellite position in ECEF cartesian coordinates (meters)
02124       // @param tt  Time (CommonTime object).
02125       // This function is deprecated; use the Position version
02126    double MOPSTropModel::correction( const Xvt& RX,
02127                                      const Xvt& SV,
02128                                      const CommonTime& tt )
02129       throw(TropModel::InvalidTropModel)
02130    {
02131       setDayOfYear(tt);
02132       Position R(RX),S(SV);
02133 
02134       return MOPSTropModel::correction(R,S);
02135    }  // end MOPSTropModel::correction(RX,SV,tt)
02136 
02137 
02138       // deprecated
02139       // Compute and return the full tropospheric delay, given the positions of
02140       // receiver and satellite and the day of the year. This version is most
02141       // useful within positioning algorithms, where the receiver position may
02142       // vary; it computes the elevation (and other receiver location
02143       // information as height and latitude) and passes them to appropriate
02144       // methods.
02145       // @param RX  Receiver position in ECEF cartesian coordinates (meters)
02146       // @param SV  Satellite position in ECEF cartesian coordinates (meters)
02147       // @param doy Day of year.
02148       // This function is deprecated; use the Position version
02149    double MOPSTropModel::correction( const Xvt& RX,
02150                                      const Xvt& SV,
02151                                      const int& doy )
02152       throw(TropModel::InvalidTropModel)
02153    {
02154       setDayOfYear(doy);
02155       Position R(RX),S(SV);
02156 
02157       return MOPSTropModel::correction(R,S);
02158    }  // end MOPSTropModel::correction(RX,SV,doy)
02159 
02160 
02161       // Compute and return the zenith delay for the dry component of the
02162       // troposphere
02163    double MOPSTropModel::dry_zenith_delay(void) const
02164       throw(TropModel::InvalidTropModel)
02165    {
02166 
02167       if(!valid) throw InvalidTropModel("Invalid model");
02168 
02169       double ddry, zh_dry, exponent;
02170 
02171          // Set the extra parameters
02172       double P = MOPSParameters(0);
02173       double T = MOPSParameters(1);
02174       double beta = MOPSParameters(3);
02175 
02176          // Zero-altitude dry zenith delay:
02177       zh_dry = 0.000001*(MOPSk1*MOPSRd)*P/MOPSgm;
02178 
02179          // Zenith delay terms at MOPSHeight meters of height above mean sea
02180          // level
02181       exponent = MOPSg/MOPSRd/beta;
02182       ddry = zh_dry * std::pow( (1.0 - beta*MOPSHeight/T), exponent );
02183 
02184       return ddry;
02185 
02186    }  // end MOPSTropModel::dry_zenith_delay()
02187 
02188 
02189       // Compute and return the zenith delay for the wet component of the
02190       // troposphere
02191    double MOPSTropModel::wet_zenith_delay(void) const
02192       throw(TropModel::InvalidTropModel)
02193    {
02194 
02195       if(!valid) throw InvalidTropModel("Invalid model");
02196 
02197          double dwet, zh_wet, exponent;
02198 
02199          // Set the extra parameters
02200       double T = MOPSParameters(1);
02201       double e = MOPSParameters(2);
02202       double beta = MOPSParameters(3);
02203       double lambda = MOPSParameters(4);
02204 
02205          // Zero-altitude wet zenith delay:
02206       zh_wet = (0.000001*MOPSk2)*MOPSRd/(MOPSgm*(lambda+1.0)-beta*MOPSRd)*e/T;
02207 
02208          // Zenith delay terms at MOPSHeight meters of height above mean sea
02209          // level
02210       exponent = ( (lambda+1.0)*MOPSg/MOPSRd/beta)-1.0;
02211       dwet= zh_wet * std::pow( (1.0 - beta*MOPSHeight/T), exponent );
02212 
02213       return dwet;
02214 
02215    }  // end MOPSTropModel::wet_zenith_delay()
02216 
02217 
02218       // This method configure the model to estimate the weather using height,
02219       // latitude and day of year (DOY). It is called automatically when
02220       // setting those parameters.
02221    void MOPSTropModel::setWeather()
02222       throw(TropModel::InvalidTropModel)
02223    {
02224 
02225       if(!validLat)
02226       {
02227          valid = false;
02228          throw InvalidTropModel(
02229             "MOPSTropModel must have Rx latitude before computing weather");
02230       }
02231 
02232       if(!validTime)
02233       {
02234          valid = false;
02235          throw InvalidTropModel(
02236             "MOPSTropModel must have day of year before computing weather");
02237       }
02238 
02239          // In order to compute tropospheric delay we need to compute some
02240          // extra parameters
02241       try
02242       {
02243          prepareParameters();
02244       }
02245       catch(InvalidTropModel& e)
02246       {
02247          GPSTK_RETHROW(e);
02248       }
02249 
02250       valid = validHeight && validLat && validTime;
02251    }
02252 
02253 
02254       /* Define the receiver height; this is required before calling
02255        *  correction() or any of the zenith_delay routines.
02256        *
02257        * @param ht   Height of the receiver above mean sea level, in meters.
02258        */
02259    void MOPSTropModel::setReceiverHeight(const double& ht)
02260    {
02261       MOPSHeight = ht;
02262       validHeight = true;
02263 
02264          // Change the value of field "valid" if everything is already set
02265       valid = validHeight && validLat && validTime;
02266 
02267          // If model is valid, set the appropriate parameters
02268       if (valid) setWeather();
02269 
02270    }
02271 
02272 
02273       /* Define the receiver latitude; this is required before calling
02274        *  correction() or any of the zenith_delay routines.
02275        *
02276        * @param lat  Latitude of receiver, in degrees.
02277        */
02278    void MOPSTropModel::setReceiverLatitude(const double& lat)
02279    {
02280       MOPSLat = lat;
02281       validLat = true;
02282 
02283          // Change the value of field "valid" if everything is already set
02284       valid = validHeight && validLat && validTime;
02285 
02286          // If model is valid, set the appropriate parameters
02287       if (valid) setWeather();
02288 
02289    }
02290 
02291 
02292       /* Set the time when tropospheric correction will be computed for, in
02293        *  days of the year.
02294        *
02295        * @param doy  Day of the year.
02296        */
02297    void MOPSTropModel::setDayOfYear(const int& doy)
02298    {
02299 
02300       if ( (doy>=1) && (doy<=366))
02301       {
02302          validTime = true;
02303       }
02304       else
02305       {
02306          validTime = false;
02307       }
02308 
02309       MOPSTime = doy;
02310 
02311          // Change the value of field "valid" if everything is already set
02312       valid = validHeight && validLat && validTime;
02313 
02314          // If model is valid, set the appropriate parameters
02315       if (valid) setWeather();
02316    }
02317 
02318 
02319       /* Set the time when tropospheric correction will be computed for, in
02320        *  days of the year.
02321        *
02322        * @param time  Time object.
02323        */
02324    void MOPSTropModel::setDayOfYear(const CommonTime& time)
02325    {
02326       MOPSTime = (int)(static_cast<YDSTime>(time)).doy;
02327       validTime = true;
02328 
02329          // Change the value of field "valid" if everything is already set
02330       valid = validHeight && validLat && validTime;
02331 
02332          // If model is valid, set the appropriate parameters
02333       if (valid) setWeather();
02334 
02335    }
02336 
02337 
02338       /* Convenient method to set all model parameters in one pass.
02339        *
02340        * @param time  Time object.
02341        * @param rxPos Receiver position object.
02342        */
02343    void MOPSTropModel::setAllParameters( const CommonTime& time,
02344                                          const Position& rxPos )
02345    {
02346 
02347       MOPSTime = (int)(static_cast<YDSTime>(time)).doy;
02348       validTime = true;
02349       MOPSLat = rxPos.getGeodeticLatitude();
02350       validHeight = true;
02351       MOPSLat = rxPos.getHeight();
02352       validLat = true;
02353 
02354          // Change the value of field "valid" if everything is already set
02355       valid = validHeight && validLat && validTime;
02356 
02357          // If model is valid, set the appropriate parameters
02358       if (valid) setWeather();
02359 
02360    }
02361 
02362 
02363       // Compute and return the sigma-squared value of tropospheric delay
02364       // residual error (meters^2)
02365       // @param elevation  Elevation of satellite as seen at receiver,
02366       //                   in degrees
02367    double MOPSTropModel::MOPSsigma2(double elevation)
02368       throw(TropModel::InvalidTropModel)
02369    {
02370 
02371       double map_f;
02372 
02373          // If elevation is below bounds, fail in a sensible way returning a
02374          // very big sigma value
02375       if(elevation < 5.0)
02376       {
02377          return 9.9e9;
02378       }
02379       else
02380       {
02381          map_f = MOPSTropModel::mapping_function(elevation);
02382       }
02383 
02384          // Compute residual error for tropospheric delay
02385       double MOPSsigma2trop = (0.12*map_f)*(0.12*map_f);
02386 
02387       return MOPSsigma2trop;
02388 
02389    }  // end MOPSTropModel::MOPSsigma(elevation)
02390 
02391 
02392       // The MOPS tropospheric model needs to compute several extra parameters
02393    void MOPSTropModel::prepareParameters(void)
02394       throw(TropModel::InvalidTropModel)
02395    {
02396 
02397       if(!valid) throw InvalidTropModel("Invalid model");
02398 
02399       try
02400       {
02401             // We need to read some data
02402          prepareTables();
02403 
02404             // Declare some variables
02405          int idmin, j, index;
02406          double fact, axfi;
02407          Vector<double> avr0(5);
02408          Vector<double> svr0(5);
02409 
02410             // Resize MOPSParameters as appropriate
02411          MOPSParameters.resize(5);
02412 
02413          if (MOPSLat >= 0.0)
02414          {
02415             idmin = 28;
02416          }
02417          else
02418          {
02419             idmin = 211;
02420          }
02421 
02422             // Fraction of the year in radians
02423          fact = 2.0*PI*((double)(MOPSTime-idmin))/365.25;
02424 
02425          axfi = ABS(MOPSLat);
02426 
02427          if ( axfi <= 15.0 )                    index=0;
02428          if ( (axfi > 15.0) && (axfi <= 30.0) ) index=1;
02429          if ( (axfi > 30.0) && (axfi <= 45.0) ) index=2;
02430          if ( (axfi > 45.0) && (axfi <= 60.0) ) index=3;
02431          if ( (axfi > 60.0) && (axfi <  75.0) ) index=4;
02432          if ( axfi >= 75.0 )                     index=5;
02433 
02434          for (j=0; j<5; j++)
02435          {
02436             if (index == 0) {
02437                avr0(j)=avr(index,j);
02438                svr0(j)=svr(index,j);
02439             }
02440             else
02441             {
02442                if (index < 5)
02443                {
02444                   avr0(j) = avr(index-1,j) + (avr(index,j)-avr(index-1,j)) *
02445                             (axfi-fi0(index-1))/(fi0( index)-fi0(index-1));
02446 
02447                   svr0(j) = svr(index-1,j) + (svr(index,j)-svr(index-1,j)) *
02448                             (axfi-fi0(index-1))/(fi0( index)-fi0(index-1));
02449                }
02450                else
02451                {
02452                   avr0(j) = avr(index-1,j);
02453                   svr0(j) = svr(index-1,j);
02454                }
02455             }
02456 
02457             MOPSParameters(j) = avr0(j)-svr0(j)*std::cos(fact);
02458          }
02459 
02460       } // end try
02461       catch (...)
02462       {
02463          InvalidTropModel e("Problem computing extra MOPS parameters.");
02464          GPSTK_RETHROW(e);
02465       }
02466 
02467    }  // end MOPSTropModel::prepareParameters()
02468 
02469 
02470       // The MOPS tropospheric model uses several predefined data tables
02471    void MOPSTropModel::prepareTables(void)
02472    {
02473       avr.resize(5,5);
02474       svr.resize(5,5);
02475       fi0.resize(5);
02476 
02477 
02478          // Table avr (Average):
02479 
02480       avr(0,0) = 1013.25; avr(0,1) = 299.65; avr(0,2) = 26.31;
02481          avr(0,3) = 0.0063; avr(0,4) = 2.77;
02482 
02483       avr(1,0) = 1017.25; avr(1,1) = 294.15; avr(1,2) = 21.79;
02484          avr(1,3) = 0.00605; avr(1,4) = 3.15;
02485 
02486       avr(2,0) = 1015.75; avr(2,1) = 283.15; avr(2,2) = 11.66;
02487          avr(2,3) = 0.00558; avr(2,4) = 2.57;
02488 
02489       avr(3,0) = 1011.75; avr(3,1) = 272.15; avr(3,2) = 6.78;
02490          avr(3,3) = 0.00539; avr(3,4) = 1.81;
02491 
02492       avr(4,0) = 1013.00; avr(4,1) = 263.65; avr(4,2) = 4.11;
02493          avr(4,3) = 0.00453; avr(4,4) = 1.55;
02494 
02495 
02496          // Table svr (Seasonal Variation):
02497 
02498       svr(0,0) = 0.00; svr(0,1) = 0.00; svr(0,2) = 0.00;
02499          svr(0,3) = 0.00000; svr(0,4) = 0.00;
02500 
02501       svr(1,0) = -3.75; svr(1,1) = 7.00; svr(1,2) = 8.85;
02502          svr(1,3) = 0.00025; svr(1,4) = 0.33;
02503 
02504       svr(2,0) = -2.25; svr(2,1) = 11.00; svr(2,2) = 7.24;
02505          svr(2,3) = 0.00032; svr(2,4) = 0.46;
02506 
02507       svr(3,0) = -1.75; svr(3,1) = 15.00; svr(3,2) = 5.36;
02508          svr(3,3) = 0.00081; svr(3,4) = 0.74;
02509 
02510       svr(4,0) = -0.50; svr(4,1) = 14.50; svr(4,2) = 3.39;
02511          svr(4,3) = 0.00062; svr(4,4) = 0.30;
02512 
02513 
02514          // Table fi0 (Latitude bands):
02515 
02516       fi0(0) = 15.0; fi0(1) = 30.0; fi0(2) = 45.0;
02517          fi0(3) = 60.0; fi0(4) = 75.0;
02518 
02519    }
02520 
02521 
02522    //---------------------------------------------------------------------
02523       /* Tropospheric model based in the Neill mapping functions.
02524        *
02525        * This model uses the mapping functions developed by A.E. Niell and
02526        * published in Neill, A.E., 1996, 'Global Mapping Functions for the
02527        * Atmosphere Delay of Radio Wavelengths,' J. Geophys. Res., 101,
02528        * pp. 3227-3246 (also see IERS TN 32).
02529        *
02530        * The coefficients of the hydrostatic mapping function depend on the
02531        * latitude and height above sea level of the receiver station, and on
02532        * the day of the year. On the other hand, the wet mapping function
02533        * depends only on latitude.
02534        *
02535        * This mapping is independent from surface meteorology, while having
02536        * comparable accuracy and precision to those that require such data.
02537        * This characteristic makes this model very useful, and it is
02538        * implemented in geodetic software such as JPL's Gipsy/OASIS.
02539        *
02540        * A typical way to use this model follows:
02541        *
02542        * @code
02543        *   NeillTropModel neillTM;
02544        *   neillTM.setReceiverLatitude(lat);
02545        *   neillTM.setReceiverHeight(height);
02546        *   neillTM.setDayOfYear(doy);
02547        * @endcode
02548        *
02549        * Once all the basic model parameters are set (latitude, height and
02550        * day of year), then we are able to compute the tropospheric correction
02551        * as a function of elevation:
02552        *
02553        * @code
02554        *   trop = neillTM.correction(elevation);
02555        * @endcode
02556        *
02557        * @warning The Neill mapping functions are defined for elevation
02558        * angles down to 3 degrees.
02559        *
02560        */
02561 
02562       // Constructor to create a Neill trop model providing the position
02563       // of the receiver and current time.
02564       //
02565       // @param RX   Receiver position.
02566       // @param time Time.
02567    NeillTropModel::NeillTropModel( const Position& RX,
02568                                    const CommonTime& time )
02569    {
02570       setReceiverHeight(RX.getAltitude());
02571       setReceiverLatitude(RX.getGeodeticLatitude( ));
02572       setDayOfYear(time);
02573    }
02574 
02575 
02576 
02577       // Parameters borrowed from Saastamoinen tropospheric model
02578       // Constants for wet mapping function
02579    static const double NeillWetA[5] =
02580                               { 0.00058021897, 0.00056794847, 0.00058118019,
02581                                 0.00059727542, 0.00061641693 };
02582    static const double NeillWetB[5] =
02583                               { 0.0014275268, 0.0015138625, 0.0014572752,
02584                                 0.0015007428, 0.0017599082 };
02585    static const double NeillWetC[5] =
02586                               { 0.043472961, 0.046729510, 0.043908931,
02587                                 0.044626982, 0.054736038 };
02588 
02589       // constants for dry mapping function
02590    static const double NeillDryA[5] =
02591                               { 0.0012769934, 0.0012683230, 0.0012465397,
02592                                 0.0012196049, 0.0012045996 };
02593    static const double NeillDryB[5] =
02594                               { 0.0029153695, 0.0029152299, 0.0029288445,
02595                                 0.0029022565, 0.0029024912 };
02596    static const double NeillDryC[5] =
02597                               { 0.062610505, 0.062837393, 0.063721774,
02598                                 0.063824265, 0.064258455 };
02599 
02600    static const double NeillDryA1[5] =
02601                               { 0.0, 0.000012709626, 0.000026523662,
02602                                 0.000034000452, 0.000041202191 };
02603    static const double NeillDryB1[5] =
02604                               { 0.0, 0.000021414979, 0.000030160779,
02605                                 0.000072562722, 0.00011723375 };
02606    static const double NeillDryC1[5] =
02607                               { 0.0, 0.000090128400, 0.000043497037,
02608                                 0.00084795348, 0.0017037206 };
02609 
02610 
02611       // Compute and return the full tropospheric delay. The receiver height,
02612       // latitude and Day oy Year must has been set before using the
02613       // appropriate constructor or the provided methods.
02614       //
02615       // @param elevation Elevation of satellite as seen at receiver,
02616       // in degrees
02617    double NeillTropModel::correction(double elevation) const
02618       throw(TropModel::InvalidTropModel)
02619    {
02620 
02621       if(!valid)
02622       {
02623          if(!validLat)
02624          {
02625             throw InvalidTropModel("Invalid Neill trop model: Rx Latitude");
02626          }
02627 
02628          if(!validHeight)
02629          {
02630             throw InvalidTropModel("Invalid Neill trop model: Rx Height");
02631          }
02632 
02633          if(!validDOY)
02634          {
02635             throw InvalidTropModel("Invalid Neill trop model: day of year");
02636          }
02637       }
02638 
02639          // Neill mapping functions work down to 3 degrees of elevation
02640       if(elevation < 3.0)
02641       {
02642          return 0.0;
02643       }
02644 
02645       double map_dry(NeillTropModel::dry_mapping_function(elevation));
02646 
02647       double map_wet(NeillTropModel::wet_mapping_function(elevation));
02648 
02649          // Compute tropospheric delay
02650       double tropDelay( (NeillTropModel::dry_zenith_delay() * map_dry) +
02651                         (NeillTropModel::wet_zenith_delay() * map_wet) );
02652 
02653       return tropDelay;
02654 
02655    }  // end NeillTropModel::correction(elevation)
02656 
02657 
02658       /* Compute and return the full tropospheric delay, given the
02659        * positions of receiver and satellite.
02660        *
02661        * This version is more useful within positioning algorithms, where
02662        * the receiver position may vary; it computes the elevation (and
02663        * other receiver location information as height and latitude) and
02664        * passes them to appropriate methods.
02665        *
02666        * You must set time using method setReceiverDOY() before calling
02667        * this method.
02668        *
02669        * @param RX  Receiver position.
02670        * @param SV  Satellite position.
02671        */
02672    double NeillTropModel::correction( const Position& RX,
02673                                       const Position& SV )
02674       throw(TropModel::InvalidTropModel)
02675    {
02676 
02677       try
02678       {
02679          setReceiverHeight( RX.getAltitude() );
02680          setReceiverLatitude(RX.getGeodeticLatitude());
02681          setWeather();
02682       }
02683       catch(GeometryException& e)
02684       {
02685          valid = false;
02686       }
02687 
02688       if(!valid)
02689       {
02690          throw InvalidTropModel("Invalid model");
02691       }
02692 
02693       double c;
02694       try
02695       {
02696          c = NeillTropModel::correction(RX.elevationGeodetic(SV));
02697       }
02698       catch(InvalidTropModel& e)
02699       {
02700          GPSTK_RETHROW(e);
02701       }
02702 
02703       return c;
02704 
02705    }  // end NeillTropModel::correction(RX,SV)
02706 
02707 
02708       /* Compute and return the full tropospheric delay, given the
02709        * positions of receiver and satellite and the time tag.
02710        *
02711        * This version is more useful within positioning algorithms, where
02712        * the receiver position may vary; it computes the elevation (and
02713        * other receiver location information as height and latitude), and
02714        * passes them to appropriate methods.
02715        *
02716        * @param RX  Receiver position.
02717        * @param SV  Satellite position.
02718        * @param tt  Time (CommonTime object).
02719        */
02720    double NeillTropModel::correction( const Position& RX,
02721                                       const Position& SV,
02722                                       const CommonTime& tt )
02723       throw(TropModel::InvalidTropModel)
02724    {
02725 
02726       setDayOfYear(tt);
02727 
02728       return NeillTropModel::correction(RX,SV);
02729 
02730    }  // end NeillTropModel::correction(RX,SV,TT)
02731 
02732 
02733       /* Compute and return the full tropospheric delay, given the
02734        * positions of receiver and satellite and the day of the year.
02735        *
02736        * This version is more useful within positioning algorithms, where
02737        * the receiver position may vary; it computes the elevation (and
02738        * other receiver location information as height and latitude), and
02739        * passes them to appropriate methods.
02740        *
02741        * @param RX  Receiver position.
02742        * @param SV  Satellite position.
02743        * @param doy Day of year.
02744        */
02745    double NeillTropModel::correction( const Position& RX,
02746                                       const Position& SV,
02747                                       const int& doy )
02748       throw(TropModel::InvalidTropModel)
02749    {
02750 
02751       setDayOfYear(doy);
02752 
02753       return NeillTropModel::correction(RX,SV);
02754 
02755    }  // end NeillTropModel::correction(RX,SV,doy)
02756 
02757 
02758 
02759       // deprecated
02760       // Compute and return the full tropospheric delay, given the positions
02761       // of receiver and satellite. . You must set time using method
02762       // setDayOfYear() before calling this method.
02763       //
02764       // @param RX  Receiver position in ECEF cartesian coordinates (meters)
02765       // @param SV  Satellite position in ECEF cartesian coordinates (meters)
02766       // This function is deprecated; use the Position version
02767    double NeillTropModel::correction( const Xvt& RX,
02768                                       const Xvt& SV )
02769       throw(TropModel::InvalidTropModel)
02770    {
02771       Position R(RX),S(SV);
02772 
02773       return NeillTropModel::correction(R,S);
02774 
02775    }  // end NeillTropModel::correction(RX,SV)
02776 
02777 
02778       // deprecated
02779       // Compute and return the full tropospheric delay, given the positions of
02780       // receiver and satellite and the time tag. This version is more useful
02781       // within positioning algorithms, where the receiver position may vary;
02782       // it computes the elevation (and other receiver location information as
02783       // height and latitude) and passes them to appropriate methods.
02784       //
02785       // @param RX  Receiver position in ECEF cartesian coordinates (meters)
02786       // @param SV  Satellite position in ECEF cartesian coordinates (meters)
02787       // @param tt  Time (CommonTime object).
02788       // This function is deprecated; use the Position version
02789    double NeillTropModel::correction( const Xvt& RX,
02790                                       const Xvt& SV,
02791                                       const CommonTime& tt )
02792       throw(TropModel::InvalidTropModel)
02793    {
02794       setDayOfYear(tt);
02795       Position R(RX),S(SV);
02796 
02797       return NeillTropModel::correction(R,S);
02798 
02799    }  // end NeillTropModel::correction(RX,SV,tt)
02800 
02801 
02802       // deprecated
02803       // Compute and return the full tropospheric delay, given the positions of
02804       // receiver and satellite and the day of the year. This version is more
02805       // within positioning algorithms, where the receiver position may vary;
02806       // it computes the elevation (and other receiver location information as
02807       // height and latitude) and passes them to appropriate methods.
02808       //
02809       // @param RX  Receiver position in ECEF cartesian coordinates (meters)
02810       // @param SV  Satellite position in ECEF cartesian coordinates (meters)
02811       // @param doy Day of year.
02812       // This function is deprecated; use the Position version
02813    double NeillTropModel::correction( const Xvt& RX,
02814                                       const Xvt& SV,
02815                                       const int& doy )
02816       throw(TropModel::InvalidTropModel)
02817    {
02818       setDayOfYear(doy);
02819       Position R(RX),S(SV);
02820 
02821       return NeillTropModel::correction(R,S);
02822 
02823    }  // end NeillTropModel::correction(RX,SV,doy)
02824 
02825 
02826       // Compute and return the zenith delay for the dry component of
02827       // the troposphere.
02828    double NeillTropModel::dry_zenith_delay(void) const
02829       throw(TropModel::InvalidTropModel)
02830    {
02831 
02832       if( !valid )
02833       {
02834          throw InvalidTropModel("Invalid model");
02835       }
02836 
02837          // Note: 1.013*2.27 = 2.29951
02838       double ddry( 2.29951*std::exp( (-0.000116 * NeillHeight) ) );
02839 
02840       return ddry;
02841 
02842    }  // end NeillTropModel::dry_zenith_delay()
02843 
02844 
02845       // Compute and return the mapping function for dry component of
02846       // the troposphere.
02847       //
02848       // @param elevation Elevation of satellite as seen at receiver, in
02849       //                  degrees
02850    double NeillTropModel::dry_mapping_function(double elevation) const
02851       throw(TropModel::InvalidTropModel)
02852    {
02853       if(!valid)
02854       {
02855          if(!validLat)
02856          {
02857             GPSTK_THROW( InvalidTropModel( "Invalid Neill trop model: Rx \
02858                                             Latitude" ) );
02859          }
02860 
02861          if(!validHeight)
02862          {
02863             GPSTK_THROW( InvalidTropModel( "Invalid Neill trop model: Rx \
02864                                             Height" ) );
02865          }
02866 
02867          if(!validDOY)
02868          {
02869             GPSTK_THROW( InvalidTropModel( "Invalid Neill trop model: day \
02870                                             of year" ) );
02871          }
02872 
02873          GPSTK_THROW(InvalidTropModel( "Valid flag corrupted in Neill trop \
02874                                         model" ) );
02875       }
02876 
02877       if(elevation < 3.0)
02878       {
02879          return 0.0;
02880       }
02881 
02882       double lat, t, ct;
02883       lat = fabs(NeillLat);         // degrees
02884       t = static_cast<double>(NeillDOY) - 28.0;  // mid-winter
02885 
02886       if(NeillLat < 0.0)              // southern hemisphere
02887       {
02888          t += 365.25/2.;
02889       }
02890 
02891       t *= 360.0/365.25;            // convert to degrees
02892       ct = ::cos(t*DEG_TO_RAD);
02893 
02894       double a, b, c;
02895       if(lat < 15.0)
02896       {
02897          a = NeillDryA[0];
02898          b = NeillDryB[0];
02899          c = NeillDryC[0];
02900       }
02901       else if(lat < 75.)      // coefficients are for 15,30,45,60,75 deg
02902       {
02903          int i=int(lat/15.0)-1;
02904          double frac=(lat-15.*(i+1))/15.;
02905          a = NeillDryA[i] + frac*(NeillDryA[i+1]-NeillDryA[i]);
02906          b = NeillDryB[i] + frac*(NeillDryB[i+1]-NeillDryB[i]);
02907          c = NeillDryC[i] + frac*(NeillDryC[i+1]-NeillDryC[i]);
02908 
02909          a -= ct * (NeillDryA1[i] + frac*(NeillDryA1[i+1]-NeillDryA1[i]));
02910          b -= ct * (NeillDryB1[i] + frac*(NeillDryB1[i+1]-NeillDryB1[i]));
02911          c -= ct * (NeillDryC1[i] + frac*(NeillDryC1[i+1]-NeillDryC1[i]));
02912       }
02913       else
02914       {
02915          a = NeillDryA[4] - ct * NeillDryA1[4];
02916          b = NeillDryB[4] - ct * NeillDryB1[4];
02917          c = NeillDryC[4] - ct * NeillDryC1[4];
02918       }
02919 
02920       double se = ::sin(elevation*DEG_TO_RAD);
02921       double map = (1.+a/(1.+b/(1.+c)))/(se+a/(se+b/(se+c)));
02922 
02923       a = 0.0000253;
02924       b = 0.00549;
02925       c = 0.00114;
02926       map += ( NeillHeight/1000.0 ) *
02927              ( 1./se - ( (1.+a/(1.+b/(1.+c))) / (se+a/(se+b/(se+c))) ) );
02928 
02929       return map;
02930 
02931    }  // end NeillTropModel::dry_mapping_function()
02932 
02933 
02934       // Compute and return the mapping function for wet component of the
02935       // troposphere.
02936       //
02937       // @param elevation Elevation of satellite as seen at receiver,
02938       //                  in degrees.
02939    double NeillTropModel::wet_mapping_function(double elevation) const
02940       throw(TropModel::InvalidTropModel)
02941    {
02942       if(!valid)
02943       {
02944          if(!validLat)
02945          {
02946             GPSTK_THROW( InvalidTropModel( "Invalid Neill trop model: Rx \
02947                                             Latitude" ) );
02948          }
02949          if(!validHeight)
02950          {
02951             GPSTK_THROW( InvalidTropModel( "Invalid Neill trop model: Rx \
02952                                             Height" ) );
02953          }
02954          if(!validDOY)
02955          {
02956             GPSTK_THROW( InvalidTropModel(" Invalid Neill trop model: day \
02957                                             of year" ) );
02958          }
02959 
02960          GPSTK_THROW(InvalidTropModel( "Valid flag corrupted in Neill trop \
02961                                         model" ) );
02962       }
02963 
02964       if(elevation < 3.0)
02965       {
02966          return 0.0;
02967       }
02968 
02969       double a,b,c,lat;
02970       lat = fabs(NeillLat);         // degrees
02971       if(lat < 15.0)
02972       {
02973          a = NeillWetA[0];
02974          b = NeillWetB[0];
02975          c = NeillWetC[0];
02976       }
02977       else if(lat < 75.)          // coefficients are for 15,30,45,60,75 deg
02978       {
02979          int i=int(lat/15.0)-1;
02980          double frac=(lat-15.*(i+1))/15.;
02981          a = NeillWetA[i] + frac*(NeillWetA[i+1]-NeillWetA[i]);
02982          b = NeillWetB[i] + frac*(NeillWetB[i+1]-NeillWetB[i]);
02983          c = NeillWetC[i] + frac*(NeillWetC[i+1]-NeillWetC[i]);
02984       }
02985       else
02986       {
02987          a = NeillWetA[4];
02988          b = NeillWetB[4];
02989          c = NeillWetC[4];
02990       }
02991 
02992       double se = ::sin(elevation*DEG_TO_RAD);
02993       double map = ( 1.+ a/ (1.+ b/(1.+c) ) ) / (se + a/(se + b/(se+c) ) );
02994 
02995       return map;
02996 
02997    }  // end NeillTropModel::wet_mapping_function()
02998 
02999 
03000       // This method configure the model to estimate the weather using height,
03001       // latitude and day of year (DOY). It is called automatically when
03002       // setting those parameters.
03003    void NeillTropModel::setWeather()
03004       throw(TropModel::InvalidTropModel)
03005    {
03006 
03007       if(!validLat)
03008       {
03009          valid = false;
03010          throw InvalidTropModel( "NeillTropModel must have Rx latitude \
03011                                   before computing weather ");
03012       }
03013       if(!validDOY)
03014       {
03015          valid = false;
03016          throw InvalidTropModel( "NeillTropModel must have day of year \
03017                                   before computing weather" );
03018       }
03019 
03020       valid = validHeight && validLat && validDOY;
03021 
03022    }
03023 
03024 
03025       // Define the receiver height; this is required before calling
03026       // correction() or any of the zenith_delay routines.
03027       //
03028       // @param ht   Height of the receiver above mean sea level,
03029       //             in meters.
03030    void NeillTropModel::setReceiverHeight(const double& ht)
03031    {
03032       NeillHeight = ht;
03033       validHeight = true;
03034 
03035          // Change the value of field "valid" if everything is already set
03036       valid = validHeight && validLat && validDOY;
03037 
03038          // If model is valid, set the appropriate parameters
03039       if (valid) setWeather();
03040 
03041    }
03042 
03043 
03044       // Define the receiver latitude; this is required before calling
03045       // correction() or any of the zenith_delay routines.
03046       //
03047       // @param lat  Latitude of receiver, in degrees.
03048    void NeillTropModel::setReceiverLatitude(const double& lat)
03049    {
03050       NeillLat = lat;
03051       validLat = true;
03052 
03053          // Change the value of field "valid" if everything is already set
03054       valid = validHeight && validLat && validDOY;
03055 
03056          // If model is valid, set the appropriate parameters
03057       if (valid) setWeather();
03058 
03059    }
03060 
03061 
03062       // Set the time when tropospheric correction will be computed for,
03063       // in days of the year.
03064       //
03065       // @param doy  Day of the year.
03066    void NeillTropModel::setDayOfYear(const int& doy)
03067    {
03068 
03069       if( (doy>=1) && (doy<=366) )
03070       {
03071          validDOY = true;
03072       }
03073       else
03074       {
03075          validDOY = false;
03076       }
03077 
03078       NeillDOY = doy;
03079 
03080          // Change the value of field "valid" if everything is already set
03081       valid = validHeight && validLat && validDOY;
03082 
03083          // If model is valid, set the appropriate parameters
03084       if (valid) setWeather();
03085 
03086    }
03087 
03088 
03089       // Set the time when tropospheric correction will be computed for,
03090       // in days of the year.
03091       //
03092       // @param time  Time object.
03093    void NeillTropModel::setDayOfYear(const CommonTime& time)
03094    {
03095 
03096       NeillDOY = static_cast<int>((static_cast<YDSTime>(time)).doy);
03097       validDOY = true;
03098 
03099          // Change the value of field "valid" if everything is already set
03100       valid = validHeight && validLat && validDOY;
03101 
03102          // If model is valid, set the appropriate parameters
03103       if (valid) setWeather();
03104 
03105    }
03106 
03107 
03108       /* Convenient method to set all model parameters in one pass.
03109        *
03110        * @param time  Time object.
03111        * @param rxPos Receiver position object.
03112        */
03113    void NeillTropModel::setAllParameters( const CommonTime& time,
03114                                           const Position& rxPos )
03115    {
03116         YDSTime ydst = static_cast<YDSTime>(time);
03117       NeillDOY = static_cast<int>(ydst.doy);
03118       validDOY = true;
03119       NeillLat = rxPos.getGeodeticLatitude();
03120       validHeight = true;
03121       NeillLat = rxPos.getHeight();
03122       validLat = true;
03123 
03124          // Change the value of field "valid" if everything is already set
03125       valid = validHeight && validLat && validDOY;
03126 
03127          // If model is valid, set the appropriate parameters
03128       if (valid) setWeather();
03129 
03130    }
03131 
03132 
03133 } // end namespace gpstk

Generated on Sun May 26 03:31:13 2013 for GPS ToolKit Software Library by  doxygen 1.3.9.1