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