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