example-a.cpp

GPS Data processing: Code and Phase

// Example A
// This program shows how to use some high-level GPSTk classes
// Dagoberto Salazar - gAGE. 2006

// Basic input/output C++ class
#include <iostream>

// Classes for handling observations RINEX files (data)
#include "RinexObsBase.hpp"
#include "RinexObsData.hpp"
#include "RinexObsStream.hpp"

// Class to extract C1 data from RinexObsData objects
#include "ExtractC1.hpp"

// Classes for handling satellite navigation parameters RINEX files (Broadcast ephemerides)
#include "RinexNavBase.hpp"
#include "RinexNavHeader.hpp"
#include "RinexNavData.hpp"
#include "RinexNavStream.hpp"

// Class to store satellite broadcast navigation data
#include "BCEphemerisStore.hpp"

// Class to model GPS data for a mobile receiver
#include "ModeledPR.hpp"

// Class to model the tropospheric delays
#include "TropModel.hpp"

// Classes to model ans store ionospheric delays
#include "IonoModel.hpp"
#include "IonoModelStore.hpp"

// Class to solve the equations system using a Least Mean Square method
#include "SolverLMS.hpp"

// Basic framework for programs in the GPSTk. process() method MUST be implemented
#include "BasicFramework.hpp"

#include "geometry.hpp"                   // DEG_TO_RAD



using namespace std;
using namespace gpstk;


// A new class is declared that will handle program behaviour
// This class inherits from BasicFramework
class examplea : public BasicFramework
{
public:
    // Constructor declaration
    examplea(char* arg0);

protected:
    // Method that will take care of processing
    virtual void process();

    // Method that hold code to be run BEFORE processing
    virtual void spinUp();

private:
    // These field represent options at command line interface (CLI)
    CommandOptionWithArg dataFile;
    CommandOptionWithArg navFile;

    // If you want to share objects and variables among methods, you'd better declare them here
    RinexObsStream rObsFile;            // Object to read Rinex observation data files
    RinexObsData rData;                 // Object to store Rinex observation data
    RinexNavStream rNavFile;            // Object to read Rinex navigation data files
    RinexNavData rNavData;              // Object to store Rinex navigation data
    RinexNavHeader rNavHeader;          // Object to read the header of Rinex navigation data files
    IonoModelStore ionoStore;           // Object to store ionospheric models
    BCEphemerisStore bceStore;          // Object to store ephemeris
    ModeledPR modelPR;                  // Declare a ModeledReferencePR object
    GCATTropModel gcatTM;               // Declare a GCATTropModel object
//    MOPSTropModel mopsTM;               // Declare a MOPSTropModel object
    ExtractC1 obsC1;                    // Declare an ExtractData object
    bool useFormerPos;                  // Flag indicating if we have an a priori position
    Position formerPosition;            // Object to store the former position
    IonoModel ioModel;                  // Declare a Ionospheric Model object
    SolverLMS solver;                   // Declare an object to apply LMS method

};


// Let's implement constructor details
examplea::examplea(char* arg0)
      : BasicFramework(arg0, "\nexample-a: Program to print the position solution with respect to a nominal solution,"
      " given a RINEX observations file and a RINEX broadcast navigation file.\n\n"
      "The output is: \n"
      "  Time(seconds of day)\t dLon(m)\t dLat(m)\t dAlt(m)\n\n"
      "This program uses C1, Klobuchar ionospheric model, "
      "a simple tropospheric model and a Least Mean Squares solver.\n"),
        // Option initialization. "true" means a mandatory option
        dataFile(CommandOption::stdType, 'i', "datainput", 
                        " [-i|--datainput]      Name of RINEX observations file.", true), 
        navFile(CommandOption::stdType, 'n', "navinput", 
                        " [-n|--navinput]      Name of RINEX broadcast navigation file.", true)
{
    // These options may appear just once at CLI
    dataFile.setMaxCount(1);
    navFile.setMaxCount(1);
} // End of constructor details


// Method that will be executed AFTER initialization but BEFORE processing
void examplea::spinUp()
{

    // From now on, some parts may look similar to example3.cpp and example4.cpp
    // Activate failbit to enable exceptions
    rObsFile.exceptions(ios::failbit);
    // First, data RINEX reading object
    try
    {
        rObsFile.open(dataFile.getValue()[0].c_str(), std::ios::in);
    }
    catch(...)
    {
        cerr << "Problem opening file " << dataFile.getValue()[0].c_str() << endl;
        cerr << "Maybe it doesn't exist or you don't have proper read permissions." << endl;
        exit (-1);
    }


    // Activate failbit to enable exceptions
    rNavFile.exceptions(ios::failbit);
    // Read nav file and store unique list of ephemerides
    try
    {
        rNavFile.open(navFile.getValue()[0].c_str(), std::ios::in);
    }
    catch(...)
    {
        cerr << "Problem opening file " << navFile.getValue()[0].c_str() << endl;
        cerr << "Maybe it doesn't exist or you don't have proper read permissions." << endl;
        exit (-1);
    }


    // We will need to read ionospheric parameters (Klobuchar model) from header
    rNavFile >> rNavHeader;

    // Let's feed the ionospheric model (Klobuchar type) from data in the Navigation file header
    ioModel.setModel(rNavHeader.ionAlpha, rNavHeader.ionBeta);
    // WARNING-WARNING-WARNING: In this case, the same model will be used for the full data span
    ionoStore.addIonoModel(DayTime::BEGINNING_OF_TIME, ioModel);

    // Storing the ephemeris in "bceStore"
    while (rNavFile >> rNavData) bceStore.addEphemeris(rNavData);
    // Setting the criteria for looking up ephemeris
    bceStore.SearchPast();  // This is the default

    // This is set to true if the former computed positon will be used as a priori position
    useFormerPos = false;   // At first, we don't have an a priori position

    // Prepare for printing later on
    cout << fixed << setprecision(8);

}// End of examplea::spinUp()


