7c02d1ed6de67e0c793fd9051a1a1e3e243df80d
[gps] / gpstk-solution6.cpp
1 /*
2   Reference : CREI
3   Mobile    : MLVL
4   Date : 2011-12-20 0h-1h/1s UTC
5 */
6 #include <iostream>
7 #include <iomanip>
8
9 #include <math.h>
10
11 // Class for handling satellite observation parameters RINEX files
12 #include "RinexObsStream.hpp"
13
14 // Classes for handling RINEX Broadcast ephemeris files
15 #include "RinexNavStream.hpp"
16 #include "RinexNavHeader.hpp"
17 #include "RinexNavData.hpp"
18
19 // Class in charge of the GPS signal modelling
20 #include "ModelObs.hpp"
21
22 // Class to store satellite broadcast navigation data
23 #include "GPSEphemerisStore.hpp"
24 #include "SP3EphemerisStore.hpp"
25
26 // Class to model the tropospheric delays
27 #include "TropModel.hpp"
28
29 // Classes to model and store ionospheric delays
30 #include "IonoModel.hpp"
31 #include "IonoModelStore.hpp"
32
33 // Class to solve the equation system using Least Mean Squares
34 #include "SolverLMS.hpp"
35
36 // Class to solve the equation system using Weighted-Least Mean Squares
37 #include "SolverWMS.hpp"
38
39 // Class to solve equations systems using a simple code-based Kalman filter
40 #include "CodeKalmanSolver.hpp"
41
42 // Class defining the GNSS data structures
43 #include "DataStructures.hpp"
44
45 // Class to filter out observables grossly out of limits
46 #include "SimpleFilter.hpp"
47
48 // Class for easily changing reference base from ECEF to NEU
49 #include "XYZ2NEU.hpp"
50
51 // Class to detect cycle slips using just one frequency
52 #include "OneFreqCSDetector.hpp"
53
54 // Class to detect cycle slips using LI combination
55 #include "LICSDetector.hpp"
56
57 // Class to detect cycle slips using the Melbourne-Wubbena combination
58 #include "MWCSDetector.hpp"
59
60 // Class to compute weights according to Appendix J of MOPS C (RTCA/DO-229C)
61 #include "ComputeMOPSWeights.hpp"
62
63 // Class to compute single differences between receiver stations
64 #include "DeltaOp.hpp"
65
66 // Class to synchronize two GNSS Data Structures data streams.
67 #include "Synchronize.hpp"
68
69 #include "geometry.hpp"                   // DEG_TO_RAD
70
71
72 using namespace std;
73 using namespace gpstk;
74
75 // various global initialisation
76
77 // mobile station nominal position
78 //Position nominalPos_true(4201577.209, 189859.856,     4779064.567);
79 //Position nominalPos_true(48.803125, 2.263235, 20, Position::Geodetic);
80 Position nominalPos_true(48.803403, 2.263000, 133.0, Position::Geodetic);
81
82 // class for weight 
83
84 class ComputeElevationWeights : public WeightBase, public ProcessingClass
85 {
86 public:
87         ComputeElevationWeights() {};
88
89         double getWeight(const SatID& sat,
90                          typeValueMap& tvMap) {
91                 double elevation = tvMap(TypeID::elevation);
92                 if (elevation < 45.0)
93                         return elevation;
94                 return 90.0;
95         };
96
97         virtual satTypeValueMap& Process(const DayTime& time,
98                                          satTypeValueMap& gData) {
99                 
100                 // Loop through all the satellites
101                 satTypeValueMap::iterator it;
102                 for (it = gData.begin(); it != gData.end(); ++it ) {
103                         (*it).second[TypeID::weight] =
104                                 getWeight((*it).first, (*it).second);
105                 }
106
107                 return gData;
108         };
109
110         virtual gnssSatTypeValue& Process(gnssSatTypeValue& gData) {
111                 Process(gData.header.epoch, gData.body); return gData;
112         };
113
114         virtual gnssRinex& Process(gnssRinex& gData) {
115                 Process(gData.header.epoch, gData.body); return gData;
116         };
117
118         virtual int getIndex(void) const {
119                 return 6; /* FIXME */
120         };
121
122         virtual string getClassName(void) const {
123                 return "ComputeElevationWeights";
124         };
125 };
126
127 class ComputeSNRWeights : public WeightBase, public ProcessingClass
128 {
129 public:
130         ComputeSNRWeights() {};
131
132         double getWeight(const SatID& sat,
133                          typeValueMap& tvMap) {
134                 double snr = tvMap(TypeID::S1);
135                 return snr;
136         };
137
138         virtual satTypeValueMap& Process(const DayTime& time,
139                                          satTypeValueMap& gData) {
140                 
141                 // Loop through all the satellites
142                 satTypeValueMap::iterator it;
143                 for (it = gData.begin(); it != gData.end(); ++it ) {
144                         (*it).second[TypeID::weight] =
145                                 getWeight((*it).first, (*it).second);
146                 }
147
148                 return gData;
149         };
150
151         virtual gnssSatTypeValue& Process(gnssSatTypeValue& gData) {
152                 Process(gData.header.epoch, gData.body); return gData;
153         };
154
155         virtual gnssRinex& Process(gnssRinex& gData) {
156                 Process(gData.header.epoch, gData.body); return gData;
157         };
158
159         virtual int getIndex(void) const {
160                 return 7; /* FIXME */
161         };
162
163         virtual string getClassName(void) const {
164                 return "ComputeSNRWeights";
165         };
166 };
167
168 void display_rinex(const gnssRinex& gRin) {
169         //      cout << gRin.header.epoch
170         //           << " " << gRin.header.epoch.second() << endl;
171         cout << gRin.header.epoch
172              << " GPS week " << gRin.header.epoch.GPSfullweek()
173              << setprecision(9)
174              << " GPS seconds " << gRin.header.epoch.GPSsow()
175              << setprecision(6)
176              << endl;
177         gRin.body.dump(cout, 1);
178 }
179
180 void display_position(const Position& p,
181                       const Position& nominalPos) {
182         Position pos = p;
183
184         double dx = pos.X() - nominalPos.X();
185         double dy = pos.Y() - nominalPos.Y();
186         double dz = pos.Z() - nominalPos.Z();
187         double distance = sqrt(dx*dx+dy*dy+dz*dz);
188
189         cout << " X = " << pos.asECEF().X()
190              << " Y = " << pos.asECEF().Y()
191              << " Z = " << pos.asECEF().Z()
192              << " lon = " << pos.asGeodetic().getLongitude()
193              << " lat = " << pos.asGeodetic().geodeticLatitude()
194              << " alt = " << pos.asGeodetic().getAltitude()
195              << " d = " << distance
196              << endl;
197 }
198
199 void display_solution(const gnssRinex& gData,
200                       const ModelObsFixedStation& model,
201                       const SolverLMS& solver) {
202
203         static double sumX = 0;
204         static double sumY = 0;
205         static double sumZ = 0;
206         static int count   = 0;
207
208         Position pos(model.rxPos.X()+solver.getSolution(TypeID::dx),
209                      model.rxPos.Y()+solver.getSolution(TypeID::dy),
210                      model.rxPos.Z()+solver.getSolution(TypeID::dz));
211
212         DayTime T = gData.header.epoch - solver.getSolution(TypeID::cdt) / C_GPS_M;
213
214         double dx = pos.X() - nominalPos_true.X();
215         double dy = pos.Y() - nominalPos_true.Y();
216         double dz = pos.Z() - nominalPos_true.Z();
217         double distance = sqrt(dx*dx+dy*dy+dz*dz);
218         double altitude = pos.asGeodetic().getAltitude();
219
220         cout << setprecision(1)
221              << " X= " << pos.X()
222              << " Y= " << pos.Y()
223              << " Z= " << pos.Z()
224              << setprecision(9)
225              << " T= " << T.GPSfullweek() << "/" << T.GPSsecond()
226              << setprecision(6)
227              << " lon= " << pos.asGeodetic().getLongitude()
228              << " lat= " << pos.asGeodetic().geodeticLatitude()
229              << setprecision(0)
230              << " alt= " << altitude
231              << " cdt= " << solver.getSolution(TypeID::cdt)
232              << setprecision(9)
233              << " (" << solver.getSolution(TypeID::cdt) / C_GPS_M << ")"
234              << setprecision(6)
235              << " d= " << distance              
236              << endl;
237
238
239         if ((0<=altitude) && (altitude<=300)) {
240                 sumX += pos.X();
241                 sumY += pos.Y();
242                 sumZ += pos.Z();
243                 count ++;
244
245                 if (count > 10) {
246                         cout << "adjusting nominalPos_true" << endl;
247
248                         nominalPos_true = Position(sumX/count,
249                                                    sumY/count,
250                                                    sumZ/count,
251                                                    Position::Cartesian);
252                 }
253         }
254 }
255
256 /* keep only GPS satellites (removing Glonass for instance, ... )*/
257 gnssRinex keepOnlySatGPS(const gnssRinex& rin) {
258         gnssRinex result;
259         SatID::SatelliteSystem gps = SatID::systemGPS;
260
261         result.header = rin.header;
262
263         for (satTypeValueMap::const_iterator it = rin.body.begin();
264              it != rin.body.end(); it++) {
265                 if ((*it).first.system == gps)
266                         result.body[(*it).first] = (*it).second;
267         }
268
269         return result;
270 }
271
272 int main(int argc, const char *argv[])
273 {
274         // first set
275         //const char obsReference[] = "crei354a.11o";
276         //const char obsMobile[] = "mlvl354a.11o";
277         //const char navReference[] = "crei354a.11n";
278
279         // second set : from rgpdata.ign.fr:/pub/data/2011/361/data_1
280         const char obsReference[] = "crei361o.11o";
281         const char obsMobile[] = "out7.ooo";
282         const char navReference[] = "brdc363z.11n";
283         const char sp3File[] = "igu16682_12.sp3";
284         
285         // reference station nominal position
286         Position nominalPosRef  (4166409.705,   182770.997,     4809802.953);
287         // initial set to nominalPos
288         Position lastPosition = nominalPosRef;
289         bool lastPositionFound = false;
290
291         ////////// Initialization phase //////////
292
293         //////////// COMMON OBJECTS //////////////
294
295         RinexNavData rNavData;         // Object to store Rinex navigation data
296         GPSEphemerisStore bceStore;    // Object to store satellites ephemeris
297         SP3EphemerisStore sp3Store;
298         RinexNavHeader rNavHeader;     // Object to read the header of Rinex
299                                        // navigation data files
300         IonoModelStore ionoStore;      // Object to store ionospheric models
301         IonoModel ioModel;             // Declare a Ionospheric Model object
302
303         cout << "reading navigation file :\t" << navReference << endl;
304
305         // Create the input navigation file stream
306         RinexNavStream rnavin(navReference);
307
308         // We need to read ionospheric parameters (Klobuchar model) from header
309         rnavin >> rNavHeader;
310
311         // Let's feed the ionospheric model (Klobuchar type) from data in the
312         // Navigation file header
313         ioModel.setModel(rNavHeader.ionAlpha, rNavHeader.ionBeta);
314
315         for (int i = 0; i < 4; i++)
316                 cout << "ionAlpha = " << rNavHeader.ionAlpha[i] << endl;
317         for (int i = 0; i < 4; i++)
318                 cout << "ionBeta  = " << rNavHeader.ionBeta[i]  << endl;
319
320         // Beware: In this case, the same model will be used for the
321         // full data span
322         ionoStore.addIonoModel(DayTime::BEGINNING_OF_TIME, ioModel);
323
324         // Storing the ephemeris in "bceStore"
325         while (rnavin >> rNavData) {
326                 bceStore.addEphemeris(rNavData);
327         }
328
329         sp3Store.loadFile(sp3File);
330
331         bceStore.SearchPast();  // This is the default
332
333         // Declare a MOPSTropModel object, setting the defaults
334         MOPSTropModel mopsTM;
335
336         // Declare a simple filter object. By default, it filters C1 with
337         // default limits
338         SimpleFilter filter;
339         filter.setMaxLimit(100000000.0);
340
341         // This is the GNSS data structure that will hold all the
342         // GNSS-related information
343         gnssRinex gOriginal;
344
345         //////////// CASE #1 OBJECTS ////////////
346
347         // Declare a SolverLMS object
348         SolverLMS solver;
349
350         //////////// CASE #3 OBJECTS ////////////
351
352         // This object will compute the appropriate MOPS weights
353         ComputeMOPSWeights mopsW;
354    
355         // Declare a solver object using Weighted-Least-Mean-Squares and
356         // a topocentric reference system (NEU)
357         SolverWMS solverWMS;
358  
359
360       //////////// CASE #5 OBJECTS ////////////
361
362       // The core of this case is to add a new equation to the equation system
363       // Such equation states that there are NO changes in height for the
364       // rover.
365       //                          dH = 0
366       //
367       // Add a "fake" satellite to identify the new equation: Sat #1 of
368       // system "UserDefined"
369    SatID satEq(1,SatID::systemUserDefined);
370
371       // Declare and fill a "typeValueMap" object that will hold
372       // the equation data
373    typeValueMap equTVMap;
374    equTVMap[TypeID::prefitC] = 0.0;    // Code prefit residual is zero
375    equTVMap[TypeID::dLat]    = 0.0;    // Geometry matrix dLat coefficient is
376                                        // zero
377    equTVMap[TypeID::dLon]    = 0.0;    // Geometry matrix dLon coefficient is
378                                        // zero
379    equTVMap[TypeID::dH]      = 1.0;    // Geometry matrix dH coefficient is
380                                        // 1.0!!!
381    equTVMap[TypeID::cdt]     = 0.0;    // Geometry matrix cdt coefficient is
382                                        // zero
383
384       // Assign a relatively high weight to this information (typical
385       // MOPS weights range from 0.01 to 0.04)
386       // This means that this equation is very important for us, but it is
387       // NOT ABSOLUTELY TRUE. Some variation is allowed
388       // Given that weights are indeed (1/variances), if we assign to our new
389       // equation a confidence of 0.5 m of sigma, it means that we should use
390       // a weight of (1/(0.5^2)) = 4 m^(-2)
391    equTVMap[TypeID::weight]  = 4.0;
392
393       ////////////////////////////////////////
394
395       //////////// CASE #10 OBJECTS ////////////
396
397       // This is the GNSS data structure that will hold the
398       // reference station data
399    gnssRinex gRef;
400
401    // Create the input observation file stream for REFERENCE STATION
402    RinexObsStream rinRef(obsReference);
403    cout << "reading observation file (reference) :\t" << obsReference << endl;
404
405       // Declare a MOPSTropModel object for the reference station, setting
406       // the defaults
407    MOPSTropModel mopsTMRef;
408
409       // Declare the modeler object, setting all the parameters in one pass
410       // Given that in this example we are using a fixed GPS station with known
411       // coordinates, you could have used the "ModeledReferencePR" class, which
412       // is a little bit simpler.
413       // However, for a rover is more appropriate to use a "ModelObs" object
414       // because it allows to update the apriori position more easily (and it
415       // may automatically compute one, if needed, using Bancroft's method)
416    //ModelObs model(nominalPos_true,
417    //             ionoStore, mopsTM,
418    //             bceStore, TypeID::C1);
419
420    ModelObs model;
421    model.setDefaultIonoModel(ionoStore);
422    model.setDefaultTropoModel(mopsTM);
423    model.setDefaultEphemeris(/*sp3Store*/bceStore);
424    //model.setMinElev(1.0);
425
426       // Declare the appropriate modeler object for a reference station
427    //ModelObsFixedStation modelRef(nominalPosRef,
428    //                    ionoStore, mopsTMRef,
429    //                            bceStore, TypeID::C1);
430
431    ModelObs modelRef;
432    modelRef.setDefaultIonoModel(ionoStore);
433    modelRef.setDefaultTropoModel(mopsTMRef);
434    modelRef.setDefaultEphemeris(/*sp3Store*/ bceStore);
435    modelRef.Prepare(nominalPosRef);
436
437    // use simple model without iono & tropo corrections
438    ModelObs myModelRef;
439    ModelObs myModel;
440    
441       // Create an object to compute the single differences of prefit residuals
442    DeltaOp delta;      // By default, it will work on code prefit residuals
443
444       // Create an object to synchronize rover and reference station
445       // data streams. This object will take data out from "rinRef" until
446       // it is synchronized with data in "gOriginal". Default synchronization
447       // tolerance is 1 s.
448    Synchronize synchro(rinRef, gOriginal);
449
450       //////////////////////////////////////////////
451
452
453
454       //////////// CASE #12 OBJECTS ////////////
455
456       // Declare a new Kalman solver
457    CodeKalmanSolver solverK12;
458
459    ComputeElevationWeights elevationW;
460    ComputeSNRWeights snrW;
461
462       ////////////////////////////////////////
463
464       //////// End of initialization phase  ////////
465
466
467       //////// Processing phase ////////
468
469    int i = 0;
470
471
472    cout << fixed << setprecision(6);   // Set a proper output format
473
474    cout << "reading observation file (mobile) : " << obsMobile << endl;
475
476       // Create the input observation file stream
477       // This is a fixed station, but here it will play as "rover"
478    RinexObsStream rin(obsMobile);
479
480       // Loop over all data epochs
481    while (rin >> gOriginal) {
482
483            cout << "--- " << ++i << " ---" << endl;
484            cout << "  Ref.   :";
485            display_position(nominalPosRef, nominalPos_true);
486            cout << "  Mobile :";
487            display_position(nominalPos_true, nominalPos_true);
488            cout << "  Last " << lastPositionFound << " :";
489            display_position(lastPosition, nominalPos_true);
490            cout << "gOriginal : " << endl;
491            display_rinex(gOriginal);
492
493            // keep only GPS satellites (removing GLONASS)
494            gOriginal = keepOnlySatGPS(gOriginal);
495
496            // keep only C1 & S1 (since that what we can record now)
497            TypeIDSet set;
498            set.insert(TypeID::C1);
499            set.insert(TypeID::S1);
500            gOriginal.keepOnlyTypeID(set);
501
502            // adjust time to the nearest millisecond
503            /*
504            DayTime old_t = gOriginal.header.epoch;
505            double clock_bias = old_t.GPSsecond() - floor(1.0*old_t.GPSsecond())/1.0;
506            DayTime new_t = old_t - clock_bias;
507            gOriginal.header.epoch = new_t;
508            printf("  clock_bias %.9f\n", clock_bias);
509            */
510
511            // adjust pseudorange
512            /*
513            double clock_bias = 0.001;
514            for (satTypeValueMap::iterator it = gOriginal.body.begin();
515                 it != gOriginal.body.end(); ++it) {
516                    double old_val = (*it).second[TypeID::C1];
517                    double new_val = old_val - clock_bias * C_GPS_M;
518                    cout << (*it).first << " adjusted " << old_val << " to " << new_val << endl;
519                    (*it).second[TypeID::C1] = new_val;
520            }
521            */
522
523            // based on SimpleFilter.hpp and Bancroft.hpp
524            // pseudorange should be within [15000000 ... 30000000]
525
526            cout << "gOriginal(adjust):" << endl;
527            display_rinex(gOriginal);
528
529            // update MOPSTropModel
530            mopsTMRef.setReceiverLatitude(nominalPosRef.getGeodeticLatitude());
531            mopsTMRef.setReceiverHeight(nominalPosRef.getAltitude());
532            mopsTMRef.setDayOfYear(gOriginal.header.epoch.DOYday());
533
534            mopsTM.setReceiverLatitude(lastPosition.getGeodeticLatitude());
535            mopsTM.setReceiverHeight(lastPosition.getAltitude());
536            mopsTM.setDayOfYear(gOriginal.header.epoch.DOYday());
537
538            // update mopsW
539            mopsW.setPosition(lastPosition);
540            mopsW.setDefaultEphemeris(/*sp3Store*/bceStore);
541
542            myModelRef.setDefaultEphemeris(/*sp3Store*/ bceStore);
543            myModel.setDefaultEphemeris(/*sp3Store*/bceStore);
544    
545            myModelRef.Prepare(nominalPosRef);
546            if (lastPositionFound) {
547                    model.Prepare(lastPosition);
548                    myModel.Prepare(lastPosition);
549            } else {
550                    model.setModelPrepared(false);
551                    myModel.setModelPrepared(false);
552                    //model.Prepare(nominalPosRef);
553                    //myModel.Prepare(nominalPosRef);
554            }
555
556            // CASE #0
557            // We are using LMS + C1 + no tropo/iono model
558            gnssRinex gRin0(gOriginal);
559            try {
560                    myModel.setModelPrepared(false);
561
562                    gRin0 >> filter;
563                    //cout << "gRin0(filter) : " << endl;
564                    //display_rinex(gRin0);
565                    gRin0 >> myModel;
566                    cout << "gRin0(myModel) : " << endl;
567                    display_rinex(gRin0);
568                    gRin0 >> solver;
569                    //cout << "gRin0(solver) : " << endl;
570                    //display_rinex(gRin0);
571
572            Position pos(myModel.rxPos.X()+solver.getSolution(TypeID::dx),
573                         myModel.rxPos.Y()+solver.getSolution(TypeID::dy),
574                         myModel.rxPos.Z()+solver.getSolution(TypeID::dz));
575               
576            //cout << "  Case 0 :";
577                    //display_solution(gRin0, myModel, solver);
578
579                    double altitude = pos.asGeodetic().getAltitude();
580                    if ((0.0 < altitude) && (altitude < 300.0))  {
581                            lastPosition = pos;
582                            lastPositionFound = true;
583                    }
584            }
585            catch (Exception e) {
586                    cout << "Case 0. Exception : " << e.what() << endl;
587            }
588
589    //////////////////////////// CASE #1  ////////////////////////////
590
591          // This case is a common C1 + Least Mean Squares solver
592          // (LMS) processing
593
594          // Let's make a working copy
595       gnssRinex gRin1(gOriginal);
596
597       try {
598             // This is the line that will process all the GPS data
599               gRin1 >> filter >> model >> solver;
600
601             // - First, a basic filter to screen out very bad observables
602             // - Second, apply a model to the observables (ionosphere,
603             //   troposphere, relativity, etc.)
604             // - Third, solve the equations using a simple Least-Mean-Squares
605             //   solver
606
607               // Get your results out of the solver object. In ECEF system
608               // by default
609               Position pos(model.rxPos.X() + solver.getSolution(TypeID::dx),
610                            model.rxPos.Y() + solver.getSolution(TypeID::dy),
611                            model.rxPos.Z() + solver.getSolution(TypeID::dz));
612               
613               //cout << "  Case 1 :";
614               //display_solution(gRin1, model, solver);
615       }
616       catch(...)
617       {
618          cerr << "Case 1. Exception at epoch: " << gRin1.header.epoch << endl;
619       }
620
621    ////////////////////////// END OF CASE #1  //////////////////////////
622
623
624
625    //////////////////////////// CASE #3  ////////////////////////////
626
627          // In this case we process data using C1 + Weighted Least Mean Squares
628          // solver (WMS)
629
630          // Let's make a working copy
631       gnssRinex gRin3(gOriginal);
632
633       try  {
634
635               gRin3 >> filter >> model >> mopsW >> solverWMS;
636             // The "mopsW" object computes weights based on MOPS algorithm
637             // The "solverWMS" object solves the system using Weighted Least
638             // Mean Squares. It is already configured to work with NEU system.
639
640               Position pos((model.rxPos.X()+solverWMS.getSolution(TypeID::dx)),
641                            (model.rxPos.Y()+solverWMS.getSolution(TypeID::dy)),
642                            (model.rxPos.Z()+solverWMS.getSolution(TypeID::dz)));
643               
644               cout << "  Case 3 :";
645               display_solution(gRin3, model, solverWMS);
646               double altitude = pos.asGeodetic().getAltitude();
647               if ((0.0 < altitude) && (altitude < 300.0))  {
648                       lastPosition = pos;
649                       lastPositionFound = true;
650               }
651       }
652       catch (Exception e)
653       {
654               cout << "Exception : " << e.what() << endl;
655       }
656
657    ////////////////////////// END OF CASE #3  //////////////////////////
658
659       // Case #4 : use elevation as weight
660          // Let's make a working copy
661       gnssRinex gRin4(gOriginal);
662
663       try  {
664               gRin4 >> filter >> model >> elevationW >> solverWMS;
665
666               Position pos((model.rxPos.X()+solverWMS.getSolution(TypeID::dx)),
667                            (model.rxPos.Y()+solverWMS.getSolution(TypeID::dy)),
668                            (model.rxPos.Z()+solverWMS.getSolution(TypeID::dz)));
669               
670               //cout << "  Case 4 :";
671               //display_solution(gRin4, model, solverWMS);
672       }
673       catch (Exception e)
674       {
675               cout << "Exception : " << e.what() << endl;
676       }
677
678
679       // Case #5 : use SNR as weight
680       gnssRinex gRin5(gOriginal);
681
682       try  {
683               gRin5 >> filter >> model >> snrW >> solverWMS;
684               //cout << "  Case 5 :";
685               //display_solution(gRin5, model, solverWMS);
686       }
687       catch (Exception e)
688       {
689               cout << "Exception : " << e.what() << endl;
690       }
691
692       continue;
693
694
695    //////////////////////////// CASE #10  ////////////////////////////
696
697
698          // This is like cases #1 and #2, but using DGPS techniques instead.
699
700          // Let's make a working copy of rover data
701       gnssRinex gRin10(gOriginal);
702
703          // First, let's synchronize and process reference station data
704       try
705       {
706               gRef >> synchro >> filter >> myModelRef;
707             // Please note that the FIRST STEP is to synchronize "gRef", the
708             // reference station data stream, with "gOriginal" (or with gRin10,
709             // which is the same), the rover receiver data stream.
710             //
711             // Also, remember that in simple DGPS the differences are computed
712             // on code prefit residuals, so "modelRef" object is mandatory.
713
714             // The "delta" object will take care of proper differencing.
715             // We must tell it which GNSS data structure will be used
716             // as reference
717               delta.setRefData(gRef.body);
718
719       }
720       catch(SynchronizeException& e)   // THIS IS VERY IMPORTANT IN ORDER TO
721       {                                // MANAGE A POSSIBLE DESYNCHRONIZATION!!
722          cout << endl;
723          continue;
724       }
725       catch(...)
726       {
727          cerr << "Case 10. Exception when processing reference station data \
728 at epoch: " << gRef.header.epoch << endl;
729       }
730
731
732          // Rover data processing is done here:
733       try
734       {
735
736               gRin10 >> filter >> myModel >> delta >> solver;
737             // This is very similar to cases #1 and #2, but we insert a "delta"
738             // object that will adjust code prefit residuals BEFORE solving the
739             // system of equations.
740
741               Position pos((myModel.rxPos.X()+solver.getSolution(TypeID::dx)),
742                            (myModel.rxPos.Y()+solver.getSolution(TypeID::dy)),
743                            (myModel.rxPos.Z()+solver.getSolution(TypeID::dz)));
744               
745               cout << "  Case 10:";
746               display_solution(gRin10, myModel, solver);
747               lastPosition = pos;
748               lastPositionFound = true;
749       }
750       catch(...)
751       {
752          cerr << "Case 10. Exception at epoch: " << gRin10.header.epoch
753               << endl;
754       }
755
756    ////////////////////////// END OF CASE #10  //////////////////////////
757
758
759
760    //////////////////////////// CASE #11  ////////////////////////////
761
762
763          // This is like case #10 (DGPS), but now let's apply a WMS solver
764          // on data
765
766          // Let's make a working copy
767          gnssRinex gRin11(gOriginal);
768
769
770          // Please note that data streams are already synchronized, and
771          // "delta" object may be reused with the same reference data
772          // obtained from Case #10.
773
774     try {
775             // Like case #10, but now with "mopsW" and "solverWMS"
776             gRin11 >> filter >> myModel >> delta >> mopsW  >> solverWMS;
777             cout << "  Case 11:";
778             //display_position(pos, nominalPos_true);
779             display_solution(gRin11, myModel, solverWMS);
780       }
781       catch(...)
782       {
783          cerr << "Case 11. Exception at epoch: " << gRin11.header.epoch
784               << endl;
785       }
786
787    ////////////////////////// END OF CASE #11  //////////////////////////
788
789
790
791    //////////////////////////// CASE #12  ////////////////////////////
792
793
794          // This is like case #11 (DGPS), but now let's apply a simple
795          // Kalman filter on data
796
797          // Let's make a working copy
798          gnssRinex gRin12(gOriginal);
799
800
801          // Please note that data streams are already synchronized, and
802          // "delta" object may be reused with the same reference data
803          // obtained from Case #10.
804
805       try
806       {
807               if (i > 1) {
808                       gRin12 >> filter >> myModel >> delta
809                              >> mopsW  >> solverK12;
810                       cout << "  Case 12:";
811                       display_solution(gRin12, myModel, solverK12);
812               }
813       }
814       catch(...)
815       {
816          cerr << "Case 12. Exception at epoch: " << gRin12.header.epoch
817               << endl;
818       }
819
820    ////////////////////////// END OF CASE #12  //////////////////////////
821
822
823       /// Case #13
824       // Let's make a working copy
825       gnssRinex gRin13(gOriginal);
826
827       try {
828               gRin13 >> filter >> myModel >> delta >> elevationW >> solverWMS;
829               cout << "  Case 13:";
830               display_solution(gRin13, myModel, solverWMS);
831       }
832       catch(...)
833       {
834               cout << "Exception #13 " << endl;
835       }
836
837       /// Case #14 : DGPS + iono/trop + WMS solver
838       // Let's make a working copy
839       gnssRinex gRin14(gOriginal);
840
841       try {
842               gRin14 >> filter >> model >> delta >> mopsW >> solverWMS;
843               cout << "  Case 14:";
844               display_solution(gRin14, model, solverWMS);
845       }
846       catch(...)
847       {
848               cout << "Exception #14 " << endl;
849       }
850    }
851
852    return 0;
853 }