diff --git a/core/lib/ClockModel/ObsRngDev.cpp b/core/lib/ClockModel/ObsRngDev.cpp index cb030b496..297d5cc1a 100644 --- a/core/lib/ClockModel/ObsRngDev.cpp +++ b/core/lib/ClockModel/ObsRngDev.cpp @@ -83,7 +83,7 @@ namespace gpstk const XvtStore& 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(time).doy); computeTrop(nb); - iono = ion.getCorrection(time, gx, elevation, azimuth, fq); + iono = ion.getCorrection(time, gx, elevation, azimuth, band); ord -= iono; } @@ -122,7 +122,7 @@ 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) { @@ -130,7 +130,7 @@ namespace gpstk 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; diff --git a/core/lib/ClockModel/ObsRngDev.hpp b/core/lib/ClockModel/ObsRngDev.hpp index e79b6d0b1..4f1a40c0c 100644 --- a/core/lib/ClockModel/ObsRngDev.hpp +++ b/core/lib/ClockModel/ObsRngDev.hpp @@ -47,6 +47,7 @@ #include #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& 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); /** diff --git a/core/lib/GNSSCore/GNSSconstants.hpp b/core/lib/GNSSCore/GNSSconstants.hpp index 47cc60994..19fe0ac5d 100644 --- a/core/lib/GNSSCore/GNSSconstants.hpp +++ b/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 diff --git a/core/lib/GNSSCore/IonoModel.cpp b/core/lib/GNSSCore/IonoModel.cpp index 796f65a4a..a39abffd8 100644 --- a/core/lib/GNSSCore/IonoModel.cpp +++ b/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; } diff --git a/core/lib/GNSSCore/IonoModel.hpp b/core/lib/GNSSCore/IonoModel.hpp index f992f74e6..dcfa1cb9f 100644 --- a/core/lib/GNSSCore/IonoModel.hpp +++ b/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,7 +119,7 @@ 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 */ @@ -133,7 +127,7 @@ namespace gpstk 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(); diff --git a/core/lib/GNSSCore/IonoModelStore.cpp b/core/lib/GNSSCore/IonoModelStore.cpp index 9434e4756..952f8e301 100644 --- a/core/lib/GNSSCore/IonoModelStore.cpp +++ b/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 { diff --git a/core/lib/GNSSCore/IonoModelStore.hpp b/core/lib/GNSSCore/IonoModelStore.hpp index 39407f946..3879d4a1d 100644 --- a/core/lib/GNSSCore/IonoModelStore.hpp +++ b/core/lib/GNSSCore/IonoModelStore.hpp @@ -47,6 +47,7 @@ #include #include "CommonTime.hpp" +#include "CarrierBand.hpp" #include "IonoModel.hpp" namespace gpstk @@ -85,7 +86,7 @@ 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 */ @@ -93,7 +94,7 @@ namespace gpstk const Position& rxgeo, double svel, double svaz, - IonoModel::Frequency freq = IonoModel::L1) + CarrierBand band = CarrierBand::L1) const; diff --git a/core/lib/ORD/ord.cpp b/core/lib/ORD/ord.cpp index 4cfd830e9..9eb2aaacb 100644 --- a/core/lib/ORD/ord.cpp +++ b/core/lib/ORD/ord.cpp @@ -99,14 +99,14 @@ double IonosphereFreeRange(const std::vector& 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); @@ -114,11 +114,7 @@ double IonosphereModelCorrection(const gpstk::IonoModelStore& ionoModel, 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& frequencies, - const std::vector& 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& bands, + const std::vector& 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& ephemeris, int range_method) { - double ps_range = IonosphereFreeRange(frequencies, pseudoranges); + const gpstk::XvtStore& 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& 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 diff --git a/core/lib/ORD/ord.hpp b/core/lib/ORD/ord.hpp index 1abf542b2..1c719b057 100644 --- a/core/lib/ORD/ord.hpp +++ b/core/lib/ORD/ord.hpp @@ -72,9 +72,6 @@ double IonosphereFreeRange(const std::vector& frequencies, /// Given an ionosphere model, and locations of receiver and satellite, /// range correction due to ionospheric effects. -/// TODO(someone): IonoModel assumes only L1 and L2 frequencies, this -/// should be updated to work with an arbitrary frequency. Currently -/// This call assumes frequency is L1. /// @param ionoModel Class that encapsulates ionospheric models /// @params time The time of interest. /// @params frequency Frequency of interest - see note above. @@ -82,7 +79,7 @@ double IonosphereFreeRange(const std::vector& frequencies, /// @param sv_loc The location of the satellite at time of interest. /// @return Range correction (delta) in meters double IonosphereModelCorrection(const gpstk::IonoModelStore& ionoModel, - const gpstk::CommonTime& time, double frequency, + const gpstk::CommonTime& time, CarrierBand band, const gpstk::Position& rxLoc, const gpstk::Xvt& svXvt); /// Given a satellite id, a time, and an ephemeris store, retrieves the @@ -164,28 +161,36 @@ double SvRelativityCorrection(gpstk::Xvt& svXvt); double TroposphereCorrection(const gpstk::TropModel& trop_model, const gpstk::Position& rx_loc, const gpstk::Xvt& sv_xvt); -/// Example method that applies _all_ corrections to generate the -/// Observed Range Deviation. +/// Example method that applies _all_ corrections to generate an Observed Range Deviation. /// This is intended to be a sample showing how the above methods will be used. +/// The example is not fully developed, just a general sketch of a generic approach. +/// E.g., if dual-band, there should be no additional ionosphere correction applied. +/// Users should construct an ORD based on their data and use case. /// Parameters: -/// @params frequencies Signal frequencies. -/// @params pseudoranges Pseudorange values, corresponding to frequency array. -/// @params trop_model Class that encapsulates ionospheric models -/// @params rx_loc The location of the receiver. -/// @params sat_id Identifier for the satellite -/// @params transmit_time The transmit time reported by satellite. -/// @params receive_time The nominal receive time. -/// @params iono_model Class that encapsulates ionospheric models -/// @params trop_model Class that encapsulates troposphere models -/// @params ephemeris The ephemeris to query against. -/// @returns Observed range deviation from 1st pseudorange -double calculate_ord(const std::vector& frequencies, - const std::vector& pseudoranges, const gpstk::Position& rx_loc, - const gpstk::SatID& sat_id, const gpstk::CommonTime& transmit_time, +/// @param[in] bands -- Signal bands (one or two enums, for single-band/dual-band). +/// @param[in] pseudoranges -- Pseudorange values, corresponding to frequency array (one or two). +/// @param[in] trop_model -- Class that encapsulates ionospheric models. +/// @param[in] rx_loc -- The location of the receiver. +/// @param[in] sat_id -- Identifier for the satellite. +/// @param[in] transmit_time -- The transmit time reported by satellite. +/// @param[in] receive_time -- The nominal receive time. +/// @param[in] iono_model -- Class that encapsulates ionospheric models. +/// @param[in] trop_model -- Class that encapsulates troposphere models. +/// @param[in] ephemeris -- The ephemeris to query against. +/// @param[in] range_method -- One of four raw range methods, depending on what data is available. +/// @returns Observed range deviation from 1st pseudorange. +/* +double calculate_ord(const std::vector& bands, + const std::vector& 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& ephemeris, int range_method); + const gpstk::XvtStore& ephemeris, + int range_method); + */ } // namespace ord } // namespace gpstk diff --git a/core/tests/ClockModel/ObsRngDev_T.cpp b/core/tests/ClockModel/ObsRngDev_T.cpp index 9ff69be34..3f4b889f5 100644 --- a/core/tests/ClockModel/ObsRngDev_T.cpp +++ b/core/tests/ClockModel/ObsRngDev_T.cpp @@ -141,7 +141,7 @@ class ObsRngDev_T TUDEF("ObsRngDev", "IonosphericConstructor"); gpstk::IonoModelStore ims = ionoModelStoreGen(cTimeVec); - gpstk::IonoModel::Frequency L1 = gpstk::IonoModel::L1; + gpstk::CarrierBand L1 = gpstk::CarrierBand::L1; try { @@ -229,7 +229,7 @@ class ObsRngDev_T gpstk::SimpleTropModel stm(18.8889, 1021.2176, 77.7777); gpstk::IonoModelStore ims = ionoModelStoreGen(cTimeVec); - gpstk::IonoModel::Frequency L1 = gpstk::IonoModel::L1; + gpstk::CarrierBand L1 = gpstk::CarrierBand::L1; try { diff --git a/core/tests/GNSSCore/IonoModel_T.cpp b/core/tests/GNSSCore/IonoModel_T.cpp index ed99eb364..83107daa5 100644 --- a/core/tests/GNSSCore/IonoModel_T.cpp +++ b/core/tests/GNSSCore/IonoModel_T.cpp @@ -211,7 +211,7 @@ int IonoModel_T :: exceptionTest( void ) // //---------------------------------------- - // // What the hell is this? Commenting it out until someone else figures it out. + // // What is this? Commenting it out until someone else figures it out. // //---------------------------------------- // // try @@ -227,7 +227,7 @@ int IonoModel_T :: exceptionTest( void ) //---------------------------------------- try { - Model.getCorrection( commonTime, rxgeo,svel, svaz, Model.L1 ); + Model.getCorrection( commonTime, rxgeo,svel, svaz, CarrierBand::L1 ); assert_message = "getCorrection(), This test should have thrown an InvalidIonoModel exception"; test4.assert( false, assert_message, __LINE__ ); } @@ -245,7 +245,7 @@ int IonoModel_T :: exceptionTest( void ) //---------------------------------------- try { - goodModel.getCorrection( commonTime, rxgeo, svel, svaz, Model.L1 ); + goodModel.getCorrection( commonTime, rxgeo, svel, svaz, CarrierBand::L1 ); assert_message = "getCorrection( L1 ), This test should NOT throw an exception"; test4.assert( true, assert_message, __LINE__ ); } @@ -265,7 +265,7 @@ int IonoModel_T :: exceptionTest( void ) //---------------------------------------- try { - goodModel.getCorrection( commonTime, rxgeo, svel, svaz, Model.L2 ); + goodModel.getCorrection( commonTime, rxgeo, svel, svaz, CarrierBand::L2 ); assert_message = "getCorrection( L2 ), This test should NOT throw an exception"; test4.assert( true, assert_message, __LINE__ ); } @@ -285,7 +285,7 @@ int IonoModel_T :: exceptionTest( void ) //---------------------------------------- try { - goodModel.getCorrection(commonTime,rxgeo,72.,45.,Model.L1); + goodModel.getCorrection( commonTime, rxgeo, 72., 45., CarrierBand::L1 ); assert_message = "getCorrection( commonTime,rxgeo,72.,45.,Model.L1 ), This test should NOT throw an exception"; test4.assert( true, assert_message, __LINE__ ); } diff --git a/swig/src/pythonfunctions.i b/swig/src/pythonfunctions.i index bc2d05156..d2458b7bd 100644 --- a/swig/src/pythonfunctions.i +++ b/swig/src/pythonfunctions.i @@ -2,6 +2,8 @@ // Python stuff /////////////////////////////////////////////// %pythoncode %{ + + def now(timeSystem=TimeSystem.UTC): """ Returns the current time (defined by what SystemTime() returns) @@ -45,6 +47,13 @@ def times(starttime, endtime, seconds=0.0, days=0): t.addDays(days) +def klobuchar_correction(iono, com_time, rx_pos, elev, azim, band=CarrierBand.L1): + """ + Returns the Klobuchar model ionospheric correction for a time, position, SV elev/azim, and GPS band (L1,L2,L5). + """ + return iono.getCorrection(com_time, rx_pos, float(elev), float(azim), band) + + def moonPosition(time): """ Returns the current position (A gpstk.Triple) of the moon. @@ -69,6 +78,7 @@ def poleTides(time, position, x, y): """ return PoleTides().getPoleTide(time, position, x, y) + def solidTides(time, position): """ Returns the effect (a gpstk.Triple) of solid Earth tides (meters)