Added some gpstk examples (modified)
[gps] / gpstk-example8.cpp
1 #pragma ident "$Id: example8.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 ). 2008, 2009
22 //
23 //============================================================================
24
25 // Example program Nro 8 for GPSTk
26 //
27 // This program shows how to use GNSS Data Structures (GDS) to obtain
28 // "Precise Point Positioning" (PPP).
29 //
30 // For details on the PPP algorithm please consult:
31 //
32 //    Kouba, J. and P. Heroux. "Precise Point Positioning using IGS Orbit
33 //       and Clock Products". GPS Solutions, vol 5, pp 2-28. October, 2001.
34 //
35
36
37 #include <iostream>
38 #include <iomanip>
39
40
41    // Class for handling satellite observation parameters RINEX files
42 #include "RinexObsStream.hpp"
43
44    // Class to store satellite precise navigation data
45 #include "SP3EphemerisStore.hpp"
46
47    // Class in charge of basic GNSS signal modelling
48 #include "BasicModel.hpp"
49
50    // Class to model the tropospheric delays
51 #include "TropModel.hpp"
52
53    // Class defining the GNSS data structures
54 #include "DataStructures.hpp"
55
56    // Class to filter out satellites without required observables
57 #include "RequireObservables.hpp"
58
59    // Class to filter out observables grossly out of limits
60 #include "SimpleFilter.hpp"
61
62    // Class for easily changing reference base from ECEF to NEU
63 #include "XYZ2NEU.hpp"
64
65    // Class to detect cycle slips using LI combination
66 #include "LICSDetector2.hpp"
67
68    // Class to detect cycle slips using the Melbourne-Wubbena combination
69 #include "MWCSDetector.hpp"
70
71    // Class to compute the effect of solid tides
72 #include "SolidTides.hpp"
73
74    // Class to compute the effect of ocean loading
75 #include "OceanLoading.hpp"
76
77    // Class to compute the effect of pole tides
78 #include "PoleTides.hpp"
79
80    // Class to correct observables
81 #include "CorrectObservables.hpp"
82
83    // Class to compute the effect of wind-up
84 #include "ComputeWindUp.hpp"
85
86    // Class to compute the effect of satellite antenna phase center
87 #include "ComputeSatPCenter.hpp"
88
89    // Class to compute the tropospheric data
90 #include "ComputeTropModel.hpp"
91
92    // Class to compute linear combinations
93 #include "ComputeLinear.hpp"
94
95    // This class pre-defines several handy linear combinations
96 #include "LinearCombinations.hpp"
97
98    // Class to compute Dilution Of Precision values
99 #include "ComputeDOP.hpp"
100
101    // Class to keep track of satellite arcs
102 #include "SatArcMarker.hpp"
103
104    // Class to compute gravitational delays
105 #include "GravitationalDelay.hpp"
106
107    // Class to align phases with code measurements
108 #include "PhaseCodeAlignment.hpp"
109
110    // Compute statistical data
111 #include "PowerSum.hpp"
112
113    // Used to delete satellites in eclipse
114 #include "EclipsedSatFilter.hpp"
115
116    // Used to decimate data. This is important because RINEX observation
117    // data is provided with a 30 s sample rate, whereas SP3 files provide
118    // satellite clock information with a 900 s sample rate.
119 #include "Decimate.hpp"
120
121    // Class to compute the Precise Point Positioning (PPP) solution
122 #include "SolverPPP.hpp"
123
124 #include "geometry.hpp"                   // DEG_TO_RAD
125
126
127 using namespace std;
128 using namespace gpstk;
129
130
131 int main(void)
132 {
133
134       /////////////////// INITIALIZATION PART /////////////////////
135
136    cout << fixed << setprecision(3);   // Set a proper output format
137
138
139       // Create the input observation file stream
140    RinexObsStream rin("onsa2240.05o");
141
142
143       // Declare a "SP3EphemerisStore" object to handle precise ephemeris
144    SP3EphemerisStore SP3EphList;
145
146       // Set flags to reject satellites with bad or absent positional
147       // values or clocks
148    SP3EphList.rejectBadPositions(true);
149    SP3EphList.rejectBadClocks(true);
150
151       // Set flags to check for data gaps and too wide interpolation intervals.
152       // Default values for "gapInterval" (901.0 s) and "maxInterval"
153       // (8105.0 s) will be used.
154    SP3EphList.enableDataGapCheck();
155    SP3EphList.enableIntervalCheck();
156
157       // Load all the SP3 ephemerides files
158    SP3EphList.loadFile("igs13354.sp3");
159    SP3EphList.loadFile("igs13355.sp3");
160    SP3EphList.loadFile("igs13356.sp3");
161
162
163       // ONSA station nominal position
164    Position nominalPos(3370658.5419, 711877.1496, 5349786.9542);
165
166
167       // Declare a NeillTropModel object, setting the defaults
168    NeillTropModel neillTM( nominalPos.getAltitude(),
169                            nominalPos.getGeodeticLatitude(), 224);
170
171
172       // This is the GNSS data structure that will hold all the
173       // GNSS-related information
174    gnssRinex gRin;
175
176
177       // Declare a base-changing object: From ECEF to North-East-Up (NEU)
178    XYZ2NEU baseChange(nominalPos);
179
180
181       // Declare an object to check that all required observables are present
182    RequireObservables requireObs(TypeID::P1);
183    requireObs.addRequiredType(TypeID::P2);
184    requireObs.addRequiredType(TypeID::L1);
185    requireObs.addRequiredType(TypeID::L2);
186
187
188       // Declare a simple filter object to screen PC
189    SimpleFilter pcFilter;
190    pcFilter.setFilteredType(TypeID::PC);
191
192
193       // Declare a basic modeler
194    BasicModel basic(nominalPos, SP3EphList);
195
196
197       // Objects to mark cycle slips
198    LICSDetector2 markCSLI;     // Checks LI cycle slips
199    MWCSDetector markCSMW;      // Checks Merbourne-Wubbena cycle slips
200
201
202       // Object to compute tidal effects
203    SolidTides  solid;
204
205       // Ocean loading model
206    OceanLoading ocean("OCEAN-GOT00.dat");
207
208       // Numerical values are x,y pole displacements for Aug/12/2005 (arcsec).
209    PoleTides   pole(0.02094, 0.42728);
210
211
212       // Vector from ONSA antenna ARP to L1 phase center [UEN] (AOAD/M_B)
213    Triple offsetL1(0.0780, 0.000, 0.000);   // Units in meters
214
215       // Vector from ONSA antenna ARP to L2 phase center [UEN] (AOAD/M_B)
216    Triple offsetL2(0.096, 0.0000, 0.000);    // Units in meters
217
218       // Vector from monument to antenna ARP [UEN] for ONSA station
219    Triple offsetARP(0.9950, 0.0, 0.0);    // Units in meters
220
221
222       // Declare an object to correct observables to monument
223    CorrectObservables corr(SP3EphList);
224    ((corr.setNominalPosition(nominalPos)).setL1pc(offsetL1)).setL2pc(offsetL2);
225    corr.setMonument(offsetARP);
226
227
228       // Object to compute wind-up effect
229    ComputeWindUp windup(SP3EphList, nominalPos, "PRN_GPS");
230
231
232       // Object to compute satellite antenna phase center effect
233    ComputeSatPCenter svPcenter(nominalPos);
234
235
236       // Object to compute the tropospheric data
237    ComputeTropModel computeTropo(neillTM);
238
239
240       // This object defines several handy linear combinations
241    LinearCombinations comb;
242
243
244       // Object to compute linear combinations for cycle slip detection
245    ComputeLinear linear1(comb.pdeltaCombination);
246    linear1.addLinear(comb.ldeltaCombination);
247    linear1.addLinear(comb.mwubbenaCombination);
248    linear1.addLinear(comb.liCombination);
249
250
251       // This object computes the ionosphere-free combinations to be used
252       // as observables in the PPP processing
253    ComputeLinear linear2(comb.pcCombination);
254    linear2.addLinear(comb.lcCombination);
255
256
257       // Let's use a different object to compute prefit residuals
258    ComputeLinear linear3(comb.pcPrefit);
259    linear3.addLinear(comb.lcPrefit);
260
261
262       // Declare an object to process the data using PPP. It is set
263       // to use a NEU system
264    SolverPPP pppSolver(true);
265
266       // The real test for a PPP processing program is to handle coordinates
267       // as white noise. In such case, position error should be about 0.25 m or
268       // better. Uncomment the following couple of lines to test this.
269 //   WhiteNoiseModel wnM(100.0);            // 100 m of sigma
270 //   pppSolver.setCoordinatesModel(&wnM);
271
272
273       // Object to keep track of satellite arcs
274    SatArcMarker markArc;
275    markArc.setDeleteUnstableSats(true);
276    markArc.setUnstablePeriod(151.0);
277
278       // Object to compute gravitational delay effects
279    GravitationalDelay grDelay(nominalPos);
280
281       // Object to align phase with code measurements
282    PhaseCodeAlignment phaseAlign;
283
284       // Object to compute DOP values
285    ComputeDOP cDOP;
286
287       // Object to remove eclipsed satellites
288    EclipsedSatFilter eclipsedSV;
289
290       // Object to decimate data. This is important because RINEX observation
291       // data is provided with a 30 s sample rate, whereas SP3 files provide
292       // satellite clock information with a 900 s sample rate.
293    Decimate decimateData(900.0, 5.0, SP3EphList.getInitialTime());
294
295       // When you are printing the model, you may want to comment the previous
296       // line and uncomment the following one, generating a 30 s model
297 //   Decimate decimateData(30.0, 1.0, SP3EphList.getInitialTime());
298
299       // Statistical summary objects
300    PowerSum errorVectorStats;
301
302
303
304       /////////////////// PROCESING PART /////////////////////
305
306
307       // Use this variable to select between position printing or model printing
308    bool printPosition(true);     // By default, print position and associated
309                                  // parameters
310
311
312
313       // Loop over all data epochs
314    while(rin >> gRin)
315    {
316
317       DayTime time(gRin.header.epoch);
318
319          // Compute the effect of solid, oceanic and pole tides
320       Triple tides( solid.getSolidTide(time, nominalPos) +
321                     ocean.getOceanLoading("ONSA", time)  +
322                     pole.getPoleTide(time, nominalPos)     );
323
324          // Update observable correction object with tides information
325       corr.setExtraBiases(tides);
326
327       try
328       {
329
330             // The following lines are indeed just one line
331          gRin >> requireObs      // Check if required observations are present
332               >> linear1         // Compute linear combinations to detect CS
333               >> markCSLI        // Mark cycle slips: LI algorithm
334               >> markCSMW        // Mark cycle slips: Melbourne-Wubbena
335               >> markArc         // Keep track of satellite arcs
336               >> decimateData    // If not a multiple of 900 s, decimate
337               >> basic           // Compute the basic components of model
338               >> eclipsedSV      // Remove satellites in eclipse
339               >> grDelay         // Compute gravitational delay
340               >> svPcenter       // Compute the effect of satellite phase center
341               >> corr            // Correct observables from tides, etc.
342               >> windup          // Compute wind-up effect
343               >> computeTropo    // Compute tropospheric effect
344               >> linear2         // Compute ionosphere-free combinations
345               >> pcFilter        // Filter out spurious data
346               >> phaseAlign      // Align phases with codes
347               >> linear3         // Compute prefit residuals
348               >> baseChange      // Prepare to use North-East-UP reference frame
349               >> cDOP            // Compute DOP figures
350               >> pppSolver;      // Solve equations with a Kalman filter
351
352       }
353       catch(DecimateEpoch& d)
354       {
355             // If 'decimateData' object detects that this epoch is not a
356             // multiple of 900 s, it issues a "DecimateEpoch" exception. Here
357             // we catch such exception, and just continue to process next epoch.
358          continue;
359       }
360       catch(Exception& e)
361       {
362          cerr << "Exception at epoch: " << time << "; " << e << endl;
363          continue;
364       }
365       catch(...)
366       {
367          cerr << "Unknown exception at epoch: " << time << endl;
368          continue;
369       }
370
371
372          // Check if we want to print model or position
373       if(printPosition)
374       {
375             // Print here the position results
376          cout << time.DOYsecond()      << "  ";     // Epoch - Output field #1
377
378          cout << pppSolver.getSolution(TypeID::dLat) << "  ";    // dLat  - #2
379          cout << pppSolver.getSolution(TypeID::dLon) << "  ";    // dLon  - #3
380          cout << pppSolver.getSolution(TypeID::dH) << "  ";      // dH    - #4
381          cout << pppSolver.getSolution(TypeID::wetMap) << "  ";  // Tropo - #5
382
383          cout << pppSolver.getVariance(TypeID::dLat) << "  "; // Cov dLat - #6
384          cout << pppSolver.getVariance(TypeID::dLon) << "  "; // Cov dLon - #7
385          cout << pppSolver.getVariance(TypeID::dH) << "  ";   // Cov dH   - #8
386          cout << pppSolver.getVariance(TypeID::wetMap) << "  ";//Cov Tropo- #9
387
388          cout << gRin.numSats()        << "  ";     // Visible satellites - #10
389
390          cout << cDOP.getGDOP()        << "  ";                   // GDOP - #11
391          cout << cDOP.getPDOP()        << "  ";                   // PDOP - #12
392          cout << cDOP.getTDOP()        << "  ";                   // TDOP - #13
393          cout << cDOP.getHDOP()        << "  ";                   // HDOP - #14
394          cout << cDOP.getVDOP()        << "  ";                   // VDOP - #15
395
396          cout << endl;
397
398             // For statistical purposes we discard the first two hours of data
399          if (time.DOYsecond() > 7200.0)
400          {
401                // Statistical summary
402             double errorV( pppSolver.solution[1]*pppSolver.solution[1] +
403                            pppSolver.solution[2]*pppSolver.solution[2] +
404                            pppSolver.solution[3]*pppSolver.solution[3] );
405
406                // Get module of position error vector
407             errorV = std::sqrt(errorV);
408
409                // Add to statistical summary object
410             errorVectorStats.add(errorV);
411          }
412
413       }  // End of position printing
414       else
415       {
416             // Print here the model results
417             // First, define types we want to keep
418          TypeIDSet types;
419          types.insert(TypeID::L1);
420          types.insert(TypeID::L2);
421          types.insert(TypeID::P1);
422          types.insert(TypeID::P2);
423          types.insert(TypeID::PC);
424          types.insert(TypeID::LC);
425          types.insert(TypeID::rho);
426          types.insert(TypeID::dtSat);
427          types.insert(TypeID::rel);
428          types.insert(TypeID::gravDelay);
429          types.insert(TypeID::tropo);
430          types.insert(TypeID::dryTropo);
431          types.insert(TypeID::dryMap);
432          types.insert(TypeID::wetTropo);
433          types.insert(TypeID::wetMap);
434          types.insert(TypeID::tropoSlant);
435          types.insert(TypeID::windUp);
436          types.insert(TypeID::satPCenter);
437          types.insert(TypeID::satX);
438          types.insert(TypeID::satY);
439          types.insert(TypeID::satZ);
440          types.insert(TypeID::elevation);
441          types.insert(TypeID::azimuth);
442          types.insert(TypeID::satArc);
443          types.insert(TypeID::prefitC);
444          types.insert(TypeID::prefitL);
445          types.insert(TypeID::dx);
446          types.insert(TypeID::dy);
447          types.insert(TypeID::dz);
448          types.insert(TypeID::dLat);
449          types.insert(TypeID::dLon);
450          types.insert(TypeID::dH);
451          types.insert(TypeID::cdt);
452
453          gRin.keepOnlyTypeID(types);   // Delete the types not in 'types'
454
455             // Iterate through the GNSS Data Structure
456          satTypeValueMap::const_iterator it;
457          for (it = gRin.body.begin(); it!= gRin.body.end(); it++) 
458          {
459
460                // Print epoch
461             cout << time.year()        << " ";
462             cout << time.DOY()         << " ";
463             cout << time.DOYsecond()   << " ";
464
465             cout << cDOP.getGDOP()        << "  ";  // GDOP #4
466             cout << cDOP.getPDOP()        << "  ";  // PDOP #5
467             cout << cDOP.getTDOP()        << "  ";  // TDOP #6
468             cout << cDOP.getHDOP()        << "  ";  // HDOP #7
469             cout << cDOP.getVDOP()        << "  ";  // VDOP #8
470
471                // Print satellite information (system and PRN)
472             cout << (*it).first << " ";
473
474                // Print model values
475             typeValueMap::const_iterator itObs;
476             for( itObs  = (*it).second.begin(); 
477                  itObs != (*it).second.end();
478                  itObs++ )
479             {
480                bool printNames(true);  // Whether print types' names or not
481                if (printNames)
482                {
483                   cout << (*itObs).first << " ";
484                }
485
486                cout << (*itObs).second << " ";
487
488             }  // End for( itObs = ... )
489
490             cout << endl;
491
492          }  // End for (it = gRin.body.begin(); ... )
493
494       }  // End of model printing
495
496    }  // End of 'while(rin >> gRin)...'
497
498
499
500       // Print statistical summary in cerr
501    if(printPosition)
502    {
503       cerr << "Module of error vector: Average = "
504            << errorVectorStats.average() << " m    Std. dev. = "
505            << std::sqrt(errorVectorStats.variance()) << " m" << endl;
506    }
507
508
509
510    exit(0);       // End of program
511
512 }