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