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