source: ntrip/trunk/BNC/src/PPP/pppParlist.cpp@ 9241

Last change on this file since 9241 was 8993, checked in by stuerze, 4 years ago

minor changes in PPP
  • Property svn:keywords set to Author Date Id Rev URL;svn:eol-style=native
  • Property svn:mime-type set to text/plain
File size: 15.9 KB
Line 
1/* -------------------------------------------------------------------------
2 * BKG NTRIP Client
3 * -------------------------------------------------------------------------
4 *
5 * Class: t_pppParlist
6 *
7 * Purpose: List of estimated parameters
8 *
9 * Author: L. Mervart
10 *
11 * Created: 29-Jul-2014
12 *
13 * Changes:
14 *
15 * -----------------------------------------------------------------------*/
16
17#include <cmath>
18#include <iostream>
19#include <sstream>
20#include <iomanip>
21#include <algorithm>
22#include <newmatio.h>
23
24#include "pppParlist.h"
25#include "pppSatObs.h"
26#include "pppStation.h"
27#include "bncutils.h"
28#include "bncconst.h"
29#include "pppClient.h"
30
31using namespace BNC_PPP;
32using namespace std;
33
34// Constructor
35////////////////////////////////////////////////////////////////////////////
36t_pppParam::t_pppParam(e_type type, const t_prn& prn, t_lc::type tLC,
37 const vector<t_pppSatObs*>* obsVector) {
38
39 _type = type;
40 _prn = prn;
41 _tLC = tLC;
42 _x0 = 0.0;
43 _indexOld = -1;
44 _indexNew = -1;
45 _noise = 0.0;
46 _ambInfo = 0;
47
48 switch (_type) {
49 case crdX:
50 _epoSpec = false;
51 _sigma0 = OPT->_aprSigCrd[0];
52 _noise = OPT->_noiseCrd[0];
53 break;
54 case crdY:
55 _epoSpec = false;
56 _sigma0 = OPT->_aprSigCrd[1];
57 _noise = OPT->_noiseCrd[1];
58 break;
59 case crdZ:
60 _epoSpec = false;
61 _sigma0 = OPT->_aprSigCrd[2];
62 _noise = OPT->_noiseCrd[2];
63 break;
64 case clkR:
65 _epoSpec = true;
66 _sigma0 = OPT->_noiseClk;
67 break;
68 case amb:
69 _epoSpec = false;
70 _sigma0 = OPT->_aprSigAmb;
71 _ambInfo = new t_ambInfo();
72 if (obsVector) {
73 for (unsigned ii = 0; ii < obsVector->size(); ii++) {
74 const t_pppSatObs* obs = obsVector->at(ii);
75 if (obs->prn() == _prn) {
76 double offGG = 0;
77 if (_prn.system() == 'R' && tLC != t_lc::MW) {
78 offGG = PPP_CLIENT->offGG();
79 }
80 double offGB = 0;
81 if (_prn.system() == 'C' && tLC != t_lc::MW) {
82 offGB = PPP_CLIENT->offGB();
83 }
84 _x0 = floor((obs->obsValue(tLC) - offGG - offGB - obs->cmpValue(tLC)) / obs->lambda(tLC) + 0.5);
85 break;
86 }
87 }
88 }
89 break;
90 case offGG:
91 _epoSpec = true;
92 _sigma0 = 1000.0;
93 _x0 = PPP_CLIENT->offGG();
94 break;
95 case offGB:
96 _epoSpec = true;
97 _sigma0 = 1000.0;
98 _x0 = PPP_CLIENT->offGB();
99 break;
100 case trp:
101 _epoSpec = false;
102 _sigma0 = OPT->_aprSigTrp;
103 _noise = OPT->_noiseTrp;
104 break;
105 case ion:
106 _epoSpec = false;
107 _sigma0 = OPT->_aprSigIon;
108 _noise = OPT->_noiseIon;
109 break;
110 case cBias1:
111 case cBias2:
112 _epoSpec = false;
113 _sigma0 = OPT->_aprSigCodeBias;
114 _noise = OPT->_noiseCodeBias;
115 break;
116 case pBias1:
117 case pBias2:
118 _epoSpec = false;
119 _sigma0 = OPT->_aprSigPhaseBias;
120 _noise = OPT->_noisePhaseBias;
121 break;
122 }
123}
124
125// Destructor
126////////////////////////////////////////////////////////////////////////////
127t_pppParam::~t_pppParam() {
128 delete _ambInfo;
129}
130
131//
132////////////////////////////////////////////////////////////////////////////
133double t_pppParam::partial(const bncTime& /* epoTime */, const t_pppSatObs* obs,
134 const t_lc::type& tLC, const t_prn refPrn) const {
135
136 // Special Case - Melbourne-Wuebbena
137 // ---------------------------------
138 if (tLC == t_lc::MW && _type != amb) {
139 return 0.0;
140 }
141
142 const t_pppStation* sta = PPP_CLIENT->staRover();
143 ColumnVector rhoV = sta->xyzApr() - obs->xc().Rows(1,3);
144
145 map<t_frequency::type, double> codeCoeff;
146 map<t_frequency::type, double> phaseCoeff;
147 map<t_frequency::type, double> ionoCoeff;
148 obs->lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
149
150 switch (_type) {
151 case crdX:
152 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
153 return (sta->xyzApr()[0] - obs->xc()[0]) / rhoV.NormFrobenius();
154 case crdY:
155 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
156 return (sta->xyzApr()[1] - obs->xc()[1]) / rhoV.NormFrobenius();
157 case crdZ:
158 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
159 return (sta->xyzApr()[2] - obs->xc()[2]) / rhoV.NormFrobenius();
160 case clkR:
161 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
162 return 1.0;
163 case offGG:
164 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
165 return (obs->prn().system() == 'R') ? 1.0 : 0.0;
166 case offGB:
167 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
168 return (obs->prn().system() == 'C') ? 1.0 : 0.0;
169 case amb:
170 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
171 else if ((OPT->_obsModelType == OPT->IF) ||
172 (OPT->_obsModelType == OPT->PPPRTK) ||
173 (OPT->_obsModelType == OPT->UncombPPP) ||
174 (OPT->_obsModelType == OPT->DCMcodeBias && !obs->isReference()) ||
175 (OPT->_obsModelType == OPT->DCMphaseBias && !obs->isReference()) ) {
176 if (obs->prn() == _prn) {
177 if (tLC == _tLC) {
178 return (obs->lambda(tLC));
179 }
180 else if (tLC == t_lc::lIF && _tLC == t_lc::MW) {
181 return obs->lambda(t_lc::lIF) * obs->lambda(t_lc::MW) / obs->lambda(t_lc::l2);
182 }
183 else {
184 if (_tLC == t_lc::l1) {
185 return obs->lambda(t_lc::l1) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
186 }
187 else if (_tLC == t_lc::l2) {
188 return obs->lambda(t_lc::l2) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
189 }
190 }
191 }
192 }
193 break;
194 case trp:
195 if (tLC == t_lc::GIM) {
196 return 0.0;
197 }
198 else if (tLC == t_lc::Tz0) {
199 return 1.0;
200 }
201 else {
202 return 1.0 / sin(obs->eleSat());
203 }
204 case ion:
205 if (obs->prn() == _prn) {
206 if (tLC == t_lc::c1) {
207 return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::c1)];
208 }
209 else if (tLC == t_lc::c2) {
210 return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::c2)];
211 }
212 else if (tLC == t_lc::l1) {
213 return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
214 }
215 else if (tLC == t_lc::l2) {
216 return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
217 }
218 else if (tLC == t_lc::GIM) {
219 return -1.0;
220 }
221 }
222 if (tLC == t_lc::GIM && _prn == refPrn) {
223 return 1.0;
224 }
225 break;
226 case cBias1:
227 if (tLC == t_lc::c1) {
228 return 1.0;
229 }
230 else {
231 return 0.0;
232 }
233 break;
234 case cBias2:
235 if (tLC == t_lc::c2) {
236 return 1.0;
237 }
238 else {
239 return 0.0;
240 }
241 break;
242 case pBias1:
243 if (tLC == t_lc::l1) {
244 return 1.0;
245 }
246 else {
247 return 0.0;
248 }
249 break;
250 case pBias2:
251 if (tLC == t_lc::l2) {
252 return 1.0;
253 }
254 else {
255 return 0.0;
256 }
257 }
258 return 0.0;
259}
260
261//
262////////////////////////////////////////////////////////////////////////////
263string t_pppParam::toString() const {
264 stringstream ss;
265 switch (_type) {
266 case crdX:
267 ss << "CRD_X";
268 break;
269 case crdY:
270 ss << "CRD_Y";
271 break;
272 case crdZ:
273 ss << "CRD_Z";
274 break;
275 case clkR:
276 ss << "REC_CLK ";
277 break;
278 case offGG:
279 ss << "OGG ";
280 break;
281 case offGB:
282 ss << "OGB ";
283 break;
284 case trp:
285 ss << "TRP ";
286 break;
287 case amb:
288 ss << "AMB " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
289 break;
290 case ion:
291 ss << "ION " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
292 break;
293 case cBias1:
294 case cBias2:
295 case pBias1:
296 case pBias2:
297 ss << "BIAS " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << "REC";
298 break;
299 }
300 return ss.str();
301}
302
303// Constructor
304////////////////////////////////////////////////////////////////////////////
305t_pppParlist::t_pppParlist() {
306}
307
308// Destructor
309////////////////////////////////////////////////////////////////////////////
310t_pppParlist::~t_pppParlist() {
311 for (unsigned ii = 0; ii < _params.size(); ii++) {
312 delete _params[ii];
313 }
314}
315
316//
317////////////////////////////////////////////////////////////////////////////
318t_irc t_pppParlist::set(const bncTime& epoTime,
319 const std::vector<t_pppSatObs*>& obsVector,
320 const QMap<char, t_pppRefSat*>& refSatMap) {
321
322 // Remove some Parameters
323 // ----------------------
324 vector<t_pppParam*>::iterator it = _params.begin();
325 while (it != _params.end()) {
326 t_pppParam* par = *it;
327
328 bool remove = false;
329
330 if (par->epoSpec()) {
331 remove = true;
332 }
333
334 else if (par->type() == t_pppParam::crdX ||
335 par->type() == t_pppParam::crdY ||
336 par->type() == t_pppParam::crdZ) {
337 if (par->lastObsTime().valid() && (epoTime - par->lastObsTime() > 60.0)) {
338 remove = true;
339 }
340 }
341 else if (par->type() == t_pppParam::amb) {
342 char system = par->prn().system();
343 t_prn refPrn = t_prn();
344 if (OPT->_refSatRequired) {
345 refPrn = (refSatMap[system])->prn();
346 }
347 if ((par->lastObsTime().valid() && (epoTime - par->lastObsTime() > 60.0)) ||
348 ( refPrn == par->prn())) {
349 remove = true;
350 }
351 }
352
353 if (remove) {
354 delete par;
355 it = _params.erase(it);
356 }
357 else {
358 ++it;
359 }
360 }
361
362 // Check whether parameters have observations
363 // ------------------------------------------
364 for (unsigned ii = 0; ii < _params.size(); ii++) {
365 t_pppParam* par = _params[ii];
366 if (par->prn() == 0) {
367 par->setLastObsTime(epoTime);
368 if (par->firstObsTime().undef()) {
369 par->setFirstObsTime(epoTime);
370 }
371 }
372 else {
373 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
374 const t_pppSatObs* satObs = obsVector[jj];
375 if (satObs->prn() == par->prn()) {
376 par->setLastObsTime(epoTime);
377 if (par->firstObsTime().undef()) {
378 par->setFirstObsTime(epoTime);
379 }
380 break;
381 }
382 }
383 }
384 }
385
386 // Required Set of Parameters
387 // --------------------------
388 vector<t_pppParam*> required;
389
390 // Coordinates
391 // -----------
392 required.push_back(new t_pppParam(t_pppParam::crdX, t_prn(), t_lc::dummy));
393 required.push_back(new t_pppParam(t_pppParam::crdY, t_prn(), t_lc::dummy));
394 required.push_back(new t_pppParam(t_pppParam::crdZ, t_prn(), t_lc::dummy));
395
396 // Receiver Clock
397 // --------------
398 required.push_back(new t_pppParam(t_pppParam::clkR, t_prn(), t_lc::dummy));
399
400 // GPS-GLONASS Clock Offset
401 // ------------------------
402 if (OPT->useSystem('R')) {
403 required.push_back(new t_pppParam(t_pppParam::offGG, t_prn(), t_lc::dummy));
404 }
405
406 // GPS-BDS Clock Offset
407 // ------------------------
408 if (OPT->useSystem('C')) {
409 required.push_back(new t_pppParam(t_pppParam::offGB, t_prn(), t_lc::dummy));
410 }
411
412
413 // Troposphere
414 // -----------
415 if (OPT->estTrp()) {
416 required.push_back(new t_pppParam(t_pppParam::trp, t_prn(), t_lc::dummy));
417 }
418
419 // Ionosphere
420 // ----------
421 if (OPT->_obsModelType == OPT->UncombPPP ||
422 OPT->_obsModelType == OPT->DCMcodeBias ||
423 OPT->_obsModelType == OPT->DCMphaseBias ) {
424 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
425 const t_pppSatObs* satObs = obsVector[jj];
426 required.push_back(new t_pppParam(t_pppParam::ion, satObs->prn(), t_lc::dummy));
427 }
428 }
429
430 // Ambiguities
431 // -----------
432 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
433 const t_pppSatObs* satObs = obsVector[jj];
434 if ((OPT->_obsModelType == OPT->IF) ||
435 (OPT->_obsModelType == OPT->PPPRTK) ||
436 (OPT->_obsModelType == OPT->UncombPPP) ||
437 (OPT->_obsModelType == OPT->DCMcodeBias && !satObs->isReference()) ||
438 (OPT->_obsModelType == OPT->DCMphaseBias && !satObs->isReference()) ) {
439 const vector<t_lc::type>& ambLCs = OPT->ambLCs(satObs->prn().system());
440 for (unsigned ii = 0; ii < ambLCs.size(); ii++) {
441 required.push_back(new t_pppParam(t_pppParam::amb, satObs->prn(), ambLCs[ii], &obsVector));
442 }
443 }
444 }
445
446 // Receiver Code Biases
447 // --------------------
448 if (OPT->_obsModelType == OPT->DCMcodeBias) {
449 required.push_back(new t_pppParam(t_pppParam::cBias1, t_prn(), t_lc::c1));
450 required.push_back(new t_pppParam(t_pppParam::cBias2, t_prn(), t_lc::c2));
451 }
452
453 // Receiver Phase Biases
454 // ---------------------
455 if ((OPT->_obsModelType == OPT->DCMphaseBias) ||
456 (OPT->_obsModelType == OPT->PPPRTK) ) {
457 required.push_back(new t_pppParam(t_pppParam::pBias1, t_prn(), t_lc::l1));
458 required.push_back(new t_pppParam(t_pppParam::pBias2, t_prn(), t_lc::l2));
459 }
460
461 // Check if all required parameters are present
462 // --------------------------------------------
463 for (unsigned ii = 0; ii < required.size(); ii++) {
464 t_pppParam* parReq = required[ii];
465
466 bool found = false;
467 for (unsigned jj = 0; jj < _params.size(); jj++) {
468 t_pppParam* parOld = _params[jj];
469 if (parOld->isEqual(parReq)) {
470 found = true;
471 break;
472 }
473 }
474 if (found) {
475 delete parReq;
476 }
477 else {
478 _params.push_back(parReq);
479 }
480 }
481
482 // Set Parameter Indices
483 // ---------------------
484 sort(_params.begin(), _params.end(), t_pppParam::sortFunction);
485
486 for (unsigned ii = 0; ii < _params.size(); ii++) {
487 t_pppParam* par = _params[ii];
488 par->setIndex(ii);
489 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
490 const t_pppSatObs* satObs = obsVector[jj];
491 if (satObs->prn() == par->prn()) {
492 par->setAmbEleSat(satObs->eleSat());
493 par->stepAmbNumEpo();
494 }
495 }
496 }
497
498 return success;
499}
500
501//
502////////////////////////////////////////////////////////////////////////////
503void t_pppParlist::printResult(const bncTime& epoTime, const SymmetricMatrix& QQ,
504 const ColumnVector& xx) const {
505
506 string epoTimeStr = string(epoTime);
507 const t_pppStation* sta = PPP_CLIENT->staRover();
508
509 LOG << endl;
510
511 t_pppParam* parX = 0;
512 t_pppParam* parY = 0;
513 t_pppParam* parZ = 0;
514 for (unsigned ii = 0; ii < _params.size(); ii++) {
515 t_pppParam* par = _params[ii];
516 if (par->type() == t_pppParam::crdX) {
517 parX = par;
518 }
519 else if (par->type() == t_pppParam::crdY) {
520 parY = par;
521 }
522 else if (par->type() == t_pppParam::crdZ) {
523 parZ = par;
524 }
525 else {
526 int ind = par->indexNew();
527 double apr = (par->type() == t_pppParam::trp) ?
528 t_tropo::delay_saast(sta->xyzApr(), M_PI/2.0) : par->x0();
529 LOG << epoTimeStr << ' ' << par->toString() << ' '
530 << setw(10) << setprecision(4) << apr << ' '
531 << showpos << setw(10) << setprecision(4) << xx[ind] << noshowpos << " +- "
532 << setw(8) << setprecision(4) << sqrt(QQ[ind][ind]);
533 if (par->type() == t_pppParam::amb) {
534 LOG << " el = " << setw(6) << setprecision(2) << par->ambEleSat() * 180.0 / M_PI
535 << " epo = " << setw(4) << par->ambNumEpo();
536 }
537 LOG << endl;
538 }
539 }
540
541 if (parX && parY && parZ) {
542
543 ColumnVector xyz(3);
544 xyz[0] = xx[parX->indexNew()];
545 xyz[1] = xx[parY->indexNew()];
546 xyz[2] = xx[parZ->indexNew()];
547
548 ColumnVector neu(3);
549 xyz2neu(sta->ellApr().data(), xyz.data(), neu.data());
550
551 SymmetricMatrix QQxyz = QQ.SymSubMatrix(1,3);
552
553 SymmetricMatrix QQneu(3);
554 covariXYZ_NEU(QQxyz, sta->ellApr().data(), QQneu);
555
556 LOG << epoTimeStr << ' ' << sta->name()
557 << " X = " << setprecision(4) << sta->xyzApr()[0] + xyz[0] << " +- "
558 << setprecision(4) << sqrt(QQxyz[0][0])
559
560 << " Y = " << setprecision(4) << sta->xyzApr()[1] + xyz[1] << " +- "
561 << setprecision(4) << sqrt(QQxyz[1][1])
562
563 << " Z = " << setprecision(4) << sta->xyzApr()[2] + xyz[2] << " +- "
564 << setprecision(4) << sqrt(QQxyz[2][2])
565
566 << " dN = " << setprecision(4) << neu[0] << " +- "
567 << setprecision(4) << sqrt(QQneu[0][0])
568
569 << " dE = " << setprecision(4) << neu[1] << " +- "
570 << setprecision(4) << sqrt(QQneu[1][1])
571
572 << " dU = " << setprecision(4) << neu[2] << " +- "
573 << setprecision(4) << sqrt(QQneu[2][2])
574
575 << endl;
576 }
577}
578
Note: See TracBrowser for help on using the repository browser.