Merge branch 'update_gamma_gps' into 'master'

Update gamma GPS

See merge request sgl/gpstk!470

core/lib/ClockModel/ObsRngDev.cpp

@@ -83,7 +83,7 @@ namespace gpstk
       const XvtStore<SatID>& eph,
       EllipsoidModel& em,
       const IonoModelStore& ion,
-      IonoModel::Frequency fq,
+      CarrierBand band,
       bool svTime)
          : obstime(time), svid(svid), ord(0), wonky(false)
    {
@@ -94,7 +94,7 @@ namespace gpstk
                      gx.getGeodeticLatitude(),
                      static_cast<YDSTime>(time).doy);
       computeTrop(nb);
-      iono = ion.getCorrection(time, gx, elevation, azimuth, fq);
+      iono = ion.getCorrection(time, gx, elevation, azimuth, band);
       ord -= iono;
    }
 
@@ -122,15 +122,15 @@ namespace gpstk
       EllipsoidModel& em,
       const TropModel& tm,
       const IonoModelStore& ion,
-      IonoModel::Frequency fq,
+      CarrierBand band,
       bool svTime)
          : obstime(time), svid(svid), ord(0), wonky(false)
    {
       computeOrd(prange, rxpos, eph, em, svTime);
       computeTrop(tm);
       Position gx(rxpos);
       gx.asGeodetic(&em);
-      iono = ion.getCorrection(time, gx, elevation, azimuth, fq);
+      iono = ion.getCorrection(time, gx, elevation, azimuth, band);
       ord -= iono;
    }
 
