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

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

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