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