@@ -147,7 +147,7 @@ namespace gpstk
       double gamma)
          : obstime(time), svid(svid), ord(0), wonky(false)
    {
-      // for dual frequency see IS-GPS-200, section 20.3.3.3.3.3
+      // for dual-frequency see IS-GPS-200, section 20.3.3.3.3.3
       double icpr = (prange2 - gamma * prange1)/(1-gamma);
       iono = prange1 - icpr;
 
@@ -174,7 +174,7 @@ namespace gpstk
       double gamma)
          : obstime(time), svid(svid), ord(0), wonky(false)
    {
-      // for dual frequency see IS-GPS-200, section 20.3.3.3.3.3
+      // for dual-frequency see IS-GPS-200, section 20.3.3.3.3.3
       double icpr = (prange2 - gamma * prange1)/(1-gamma);
       iono = prange1 - icpr;
 

core/lib/ClockModel/ObsRngDev.hpp

@@ -47,6 +47,7 @@
 #include <ostream>
 
 #include "CommonTime.hpp"
+#include "CarrierBand.hpp"
 #include "XvtStore.hpp"
 #include "Exception.hpp"
 #include "GPSEllipsoid.hpp"
@@ -90,7 +91,7 @@ namespace gpstk
           * @param eph a store of either broadcast or precise ephemerides
           * @param em an EllipsoidModel for performing range calculations
           * @param ion a store of nav based ionospheric models
-          * @param fq the GPS frequency (L1 or L2) from which the obs was made
+          * @param fq the GPS band (L1, L2, L5) from which the obs was made
           * @param svTime true if prange is in SV time, false for RX time.
           */ 
       ObsRngDev(const double prange,
@@ -111,7 +112,7 @@ namespace gpstk
           * @param eph a store of either broadcast or precise ephemerides
           * @param em an EllipsoidModel for performing range calculations
           * @param ion a store of nav based ionospheric models
-          * @param fq the GPS frequency (L1 or L2) from which the obs was made
+          * @param fq the GPS band (L1, L2, L5) from which the obs was made
           * @param svTime true if prange is in SV time, false for RX time.
           */ 
       ObsRngDev(const double prange,
@@ -121,7 +122,7 @@ namespace gpstk
                 const XvtStore<SatID>& eph,
                 EllipsoidModel& em,
                 const IonoModelStore& ion,
-                IonoModel::Frequency fq,
+                CarrierBand band,
                 bool svTime = false);
          /**
           * constructor.
@@ -135,7 +136,7 @@ namespace gpstk
           * @param em an EllipsoidModel for performing range calculations
           * @param tm a TropModel for performing trop calculation
           * @param ion a store of nav based ionospheric models
-          * @param fq the GPS frequency (L1 or L2) from which the obs was made
+          * @param fq the GPS band (L1, L2, L5) from which the obs was made
           * @param svTime true if prange is in SV time, false for RX time.
           */ 
       ObsRngDev(const double prange,
@@ -159,7 +160,7 @@ namespace gpstk
           * @param em an EllipsoidModel for performing range calculations
           * @param tm a TropModel for performing trop calculation
           * @param ion a store of nav based ionospheric models
-          * @param fq the GPS frequency (L1 or L2) from which the obs was made
+          * @param fq the GPS band (L1, L2, L5) from which the obs was made
           * @param svTime true if prange is in SV time, false for RX time.
           */ 
       ObsRngDev(const double prange,
@@ -170,7 +171,7 @@ namespace gpstk
                 EllipsoidModel& em,
                 const TropModel& tm,
                 const IonoModelStore& ion,
-                IonoModel::Frequency fq,
+                CarrierBand band,
                 bool svTime = false);
 
          /**

core/lib/GNSSCore/GNSSconstants.hpp

@@ -88,8 +88,10 @@ namespace gpstk
    const double L2_MULT_GPS   = 120.0;
       /// GPS L5 frequency in units of oscillator frequency.
    const double L5_MULT_GPS   = 115.0;
-      /// GPS Gamma constant
-   const double GAMMA_GPS = 1.646944444;
+      /// GPS Gamma constants
+   const double GAMMA_GPS_12 = (L1_MULT_GPS/L2_MULT_GPS) * (L1_MULT_GPS/L2_MULT_GPS);
+   const double GAMMA_GPS_15 = (L1_MULT_GPS/L5_MULT_GPS) * (L1_MULT_GPS/L5_MULT_GPS);
+   const double GAMMA_GPS = 1.646944444;  // legacy notation and hard-coded value; wherefore the extra "44"??
       /// Reference Semi-major axis. From IS-GPS-800 Table 3.5-2 in meters.
    const double A_REF_GPS = 26559710.0;
       /// Omega reference value from Table 30-I converted to radians

core/lib/GNSSCore/IonoModel.cpp

@@ -84,7 +84,7 @@ namespace gpstk
                                    const Position& rxgeo,
                                    double svel,
                                    double svaz,
-                                   Frequency freq) const
+                                   CarrierBand band) const
    {
 
       if (!valid)
@@ -141,13 +141,21 @@ namespace gpstk
       else
          t_iono = iF * 5.0e-9;
 
-      if (freq == L2)
+      // Correction factor for GPS band; see ICD-GPS-200 20.3.3.3.3.2.
+      if (band == CarrierBand::L2)
       {
-            // see ICD-GPS-200 20.3.3.3.3.2
-         t_iono *= GAMMA_GPS;  //  GAMMA_GPS = (fL1 / fL2)^2
+         t_iono *= GAMMA_GPS_12;  //  GAMMA_GPS = (fL1 / fL2)^2
+      }
+      else if (band == CarrierBand::L5)
+      {
+         t_iono *= GAMMA_GPS_15;  //  GAMMA_GPS = (fL1 / fL5)^2
+      }
+      else if (band != CarrierBand::L1)
+      {
+         throw InvalidIonoModel("Invalid CarrierBand, not one of L1,L2,L5.");
       }
 
-      double correction = t_iono * C_MPS;
+      double correction = t_iono * C_MPS;  // return correction in [m]
 
       return correction;
    }

core/lib/GNSSCore/IonoModel.hpp

@@ -45,6 +45,7 @@
 #define GPSTK_IONOMODEL_HPP
 
 #include "CommonTime.hpp"
+#include "CarrierBand.hpp"
 #include "EngAlmanac.hpp"
 #include "Position.hpp"
 
@@ -77,13 +78,6 @@ namespace gpstk
          /// @ingroup exceptiongroup
       NEW_EXCEPTION_CLASS(InvalidIonoModel, gpstk::Exception);
 
-
-      enum Frequency
-      {
-         L1,  ///< L1 frequency (1575.42 MHz)
-         L2   ///< L2 frequency (1227.60 MHz)
-      };
-
          /// default constructor, creates an invalid model
       IonoModel() throw() : valid(false) {}
 
@@ -125,15 +119,15 @@ namespace gpstk
           * @param rxgeo the WGS84 geodetic position of the receiver
           * @param svel the elevation angle between the rx and SV (degrees)
           * @param svaz the azimuth angle between the rx and SV (degrees)
-          * @param freq the GPS frequency the observation was made from
+          * @param band the GPS frequency band the observation was made from
           * @return the ionospheric correction (meters)
           * @throw InvalidIonoModel
           */
       double getCorrection(const CommonTime& time,
                            const Position& rxgeo,
                            double svel,
                            double svaz,
-                           Frequency freq = L1) const;
+                           CarrierBand band = CarrierBand::L1) const;
 
          /// equality operator
       bool operator==(const IonoModel& right) const throw();

core/lib/GNSSCore/IonoModelStore.cpp

@@ -54,21 +54,21 @@ namespace gpstk
        * \param rxgeo the WGS84 geodetic position of the receiver
        * \param svel the elevation angle between the rx and SV (degrees)
        * \param svaz the azimuth angle between the rx and SV (degrees)
-       * \param freq the GPS frequency the observation was made from
+       * \param band the GPS band the observation was made from
        * \return the ionospheric correction (meters)
        */
    double IonoModelStore::getCorrection(const CommonTime& time,
                                         const Position& rxgeo,
                                         double svel,
                                         double svaz,
-                                        IonoModel::Frequency freq) const
+                                        CarrierBand band) const
    {
 
       IonoModelMap::const_iterator i = ims.upper_bound(time);
       if (!ims.empty() && i != ims.begin())
       {
          i--;
-         return i->second.getCorrection(time, rxgeo, svel, svaz, freq);
+         return i->second.getCorrection(time, rxgeo, svel, svaz, band);
       }
       else
       {

core/lib/GNSSCore/IonoModelStore.hpp

@@ -47,6 +47,7 @@
 
 #include <map>
 #include "CommonTime.hpp"
+#include "CarrierBand.hpp"
 #include "IonoModel.hpp"
 
 namespace gpstk
@@ -85,15 +86,15 @@ namespace gpstk
           * @param rxgeo the WGS84 geodetic position of the receiver
           * @param svel the elevation angle between the rx and SV (degrees)
           * @param svaz the azimuth angle between the rx and SV (degrees)
-          * @param freq the GPS frequency the observation was made from
+          * @param freq the GPS band the observation was made from
           * @return the ionospheric correction (meters)
           * @throw NoIonoModelFound
           */
       virtual double getCorrection(const CommonTime& time,
                                    const Position& rxgeo,
                                    double svel,
                                    double svaz,
-                                   IonoModel::Frequency freq = IonoModel::L1)
+                                   CarrierBand band = CarrierBand::L1)
          const;
 
 

core/lib/ORD/ord.cpp

@@ -99,26 +99,22 @@ double IonosphereFreeRange(const std::vector<double>& frequencies,
     const double gamma = (frequencies[0]/frequencies[1]) *
                          (frequencies[0]/frequencies[1]);
 
-    // for dual frequency see IS-GPS-200, section 20.3.3.3.3.3
+    // for dual-frequency see IS-GPS-200, section 20.3.3.3.3.3
     double icpr = (pseudoranges[1] - gamma * pseudoranges[0])/(1-gamma);
 
     return icpr;
 }
 
 double IonosphereModelCorrection(const gpstk::IonoModelStore& ionoModel,
-        const gpstk::CommonTime& time, double freq,
+        const gpstk::CommonTime& time, CarrierBand band,
         const gpstk::Position& rxLoc, const gpstk::Xvt& svXvt) {
     Position trx(rxLoc);
     Position svPos(svXvt);
 
     double elevation = trx.elevation(svPos);
     double azimuth = trx.azimuth(svPos);
 
-    // TODO(someone): IonoModel assumes only L1 and L2 frequencies, this
-    // should be updated to work with an arbitrary frequency.
-    double iono = ionoModel.getCorrection(time, trx,
-                                          elevation, azimuth,
-                                          IonoModel::L1);
+    double iono = ionoModel.getCorrection(time, trx, elevation, azimuth, band);
     return -iono;
 }
 
@@ -276,14 +272,22 @@ double TroposphereCorrection(const gpstk::TropModel& tropModel,
     return trop;
 }
 
-double calculate_ord(const std::vector<double>& frequencies,
-                     const std::vector<double>& pseudoranges, const gpstk::Position& rx_loc,
-                     const gpstk::SatID& sat_id, const gpstk::CommonTime& transmit_time,
+/*
+ * Example not fully fleshed-out.  If dual-band data given, for example,
+ * then the last IonosphereModelCorrection call must not be made.
+ * Users should construct an ORD algorithm based on their data and use case.
+ *
+double calculate_ord(const std::vector<CarrierBand>& bands,
+                     const std::vector<double>& pseudoranges,
+                     const gpstk::Position& rx_loc,
+                     const gpstk::SatID& sat_id,
+                     const gpstk::CommonTime& transmit_time,
                      const gpstk::CommonTime& receive_time,
                      const gpstk::IonoModelStore& iono_model,
                      const gpstk::TropModel& trop_model,
-                     const gpstk::XvtStore<gpstk::SatID>& ephemeris, int range_method) {
-    double ps_range = IonosphereFreeRange(frequencies, pseudoranges);
+                     const gpstk::XvtStore<gpstk::SatID>& ephemeris,
+                     int range_method) {
+    double ps_range = IonosphereFreeRange(bands, pseudoranges);
 
     gpstk::Xvt sv_xvt;
     // find raw_range
@@ -314,13 +318,14 @@ double calculate_ord(const std::vector<double>& frequencies,
     // apply troposphere model correction
     range += TroposphereCorrection(trop_model, rx_loc, sv_xvt);
 
-    // apply ionosphere model correction
-    range += IonosphereModelCorrection(iono_model, receive_time,
-                                       frequencies[0],
+    // apply ionosphere model correction -- GPS only (this is Klobuchar)
+    range += IonosphereModelCorrection(iono_model,
+                                       receive_time, bands[0],
                                        rx_loc, sv_xvt);
 
     return ps_range - range;
 }
+ */
 
 }  // namespace ord
 }  // namespace gpstk