Backup
[gps] / gpstk-solution6.cpp
index 96564ec..9180f88 100644 (file)
@@ -6,73 +6,62 @@
 #include <iostream>
 #include <iomanip>
 
-   // Class for handling satellite observation parameters RINEX files
+// Class for handling satellite observation parameters RINEX files
 #include "RinexObsStream.hpp"
 
-   // Classes for handling RINEX Broadcast ephemeris files
+// Classes for handling RINEX Broadcast ephemeris files
 #include "RinexNavStream.hpp"
 #include "RinexNavHeader.hpp"
 #include "RinexNavData.hpp"
 
-   // Class in charge of the GPS signal modelling
+// Class in charge of the GPS signal modelling
 #include "ModelObs.hpp"
 
-   // Class to store satellite broadcast navigation data
+// Class to store satellite broadcast navigation data
 #include "GPSEphemerisStore.hpp"
+#include "SP3EphemerisStore.hpp"
 
-   // Class to model the tropospheric delays
+// Class to model the tropospheric delays
 #include "TropModel.hpp"
 
-   // Classes to model and store ionospheric delays
+// Classes to model and store ionospheric delays
 #include "IonoModel.hpp"
 #include "IonoModelStore.hpp"
 
-   // Class to solve the equation system using Least Mean Squares
+// Class to solve the equation system using Least Mean Squares
 #include "SolverLMS.hpp"
 
-   // Class to solve the equation system using Weighted-Least Mean Squares
+// 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
+// Class to solve equations systems using a simple code-based Kalman filter
 #include "CodeKalmanSolver.hpp"
 
-   // Class defining the GNSS data structures
+// Class defining the GNSS data structures
 #include "DataStructures.hpp"
 
-   // Class to filter out observables grossly out of limits
+// Class to filter out observables grossly out of limits
 #include "SimpleFilter.hpp"
 
-   // Class for easily changing reference base from ECEF to NEU
+// Class for easily changing reference base from ECEF to NEU
 #include "XYZ2NEU.hpp"
 
-   // Class to detect cycle slips using just one frequency
+// Class to detect cycle slips using just one frequency
 #include "OneFreqCSDetector.hpp"
 
-   // Class to detect cycle slips using LI combination
+// Class to detect cycle slips using LI combination
 #include "LICSDetector.hpp"
 
-   // Class to detect cycle slips using the Melbourne-Wubbena combination
+// 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)
+// 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
+// Class to compute single differences between receiver stations
 #include "DeltaOp.hpp"
 
-   // Class to synchronize two GNSS Data Structures data streams.
+// Class to synchronize two GNSS Data Structures data streams.
 #include "Synchronize.hpp"
 
 #include "geometry.hpp"                   // DEG_TO_RAD
 using namespace std;
 using namespace gpstk;
 