// Method that will really process information
void examplea::process()
{

    // Let's read the observations RINEX, epoch by epoch
    while (rObsFile >> rData)
    {
        if ( (rData.epochFlag == 0 || rData.epochFlag == 1) && (rData.numSvs > 3) )  // Begin usable data with enough number of satellites
        {
            int validSats = 0;   // Number of satellites with valid data in this epoch
            int prepareResult;
            double rxAltitude;  // Receiver altitude for tropospheric model
//            double rxLatitude;  // Receiver latitude for tropospheric model

            // We need extract C1 data from this epoch. Skip epoch if not enough data (4 SV at least)
            if ( obsC1.getData(rData) < 4) {
                useFormerPos = false;  // The former position will not be valid next time
                continue;
            }


            if (useFormerPos) {  // If possible, use former position as a priori
                prepareResult = modelPR.Prepare(formerPosition);
                // We need to seed this kind of tropospheric model with receiver altitude
                rxAltitude = formerPosition.getAltitude();
//                rxLatitude = formerPosition.getGeodeticLatitude();
            }
            else {  // Use Bancroft method is no a priori position is available
                cerr << "Bancroft method was used at epoch " << rData.time.DOYsecond() << endl;
                prepareResult = modelPR.Prepare(rData.time, obsC1.availableSV, obsC1.obsData, bceStore);
                // We need to seed this kind of tropospheric model with receiver altitude
                rxAltitude = modelPR.rxPos.getAltitude();
//                rxLatitude = modelPR.rxPos.getGeodeticLatitude();
            }

            if (prepareResult) {  // If there were problems with Prepare(), skip this epoch
                useFormerPos = false;  // The former position will not be valid next time
                continue;
            }
            // If there were no problems, let's feed the tropospheric model
            gcatTM.setReceiverHeight(rxAltitude);
//            mopsTM.setReceiverLatitude(rxLatitude);
//            mopsTM.setDayOfYear(rData.time.DOY());


            // Now, let's compute the GPS model for our observable (C1)
            validSats = modelPR.Compute(rData.time, obsC1.availableSV, obsC1.obsData, bceStore, &gcatTM, &ionoStore);


            // Only get into further computations if there are enough satellites
            if (validSats >= 4) {

                // Now, let's solve the navigation equations using the WMS method
                try {
                    // Then, solve the system
                    solver.Compute(modelPR.prefitResiduals, modelPR.geoMatrix);
                }
                catch(InvalidSolver& e) {
                    cerr << "Couldn't solve equations system at epoch " << rData.time.DOYsecond() << endl;
                    cerr << e << endl;
                    useFormerPos = false;  // The former position will not be valid next time
                    continue;
                }

                // With "solver", we got the difference vector between a priori position and
                // real position. Then, let's convert the solution to a Position object
                Position solPos( (modelPR.rxPos.X() + solver.solution[0]), (modelPR.rxPos.Y() + solver.solution[1]), (modelPR.rxPos.Z() + solver.solution[2]) );

                // Known position
                Position nominalPos(4833520.2269, 41537.00768, 4147461.489);
                // Difference in position
                Position diffPos;
                diffPos = solPos - nominalPos;
                double azimuth = nominalPos.azimuthGeodetic(solPos);
                double elev = nominalPos.elevationGeodetic(solPos);
                double magnitude = RSS(diffPos.X(), diffPos.Y(), diffPos.Z());

                // Print results
                cout << rData.time.DOYsecond()  << "   ";   // Output field #1
                // Longitude change
                cout << magnitude*sin(azimuth*DEG_TO_RAD)*cos(elev*DEG_TO_RAD) << "   ";
                // Latitude change
                cout << magnitude*cos(azimuth*DEG_TO_RAD)*cos(elev*DEG_TO_RAD) << "   ";
                // Altitude change
                cout << magnitude*sin(elev*DEG_TO_RAD) << "   ";
                cout << endl;

                formerPosition = solPos;
                useFormerPos = true;  // Next time, former position will be used as a priori

            }  // End of if (validSats...
            else {
                useFormerPos = false;  // The former position will not be valid next time
            }
        } // End of "if" for usable data
        else {
            useFormerPos = false;  // The former position will not be valid next time
        }
    } // End of while

    return;

}// End of examplea::process()


// Main function
int main(int argc, char* argv[])
{
   try
   {
      examplea program(argv[0]);
      if (!program.initialize(argc, argv))
         return 0;
      if (!program.run())
         return 1;

      return 0;
   }
   catch(Exception& e)
   {
      cout << "Problem: " << e << endl;
      return 1;
   }
   catch(...)
   {
      cout << "Unknown error." << endl;
      return 1;
   }
   return 0;
}