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

Last change on this file since 10938 was 10938, checked in by stuerze, 11 days ago

bugfix regarding ionospheric constraints

  • Property svn:keywords set to Author Date Id Rev URL;svn:eol-style=native
  • Property svn:mime-type set to text/plain
File size: 14.6 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 LC,
37 const vector<t_pppSatObs*>* obsVector) {
38
39 _type = type;
40 _prn = prn;
41 _sys = '?';
42 _LC = LC;
43 _x0 = 0.0;
44 _indexOld = -1;
45 _indexNew = -1;
46 _noise = 0.0;
47 _ambInfo = new t_ambInfo();
48
49 switch (_type) {
50 case crdX:
51 _epoSpec = false;
52 _sigma0 = OPT->_aprSigCrd[0];
53 _noise = OPT->_noiseCrd[0];
54 break;
55 case crdY:
56 _epoSpec = false;
57 _sigma0 = OPT->_aprSigCrd[1];
58 _noise = OPT->_noiseCrd[1];
59 break;
60 case crdZ:
61 _epoSpec = false;
62 _sigma0 = OPT->_aprSigCrd[2];
63 _noise = OPT->_noiseCrd[2];
64 break;
65 case rClk:
66 _epoSpec = true;
67 _sigma0 = OPT->_aprSigClk;
68 break;
69 case amb:
70 _epoSpec = false;
71 _sigma0 = OPT->_aprSigAmb;
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 _x0 = floor((obs->obsValue(LC) - obs->cmpValue(LC)) / obs->lambda(LC) + 0.5);
77 break;
78 }
79 }
80 }
81 break;
82 case trp:
83 _epoSpec = false;
84 _sigma0 = OPT->_aprSigTrp;
85 _noise = OPT->_noiseTrp;
86 break;
87 case ion:
88 _epoSpec = true;
89 _sigma0 = OPT->_aprSigIon;
90 break;
91 case bias:
92 _epoSpec = true;
93 _sigma0 = OPT->_aprSigCodeBias;
94 break;
95 }
96}
97
98// Destructor
99////////////////////////////////////////////////////////////////////////////
100t_pppParam::~t_pppParam() {
101 delete _ambInfo;
102}
103//
104////////////////////////////////////////////////////////////////////////////
105double t_pppParam::partial(const bncTime& /* epoTime */, const t_pppSatObs* obs,
106 const t_lc& obsLC) const {
107
108 // Special Case - Melbourne-Wuebbena
109 // ---------------------------------
110 if (obsLC._type == t_lc::MW && _type != amb) {
111 return 0.0;
112 }
113
114 // Special Case - GIM (single-differenced: A[ION_ref]=+1, A[ION_i]=-1)
115 // -------------------------------------------------------------------
116 if (obsLC._type == t_lc::GIM) {
117 if (_type == ion) {
118 if (obs->prn() == _prn) {
119 return -1.0;
120 }
121 if (_prn == _refPrn && _prn.system() == obs->prn().system()) {
122 return 1.0;
123 }
124 }
125 return 0.0;
126 }
127
128 const t_pppStation* sta = PPP_CLIENT->staRover();
129 ColumnVector rhoV = sta->xyzApr() - obs->xc().Rows(1,3);
130
131 map<t_frequency::type, double> codeCoeff;
132 map<t_frequency::type, double> phaseCoeff;
133 map<t_frequency::type, double> ionoCoeff;
134 obs->lcCoeff(obsLC, codeCoeff, phaseCoeff, ionoCoeff);
135
136 switch (_type) {
137 case crdX:
138 return (sta->xyzApr()[0] - obs->xc()[0]) / rhoV.NormFrobenius();
139 case crdY:
140 return (sta->xyzApr()[1] - obs->xc()[1]) / rhoV.NormFrobenius();
141 case crdZ:
142 return (sta->xyzApr()[2] - obs->xc()[2]) / rhoV.NormFrobenius();
143 case rClk:
144 return (_sys != obsLC.system() || obsLC.isGeometryFree()) ? 0.0 : 1.0;
145 case amb:
146 if (obs->prn() == _prn) {
147 if (obsLC == _LC) {
148 return (obs->lambda(obsLC));
149 }
150 else if (_LC._type == t_lc::phase) {
151 return obs->lambda(_LC) * phaseCoeff[_LC._frq1];
152 }
153 else if (obsLC._type == t_lc::phaseIF && _LC._type == t_lc::MW) {
154 return obs->lambda(obsLC) * obs->lambda(_LC) / obs->lambda(t_lc(t_lc::phase, _LC._frq2));
155 }
156 }
157 break;
158 case trp:
159 return 1.0 / sin(obs->eleSat());
160 case ion:
161 if (obs->prn() == _prn) {
162 if (obsLC._type == t_lc::code || obsLC._type == t_lc::phase) {
163 return ionoCoeff[obsLC._frq1];
164 }
165 }
166 break;
167 case bias:
168 return (_LC == obsLC ? 1.0 : 0.0);
169 }
170 return 0.0;
171}
172
173//
174////////////////////////////////////////////////////////////////////////////
175string t_pppParam::toString() const {
176 stringstream ss;
177 switch (_type) {
178 case crdX:
179 ss << "CRD_X";
180 break;
181 case crdY:
182 ss << "CRD_Y";
183 break;
184 case crdZ:
185 ss << "CRD_Z";
186 break;
187 case rClk:
188 ss << "REC_CLK " << _sys << " ";
189 break;
190 case trp:
191 ss << "TRP ";
192 break;
193 case amb:
194 ss << "AMB " << left << setw(3) << _LC.toString() << right << ' ' << _prn.toString();
195 break;
196 case ion:
197 ss << "ION " << left << setw(3) << _LC.toString() << right << ' ' << _prn.toString();
198 break;
199 case bias:
200 char sys = t_frequency::toSystem(_LC._frq1);
201 ss << "BIA " << left << setw(3) << _LC.toString() << right << ' ' << sys << " ";
202 break;
203 }
204 return ss.str();
205}
206
207// Constructor
208////////////////////////////////////////////////////////////////////////////
209t_pppParlist::t_pppParlist() {
210}
211
212// Destructor
213////////////////////////////////////////////////////////////////////////////
214t_pppParlist::~t_pppParlist() {
215
216 for (unsigned ii = 0; ii < _params.size(); ii++) {
217 delete _params[ii];
218 }
219}
220
221//
222////////////////////////////////////////////////////////////////////////////
223t_irc t_pppParlist::set(const bncTime& epoTime, const std::vector<t_pppSatObs*>& obsVector) {
224
225 // Remove some Parameters
226 // ----------------------
227 vector<t_pppParam*>::iterator it = _params.begin();
228 while (it != _params.end()) {
229 t_pppParam* par = *it;
230
231 bool remove = false;
232
233 if (par->epoSpec()) {
234 remove = true;
235 }
236
237 else if (par->type() == t_pppParam::amb ||
238 par->type() == t_pppParam::crdX ||
239 par->type() == t_pppParam::crdY ||
240 par->type() == t_pppParam::crdZ ||
241 par->type() == t_pppParam::ion ) {
242 if (par->lastObsTime().valid() && (epoTime - par->lastObsTime() > 60.0)) {
243 remove = true;
244 }
245 }
246
247 if (remove) {
248#ifdef BNC_DEBUG_PPP
249// LOG << "remove0 " << par->toString() << std::endl;
250#endif
251 delete par;
252 it = _params.erase(it);
253 }
254 else {
255 ++it;
256 }
257 }
258
259 // Check which systems have observations
260 // -------------------------------------
261 vector<char> systems = OPT->systems();
262 for (char sys : systems) {
263 _usedSystems[sys] = 0;
264 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
265 const t_pppSatObs* satObs = obsVector[jj];
266 if (satObs->prn().system() == sys) {
267 _usedSystems[sys]++;
268 }
269 }
270 };
271
272 // Check whether parameters have observations
273 // ------------------------------------------
274 for (unsigned ii = 0; ii < _params.size(); ii++) {
275 t_pppParam* par = _params[ii];
276 if (par->prn() == 0) {
277 par->setLastObsTime(epoTime);
278 if (par->firstObsTime().undef()) {
279 par->setFirstObsTime(epoTime);
280 }
281 }
282 else {
283 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
284 const t_pppSatObs* satObs = obsVector[jj];
285 if (satObs->prn() == par->prn()) {
286 par->setLastObsTime(epoTime);
287 if (par->firstObsTime().undef()) {
288 par->setFirstObsTime(epoTime);
289 }
290 break;
291 }
292 }
293 }
294 }
295
296 // Required Set of Parameters
297 // --------------------------
298 vector<t_pppParam*> required;
299
300 // Coordinates
301 // -----------
302 required.push_back(new t_pppParam(t_pppParam::crdX, t_prn()));
303 required.push_back(new t_pppParam(t_pppParam::crdY, t_prn()));
304 required.push_back(new t_pppParam(t_pppParam::crdZ, t_prn()));
305
306 // Receiver Clocks
307 // ---------------
308 for (const auto& [sys, numObs] : _usedSystems) {
309 if (numObs > 0) {
310 t_pppParam* clk = new t_pppParam(t_pppParam::rClk, t_prn());
311 clk->setSystem(sys);
312 required.push_back(clk);
313 }
314 }
315
316 // Troposphere
317 // -----------
318 if (OPT->estTrp()) {
319 required.push_back(new t_pppParam(t_pppParam::trp, t_prn()));
320 }
321
322 // Ionosphere
323 // ----------
324 if (OPT->_ionoModelType == OPT->est) {
325 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
326 const t_pppSatObs* satObs = obsVector[jj];
327 std::vector<t_lc> LCs = OPT->LCs(satObs->prn().system());
328 for (auto it = LCs.begin(); it != LCs.end(); ++it) {
329 const t_lc& lc = *it;
330 if (!lc.isIonoFree()) {
331 required.push_back(new t_pppParam(t_pppParam::ion, satObs->prn()));
332 break;
333 }
334 }
335 }
336 }
337
338 // Reference satellite PRN per system (from satellite marked by preparePseudoObs)
339 // -------------------------------------------------------------------------------
340 if (OPT->_pseudoObsIono) {
341 map<char, t_prn> refPrnMap;
342 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
343 if (obsVector[jj]->isReference()) {
344 refPrnMap[obsVector[jj]->prn().system()] = obsVector[jj]->prn();
345 }
346 }
347 for (unsigned ii = 0; ii < required.size(); ii++) {
348 if (required[ii]->type() == t_pppParam::ion) {
349 char sys = required[ii]->prn().system();
350 if (refPrnMap.count(sys)) {
351 required[ii]->setRefPrn(refPrnMap[sys]);
352 }
353 }
354 }
355 }
356
357 // Ambiguities
358 // -----------
359 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
360 const t_pppSatObs* satObs = obsVector[jj];
361 char sys = satObs->prn().system();
362 const vector<t_lc>& ambLCs = OPT->ambLCs(sys);
363 for (unsigned ii = 0; ii < ambLCs.size(); ii++) {
364 if (ambLCs[ii]._frq1 == t_frequency::G5 && !satObs->isValid(ambLCs[ii])) {
365 continue;
366 }
367 required.push_back(new t_pppParam(t_pppParam::amb, satObs->prn(), ambLCs[ii], &obsVector));
368 }
369 }
370
371 // Biases
372 // ------
373 int maxSkip = (OPT->_pseudoObsIono ? 1 : 2);
374 for (const auto& [sys, numObs] : _usedSystems) {
375 if (numObs > 0) {
376 bool ar = OPT->arSystem(sys);
377 int skip = 0;
378 vector<t_lc> LCs = OPT->LCs(sys);
379 for (const t_lc& lc : LCs) {
380 if (ar) {
381 if (skip < maxSkip && lc.includesPhase()) {
382 skip += 1;
383 }
384 else {
385 required.push_back(new t_pppParam(t_pppParam::bias, t_prn(), lc));
386 }
387 }
388 else {
389 if (lc.includesCode()) {
390 if (skip < maxSkip) {
391 skip += 1;
392 }
393 else {
394 required.push_back(new t_pppParam(t_pppParam::bias, t_prn(), lc));
395 }
396 }
397 }
398 }
399 }
400 }
401
402 // Check if all required parameters are present
403 // --------------------------------------------
404 for (unsigned ii = 0; ii < required.size(); ii++) {
405 t_pppParam* parReq = required[ii];
406
407 bool found = false;
408 for (unsigned jj = 0; jj < _params.size(); jj++) {
409 t_pppParam* parOld = _params[jj];
410 if (parOld->isEqual(parReq)) {
411 found = true;
412 break;
413 }
414 }
415 if (found) {
416 delete parReq;
417 }
418 else {
419 _params.push_back(parReq);
420 }
421 }
422
423 // Set Parameter Indices
424 // ---------------------
425 sort(_params.begin(), _params.end(), t_pppParam::sortFunction);
426
427 for (unsigned ii = 0; ii < _params.size(); ii++) {
428 t_pppParam* par = _params[ii];
429 par->setIndex(ii);
430 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
431 const t_pppSatObs* satObs = obsVector[jj];
432 if (satObs->prn() == par->prn()) {
433 par->setAmbEleSat(satObs->eleSat());
434 par->stepAmbNumEpo();
435 }
436 }
437 }
438
439 return success;
440}
441
442//
443////////////////////////////////////////////////////////////////////////////
444void t_pppParlist::printResult(const bncTime& epoTime, const SymmetricMatrix& QQ,
445 const ColumnVector& xx, double fixRatio) const {
446
447 string epoTimeStr = string(epoTime);
448 const t_pppStation* sta = PPP_CLIENT->staRover();
449
450 LOG << endl;
451
452 t_pppParam* parX = 0;
453 t_pppParam* parY = 0;
454 t_pppParam* parZ = 0;
455 for (unsigned ii = 0; ii < _params.size(); ii++) {
456 t_pppParam* par = _params[ii];
457 if (par->type() == t_pppParam::crdX) {
458 parX = par;
459 }
460 else if (par->type() == t_pppParam::crdY) {
461 parY = par;
462 }
463 else if (par->type() == t_pppParam::crdZ) {
464 parZ = par;
465 }
466 else {
467 int ind = par->indexNew();
468 double apr = (par->type() == t_pppParam::trp) ?
469 t_tropo::delay_saast(sta->xyzApr(), M_PI/2.0) : par->x0();
470 LOG << epoTimeStr << ' ' << par->toString() << ' '
471 << setw(10) << setprecision(4) << apr << ' '
472 << showpos << setw(10) << setprecision(4) << xx[ind] << noshowpos << " +- "
473 << setw(8) << setprecision(4) << sqrt(QQ[ind][ind]);
474 if (par->type() == t_pppParam::amb) {
475 LOG << " el = " << setw(6) << setprecision(2) << par->ambEleSat() * 180.0 / M_PI
476 << " epo = " << setw(4) << par->ambNumEpo();
477 }
478 LOG << endl;
479 }
480 }
481
482 if (parX && parY && parZ) {
483
484 ColumnVector xyz(3);
485 xyz[0] = xx[parX->indexNew()];
486 xyz[1] = xx[parY->indexNew()];
487 xyz[2] = xx[parZ->indexNew()];
488
489 ColumnVector neu(3);
490 xyz2neu(sta->ellApr().data(), xyz.data(), neu.data());
491
492 SymmetricMatrix QQxyz = QQ.SymSubMatrix(1,3);
493
494 SymmetricMatrix QQneu(3);
495 covariXYZ_NEU(QQxyz, sta->ellApr().data(), QQneu);
496
497 LOG << epoTimeStr << ' ' << sta->name()
498 << " X = " << setprecision(4) << sta->xyzApr()[0] + xyz[0] << " +- "
499 << setprecision(4) << sqrt(QQxyz[0][0])
500
501 << " Y = " << setprecision(4) << sta->xyzApr()[1] + xyz[1] << " +- "
502 << setprecision(4) << sqrt(QQxyz[1][1])
503
504 << " Z = " << setprecision(4) << sta->xyzApr()[2] + xyz[2] << " +- "
505 << setprecision(4) << sqrt(QQxyz[2][2])
506
507 << " dN = " << setprecision(4) << neu[0] << " +- "
508 << setprecision(4) << sqrt(QQneu[0][0])
509
510 << " dE = " << setprecision(4) << neu[1] << " +- "
511 << setprecision(4) << sqrt(QQneu[1][1])
512
513 << " dU = " << setprecision(4) << neu[2] << " +- "
514 << setprecision(4) << sqrt(QQneu[2][2]);
515
516 if (fixRatio > 0.0) {
517 LOG << " fix " << int(100*fixRatio) << " %";
518 }
519 else {
520 LOG << " flt ";
521 }
522
523 LOG << endl;
524 }
525 return;
526}
527
528//
529////////////////////////////////////////////////////////////////////////////
530void t_pppParlist::printParams(const bncTime& epoTime) {
531
532 for (unsigned iPar = 0; iPar < _params.size(); iPar++) {
533 LOG << _params[iPar]->toString()
534 << "\t lastObsTime().valid() \t" << _params[iPar]->lastObsTime().valid()
535 << "\t epoTime - par->lastObsTime() \t" << (epoTime - _params[iPar]->lastObsTime())
536 << endl;
537 }
538}
539
Note: See TracBrowser for help on using the repository browser.