+// class for weight 
+
+class ComputeElevationWeights : public WeightBase, public ProcessingClass
+{
+public:
+       ComputeElevationWeights() {};
+
+       double getWeight(const SatID& sat,
+                        typeValueMap& tvMap) {
+               double elevation = tvMap(TypeID::elevation);
+               if (elevation < 45.0)
+                       return elevation;
+               return 90.0;
+       };
+
+       virtual satTypeValueMap& Process(const DayTime& time,
+                                        satTypeValueMap& gData) {
+               
+               // Loop through all the satellites
+               satTypeValueMap::iterator it;
+               for (it = gData.begin(); it != gData.end(); ++it ) {
+                       (*it).second[TypeID::weight] =
+                               getWeight((*it).first, (*it).second);
+               }
+
+               return gData;
+       };
+
+       virtual gnssSatTypeValue& Process(gnssSatTypeValue& gData) {
+               Process(gData.header.epoch, gData.body); return gData;
+       };
+
+       virtual gnssRinex& Process(gnssRinex& gData) {
+               Process(gData.header.epoch, gData.body); return gData;
+       };
+
+       virtual int getIndex(void) const {
+               return 6; /* FIXME */
+       };
+
+       virtual string getClassName(void) const {
+               return "ComputeElevationWeights";
+       };
+};
+
+void display_rinex(const gnssRinex& gRin) {
+       cout << gRin.header.epoch
+            << " " << gRin.header.epoch.second() << endl;
+       gRin.body.dump(cout, 1);
+}
+
 void display_position(const Position& p,
                      const Position& nominalPos) {
        Position pos = p;
@@ -116,11 +156,18 @@ gnssRinex keepOnlySatGPS(const gnssRinex& rin) {
        return result;
 }
 
-int main(void)
+int main(int argc, const char *argv[])
 {
-       const char obsReference[] = "crei354a.11o";
-       const char obsMobile[] = "mlvl354a.11o";
-       const char navReference[] = "crei354a.11n";
+       // first set
+       //const char obsReference[] = "crei354a.11o";
+       //const char obsMobile[] = "mlvl354a.11o";
+       //const char navReference[] = "crei354a.11n";
+
+       // second set : from rgpdata.ign.fr:/pub/data/2011/361/data_1
+       const char obsReference[] = "crei361o.11o";
+       const char obsMobile[] = "4.ooo";
+       const char navReference[] = "crei361o.11n";
+       const char sp3File[] = "igu16682_12.sp3";
        
        // reference station nominal position
        Position nominalPosRef  (4166409.705,   182770.997,     4809802.953);
@@ -128,107 +175,79 @@ int main(void)
        Position nominalPos_true(4201577.209,   189859.856,     4779064.567);   
        Position nominalPos     (4201000.000,   189000.000,     4770000.000);   
 
-       Position lastPosition; /* only valid if i > 1 */
+       // initial set to nominalPos
+       Position lastPosition = nominalPos;
 
-      ////////// Initialization phase //////////
+       ////////// Initialization phase //////////
 
-      //////////// COMMON OBJECTS //////////////
+       //////////// COMMON OBJECTS //////////////
 
-   RinexNavData rNavData;              // Object to store Rinex navigation data
-   GPSEphemerisStore bceStore;         // Object to store satellites ephemeris
-   RinexNavHeader rNavHeader;          // Object to read the header of Rinex
+       RinexNavData rNavData;         // Object to store Rinex navigation data
+       GPSEphemerisStore bceStore;    // Object to store satellites ephemeris
+       SP3EphemerisStore sp3Store;
+       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
-
-   cout << "reading navigation file : " << navReference << endl;
+       IonoModelStore ionoStore;      // Object to store ionospheric models
+       IonoModel ioModel;             // Declare a Ionospheric Model object
 
-      // Create the input navigation file stream
-   RinexNavStream rnavin(navReference);
+       cout << "reading navigation file :\t" << navReference << endl;
 
-      // We need to read ionospheric parameters (Klobuchar model) from header
-   rnavin >> rNavHeader;
+       // Create the input navigation file stream
+       RinexNavStream rnavin(navReference);
 
-      // Let's feed the ionospheric model (Klobuchar type) from data in the
-      // Navigation file header
-   ioModel.setModel(rNavHeader.ionAlpha, rNavHeader.ionBeta);
+       // We need to read ionospheric parameters (Klobuchar model) from header
+       rnavin >> rNavHeader;
 
-   for (int i = 0; i < 4; i++)
-          cout << "ionAlpha = " << rNavHeader.ionAlpha[i] << endl;
-   for (int i = 0; i < 4; i++)
-          cout << "ionBeta  = " << rNavHeader.ionBeta[i]  << endl;
+       // Let's feed the ionospheric model (Klobuchar type) from data in the
+       // Navigation file header
+       ioModel.setModel(rNavHeader.ionAlpha, rNavHeader.ionBeta);
 
-      // Beware: In this case, the same model will be used for the
-      // full data span
-   ionoStore.addIonoModel(DayTime::BEGINNING_OF_TIME, ioModel);
+       for (int i = 0; i < 4; i++)
+               cout << "ionAlpha = " << rNavHeader.ionAlpha[i] << endl;
+       for (int i = 0; i < 4; i++)
+               cout << "ionBeta  = " << rNavHeader.ionBeta[i]  << endl;
 
-      // Storing the ephemeris in "bceStore"
-   while (rnavin >> rNavData)
-   {
-      bceStore.addEphemeris(rNavData);
-   }
+       // Beware: In this case, the same model will be used for the
+       // full data span
+       ionoStore.addIonoModel(DayTime::BEGINNING_OF_TIME, ioModel);
 
-   bceStore.SearchPast();  // This is the default
+       // Storing the ephemeris in "bceStore"
+       while (rnavin >> rNavData) {
+               bceStore.addEphemeris(rNavData);
+       }
 
-      // Declare a MOPSTropModel object, setting the defaults
-   MOPSTropModel mopsTM( nominalPos_true.getAltitude(),
-                         nominalPos_true.getGeodeticLatitude(),
-                         30 );
+       sp3Store.loadFile(sp3File);
 
-      // Declare a simple filter object. By default, it filters C1 with
-      // default limits
-   SimpleFilter filter;
+       bceStore.SearchPast();  // This is the default
 
-      // This is the GNSS data structure that will hold all the
-      // GNSS-related information
-   gnssRinex gOriginal;
+       // Declare a MOPSTropModel object, setting the defaults
+       MOPSTropModel mopsTM;
 
-      ////////////////////////////////////////
+       // Declare a simple filter object. By default, it filters C1 with
+       // default limits
+       SimpleFilter filter;
 
+       // This is the GNSS data structure that will hold all the
+       // GNSS-related information
+       gnssRinex gOriginal;
 
+       //////////// CASE #1 OBJECTS ////////////
 
-      //////////// CASE #1 OBJECTS ////////////
+       // Declare a SolverLMS object
+       SolverLMS solver;
 
-      // Declare a SolverLMS object
-   SolverLMS solver;
+       //////////// CASE #3 OBJECTS ////////////
 
-      ////////////////////////////////////////
-
-      //////////// CASE #3 OBJECTS ////////////
-
-      // This object will compute the appropriate MOPS weights
-   ComputeMOPSWeights mopsW(nominalPos, bceStore);
+       // This object will compute the appropriate MOPS weights
+       ComputeMOPSWeights mopsW;
    
-   // Declare a solver object using Weighted-Least-Mean-Squares and
-   // a topocentric reference system (NEU)
-   SolverWMS solverWMS;
-
-      ////////////////////////////////////////
-
-
-
-      //////////// 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
-
-      ////////////////////////////////////////
-
-
+       // Declare a solver object using Weighted-Least-Mean-Squares and
+       // a topocentric reference system (NEU)
+       SolverWMS solverWMS;
 
       //////////// 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.
@@ -263,75 +282,8 @@ int main(void)
       ////////////////////////////////////////
 
 
-
-      //////////// 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
    CodeKalmanSolver solverK9;
 
@@ -345,14 +297,13 @@ int main(void)
       // reference station data
    gnssRinex gRef;
 
-      // Create the input observation file stream for REFERENCE STATION
+   // Create the input observation file stream for REFERENCE STATION
    RinexObsStream rinRef(obsReference);
+   cout << "reading observation file (reference) :\t" << obsReference << endl;
 
       // Declare a MOPSTropModel object for the reference station, setting
       // the defaults
-   MOPSTropModel mopsTMRef( nominalPosRef.getAltitude(),
-                            nominalPosRef.getGeodeticLatitude(),
-                            30 );
+   MOPSTropModel mopsTMRef;
 
       // 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
@@ -366,7 +317,10 @@ int main(void)
    //            bceStore, TypeID::C1);
 
    ModelObs model;
-   model.setDefaultEphemeris(bceStore);
+   model.setDefaultIonoModel(ionoStore);
+   model.setDefaultTropoModel(mopsTM);
+   model.setDefaultEphemeris(sp3Store /*bceStore*/);
+   model.setMinElev(1.0);
 
       // Declare the appropriate modeler object for a reference station
    //ModelObsFixedStation modelRef(nominalPosRef,
@@ -374,9 +328,15 @@ int main(void)
    //                           bceStore, TypeID::C1);
 
    ModelObs modelRef;
-   modelRef.setDefaultEphemeris(bceStore);
+   modelRef.setDefaultIonoModel(ionoStore);
+   modelRef.setDefaultTropoModel(mopsTMRef);
+   modelRef.setDefaultEphemeris(sp3Store /*bceStore*/);
    modelRef.Prepare(nominalPosRef);
 
+   // use simple model without iono & tropo corrections
+   ModelObs myModelRef;
+   ModelObs myModel;
+   
       // Create an object to compute the single differences of prefit residuals
    DeltaOp delta;      // By default, it will work on code prefit residuals
 
@@ -395,6 +355,8 @@ int main(void)
       // Declare a new Kalman solver
    CodeKalmanSolver solverK12;
 
+   ComputeElevationWeights elevationW;
+
       ////////////////////////////////////////
 
       //////// End of initialization phase  ////////
@@ -412,15 +374,21 @@ int main(void)
       // Create the input observation file stream
       // This is a fixed station, but here it will play as "rover"
    RinexObsStream rin(obsMobile);
-      // Please note that data was collected in year 2002, when the Sun
-      // was very active
 
       // Loop over all data epochs
    while (rin >> gOriginal) {
 
           // remove glonass satellite
           gOriginal = keepOnlySatGPS(gOriginal);
-
+          gOriginal.keepOnlyTypeID(TypeID::C1);
+
+          // shift pseudorange by 37900000 m
+          /*
+          for (satTypeValueMap::iterator it = gOriginal.body.begin();
+               it != gOriginal.body.end(); it++) {
+                  (*it).second[TypeID::C1] -= 37900000.0;
+          }
+          */
           cout << "--- " << ++i << " ---" << endl;
           cout << "  Ref.   :";
           display_position(nominalPosRef, nominalPos_true);
@@ -431,9 +399,49 @@ int main(void)
           cout << "  Last   :";
           display_position(lastPosition, nominalPos_true);
 
-          // prepare the model using the last position (if available)
-          if (i > 1)
-                  model.Prepare(lastPosition);
+          cout << "gOriginal : " << endl;
+          display_rinex(gOriginal);
+
+          // update MOPSTropModel
+          mopsTMRef.setReceiverLatitude(nominalPosRef.getGeodeticLatitude());
+          mopsTMRef.setReceiverHeight(nominalPosRef.getAltitude());
+          mopsTMRef.setDayOfYear(gOriginal.header.epoch.DOYday());
+
+          mopsTM.setReceiverLatitude(lastPosition.getGeodeticLatitude());
+          mopsTM.setReceiverHeight(lastPosition.getAltitude());
+          mopsTM.setDayOfYear(gOriginal.header.epoch.DOYday());
+
+          // update mopsW
+          mopsW.setPosition(lastPosition);
+          mopsW.setDefaultEphemeris(sp3Store /*bceStore*/);
+
+          myModelRef.setDefaultEphemeris(sp3Store /*bceStore*/);
+          myModel.setDefaultEphemeris(sp3Store /*bceStore*/);
+   
+          model.Prepare(lastPosition);
+          myModelRef.Prepare(nominalPosRef);
+          myModel.Prepare(lastPosition);
+
+          // CASE #0
+          // We are using LMS + C1 + no tropo/iono model
+          gnssRinex gRin0(gOriginal);
+          try {
+                  gRin0 >> filter >> myModel;
+                  cout << "gRin0 : " << endl;
+                  display_rinex(gRin0);
+                  gRin0 >> solver;
+                  display_rinex(gRin0);
+
+          Position pos(myModel.rxPos.X()+solver.getSolution(TypeID::dx),
+                       myModel.rxPos.Y()+solver.getSolution(TypeID::dy),
+                       myModel.rxPos.Z()+solver.getSolution(TypeID::dz));
+             
+                  cout << "  Case 0 :";
+                  display_position(pos, nominalPos_true);
+          }
+          catch (Exception e) {
+                  cout << "Case 0. Exception : " << e.what() << endl;
+          }
 
    //////////////////////////// CASE #1  ////////////////////////////
 
@@ -455,12 +463,12 @@ int main(void)
 
              // 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]) );
+             Position pos(model.rxPos.X() + solver.getSolution(TypeID::dx),
+                          model.rxPos.Y() + solver.getSolution(TypeID::dy),
+                          model.rxPos.Z() + solver.getSolution(TypeID::dz));
              
              cout << "  Case 1 :";
-             display_position(solPos, nominalPos_true);
+             display_position(pos, nominalPos_true);
       }
       catch(...)
       {
@@ -492,6 +500,9 @@ int main(void)
              
              cout << "  Case 3 :";
              display_position(pos, nominalPos_true);
+             double altitude = pos.asGeodetic().getAltitude();
+             if ((0.0 < altitude) && (altitude < 300.0))
+                     lastPosition = pos;
       }
       catch (Exception e)
       {
@@ -500,219 +511,26 @@ int main(void)
 
    ////////////////////////// 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
-
+      // Case #4 : use elevation as weight
          // Let's make a working copy
       gnssRinex gRin4(gOriginal);
 
-      try
-      {
-
-             gRin4 >> filter >> markCSC1 >> smoothC1 >> model >> mopsW
-                   >> 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
+      try  {
+             gRin4 >> filter >> model >> elevationW >> solverWMS;
 
              Position pos((model.rxPos.X()+solverWMS.getSolution(TypeID::dx)),
                           (model.rxPos.Y()+solverWMS.getSolution(TypeID::dy)),
                           (model.rxPos.Z()+solverWMS.getSolution(TypeID::dz)));
-
+             
              cout << "  Case 4 :";
              display_position(pos, nominalPos_true);
       }
-      catch(...)
-      {
-         cerr << "Case 4. Exception at epoch: " << gRin4.header.epoch << 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);
-      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 >> filter >> markCSC1case5 >> smoothC1case5 >> model
-               >> mopsW;
-            // 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;
-      }
-
-   ////////////////////////// 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);
-      try {
-
-         gRin6 >> getPC >> pcFilter >> modelPC >> mopsW >> 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.
-
-             Position pos((model.rxPos.X()+solverWMS.getSolution(TypeID::dx)),
-                          (model.rxPos.Y()+solverWMS.getSolution(TypeID::dy)),
-                          (model.rxPos.Z()+solverWMS.getSolution(TypeID::dz)));
-
-             cout << "  Case 6 :";
-             display_position(pos, nominalPos_true);
-      }
-      catch(...)
-      {
-         cerr << "Case 6. Exception at epoch: " << gRin6.header.epoch << 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);
-      try
-      {
-
-         gRin7 >> getPC >> getLC >> getLI >> getMW >> markCSLI >> markCSMW
-               >> smoothPC >> pcFilter >> modelPC >> mopsW >> 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
-        Position pos((modelPC.rxPos.X()+solverWMS.getSolution(TypeID::dx)),
-                     (modelPC.rxPos.Y()+solverWMS.getSolution(TypeID::dy)),
-                     (modelPC.rxPos.Z()+solverWMS.getSolution(TypeID::dz)));
-
-        cout << "  Case 7 :";
-        display_position(pos, nominalPos_true);
-      }
-      catch(...)
-      {
-         cerr << "Case 7. Exception at epoch: " << gRin7.header.epoch << 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);
-      try
-      {
-            // First, let's process data up to the change of reference frame
-         gRin8 >> getPC >> getLC >> getLI >> getMW >> markCSLIcase8
-               >> markCSMWcase8 >> smoothPCcase8 >> pcFilter >> modelPC
-               >> mopsW;
-            // 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;
-      }
-
-   ////////////////////////// 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);
-      try
-      {
-
-         gRin9 >> getPC >> getLC >> getLI >> getMW >> markCSLIcase9
-               >> markCSMWcase9 >> smoothPCcase9 >> pcFilter >> modelPC
-               >> mopsW >> solverK9;
-            // VERY IMPORTANT: Note that in this case the coordinates are
-            // handled as constants, whereas the receiver clock is modeled as
-            // white noise.
-
-        Position pos((modelPC.rxPos.X()+solverK9.getSolution(TypeID::dx)),
-                     (modelPC.rxPos.Y()+solverK9.getSolution(TypeID::dy)),
-                     (modelPC.rxPos.Z()+solverK9.getSolution(TypeID::dz)));
-
-        cout << "  Case 9 :";
-        display_position(pos, nominalPos_true);
-      }
-      catch(...)
+      catch (Exception e)
       {
-         cerr << "Case 9. Exception at epoch: " << gRin9.header.epoch << endl;
+             cout << "Exception : " << e.what() << endl;
       }
 
-   ////////////////////////// END OF CASE #9  //////////////////////////
-
+      continue;
 
 
    //////////////////////////// CASE #10  ////////////////////////////
@@ -726,7 +544,7 @@ int main(void)
          // First, let's synchronize and process reference station data
       try
       {
-             gRef >> synchro >> filter >> modelRef;
+             gRef >> synchro >> filter >> myModelRef;
             // 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.
@@ -737,7 +555,7 @@ int main(void)
             // 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);
+             delta.setRefData(gRef.body);
 
       }
       catch(SynchronizeException& e)   // THIS IS VERY IMPORTANT IN ORDER TO
@@ -756,14 +574,14 @@ at epoch: " << gRef.header.epoch << endl;
       try
       {
 
-             gRin10 >> filter >> model >> delta >> solver;
+             gRin10 >> filter >> myModel >> delta >> solver;
             // 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.
 
-             Position pos((model.rxPos.X()+solver.getSolution(TypeID::dx)),
-                          (model.rxPos.Y()+solver.getSolution(TypeID::dy)),
-                          (model.rxPos.Z()+solver.getSolution(TypeID::dz)));
+             Position pos((myModel.rxPos.X()+solver.getSolution(TypeID::dx)),
+                          (myModel.rxPos.Y()+solver.getSolution(TypeID::dy)),
+                          (myModel.rxPos.Z()+solver.getSolution(TypeID::dz)));
              
              cout << "  Case 10:";
              display_position(pos, nominalPos_true);
@@ -794,11 +612,11 @@ at epoch: " << gRef.header.epoch << endl;
          // obtained from Case #10.
 
     try {
-           gRin11 >> filter >> model >> delta >> mopsW  >> solverWMS;
+           gRin11 >> filter >> myModel >> delta >> mopsW  >> solverWMS;
            // Like case #10, but now with "mopsW" and "solverWMS"
-           Position pos((model.rxPos.X()+solverWMS.getSolution(TypeID::dx)),
-                        (model.rxPos.Y()+solverWMS.getSolution(TypeID::dy)),
-                        (model.rxPos.Z()+solverWMS.getSolution(TypeID::dz)));
+           Position pos(myModel.rxPos.X()+solverWMS.getSolution(TypeID::dx),
+                        myModel.rxPos.Y()+solverWMS.getSolution(TypeID::dy),
+                        myModel.rxPos.Z()+solverWMS.getSolution(TypeID::dz));
            
            cout << "  Case 11:";
            display_position(pos, nominalPos_true);
@@ -829,19 +647,21 @@ at epoch: " << gRef.header.epoch << endl;
 
       try
       {
-
-         gRin12 >> filter >> model >> delta >> mopsW  >> solverK12;
+             if (i > 1) {
+                     gRin12 >> filter >> myModel >> delta
+                            >> mopsW  >> 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.
 
-        Position pos((modelPC.rxPos.X()+solverK12.getSolution(TypeID::dx)),
-                     (modelPC.rxPos.Y()+solverK12.getSolution(TypeID::dy)),
-                     (modelPC.rxPos.Z()+solverK12.getSolution(TypeID::dz)));
+             Position pos(myModel.rxPos.X()+solverK12.getSolution(TypeID::dx),
+                          myModel.rxPos.Y()+solverK12.getSolution(TypeID::dy),
+                          myModel.rxPos.Z()+solverK12.getSolution(TypeID::dz));
 
-        cout << "  Case 12:";
-        display_position(pos, nominalPos_true);
+             cout << "  Case 12:";
+             display_position(pos, nominalPos_true);
+             }
       }
       catch(...)
       {
@@ -850,6 +670,45 @@ at epoch: " << gRef.header.epoch << endl;
       }
 
    ////////////////////////// END OF CASE #12  //////////////////////////
+
+
+      /// Case #13
+      // Let's make a working copy
+      gnssRinex gRin13(gOriginal);
+
+      try {
+             gRin13 >> filter >> myModel >> delta >> elevationW >> solverWMS;
+
+             Position pos(myModel.rxPos.X()+solverWMS.getSolution(TypeID::dx),
+                          myModel.rxPos.Y()+solverWMS.getSolution(TypeID::dy),
+                          myModel.rxPos.Z()+solverWMS.getSolution(TypeID::dz));
+           
+             cout << "  Case 13:";
+             display_position(pos, nominalPos_true);
+      }
+      catch(...)
+      {
+             cout << "Exception #13 " << endl;
+      }
+
+      /// Case #14 : DGPS + iono/trop + WMS solver
+      // Let's make a working copy
+      gnssRinex gRin14(gOriginal);
+
+      try {
+             gRin14 >> filter >> model >> delta >> mopsW >> solverWMS;
+
+             Position pos(model.rxPos.X()+solverWMS.getSolution(TypeID::dx),
+                          model.rxPos.Y()+solverWMS.getSolution(TypeID::dy),
+                          model.rxPos.Z()+solverWMS.getSolution(TypeID::dz));
+           
+             cout << "  Case 14:";
+             display_position(pos, nominalPos_true);
+      }
+      catch(...)
+      {
+             cout << "Exception #14 " << endl;
+      }
    }
 
    return 0;