Added some gpstk examples (modified)
authorBenoit Papillault <benoit.papillault@free.fr>
Wed, 21 Dec 2011 21:57:27 +0000 (22:57 +0100)
committerBenoit Papillault <benoit.papillault@free.fr>
Wed, 21 Dec 2011 21:57:27 +0000 (22:57 +0100)
gpstk-example1.cpp [new file with mode: 0644]
gpstk-example2.cpp [new file with mode: 0644]
gpstk-example7.cpp [new file with mode: 0644]
gpstk-example8.cpp [new file with mode: 0644]

diff --git a/gpstk-example1.cpp b/gpstk-example1.cpp
new file mode 100644 (file)
index 0000000..b709d32
--- /dev/null
@@ -0,0 +1,28 @@
+#include <iostream>
+#include "DayTime.hpp"
+
+using namespace std;
+using namespace gpstk;
+
+int main(int argc, const char *argv[])
+{
+       try {
+
+               DayTime time;
+
+               cout << "Hello world!" << endl;
+               cout << "The current GPS week is " 
+                    << time.GPSfullweek() << endl;
+               cout << "The day of the GPS week is "
+                    << time.GPSday() << endl;
+               cout << "The seconds of the GPS week is "
+                    << time.GPSsecond() << endl;
+       }
+
+       catch (Exception error) {
+               cout << error << endl;
+               exit (-1);
+       }
+
+       return 0;
+}
diff --git a/gpstk-example2.cpp b/gpstk-example2.cpp
new file mode 100644 (file)
index 0000000..088ce89
--- /dev/null
@@ -0,0 +1,37 @@
+#include <iostream>
+#include <iomanip>
+
+#include "RinexObsBase.hpp"
+#include "RinexObsHeader.hpp"
+#include "RinexObsData.hpp"
+#include "RinexObsStream.hpp"
+
+using namespace std;
+using namespace gpstk;
+
+int main(void) {
+
+       // from : ftp://rgpdata.ign.fr/pub/data/2011/350/data_30/
+       // Create the input file stream
+       RinexObsStream rin("cena0010.03d");
+       
+       // Create the output file stream
+       RinexObsStream rout("cena0010.03d.new", ios::out|ios::trunc);
+       
+       // Read the RINEX header
+        RinexObsHeader head; //RINEX header object
+       rin >> head;
+
+       cout << "versionString : " << head.versionString << endl;
+
+       rout.header = rin.header;
+       rout << rout.header;
+       
+       // Loop over all data epochs
+        RinexObsData data; //RINEX data object
+       while (rin >> data) {
+               rout << data;
+       }
+       
+       return 0;
+}
diff --git a/gpstk-example7.cpp b/gpstk-example7.cpp
new file mode 100644 (file)
index 0000000..f3422d2
--- /dev/null
@@ -0,0 +1,957 @@
+#pragma ident "$Id: example7.cpp 1897 2009-05-13 09:04:32Z architest $"
+
+//============================================================================
+//
+//  This file is part of GPSTk, the GPS Toolkit.
+//
+//  The GPSTk is free software; you can redistribute it and/or modify
+//  it under the terms of the GNU Lesser General Public License as published
+//  by the Free Software Foundation; either version 2.1 of the License, or
+//  any later version.
+//
+//  The GPSTk is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU Lesser General Public License for more details.
+//
+//  You should have received a copy of the GNU Lesser General Public
+//  License along with GPSTk; if not, write to the Free Software Foundation,
+//  Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+//
+//  Copyright Dagoberto Salazar - gAGE ( http://www.gage.es ). 2007, 2008
+//
+//============================================================================
+
+// Example program Nro 7 for GPSTk
+// This program shows several different ways to process GPS data
+// using GNSS Data Structures (DataStructures.hpp).
+
+
+#include <iostream>
+#include <iomanip>
+
+   // Class for handling satellite observation parameters RINEX files
+#include "RinexObsStream.hpp"
+
+   // Classes for handling RINEX Broadcast ephemeris files
+#include "RinexNavStream.hpp"
+#include "RinexNavHeader.hpp"
+#include "RinexNavData.hpp"
+
+   // Class in charge of the GPS signal modelling
+#include "ModelObs.hpp"
+
+   // Class to store satellite broadcast navigation data
+#include "GPSEphemerisStore.hpp"
+
+   // Class to model the tropospheric delays
+#include "TropModel.hpp"
+
+   // Classes to model and store ionospheric delays
+#include "IonoModel.hpp"
+#include "IonoModelStore.hpp"
+
+   // Class to solve the equation system using Least Mean Squares
+#include "SolverLMS.hpp"
+
+   // Class to solve the equation system using Weighted-Least Mean Squares
+#include "SolverWMS.hpp"
+
+   // Class to solve equations systems using a simple code-based Kalman filter
+#include "CodeKalmanSolver.hpp"
+
+   // Class defining the GNSS data structures
+#include "DataStructures.hpp"
+
+   // Class to filter out observables grossly out of limits
+#include "SimpleFilter.hpp"
+
+   // Class for easily changing reference base from ECEF to NEU
+#include "XYZ2NEU.hpp"
+
+   // Class to detect cycle slips using just one frequency
+#include "OneFreqCSDetector.hpp"
+
+   // Class to detect cycle slips using LI combination
+#include "LICSDetector.hpp"
+
+   // Class to detect cycle slips using the Melbourne-Wubbena combination
+#include "MWCSDetector.hpp"
+
+   // Class to compute weights according to Appendix J of MOPS C (RTCA/DO-229C)
+#include "ComputeMOPSWeights.hpp"
+
+   // Class to smooth code observables (by default, C1)
+#include "CodeSmoother.hpp"
+
+   // Class to smooth the PC combination
+#include "PCSmoother.hpp"
+
+   // Classes to compute several combinations
+#include "ComputePC.hpp"
+#include "ComputeLC.hpp"
+#include "ComputeLI.hpp"
+#include "ComputeMelbourneWubbena.hpp"
+
+   // Class to compute single differences between receiver stations
+#include "DeltaOp.hpp"
+
+   // Class to synchronize two GNSS Data Structures data streams.
+#include "Synchronize.hpp"
+
+#include "geometry.hpp"                   // DEG_TO_RAD
+
+
+using namespace std;
+using namespace gpstk;
+
+int main(void)
+{
+
+      ////////// Initialization phase //////////
+
+      //////////// COMMON OBJECTS //////////////
+
+   cout << fixed << setprecision(3);   // Set a proper output format
+
+   RinexNavData rNavData;              // Object to store Rinex navigation data
+   GPSEphemerisStore bceStore;         // Object to store satellites ephemeris
+   RinexNavHeader rNavHeader;          // Object to read the header of Rinex
+                                       // navigation data files
+   IonoModelStore ionoStore;           // Object to store ionospheric models
+   IonoModel ioModel;                  // Declare a Ionospheric Model object
+
+      // Create the input observation file stream
+      // This is a fixed station, but here it will play as "rover"
+   RinexObsStream rin("ebre030a.02o");
+      // Please note that data was collected in year 2002, when the Sun
+      // was very active
+
+      // Create the input navigation file stream
+   RinexNavStream rnavin("brdc0300.02n");
+
+      // We need to read ionospheric parameters (Klobuchar model) from header
+   rnavin >> rNavHeader;
+
+      // Let's feed the ionospheric model (Klobuchar type) from data in the
+      // Navigation file header
+   ioModel.setModel(rNavHeader.ionAlpha, rNavHeader.ionBeta);
+
+   for (int i = 0; i < 4; i++) {
+          cout << "ionAlpha = " << rNavHeader.ionAlpha[i] << endl;
+          cout << "ionBeta  = " << rNavHeader.ionBeta[i]  << endl;
+   }
+
+
+      // Beware: 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 (rnavin >> rNavData)
+   {
+      bceStore.addEphemeris(rNavData);
+   }
+
+   bceStore.SearchPast();  // This is the default
+
+      // EBRE station nominal positional
+   Position nominalPos(4833520.3800, 41536.8300, 4147461.2800);
+
+      // Declare a MOPSTropModel object, setting the defaults
+   MOPSTropModel mopsTM( nominalPos.getAltitude(),
+                         nominalPos.getGeodeticLatitude(),
+                         30 );
+
+      // Declare the modeler object, setting all the parameters in one pass
+      // Given that in this example we are using a fixed GPS station with known
+      // coordinates, you could have used the "ModeledReferencePR" class, which
+      // is a little bit simpler.
+      // However, for a rover is more appropriate to use a "ModelObs" object
+      // because it allows to update the apriori position more easily (and it
+      // may automatically compute one, if needed, using Bancroft's method)
+   ModelObs model(nominalPos, ionoStore, mopsTM, bceStore, TypeID::C1);
+
+      // On the other hand, the usual way to use "ModelObs" is setting just the
+      // models in the constructor, and calling method "Prepare()" later, like
+      // in the following lines:
+      // ModelObs model(ionoStore, mopsTM, bceStore, TypeID::C1);
+      // model.Prepare(nominalPos);       // Set the reference position
+
+      // Declare a simple filter object. By default, it filters C1 with
+      // default limits
+   SimpleFilter myFilter;
+
+      // This is the GNSS data structure that will hold all the
+      // GNSS-related information
+   gnssRinex gOriginal;
+
+      ////////////////////////////////////////
+
+
+
+      //////////// CASE #1 OBJECTS ////////////
+
+      // Declare a SolverLMS object
+   SolverLMS solver;
+
+      ////////////////////////////////////////
+
+
+
+      //////////// CASE #2 OBJECTS ////////////
+
+      // Declare a base-changing object: From ECEF to North-East-Up (NEU)
+   XYZ2NEU baseChange(nominalPos);
+
+      // For some examples we need to reconfigure the solver in order
+      // to use a NEU system
+   TypeIDSet typeSet;
+   typeSet.insert(TypeID::dLat);
+   typeSet.insert(TypeID::dLon);
+   typeSet.insert(TypeID::dH);
+   typeSet.insert(TypeID::cdt);
+
+      // This is the proper equation structure to use with a NEU system
+   gnssEquationDefinition newEq(TypeID::prefitC, typeSet);
+
+      // Declare another SolverLMS object, but configure it to use a
+      // topocentric reference system (North-East-Up: NEU)
+   SolverLMS solverNEU;
+   solverNEU.setDefaultEqDefinition(newEq);    // NEU reconfiguration
+
+      ////////////////////////////////////////
+
+
+
+      //////////// CASE #3 OBJECTS ////////////
+
+      // This object will compute the appropriate MOPS weights
+   ComputeMOPSWeights mopsW(nominalPos, bceStore);
+
+      // Declare a solver object using Weighted-Least-Mean-Squares and
+      // a topocentric reference system (NEU)
+   SolverWMS solverWMS;
+   solverWMS.setDefaultEqDefinition(newEq);    // NEU reconfiguration
+
+      ////////////////////////////////////////
+
+
+
+      //////////// CASE #4 OBJECTS ////////////
+
+      // Let's declare a cycle slip detector using just one frequency
+   OneFreqCSDetector markCSC1;
+
+      // Declare an object to smooth code (C1 by default)
+   CodeSmoother smoothC1;
+   smoothC1.setMaxWindowSize(8);   // Configure smoother for 30 s sampling data
+
+      ////////////////////////////////////////
+
+
+
+      //////////// CASE #5 OBJECTS ////////////
+
+      // Let's declare another cycle slip detector using just one frequency
+   OneFreqCSDetector markCSC1case5;
+
+      // Declare another object to smooth code (C1 by default)
+   CodeSmoother smoothC1case5;
+   smoothC1case5.setMaxWindowSize(8);   // Configure for 30 s sampling data
+
+      // The core of this case is to add a new equation to the equation system
+      // Such equation states that there are NO changes in height for the
+      // rover.
+      //                          dH = 0
+      //
+      // Add a "fake" satellite to identify the new equation: Sat #1 of
+      // system "UserDefined"
+   SatID satEq(1,SatID::systemUserDefined);
+
+      // Declare and fill a "typeValueMap" object that will hold
+      // the equation data
+   typeValueMap equTVMap;
+   equTVMap[TypeID::prefitC] = 0.0;    // Code prefit residual is zero
+   equTVMap[TypeID::dLat]    = 0.0;    // Geometry matrix dLat coefficient is
+                                       // zero
+   equTVMap[TypeID::dLon]    = 0.0;    // Geometry matrix dLon coefficient is
+                                       // zero
+   equTVMap[TypeID::dH]      = 1.0;    // Geometry matrix dH coefficient is
+                                       // 1.0!!!
+   equTVMap[TypeID::cdt]     = 0.0;    // Geometry matrix cdt coefficient is
+                                       // zero
+
+      // Assign a relatively high weight to this information (typical
+      // MOPS weights range from 0.01 to 0.04)
+      // This means that this equation is very important for us, but it is
+      // NOT ABSOLUTELY TRUE. Some variation is allowed
+      // Given that weights are indeed (1/variances), if we assign to our new
+      // equation a confidence of 0.5 m of sigma, it means that we should use
+      // a weight of (1/(0.5^2)) = 4 m^(-2)
+   equTVMap[TypeID::weight]  = 4.0;
+
+      ////////////////////////////////////////
+
+
+
+      //////////// CASE #6 OBJECTS ////////////
+
+      // Object to compute the PC (ionosphere-free) combination.
+   ComputePC getPC;
+      // Use C1 instead of P1. P1 observables are declared in available RINEX
+      // files, but often they are indeed missing (like in this case). When
+      // that happens, this step is mandatory
+   getPC.useC1();
+
+      // Declare a simple filter object to screen PC
+   SimpleFilter pcFilter;
+   pcFilter.setFilteredType(TypeID::PC);
+
+      // Declare the modeler object for PC, setting all the parameters
+      // in one pass
+   ModelObs modelPC(nominalPos, mopsTM, bceStore, TypeID::PC, false);
+      // Take notice that PC combination doesn't use ionosphere modelling, nor
+      // TGD computation.
+      // WARNING: When using C1 instead of P1 to compute PC combination, be
+      // aware that instrumental errors will NOT cancel, introducing a bias
+      // that must be taken into account by other means. This will not work out
+      // in this example.
+
+      ////////////////////////////////////////
+
+
+
+      //////////// CASE #7 OBJECTS ////////////
+
+      // Objects to compute several common combinations.
+   ComputeLC getLC;
+   ComputeLI getLI;
+   ComputeMelbourneWubbena getMW;
+   getMW.useC1();      // Use C1 instead of P1
+
+      // Objects to mark cycle slips
+   LICSDetector markCSLI;      // Checks LI cycle slips
+   MWCSDetector markCSMW;      // Checks Merbourne-Wubbena cycle slips
+
+      // Object to smooth the PC combination. Defaults are usually fine
+   PCSmoother smoothPC;
+
+      ////////////////////////////////////////
+
+
+
+      //////////// CASE #8 OBJECTS ////////////
+
+      // Objects to mark cycle slips
+   LICSDetector markCSLIcase8;      // Checks LI cycle slips
+   MWCSDetector markCSMWcase8;      // Checks Merbourne-Wubbena cycle slips
+
+      // Object to smooth the PC combination. Defaults are usually fine
+   PCSmoother smoothPCcase8;
+
+      ////////////////////////////////////////
+
+
+
+      //////////// CASE #9 OBJECTS ////////////
+
+      // Objects to mark cycle slips
+   LICSDetector markCSLIcase9;      // Checks LI cycle slips
+   MWCSDetector markCSMWcase9;      // Checks Merbourne-Wubbena cycle slips
+
+      // Object to smooth the PC combination. Defaults are usually fine
+   PCSmoother smoothPCcase9;
+
+      // Declare a new Kalman solver, already reconfigured for NEU system
+   CodeKalmanSolver solverK9(newEq);
+
+      ////////////////////////////////////////
+
+
+
+      //////////// CASE #10 OBJECTS ////////////
+
+      // This is the GNSS data structure that will hold the
+      // reference station data
+   gnssRinex gRef;
+
+      // Create the input observation file stream for REFERENCE STATION
+   RinexObsStream rinRef("bell030a.02o");
+
+      // BELL reference station nominal position
+   Position nominalPosRef(4775849.6200, 116814.1000, 4213018.7100);
+
+      // Declare a MOPSTropModel object for the reference station, setting
+      // the defaults
+   MOPSTropModel mopsTMRef( nominalPosRef.getAltitude(),
+                            nominalPosRef.getGeodeticLatitude(),
+                            30 );
+
+      // Declare the appropriate modeler object for a reference station
+   ModelObsFixedStation modelRef( nominalPosRef,
+                                  ionoStore,
+                                  mopsTMRef,
+                                  bceStore,
+                                  TypeID::C1 );
+
+      // Create an object to compute the single differences of prefit residuals
+   DeltaOp delta;      // By default, it will work on code prefit residuals
+
+      // Create an object to synchronize rover and reference station
+      // data streams. This object will take data out from "rinRef" until
+      // it is synchronized with data in "gOriginal". Default synchronization
+      // tolerance is 1 s.
+   Synchronize synchro(rinRef, gOriginal);
+
+      //////////////////////////////////////////////
+
+
+
+      //////////// CASE #12 OBJECTS ////////////
+
+      // Declare a new Kalman solver, already reconfigured for NEU system
+   CodeKalmanSolver solverK12(newEq);
+
+      ////////////////////////////////////////
+
+      //////// End of initialization phase  ////////
+
+
+      //////// Processing phase ////////
+
+   int i = 0;
+
+      // Loop over all data epochs
+   while(rin >> gOriginal)
+   {
+
+          cout << "--- " << ++i << " ---" << endl;
+          cout << "  nominalPos\t: X=" << nominalPos.X()
+               << " Y=" << nominalPos.Y()
+               << " Z=" << nominalPos.Z() << endl;
+
+         // Let's output the time stamp (in seconds of day)
+          //cout << gOriginal.header.epoch.DOYsecond() << "  ";   // Output field #1
+
+
+   //////////////////////////// CASE #1  ////////////////////////////
+
+         // This case is a common C1 + Least Mean Squares solver
+         // (LMS) processing
+
+         // Let's make a working copy
+      gnssRinex gRin1(gOriginal);
+
+      try
+      {
+
+            // This is the line that will process all the GPS data
+             gRin1 >> myFilter >> model >> solver;
+
+            // - First, a basic filter to screen out very bad observables
+            // - Second, apply a model to the observables (ionosphere,
+            //   troposphere, relativity, etc.)
+            // - Third, solve the equations using a simple Least-Mean-Squares
+            //   solver
+      }
+      catch(...)
+      {
+         cerr << "Case 1. Exception at epoch: " << gRin1.header.epoch << endl;
+      }
+
+         // Get your results out of the solver object. In ECEF system
+         // by default
+      Position solPos( (model.rxPos.X() + solver.solution[0]),
+                       (model.rxPos.Y() + solver.solution[1]),
+                       (model.rxPos.Z() + solver.solution[2]) );
+
+      cout << "  Case 1\t: X=" << solPos.X()
+          << " Y=" << solPos.Y()
+          << " Z=" << solPos.Z() << endl;
+
+
+         // Let's change results to a North-East-Up (NEU) reference frame
+         // Compute the difference regarding the nominal 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 of case #1 in a topocentrical, North-East-Up
+         // reference frame
+         // Latitude change, output field #2
+      //cout << magnitude*cos(azimuth*DEG_TO_RAD)*cos(elev*DEG_TO_RAD) << "  ";
+        // Longitude change, output field #3
+      //cout << magnitude*sin(azimuth*DEG_TO_RAD)*cos(elev*DEG_TO_RAD) << "  ";
+         // Altitude change, output field #4
+      //cout << magnitude*sin(elev*DEG_TO_RAD) << "  ";
+
+   ////////////////////////// END OF CASE #1  //////////////////////////
+
+
+
+   //////////////////////////// CASE #2  ////////////////////////////
+
+         // This is exactly the same as CASE #1, but using a nice class
+         // to change the reference frame: ECEF -> NEU
+
+         // Let's make a working copy
+      gnssRinex gRin2(gOriginal);
+
+      cout << "  Case 2\t: ";
+      try
+      {
+
+         gRin2 >> myFilter >> model >> baseChange >> solverNEU;
+
+      }
+      catch(...)
+      {
+         cerr << "Case 2. Exception at epoch: " << gRin2.header.epoch << endl;
+      }
+
+      cout << solverNEU.solution(0) << "  ";  // dLat - Output field #5
+      cout << solverNEU.solution(1) << "  ";  // dLon - Output field #6
+      cout << solverNEU.solution(2) << "  ";  // dH   - Output field #7
+      cout << endl;
+
+
+         // Quite easier with respect to CASE #1, isn't it?  ;-)
+
+         // - "baseChange" object changes reference frame from ECEF to NEU
+         // - "solverNEU" is a simple Least-Mean-Squares solver, but
+         //    reconfigured to solve the dLat, dLon, dH, cdt (NEU) system
+         //    instead of the dx, dy, dz, cdt (ECEF) system
+         // - The other steps are exactly the same as case #1, and results
+         //   MUST match
+         // - If you want to see an even easier method to report the solution,
+         //   please see Case #3.
+
+         // By the way, if you want to inspect what is inside the body of a
+         // given GNSS data structure, you may write something like:
+         //
+         //      gRin2.body.dump(cout, 1);
+
+   ////////////////////////// END OF CASE #2  //////////////////////////
+
+
+
+   //////////////////////////// CASE #3  ////////////////////////////
+
+         // In this case we process data using C1 + Weighted Least Mean Squares
+         // solver (WMS)
+
+         // Let's make a working copy
+      gnssRinex gRin3(gOriginal);
+
+      cout << "  Case 3:\t ";
+      try
+      {
+
+         gRin3 >> myFilter >> model >> mopsW >> baseChange >> solverWMS;
+            // The "mopsW" object computes weights based on MOPS algorithm
+            // The "solverWMS" object solves the system using Weighted Least
+            // Mean Squares. It is already configured to work with NEU system.
+      }
+      catch(...)
+      {
+         cerr << "Case 3. Exception at epoch: " << gRin3.header.epoch << endl;
+      }
+
+         // An alternative way to report the solution is to access it
+         // using the TypeID's defined in the "gnssEquationDefinition" object
+         // assigned to the solver.
+         // With this method we avoid the possibility of getting the wrong
+         // type of solution from the "solution" vector.
+      cout << solverWMS.getSolution(TypeID::dLat) << "  ";  // dLat - Field #8
+      cout << solverWMS.getSolution(TypeID::dLon) << "  ";  // dLon - Field #9
+      cout << solverWMS.getSolution(TypeID::dH) << "  ";    // dH   - Field #10
+      cout << endl;
+
+   ////////////////////////// END OF CASE #3  //////////////////////////
+
+
+
+   //////////////////////////// CASE #4  ////////////////////////////
+
+         // This case does about the same as a modern GPS aircraft receiver,
+         // except for SBAS corrections and RAIM: C1smoothed + MOPS weights
+         // + WMS
+
+         // Let's make a working copy
+      gnssRinex gRin4(gOriginal);
+
+      cout << "  Case 4:\t ";
+      try
+      {
+
+         gRin4 >> myFilter >> markCSC1 >> smoothC1 >> model >> mopsW
+               >> baseChange >> solverWMS;
+            // The "markCSC1" object will try to detect cycle slips using just
+            // one frequency data (C1 and L1 observables), marking the CS flags
+            // Then, "smoothC1" will use the former information to smooth C1
+            // observations using phase data (L1)
+
+            // BEWARE: Both cycle slip detectors and "smoothers" are objects
+            // that store their internal state, so you MUST NOT use the SAME
+            // object to process DIFFERENT data streams
+      }
+      catch(...)
+      {
+         cerr << "Case 4. Exception at epoch: " << gRin4.header.epoch << endl;
+      }
+
+      cout << solverWMS.getSolution(TypeID::dLat) << "  ";  // dLat - Field #11
+      cout << solverWMS.getSolution(TypeID::dLon) << "  ";  // dLon - Field #12
+      cout << solverWMS.getSolution(TypeID::dH) << "  ";    // dH   - Field #13
+      cout <<endl;
+
+   ////////////////////////// END OF CASE #4  //////////////////////////
+
+
+
+   //////////////////////////// CASE #5  ////////////////////////////
+
+         // This case is like the former, but now let's suppose that one of
+         // the unknowns is indeed known: In this case dH is constant an equal
+         // to zero (i.e.: the "rover" doesn't change altitude), and we assign
+         //  a high "weight" to this information.
+
+         // Let's make a working copy
+      gnssRinex gRin5(gOriginal);
+
+      cout << "  Case 5:\t ";
+      try
+      {
+            // First, the typical processing up to the change of reference
+            // frame. Please note that all changes are stored in gRin5 GNSS
+            // data structure
+         gRin5 >> myFilter >> markCSC1case5 >> smoothC1case5 >> model
+               >> mopsW >> baseChange;
+            // Remember that both cycle slip detectors and "smoothers" are
+            // objects that store their internal state, so you MUST NOT use
+            // the SAME object to process DIFFERENT data streams (please
+            // compare with case #4)
+
+            // Now, let's insert the new equation data, including its
+            // corresponding weight
+         gRin5.body[satEq] = equTVMap;
+
+            // Let's continue processing data as usual
+         gRin5 >> solverWMS;
+
+      }
+      catch(...)
+      {
+         cerr << "Case 5. Exception at epoch: " << gRin5.header.epoch << endl;
+      }
+
+      cout << solverWMS.getSolution(TypeID::dLat) << "  ";  // dLat - Field #14
+      cout << solverWMS.getSolution(TypeID::dLon) << "  ";  // dLon - Field #15
+      cout << solverWMS.getSolution(TypeID::dH) << "  ";    // dH   - Field #16
+      cout << endl;
+
+   ////////////////////////// END OF CASE #5  //////////////////////////
+
+
+
+   //////////////////////////// CASE #6  ////////////////////////////
+
+         // This case uses de PC combination plus a WMS solver
+
+         // Let's make a working copy
+      gnssRinex gRin6(gOriginal);
+
+      cout << "  Case 6:\t ";
+      try
+      {
+
+         gRin6 >> getPC >> pcFilter >> modelPC >> mopsW
+               >> baseChange >> solverWMS;
+            // First, we need to compute the PC combination with "getPC" and
+            // insert it into the "gRin6" data structure.
+            // Then, use "pcFilter" to take out grossly out of range
+            // results in PC.
+            // After that, use an specific model ("modelPC") for this
+            // combination. It doesn't use ionospheric model nor TGD.
+            // The remaining steps are similar to the other cases.
+      }
+      catch(...)
+      {
+         cerr << "Case 6. Exception at epoch: " << gRin6.header.epoch << endl;
+      }
+
+      cout << solverWMS.getSolution(TypeID::dLat) << "  ";  // dLat - Field #17
+      cout << solverWMS.getSolution(TypeID::dLon) << "  ";  // dLon - Field #18
+      cout << solverWMS.getSolution(TypeID::dH) << "  ";    // dH   - Field #19
+      cout << endl;
+
+   ////////////////////////// END OF CASE #6  //////////////////////////
+
+
+
+   //////////////////////////// CASE #7  ////////////////////////////
+
+         // This case uses the smoothed-PC combination plus WMS
+
+         // Let's make a working copy
+      gnssRinex gRin7(gOriginal);
+
+      cout << "  Case 7:\t ";
+      try
+      {
+
+         gRin7 >> getPC >> getLC >> getLI >> getMW >> markCSLI >> markCSMW
+               >> smoothPC >> pcFilter >> modelPC >> mopsW >> baseChange
+               >> solverWMS;
+            // In addition to PC, we will also neet LC ("getLC"), LI ("getLI")
+            // and MW ("getMW") combinations:
+            //
+            // - LC (as well as PC) is needed by "smoothPC" in order to smooth
+            //   PC data. Also, the smoother works better with cycle slip
+            //   information, and therefore:
+            //
+            //   - LI feeds "markCSLI": The LI-based cycle slip detector
+            //   - MW feeds "markCSMW": The MW-based cycle slip detector
+            //
+            // - The remaining steps are essentially the same
+      }
+      catch(...)
+      {
+         cerr << "Case 7. Exception at epoch: " << gRin7.header.epoch << endl;
+      }
+
+      cout << solverWMS.getSolution(TypeID::dLat) << "  ";  // dLat - Field #20
+      cout << solverWMS.getSolution(TypeID::dLon) << "  ";  // dLon - Field #21
+      cout << solverWMS.getSolution(TypeID::dH) << "  ";    // dH   - Field #22
+      cout << endl;
+
+   ////////////////////////// END OF CASE #7  //////////////////////////
+
+
+
+   //////////////////////////// CASE #8  ////////////////////////////
+
+         // This case uses the smoothed-PC combination + WMS + information
+         // about dH (constant and equal to zero with a high confidence).
+         // It is a mix of the former case (#7) and case #5.
+
+         // Let's make a working copy
+      gnssRinex gRin8(gOriginal);
+
+      cout << "  Case 8:\t ";
+      try
+      {
+            // First, let's process data up to the change of reference frame
+         gRin8 >> getPC >> getLC >> getLI >> getMW >> markCSLIcase8
+               >> markCSMWcase8 >> smoothPCcase8 >> pcFilter >> modelPC
+               >> mopsW >> baseChange;
+            // Remember that both cycle slip detectors and "smoothers" are
+            // objects that store their internal state, so you MUST NOT use
+            // the SAME object to process DIFFERENT data streams (compare with
+            // case #7).
+
+            // Now, let's insert the new equation data, including its weight.
+            // It is the same equation as case #5.
+         gRin8.body[satEq] = equTVMap;
+
+            // Let's continue processing data as usual
+         gRin8 >> solverWMS;
+
+      }
+      catch(...)
+      {
+         cerr << "Case 8. Exception at epoch: " << gRin8.header.epoch << endl;
+      }
+
+      cout << solverWMS.getSolution(TypeID::dLat) << "  ";  // dLat - Field #23
+      cout << solverWMS.getSolution(TypeID::dLon) << "  ";  // dLon - Field #24
+      cout << solverWMS.getSolution(TypeID::dH) << "  ";    // dH   - Field #25
+      cout << endl;
+
+   ////////////////////////// END OF CASE #8  //////////////////////////
+
+
+
+   //////////////////////////// CASE #9  ////////////////////////////
+
+         // This case uses the smoothed-PC combination, exactly like case #7,
+         // but solves the equation system using a simple Kalman filter.
+
+         // Let's make a working copy
+      gnssRinex gRin9(gOriginal);
+
+      cout << "  Case 9:\t ";
+      try
+      {
+
+         gRin9 >> getPC >> getLC >> getLI >> getMW >> markCSLIcase9
+               >> markCSMWcase9 >> smoothPCcase9 >> pcFilter >> modelPC
+               >> mopsW >> baseChange >> solverK9;
+            // VERY IMPORTANT: Note that in this case the coordinates are
+            // handled as constants, whereas the receiver clock is modeled as
+            // white noise.
+      }
+      catch(...)
+      {
+         cerr << "Case 9. Exception at epoch: " << gRin9.header.epoch << endl;
+      }
+
+      cout << solverK9.getSolution(TypeID::dLat) << "  ";  // dLat - Field #26
+      cout << solverK9.getSolution(TypeID::dLon) << "  ";  // dLon - Field #27
+      cout << solverK9.getSolution(TypeID::dH) << "  ";    // dH   - Field #28
+      cout << endl;
+
+   ////////////////////////// END OF CASE #9  //////////////////////////
+
+
+
+   //////////////////////////// CASE #10  ////////////////////////////
+
+
+         // This is like cases #1 and #2, but using DGPS techniques instead.
+
+         // Let's make a working copy of rover data
+      gnssRinex gRin10(gOriginal);
+
+
+         // First, let's synchronize and process reference station data
+      cout << "  Case 10:\t ";
+      try
+      {
+         gRef >> synchro >> myFilter >> modelRef;
+            // Please note that the FIRST STEP is to synchronize "gRef", the
+            // reference station data stream, with "gOriginal" (or with gRin10,
+            // which is the same), the rover receiver data stream.
+            //
+            // Also, remember that in simple DGPS the differences are computed
+            // on code prefit residuals, so "modelRef" object is mandatory.
+
+            // The "delta" object will take care of proper differencing.
+            // We must tell it which GNSS data structure will be used
+            // as reference
+         delta.setRefData(gRef.body);
+
+      }
+      catch(SynchronizeException& e)   // THIS IS VERY IMPORTANT IN ORDER TO
+      {                                // MANAGE A POSSIBLE DESYNCHRONIZATION!!
+         cout << endl;
+         continue;
+      }
+      catch(...)
+      {
+         cerr << "Case 10. Exception when processing reference station data \
+at epoch: " << gRef.header.epoch << endl;
+      }
+
+
+         // Rover data processing is done here:
+      try
+      {
+
+         gRin10 >> myFilter >> model >> delta >> baseChange >> solverNEU;
+            // This is very similar to cases #1 and #2, but we insert a "delta"
+            // object that will adjust code prefit residuals BEFORE solving the
+            // system of equations.
+      }
+      catch(...)
+      {
+         cerr << "Case 10. Exception at epoch: " << gRin10.header.epoch
+              << endl;
+      }
+
+      cout << solverNEU.getSolution(TypeID::dLat) << "  ";  // dLat - Field #29
+      cout << solverNEU.getSolution(TypeID::dLon) << "  ";  // dLon - Field #30
+      cout << solverNEU.getSolution(TypeID::dH) << "  ";    // dH   - Field #31
+      cout << endl;
+
+   ////////////////////////// END OF CASE #10  //////////////////////////
+
+
+
+   //////////////////////////// CASE #11  ////////////////////////////
+
+
+         // This is like case #10 (DGPS), but now let's apply a WMS solver
+         // on data
+
+         // Let's make a working copy
+         gnssRinex gRin11(gOriginal);
+
+
+         // Please note that data streams are already synchronized, and
+         // "delta" object may be reused with the same reference data
+         // obtained from Case #10.
+
+      cout << "  Case 11:\t ";
+      try
+      {
+
+         gRin11 >> myFilter >> model >> delta >> mopsW >> baseChange
+                >> solverWMS;
+            // Like case #10, but now with "mopsW" and "solverWMS"
+      }
+      catch(...)
+      {
+         cerr << "Case 11. Exception at epoch: " << gRin11.header.epoch
+              << endl;
+      }
+
+      cout << solverWMS.getSolution(TypeID::dLat) << "  ";  // dLat - Field #32
+      cout << solverWMS.getSolution(TypeID::dLon) << "  ";  // dLon - Field #33
+      cout << solverWMS.getSolution(TypeID::dH) << "  ";    // dH   - Field #34
+      cout << endl;
+
+   ////////////////////////// END OF CASE #11  //////////////////////////
+
+
+
+   //////////////////////////// CASE #12  ////////////////////////////
+
+
+         // This is like case #11 (DGPS), but now let's apply a simple
+         // Kalman filter on data
+
+         // Let's make a working copy
+         gnssRinex gRin12(gOriginal);
+
+
+         // Please note that data streams are already synchronized, and
+         // "delta" object may be reused with the same reference data
+         // obtained from Case #10.
+
+      cout << "  Case 12:\t ";
+      try
+      {
+
+         gRin12 >> myFilter >> model >> delta >> mopsW >> baseChange
+                >> solverK12;
+            // Like case #11, but now with "solverK12"
+            // VERY IMPORTANT: Note that in this case the coordinates are
+            // handled as constants, whereas the receiver clock is modeled as
+            // white noise.
+      }
+      catch(...)
+      {
+         cerr << "Case 12. Exception at epoch: " << gRin12.header.epoch
+              << endl;
+      }
+
+      cout << solverK12.getSolution(TypeID::dLat) << "  ";  // dLat - Field #35
+      cout << solverK12.getSolution(TypeID::dLon) << "  ";  // dLon - Field #36
+      cout << solverK12.getSolution(TypeID::dH) << "  ";    // dH   - Field #37
+      cout << endl;
+
+   ////////////////////////// END OF CASE #12  //////////////////////////
+
+
+         // End of data processing for this epoch
+      cout << endl;
+
+   }
+
+   exit(0);
+
+}
diff --git a/gpstk-example8.cpp b/gpstk-example8.cpp
new file mode 100644 (file)
index 0000000..dfca675
--- /dev/null
@@ -0,0 +1,512 @@
+#pragma ident "$Id: example8.cpp 1897 2009-05-13 09:04:32Z architest $"
+
+//============================================================================
+//
+//  This file is part of GPSTk, the GPS Toolkit.
+//
+//  The GPSTk is free software; you can redistribute it and/or modify
+//  it under the terms of the GNU Lesser General Public License as published
+//  by the Free Software Foundation; either version 2.1 of the License, or
+//  any later version.
+//
+//  The GPSTk is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU Lesser General Public License for more details.
+//
+//  You should have received a copy of the GNU Lesser General Public
+//  License along with GPSTk; if not, write to the Free Software Foundation,
+//  Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+//
+//  Copyright Dagoberto Salazar - gAGE ( http://www.gage.es ). 2008, 2009
+//
+//============================================================================
+
+// Example program Nro 8 for GPSTk
+//
+// This program shows how to use GNSS Data Structures (GDS) to obtain
+// "Precise Point Positioning" (PPP).
+//
+// For details on the PPP algorithm please consult:
+//
+//    Kouba, J. and P. Heroux. "Precise Point Positioning using IGS Orbit
+//       and Clock Products". GPS Solutions, vol 5, pp 2-28. October, 2001.
+//
+
+
+#include <iostream>
+#include <iomanip>
+
+
+   // Class for handling satellite observation parameters RINEX files
+#include "RinexObsStream.hpp"
+
+   // Class to store satellite precise navigation data
+#include "SP3EphemerisStore.hpp"
+
+   // Class in charge of basic GNSS signal modelling
+#include "BasicModel.hpp"
+
+   // Class to model the tropospheric delays
+#include "TropModel.hpp"
+
+   // Class defining the GNSS data structures
+#include "DataStructures.hpp"
+
+   // Class to filter out satellites without required observables
+#include "RequireObservables.hpp"
+
+   // Class to filter out observables grossly out of limits
+#include "SimpleFilter.hpp"
+
+   // Class for easily changing reference base from ECEF to NEU
+#include "XYZ2NEU.hpp"
+
+   // Class to detect cycle slips using LI combination
+#include "LICSDetector2.hpp"
+
+   // Class to detect cycle slips using the Melbourne-Wubbena combination
+#include "MWCSDetector.hpp"
+
+   // Class to compute the effect of solid tides
+#include "SolidTides.hpp"
+
+   // Class to compute the effect of ocean loading
+#include "OceanLoading.hpp"
+
+   // Class to compute the effect of pole tides
+#include "PoleTides.hpp"
+
+   // Class to correct observables
+#include "CorrectObservables.hpp"
+
+   // Class to compute the effect of wind-up
+#include "ComputeWindUp.hpp"
+
+   // Class to compute the effect of satellite antenna phase center
+#include "ComputeSatPCenter.hpp"
+
+   // Class to compute the tropospheric data
+#include "ComputeTropModel.hpp"
+
+   // Class to compute linear combinations
+#include "ComputeLinear.hpp"
+
+   // This class pre-defines several handy linear combinations
+#include "LinearCombinations.hpp"
+
+   // Class to compute Dilution Of Precision values
+#include "ComputeDOP.hpp"
+
+   // Class to keep track of satellite arcs
+#include "SatArcMarker.hpp"
+
+   // Class to compute gravitational delays
+#include "GravitationalDelay.hpp"
+
+   // Class to align phases with code measurements
+#include "PhaseCodeAlignment.hpp"
+
+   // Compute statistical data
+#include "PowerSum.hpp"
+
+   // Used to delete satellites in eclipse
+#include "EclipsedSatFilter.hpp"
+
+   // Used to decimate data. This is important because RINEX observation
+   // data is provided with a 30 s sample rate, whereas SP3 files provide
+   // satellite clock information with a 900 s sample rate.
+#include "Decimate.hpp"
+
+   // Class to compute the Precise Point Positioning (PPP) solution
+#include "SolverPPP.hpp"
+
+#include "geometry.hpp"                   // DEG_TO_RAD
+
+
+using namespace std;
+using namespace gpstk;
+
+
+int main(void)
+{
+
+      /////////////////// INITIALIZATION PART /////////////////////
+
+   cout << fixed << setprecision(3);   // Set a proper output format
+
+
+      // Create the input observation file stream
+   RinexObsStream rin("onsa2240.05o");
+
+
+      // Declare a "SP3EphemerisStore" object to handle precise ephemeris
+   SP3EphemerisStore SP3EphList;
+
+      // Set flags to reject satellites with bad or absent positional
+      // values or clocks
+   SP3EphList.rejectBadPositions(true);
+   SP3EphList.rejectBadClocks(true);
+
+      // Set flags to check for data gaps and too wide interpolation intervals.
+      // Default values for "gapInterval" (901.0 s) and "maxInterval"
+      // (8105.0 s) will be used.
+   SP3EphList.enableDataGapCheck();
+   SP3EphList.enableIntervalCheck();
+
+      // Load all the SP3 ephemerides files
+   SP3EphList.loadFile("igs13354.sp3");
+   SP3EphList.loadFile("igs13355.sp3");
+   SP3EphList.loadFile("igs13356.sp3");
+
+
+      // ONSA station nominal position
+   Position nominalPos(3370658.5419, 711877.1496, 5349786.9542);
+
+
+      // Declare a NeillTropModel object, setting the defaults
+   NeillTropModel neillTM( nominalPos.getAltitude(),
+                           nominalPos.getGeodeticLatitude(), 224);
+
+
+      // This is the GNSS data structure that will hold all the
+      // GNSS-related information
+   gnssRinex gRin;
+
+
+      // Declare a base-changing object: From ECEF to North-East-Up (NEU)
+   XYZ2NEU baseChange(nominalPos);
+
+
+      // Declare an object to check that all required observables are present
+   RequireObservables requireObs(TypeID::P1);
+   requireObs.addRequiredType(TypeID::P2);
+   requireObs.addRequiredType(TypeID::L1);
+   requireObs.addRequiredType(TypeID::L2);
+
+
+      // Declare a simple filter object to screen PC
+   SimpleFilter pcFilter;
+   pcFilter.setFilteredType(TypeID::PC);
+
+
+      // Declare a basic modeler
+   BasicModel basic(nominalPos, SP3EphList);
+
+
+      // Objects to mark cycle slips
+   LICSDetector2 markCSLI;     // Checks LI cycle slips
+   MWCSDetector markCSMW;      // Checks Merbourne-Wubbena cycle slips
+
+
+      // Object to compute tidal effects
+   SolidTides  solid;
+
+      // Ocean loading model
+   OceanLoading ocean("OCEAN-GOT00.dat");
+
+      // Numerical values are x,y pole displacements for Aug/12/2005 (arcsec).
+   PoleTides   pole(0.02094, 0.42728);
+
+
+      // Vector from ONSA antenna ARP to L1 phase center [UEN] (AOAD/M_B)
+   Triple offsetL1(0.0780, 0.000, 0.000);   // Units in meters
+
+      // Vector from ONSA antenna ARP to L2 phase center [UEN] (AOAD/M_B)
+   Triple offsetL2(0.096, 0.0000, 0.000);    // Units in meters
+
+      // Vector from monument to antenna ARP [UEN] for ONSA station
+   Triple offsetARP(0.9950, 0.0, 0.0);    // Units in meters
+
+
+      // Declare an object to correct observables to monument
+   CorrectObservables corr(SP3EphList);
+   ((corr.setNominalPosition(nominalPos)).setL1pc(offsetL1)).setL2pc(offsetL2);
+   corr.setMonument(offsetARP);
+
+
+      // Object to compute wind-up effect
+   ComputeWindUp windup(SP3EphList, nominalPos, "PRN_GPS");
+
+
+      // Object to compute satellite antenna phase center effect
+   ComputeSatPCenter svPcenter(nominalPos);
+
+
+      // Object to compute the tropospheric data
+   ComputeTropModel computeTropo(neillTM);
+
+
+      // This object defines several handy linear combinations
+   LinearCombinations comb;
+
+
+      // Object to compute linear combinations for cycle slip detection
+   ComputeLinear linear1(comb.pdeltaCombination);
+   linear1.addLinear(comb.ldeltaCombination);
+   linear1.addLinear(comb.mwubbenaCombination);
+   linear1.addLinear(comb.liCombination);
+
+
+      // This object computes the ionosphere-free combinations to be used
+      // as observables in the PPP processing
+   ComputeLinear linear2(comb.pcCombination);
+   linear2.addLinear(comb.lcCombination);
+
+
+      // Let's use a different object to compute prefit residuals
+   ComputeLinear linear3(comb.pcPrefit);
+   linear3.addLinear(comb.lcPrefit);
+
+
+      // Declare an object to process the data using PPP. It is set
+      // to use a NEU system
+   SolverPPP pppSolver(true);
+
+      // The real test for a PPP processing program is to handle coordinates
+      // as white noise. In such case, position error should be about 0.25 m or
+      // better. Uncomment the following couple of lines to test this.
+//   WhiteNoiseModel wnM(100.0);            // 100 m of sigma
+//   pppSolver.setCoordinatesModel(&wnM);
+
+
+      // Object to keep track of satellite arcs
+   SatArcMarker markArc;
+   markArc.setDeleteUnstableSats(true);
+   markArc.setUnstablePeriod(151.0);
+
+      // Object to compute gravitational delay effects
+   GravitationalDelay grDelay(nominalPos);
+
+      // Object to align phase with code measurements
+   PhaseCodeAlignment phaseAlign;
+
+      // Object to compute DOP values
+   ComputeDOP cDOP;
+
+      // Object to remove eclipsed satellites
+   EclipsedSatFilter eclipsedSV;
+
+      // Object to decimate data. This is important because RINEX observation
+      // data is provided with a 30 s sample rate, whereas SP3 files provide
+      // satellite clock information with a 900 s sample rate.
+   Decimate decimateData(900.0, 5.0, SP3EphList.getInitialTime());
+
+      // When you are printing the model, you may want to comment the previous
+      // line and uncomment the following one, generating a 30 s model
+//   Decimate decimateData(30.0, 1.0, SP3EphList.getInitialTime());
+
+      // Statistical summary objects
+   PowerSum errorVectorStats;
+
+
+
+      /////////////////// PROCESING PART /////////////////////
+
+
+      // Use this variable to select between position printing or model printing
+   bool printPosition(true);     // By default, print position and associated
+                                 // parameters
+
+
+
+      // Loop over all data epochs
+   while(rin >> gRin)
+   {
+
+      DayTime time(gRin.header.epoch);
+
+         // Compute the effect of solid, oceanic and pole tides
+      Triple tides( solid.getSolidTide(time, nominalPos) +
+                    ocean.getOceanLoading("ONSA", time)  +
+                    pole.getPoleTide(time, nominalPos)     );
+
+         // Update observable correction object with tides information
+      corr.setExtraBiases(tides);
+
+      try
+      {
+
+            // The following lines are indeed just one line
+         gRin >> requireObs      // Check if required observations are present
+              >> linear1         // Compute linear combinations to detect CS
+              >> markCSLI        // Mark cycle slips: LI algorithm
+              >> markCSMW        // Mark cycle slips: Melbourne-Wubbena
+              >> markArc         // Keep track of satellite arcs
+              >> decimateData    // If not a multiple of 900 s, decimate
+              >> basic           // Compute the basic components of model
+              >> eclipsedSV      // Remove satellites in eclipse
+              >> grDelay         // Compute gravitational delay
+              >> svPcenter       // Compute the effect of satellite phase center
+              >> corr            // Correct observables from tides, etc.
+              >> windup          // Compute wind-up effect
+              >> computeTropo    // Compute tropospheric effect
+              >> linear2         // Compute ionosphere-free combinations
+              >> pcFilter        // Filter out spurious data
+              >> phaseAlign      // Align phases with codes
+              >> linear3         // Compute prefit residuals
+              >> baseChange      // Prepare to use North-East-UP reference frame
+              >> cDOP            // Compute DOP figures
+              >> pppSolver;      // Solve equations with a Kalman filter
+
+      }
+      catch(DecimateEpoch& d)
+      {
+            // If 'decimateData' object detects that this epoch is not a
+            // multiple of 900 s, it issues a "DecimateEpoch" exception. Here
+            // we catch such exception, and just continue to process next epoch.
+         continue;
+      }
+      catch(Exception& e)
+      {
+         cerr << "Exception at epoch: " << time << "; " << e << endl;
+         continue;
+      }
+      catch(...)
+      {
+         cerr << "Unknown exception at epoch: " << time << endl;
+         continue;
+      }
+
+
+         // Check if we want to print model or position
+      if(printPosition)
+      {
+            // Print here the position results
+         cout << time.DOYsecond()      << "  ";     // Epoch - Output field #1
+
+         cout << pppSolver.getSolution(TypeID::dLat) << "  ";    // dLat  - #2
+         cout << pppSolver.getSolution(TypeID::dLon) << "  ";    // dLon  - #3
+         cout << pppSolver.getSolution(TypeID::dH) << "  ";      // dH    - #4
+         cout << pppSolver.getSolution(TypeID::wetMap) << "  ";  // Tropo - #5
+
+         cout << pppSolver.getVariance(TypeID::dLat) << "  "; // Cov dLat - #6
+         cout << pppSolver.getVariance(TypeID::dLon) << "  "; // Cov dLon - #7
+         cout << pppSolver.getVariance(TypeID::dH) << "  ";   // Cov dH   - #8
+         cout << pppSolver.getVariance(TypeID::wetMap) << "  ";//Cov Tropo- #9
+
+         cout << gRin.numSats()        << "  ";     // Visible satellites - #10
+
+         cout << cDOP.getGDOP()        << "  ";                   // GDOP - #11
+         cout << cDOP.getPDOP()        << "  ";                   // PDOP - #12
+         cout << cDOP.getTDOP()        << "  ";                   // TDOP - #13
+         cout << cDOP.getHDOP()        << "  ";                   // HDOP - #14
+         cout << cDOP.getVDOP()        << "  ";                   // VDOP - #15
+
+         cout << endl;
+
+            // For statistical purposes we discard the first two hours of data
+         if (time.DOYsecond() > 7200.0)
+         {
+               // Statistical summary
+            double errorV( pppSolver.solution[1]*pppSolver.solution[1] +
+                           pppSolver.solution[2]*pppSolver.solution[2] +
+                           pppSolver.solution[3]*pppSolver.solution[3] );
+
+               // Get module of position error vector
+            errorV = std::sqrt(errorV);
+
+               // Add to statistical summary object
+            errorVectorStats.add(errorV);
+         }
+
+      }  // End of position printing
+      else
+      {
+            // Print here the model results
+            // First, define types we want to keep
+         TypeIDSet types;
+         types.insert(TypeID::L1);
+         types.insert(TypeID::L2);
+         types.insert(TypeID::P1);
+         types.insert(TypeID::P2);
+         types.insert(TypeID::PC);
+         types.insert(TypeID::LC);
+         types.insert(TypeID::rho);
+         types.insert(TypeID::dtSat);
+         types.insert(TypeID::rel);
+         types.insert(TypeID::gravDelay);
+         types.insert(TypeID::tropo);
+         types.insert(TypeID::dryTropo);
+         types.insert(TypeID::dryMap);
+         types.insert(TypeID::wetTropo);
+         types.insert(TypeID::wetMap);
+         types.insert(TypeID::tropoSlant);
+         types.insert(TypeID::windUp);
+         types.insert(TypeID::satPCenter);
+         types.insert(TypeID::satX);
+         types.insert(TypeID::satY);
+         types.insert(TypeID::satZ);
+         types.insert(TypeID::elevation);
+         types.insert(TypeID::azimuth);
+         types.insert(TypeID::satArc);
+         types.insert(TypeID::prefitC);
+         types.insert(TypeID::prefitL);
+         types.insert(TypeID::dx);
+         types.insert(TypeID::dy);
+         types.insert(TypeID::dz);
+         types.insert(TypeID::dLat);
+         types.insert(TypeID::dLon);
+         types.insert(TypeID::dH);
+         types.insert(TypeID::cdt);
+
+         gRin.keepOnlyTypeID(types);   // Delete the types not in 'types'
+
+            // Iterate through the GNSS Data Structure
+         satTypeValueMap::const_iterator it;
+         for (it = gRin.body.begin(); it!= gRin.body.end(); it++) 
+         {
+
+               // Print epoch
+            cout << time.year()        << " ";
+            cout << time.DOY()         << " ";
+            cout << time.DOYsecond()   << " ";
+
+            cout << cDOP.getGDOP()        << "  ";  // GDOP #4
+            cout << cDOP.getPDOP()        << "  ";  // PDOP #5
+            cout << cDOP.getTDOP()        << "  ";  // TDOP #6
+            cout << cDOP.getHDOP()        << "  ";  // HDOP #7
+            cout << cDOP.getVDOP()        << "  ";  // VDOP #8
+
+               // Print satellite information (system and PRN)
+            cout << (*it).first << " ";
+
+               // Print model values
+            typeValueMap::const_iterator itObs;
+            for( itObs  = (*it).second.begin(); 
+                 itObs != (*it).second.end();
+                 itObs++ )
+            {
+               bool printNames(true);  // Whether print types' names or not
+               if (printNames)
+               {
+                  cout << (*itObs).first << " ";
+               }
+
+               cout << (*itObs).second << " ";
+
+            }  // End for( itObs = ... )
+
+            cout << endl;
+
+         }  // End for (it = gRin.body.begin(); ... )
+
+      }  // End of model printing
+
+   }  // End of 'while(rin >> gRin)...'
+
+
+
+      // Print statistical summary in cerr
+   if(printPosition)
+   {
+      cerr << "Module of error vector: Average = "
+           << errorVectorStats.average() << " m    Std. dev. = "
+           << std::sqrt(errorVectorStats.variance()) << " m" << endl;
+   }
+
+
+
+   exit(0);       // End of program
+
+}