Merge branch 'master' of ssh://git.popipo.fr/~benoit/gps
authorBenoit Papillault <benoit.papillault@free.fr>
Mon, 26 Dec 2011 18:26:07 +0000 (19:26 +0100)
committerBenoit Papillault <benoit.papillault@free.fr>
Mon, 26 Dec 2011 18:26:07 +0000 (19:26 +0100)
Conflicts:
rnx2crx.c

1  2 
gpstk-example2.cpp
gpstk-example7.cpp
gpstk-solution4.cpp

@@@ -10,19 -10,19 +10,27 @@@ using namespace std
  using namespace gpstk;
  
  int main(void) {
++      int count = 0;
  
        // from : ftp://rgpdata.ign.fr/pub/data/2011/350/data_30/
        // Create the input file stream
--      RinexObsStream rin("cena0010.03d");
++      RinexObsStream rin("crei354k.11d");
        
        // Create the output file stream
--      RinexObsStream rout("cena0010.03d.new", ios::out|ios::trunc);
++      RinexObsStream rout("crei354k.11d.new", ios::out|ios::trunc);
        
        // Read the RINEX header
          RinexObsHeader head; //RINEX header object
--      rin >> head;
  
--      cout << "versionString : " << head.versionString << endl;
++      try {
++              rin >> head;
++      }
++
++      catch (Exception e) {
++              cout << "Exception : " << e.what() << endl;
++      }
++
++      head.dump(cout);
  
        rout.header = rin.header;
        rout << rout.header;
        // Loop over all data epochs
          RinexObsData data; //RINEX data object
        while (rin >> data) {
++              data.dump(cout);
                rout << data;
++              count ++;
        }
++
++      cout << count << " records" << endl;
        
        return 0;
  }
@@@ -821,7 -821,7 +821,11 @@@ int main(void
        cout << "  Case 10:\t ";
        try
        {
--         gRef >> synchro >> myFilter >> modelRef;
++            /*
++              gRef    : from rinRef("bell030a.02o")   - base
++              gRin10  : from rin("ebre030a.02o")      - rover
++            */
++            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.
              // The "delta" object will take care of proper differencing.
              // We must tell it which GNSS data structure will be used
              // as reference
++
++            cout << "gRin10 : " << endl;
++            cout << "  " << gRin10.header.epoch << endl;
++            gRin10.body.dump(cout,1);
++            cout << "gRef : " << endl;
++            cout << "  " << gRef.header.epoch << endl;
++            gRef.body.dump(cout,1);
++
++
           delta.setRefData(gRef.body);
  
        }
        catch(SynchronizeException& e)   // THIS IS VERY IMPORTANT IN ORDER TO
        {                                // MANAGE A POSSIBLE DESYNCHRONIZATION!!
--         cout << endl;
++            cout << "SynchronizeException: " << e.what() << endl;
           continue;
        }
        catch(...)
@@@ -856,15 -856,15 +869,19 @@@ at epoch: " << gRef.header.epoch << end
              // object that will adjust code prefit residuals BEFORE solving the
              // system of equations.
        }
--      catch(...)
++      catch(Exception e)
        {
--         cerr << "Case 10. Exception at epoch: " << gRin10.header.epoch
--              << endl;
++            cout << "Case 10 exception : " << e.what() << endl;
++            //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 << "dLat = "
++         << solverNEU.getSolution(TypeID::dLat) << "  ";  // dLat - Field #29
++      cout << "dLon = "
++         << solverNEU.getSolution(TypeID::dLon) << "  ";  // dLon - Field #30
++      cout << "dH   = "
++         << solverNEU.getSolution(TypeID::dH) << "  ";    // dH   - Field #31
        cout << endl;
  
     ////////////////////////// END OF CASE #10  //////////////////////////
index a7a728c,0000000..67328d5
mode 100644,000000..100644
--- /dev/null
@@@ -1,358 -1,0 +1,392 @@@
-       cout << fixed << setprecision(6);   // Set a proper output format
 +// Using DGPS ...
 +// reference station : MLVL (Marne la VallĂ©e)
 +// ECEF : 4201577.209, 189859.856, 4779064.567
++//
 +// rover station : CREI (Creil, Val d'Oise)
 +// ECEF : 4166409.705, 182770.997, 4809802.953
++//
 +// periode : 2011-12-20 (mardi = 2)
 +// GPS Week : 1667
 +// GPS Time : 210433
 +
 +// Files to download : 
 +// MLVL observation file :
 +// ftp://rgpdata.ensg.ign.fr:/pub/data/2011/354/data_1/mlvl354k.11d.Z
 +// ftp://rgpdata.ensg.ign.fr:/pub/data/2011/354/data_1/mlvl354k.11n.Z
 +// CREI observation file :
 +// ftp://rgpdata.ensg.ign.fr:/pub/data/2011/354/data_1/crei354k.11d.Z
 +// ftp://rgpdata.ensg.ign.fr:/pub/data/2011/354/data_1/crei354k.11n.Z
 +// ultra rapid ephemeris :
 +// ftp://rgpdata.ensg.ign.fr:/pub/products/ephemerides/1667/igr16672.sp3.Z
 +
 +#include <iostream>
 +#include <iomanip>
 +#include <math.h>
 +
 +// 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"
 +
 +using namespace std;
 +using namespace gpstk;
 +
 +void display(const gnssRinex& gRin) {
 +      cout << gRin.header.epoch
 +           << " " << gRin.header.epoch.second() << endl;
 +      gRin.body.dump(cout, 1);
 +}
 +
++Position display_solution(const ModelObs& model, const SolverLMS& solver) {
++
++      Position newPos(model.rxPos.X()+solver.getSolution(TypeID::dx),
++                      model.rxPos.Y()+solver.getSolution(TypeID::dy),
++                      model.rxPos.Z()+solver.getSolution(TypeID::dz),
++                      Position::Cartesian);
++
++      cout << "Solution is : "
++           << " X = " << newPos.X()
++           << " Y = " << newPos.Y()
++           << " Z = " << newPos.Z()
++           << endl;
++
++      return newPos;
++}
++
 +double distance(const Position& p1, const Position& p2) {
 +      double dx = p1.X() - p2.X();
 +      double dy = p1.Y() - p2.Y();
 +      double dz = p1.Z() - p2.Z();
 +
 +      return sqrt(dx*dx + dy*dy + dz*dz);
 +}
 +
 +int main(int argc, const char *argv[])
 +{
 +      DayTime t((short)643, (double)210433, (short)2011);
 +
 +      // display GPS week , GPS seconds
 +      cout << "Time : " << endl;
 +      cout << "  GPS week / GPS seconds : "   // 1667 / 210433
 +           << t.GPSfullweek() << " / " << t.GPSsow()
 +           << endl;
 +
 +      cout << "  Y-M-D h:m:s : "              // 2011-12-20 10:27:13
 +           << t.year() << "-" << t.month() << "-" << t.day() << " "
 +           << t.hour() << ":" << t.minute() << ":" << t.second()
 +           << endl;
 +
 +      cout << "  Year / Day of year : "       // 2011 / 354
 +           << t.DOYyear() << " / " << t.DOYday()
 +           << endl;
 +
 +      ////////// Initialization phase //////////
 +
 +      //////////// COMMON OBJECTS //////////////
 +      
-    RinexNavStream rnavin("mlvl354k.11n");
 +   // Read the navigation file
 +
 +   // Create the input navigation file stream
-    while (rnavin >> rNavData)
-    {
-          cout << "rNavData" << endl;
++      RinexNavStream rnavin("crei354k.11n" /*"mlvl354k.11n"*/);
 +   RinexNavHeader rNavHeader;          // Object to read the header of Rinex
 +                                       // navigation data files
 +   RinexNavData rNavData;              // Object to store Rinex navigation data
 +
 +   IonoModelStore ionoStore;           // Object to store ionospheric models
 +   IonoModel ioModel;                  // Declare a Ionospheric Model object
 +
 +   GPSEphemerisStore bceStore;         // Object to store satellites ephemeris
 +
 +   // 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);
 +
 +   // Display ionospheric parameters
 +   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;
 +
 +   // 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"
-    Position nominalPos(4166409.705, 182770.997, 4809802.953);
++   while (rnavin >> rNavData) {
 +         bceStore.addEphemeris(rNavData);
 +   }
 +
 +   bceStore.SearchPast();  // This is the default
 +
 +   // rover station nominal positional
-    SimpleFilter myFilter;
++   Position nominalPos(4166409.705, 182770.997, 4809802.953, Position::Cartesian);
 +
 +   // 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
-    Position nominalPosRef(4201577.209, 189859.856, 4779064.567);
++   SimpleFilter filter;
 +
 +   // 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("mlvl354k.11o");
 +
 +   // reference station nominal position
-          gRin.keepOnlyTypeID(TypeID::C1);
++   Position nominalPosRef(4201577.209, 189859.856, 4779064.567, Position::Cartesian);
 +
 +   // 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
 +
 +   // Declare a SolverLMS object
 +   SolverLMS solver;
 +
 +   //////// End of initialization phase  ////////
 +
 +
 +   //////// Processing phase ////////
 +
 +   // Create the input observation file stream
 +   // This is a fixed station, but here it will play as "rover"
 +   RinexObsStream rin("crei354k.11o");
 +   // Please note that data was collected in year 2002, when the Sun
 +   // was very active
 +
 +   // This is the GNSS data structure that will hold all the
 +   // GNSS-related information
 +   gnssRinex gRin;
 +
 +   // 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 "gRin". Default synchronization tolerance is
 +   // 1 s.
 +   Synchronize synchro(rinRef, gRin);
 +
 +   int i = 0;
 +
 +   // reference station / nominalPosRef
 +   // rover station / nominalPos
 +
 +   // 
 +   cout << "distance between ref and rover : "
 +      << ::distance(nominalPos, nominalPosRef) << endl;
 +
 +   // Loop over all data epochs
 +   for (;;) {
 +
++         cout << fixed << setprecision(6);   // Set a proper output format
++
 +         try {
 +                 rin >> gRin;
++                 // removed GLONASS satellites since they were producing
++                 // incorrect solution
++                 for (int i=1; i<= 30; i++) {
++                         SatID sat = SatID(i, SatID::systemGlonass);
++                         gRin.removeSatID(sat);
++                 }
 +         }
 +
 +         catch (Exception e) {
 +                 cout << "Exception : " << e.what() << endl;
 +         }
 +
 +         cout << "--- " << ++i << " ---" << endl;
 +         cout << "  nominalPos\t: X=" << nominalPos.X()
 +              << " Y=" << nominalPos.Y()
 +              << " Z=" << nominalPos.Z() << endl;
 +
-                  gRef >> synchro >> myFilter >> modelRef;
++         //gRin.keepOnlyTypeID(TypeID::C1);
 +
 +         // display gRin
 +         cout << "initial gRin : " << endl;
 +         display(gRin);
 +
++         // compute gRin position without DGPS
++         try {
++                 gRin >> filter >> model >> solver;
++                 cout << "gRin / initial Position : " << endl;
++                 display(gRin);
++                 display_solution(model, solver);
++         }
++
++         catch (Exception e) {
++                 cout << "Exception : " << e.what() << endl;
++         }
++
 +         // First, let's synchronize and process reference station data
 +         try {
-             gRin >> myFilter >> model >> delta >> solver;
++                 gRef >> synchro >> filter >> modelRef;
 +                 // Please note that the FIRST STEP is to synchronize
 +                 // "gRef", the reference station data stream, with "gRin"
 +                 // 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
 +
 +                 cout << "gRin : " << endl;
 +                 display(gRin);
 +                 cout << "gRef : " << endl;
 +                 display(gRef);
 +
 +                 delta.setRefData(gRef.body);
 +         }
 +      catch(SynchronizeException& e)   // THIS IS VERY IMPORTANT IN ORDER TO
 +      {                                // MANAGE A POSSIBLE DESYNCHRONIZATION!!
 +            cout << "SynchronizeException: " << e.what() << endl;
 +         continue;
 +      }
 +      catch(Exception e)
 +      {
 +            cout << "Exception1 : " << e.what() << endl;
 +      }
 +
 +      // Rover data processing is done here:
 +      try
 +      {
++            gRin >> filter >> model >> 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.
 +      }
 +      catch(Exception e)
 +      {
 +            cout << "Exception2 : " << e.what() << endl;
 +      }
 +
 +      cout << "final gRin : " << endl;
 +      display(gRin);
 +      cout << "final gRef : " << endl;
 +      display(gRef);
 +
 +      cout << "Solution is : " << solver.isValid() << endl;
 +      cout << "dx = " << solver.getSolution(TypeID::dx) << endl;
 +      cout << "dy = " << solver.getSolution(TypeID::dy) << endl;
 +      cout << "dz = " << solver.getSolution(TypeID::dz) << endl;
 +      cout << "cdt = " << solver.getSolution(TypeID::cdt) << endl;
 +      
 +      cout << "model.rxPos is : " <<endl;
 +      cout << "x = " << model.rxPos.X() << endl;
 +      cout << "y = " << model.rxPos.Y() << endl;
 +      cout << "z = " << model.rxPos.Z() << endl;
 +      
 +      Position newPos(model.rxPos.X()+solver.getSolution(TypeID::dx),
 +                    model.rxPos.Y()+solver.getSolution(TypeID::dy),
 +                    model.rxPos.Z()+solver.getSolution(TypeID::dz));
 +                    
 +      cout << "ECEF result : " << endl;
 +      cout << "x = " << newPos.X() << endl;
 +      cout << "y = " << newPos.Y() << endl;
 +      cout << "z = " << newPos.Z() << endl;
 +
 +      cout << "distance to nominalPos : "
 +         << ::distance(newPos, nominalPos) << endl;
 +   }
 +
 +
 +   return 0;
 +}