source: ntrip/trunk/BNC/src/combination/bnccomb.cpp@ 9710

Last change on this file since 9710 was 9710, checked in by stuerze, 7 months ago

minor changes

File size: 41.5 KB
Line 
1/* -------------------------------------------------------------------------
2 * BKG NTRIP Client
3 * -------------------------------------------------------------------------
4 *
5 * Class: bncComb
6 *
7 * Purpose: Combinations of Orbit/Clock Corrections
8 *
9 * Author: L. Mervart
10 *
11 * Created: 22-Jan-2011
12 *
13 * Changes:
14 *
15 * -----------------------------------------------------------------------*/
16
17#include <newmatio.h>
18#include <iomanip>
19#include <sstream>
20#include <map>
21
22#include "bnccomb.h"
23#include <combination/bncbiassnx.h>
24#include "bnccore.h"
25#include "upload/bncrtnetdecoder.h"
26#include "bncsettings.h"
27#include "bncutils.h"
28#include "bncsp3.h"
29#include "bncantex.h"
30#include "t_prn.h"
31
32const double sig0_offAC = 1000.0;
33const double sig0_offACSat = 100.0;
34const double sigP_offACSat = 0.01;
35const double sig0_clkSat = 100.0;
36
37const double sigObs = 0.05;
38
39using namespace std;
40
41// Constructor
42////////////////////////////////////////////////////////////////////////////
43bncComb::cmbParam::cmbParam(parType type_, int index_, const QString& ac_, const QString& prn_) {
44
45 type = type_;
46 index = index_;
47 AC = ac_;
48 prn = prn_;
49 xx = 0.0;
50 eph = 0;
51
52 if (type == offACgnss) {
53 epoSpec = true;
54 sig0 = sig0_offAC;
55 sigP = sig0;
56 }
57 else if (type == offACSat) {
58 epoSpec = false;
59 sig0 = sig0_offACSat;
60 sigP = sigP_offACSat;
61 }
62 else if (type == clkSat) {
63 epoSpec = true;
64 sig0 = sig0_clkSat;
65 sigP = sig0;
66 }
67}
68
69// Destructor
70////////////////////////////////////////////////////////////////////////////
71bncComb::cmbParam::~cmbParam() {
72}
73
74// Partial
75////////////////////////////////////////////////////////////////////////////
76double bncComb::cmbParam::partial(char sys, const QString& AC_, const QString& prn_) {
77
78 if (type == offACgnss) {
79 if (AC == AC_ && prn_[0].toLatin1() == sys) {
80 return 1.0;
81 }
82 }
83 else if (type == offACSat) {
84 if (AC == AC_ && prn == prn_) {
85 return 1.0;
86 }
87 }
88 else if (type == clkSat) {
89 if (prn == prn_) {
90 return 1.0;
91 }
92 }
93
94 return 0.0;
95}
96
97//
98////////////////////////////////////////////////////////////////////////////
99QString bncComb::cmbParam::toString(char sys) const {
100
101 QString outStr;
102
103 if (type == offACgnss) {
104 outStr = "AC Offset " + AC + " " + sys + " ";
105 }
106 else if (type == offACSat) {
107 outStr = "Sat Offset " + AC + " " + prn.mid(0,3);
108 }
109 else if (type == clkSat) {
110 outStr = "Clk Corr " + prn.mid(0,3);
111 }
112
113 return outStr;
114}
115
116// Singleton
117////////////////////////////////////////////////////////////////////////////
118bncComb* bncComb::instance() {
119 static bncComb _bncComb;
120 return &_bncComb;
121}
122
123// Constructor
124////////////////////////////////////////////////////////////////////////////
125bncComb::bncComb() : _ephUser(true) {
126
127 bncSettings settings;
128
129 QStringList cmbStreams = settings.value("cmbStreams").toStringList();
130
131 _cmbSampl = settings.value("cmbSampl").toInt();
132 if (_cmbSampl <= 0) {
133 _cmbSampl = 5;
134 }
135 _useGps = (Qt::CheckState(settings.value("cmbGps").toInt()) == Qt::Checked) ? true : false;
136 if (_useGps) {
137 _cmbSysPrn['G'] = t_prn::MAXPRN_GPS;
138 _masterMissingEpochs['G'] = 0;
139 }
140 _useGlo = (Qt::CheckState(settings.value("cmbGlo").toInt()) == Qt::Checked) ? true : false;
141 if (_useGlo) {
142 _cmbSysPrn['R'] = t_prn::MAXPRN_GLONASS;
143 _masterMissingEpochs['R'] = 0;
144 }
145 _useGal = (Qt::CheckState(settings.value("cmbGal").toInt()) == Qt::Checked) ? true : false;
146 if (_useGal) {
147 _cmbSysPrn['E'] = t_prn::MAXPRN_GALILEO;
148 _masterMissingEpochs['E'] = 0;
149 }
150 _useBds = (Qt::CheckState(settings.value("cmbBds").toInt()) == Qt::Checked) ? true : false;
151 if (_useBds) {
152 _cmbSysPrn['C'] = t_prn::MAXPRN_BDS;
153 _masterMissingEpochs['C'] = 0;
154 }
155 _useQzss = (Qt::CheckState(settings.value("cmbQzss").toInt()) == Qt::Checked) ? true : false;
156 if (_useQzss) {
157 _cmbSysPrn['J'] = t_prn::MAXPRN_QZSS;
158 _masterMissingEpochs['J'] = 0;
159 }
160 _useSbas = (Qt::CheckState(settings.value("cmbSbas").toInt()) == Qt::Checked) ? true : false;
161 if (_useSbas) {
162 _cmbSysPrn['S'] = t_prn::MAXPRN_SBAS;
163 _masterMissingEpochs['S'] = 0;
164 }
165 _useIrnss = (Qt::CheckState(settings.value("cmbIrnss").toInt()) == Qt::Checked) ? true : false;
166 if (_useIrnss) {
167 _cmbSysPrn['I'] = t_prn::MAXPRN_IRNSS;
168 _masterMissingEpochs['I'] = 0;
169 }
170
171 if (cmbStreams.size() >= 1 && !cmbStreams[0].isEmpty()) {
172 QListIterator<QString> it(cmbStreams);
173 while (it.hasNext()) {
174 QStringList hlp = it.next().split(" ");
175 cmbAC* newAC = new cmbAC();
176 newAC->mountPoint = hlp[0];
177 newAC->name = hlp[1];
178 newAC->weight = hlp[2].toDouble();
179 QMapIterator<char, unsigned> itSys(_cmbSysPrn);
180 // init
181 while (itSys.hasNext()) {
182 itSys.next();
183 char sys = itSys.key();
184 if (_masterOrbitAC.isEmpty()) {
185 _masterOrbitAC[sys] = newAC->name;
186 }
187 }
188 _ACs.append(newAC);
189 }
190 }
191
192 QString ssrFormat;
193 _ssrCorr = 0;
194 QListIterator<QString> it(settings.value("uploadMountpointsOut").toStringList());
195 while (it.hasNext()) {
196 QStringList hlp = it.next().split(",");
197 if (hlp.size() > 7) {
198 ssrFormat = hlp[7];
199 }
200 }
201 if (ssrFormat == "IGS-SSR") {
202 _ssrCorr = new SsrCorrIgs();
203 }
204 else if (ssrFormat == "RTCM-SSR") {
205 _ssrCorr = new SsrCorrRtcm();
206 }
207 else { // default
208 _ssrCorr = new SsrCorrIgs();
209 }
210
211 _rtnetDecoder = 0;
212
213 connect(this, SIGNAL(newMessage(QByteArray,bool)),
214 BNC_CORE, SLOT(slotMessage(const QByteArray,bool)));
215
216 connect(BNC_CORE, SIGNAL(providerIDChanged(QString)),
217 this, SLOT(slotProviderIDChanged(QString)));
218
219 connect(BNC_CORE, SIGNAL(newOrbCorrections(QList<t_orbCorr>)),
220 this, SLOT(slotNewOrbCorrections(QList<t_orbCorr>)));
221
222 connect(BNC_CORE, SIGNAL(newClkCorrections(QList<t_clkCorr>)),
223 this, SLOT(slotNewClkCorrections(QList<t_clkCorr>)));
224
225 connect(BNC_CORE, SIGNAL(newCodeBiases(QList<t_satCodeBias>)),
226 this, SLOT(slotNewCodeBiases(QList<t_satCodeBias>)));
227
228 // Combination Method
229 // ------------------
230 if (settings.value("cmbMethod").toString() == "Single-Epoch") {
231 _method = singleEpoch;
232 }
233 else {
234 _method = filter;
235 }
236
237 // Initialize Parameters (model: Clk_Corr = AC_Offset + Sat_Offset + Clk)
238 // ----------------------------------------------------------------------
239 if (_method == filter) {
240 // SYSTEM
241 QMapIterator<char, unsigned> itSys(_cmbSysPrn);
242 while (itSys.hasNext()) {
243 itSys.next();
244 int nextPar = 0;
245 char sys = itSys.key();
246 unsigned maxPrn = itSys.value();
247 unsigned flag = 0;
248 if (sys == 'E') {
249 flag = 1;
250 }
251 // AC
252 QListIterator<cmbAC*> itAc(_ACs);
253 while (itAc.hasNext()) {
254 cmbAC* AC = itAc.next();
255 _params[sys].push_back(new cmbParam(cmbParam::offACgnss, ++nextPar, AC->name, ""));
256 for (unsigned iGnss = 1; iGnss <= maxPrn; iGnss++) {
257 QString prn = QString("%1%2_%3").arg(sys).arg(iGnss, 2, 10, QChar('0')).arg(flag);
258 _params[sys].push_back(new cmbParam(cmbParam::offACSat, ++nextPar, AC->name, prn));
259 }
260 }
261 for (unsigned iGnss = 1; iGnss <= maxPrn; iGnss++) {
262 QString prn = QString("%1%2_%3").arg(sys).arg(iGnss, 2, 10, QChar('0')).arg(flag);
263 _params[sys].push_back(new cmbParam(cmbParam::clkSat, ++nextPar, "", prn));
264 }
265 // Initialize Variance-Covariance Matrix
266 // -------------------------------------
267 _QQ[sys].ReSize(_params[sys].size());
268 _QQ[sys] = 0.0;
269 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
270 cmbParam* pp = _params[sys][iPar-1];
271 _QQ[sys](iPar,iPar) = pp->sig0 * pp->sig0;
272 }
273 }
274 }
275
276 // ANTEX File
277 // ----------
278 _antex = 0;
279 QString antexFileName = settings.value("uploadAntexFile").toString();
280 if (!antexFileName.isEmpty()) {
281 _antex = new bncAntex();
282 if (_antex->readFile(antexFileName) != success) {
283 emit newMessage("bncCmb: wrong ANTEX file", true);
284 delete _antex;
285 _antex = 0;
286 }
287 }
288
289
290 // Bias SINEX File
291 // ---------------
292 _bsx = 0;
293 QString bsxFileName = settings.value("cmbBsxFile").toString();
294 if (!bsxFileName.isEmpty()) {
295 _bsx = new bncBiasSnx();
296 if (_bsx->readFile(bsxFileName) != success) {
297 emit newMessage("bncComb: wrong Bias SINEX file", true);
298 delete _bsx;
299 _bsx = 0;
300 }
301 }
302 if (_bsx) {
303 _ms = 6.0 * 3600 * 1000.0;
304 QTimer::singleShot(_ms, this, SLOT(slotReadBiasSnxFile()));
305 }
306
307
308 // Maximal Residuum
309 // ----------------
310 _MAXRES = settings.value("cmbMaxres").toDouble();
311 if (_MAXRES <= 0.0) {
312 _MAXRES = 999.0;
313 }
314}
315
316// Destructor
317////////////////////////////////////////////////////////////////////////////
318bncComb::~bncComb() {
319 QListIterator<cmbAC*> icAC(_ACs);
320 while (icAC.hasNext()) {
321 delete icAC.next();
322 }
323 delete _rtnetDecoder;
324 if (_ssrCorr) {
325 delete _ssrCorr;
326 }
327 delete _antex;
328 delete _bsx;
329 QMapIterator<char, unsigned> itSys(_cmbSysPrn);
330 while (itSys.hasNext()) {
331 itSys.next();
332 char sys = itSys.key();
333 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
334 delete _params[sys][iPar-1];
335 }
336 QListIterator<bncTime> itTime(_buffer[sys].keys());
337 while (itTime.hasNext()) {
338 bncTime epoTime = itTime.next();
339 _buffer[sys].remove(epoTime);
340 }
341 }
342
343}
344
345void bncComb::slotReadBiasSnxFile() {
346 bncSettings settings;
347 QString bsxFileName = settings.value("cmbBsxFile").toString();
348
349 if (bsxFileName.isEmpty()) {
350 emit newMessage("bncComb: no Bias SINEX file specified", true);
351 return;
352 }
353 if (!_bsx) {
354 _bsx = new bncBiasSnx();
355 }
356 if (_bsx->readFile(bsxFileName) != success) {
357 emit newMessage("bncComb: wrong Bias SINEX file", true);
358 delete _bsx;
359 _bsx = 0;
360 }
361 else {
362 emit newMessage("bncComb: successfully read Bias SINEX file", true);
363 }
364
365 QTimer::singleShot(_ms, this, SLOT(slotReadBiasSnxFile()));
366}
367
368
369// Remember orbit corrections
370////////////////////////////////////////////////////////////////////////////
371void bncComb::slotNewOrbCorrections(QList<t_orbCorr> orbCorrections) {
372 QMutexLocker locker(&_mutex);
373
374 for (int ii = 0; ii < orbCorrections.size(); ii++) {
375 t_orbCorr& orbCorr = orbCorrections[ii];
376 QString staID(orbCorr._staID.c_str());
377 char sys = orbCorr._prn.system();
378
379 if (!_cmbSysPrn.contains(sys)){
380 continue;
381 }
382
383 // Find/Check the AC Name
384 // ----------------------
385 QString acName;
386 QListIterator<cmbAC*> icAC(_ACs);
387 while (icAC.hasNext()) {
388 cmbAC* AC = icAC.next();
389 if (AC->mountPoint == staID) {
390 acName = AC->name;
391 break;
392 }
393 }
394 if (acName.isEmpty()) {
395 continue;
396 }
397
398 // Store the correction
399 // --------------------
400 QMap<t_prn, t_orbCorr>& storage = _orbCorrections[acName];
401 storage[orbCorr._prn] = orbCorr;
402 }
403}
404
405// Remember Satellite Code Biases
406////////////////////////////////////////////////////////////////////////////
407void bncComb::slotNewCodeBiases(QList<t_satCodeBias> satCodeBiases) {
408 QMutexLocker locker(&_mutex);
409
410 for (int ii = 0; ii < satCodeBiases.size(); ii++) {
411 t_satCodeBias& satCodeBias = satCodeBiases[ii];
412 QString staID(satCodeBias._staID.c_str());
413 char sys = satCodeBias._prn.system();
414
415 if (!_cmbSysPrn.contains(sys)){
416 continue;
417 }
418
419 // Find/Check the AC Name
420 // ----------------------
421 QString acName;
422 QListIterator<cmbAC*> icAC(_ACs);
423 while (icAC.hasNext()) {
424 cmbAC* AC = icAC.next();
425 if (AC->mountPoint == staID) {
426 acName = AC->name;
427 break;
428 }
429 }
430 if (acName.isEmpty()) {
431 continue;
432 }
433
434 // Store the correction
435 // --------------------
436 QMap<t_prn, t_satCodeBias>& storage = _satCodeBiases[acName];
437 storage[satCodeBias._prn] = satCodeBias;
438 }
439}
440
441// Process clock corrections
442////////////////////////////////////////////////////////////////////////////
443void bncComb::slotNewClkCorrections(QList<t_clkCorr> clkCorrections) {
444 QMutexLocker locker(&_mutex);
445
446 bncTime lastTime;
447
448 for (int ii = 0; ii < clkCorrections.size(); ii++) {
449 t_clkCorr& clkCorr = clkCorrections[ii];
450 QString staID(clkCorr._staID.c_str());
451 QString prn(clkCorr._prn.toInternalString().c_str());
452 char sys = clkCorr._prn.system();
453
454 if (!_cmbSysPrn.contains(sys)){
455 continue;
456 }
457
458 // Set the last time
459 // -----------------
460 if (lastTime.undef() || clkCorr._time > lastTime) {
461 lastTime = clkCorr._time;
462 }
463
464 // Find/Check the AC Name
465 // ----------------------
466 QString acName;
467 QListIterator<cmbAC*> icAC(_ACs);
468 while (icAC.hasNext()) {
469 cmbAC* AC = icAC.next();
470 if (AC->mountPoint == staID) {
471 acName = AC->name;
472 break;
473 }
474 }
475 if (acName.isEmpty()) {
476 continue;
477 }
478
479 // Check Modulo Time
480 // -----------------
481 int sec = int(nint(clkCorr._time.gpssec()*10));
482 if (sec % (_cmbSampl*10) != 0.0) {
483 continue;
484 }
485
486 // Check Correction Age
487 // --------------------
488 if (_resTime.valid() && clkCorr._time <= _resTime) {
489#ifdef BNC_DEBUG_CMB
490 emit newMessage("bncComb: old correction: " + acName.toLatin1() + " " + prn.mid(0,3).toLatin1(), true);
491#endif
492 continue;
493 }
494
495 // Create new correction
496 // ---------------------
497 cmbCorr* newCorr = new cmbCorr();
498 newCorr->_prn = prn;
499 newCorr->_time = clkCorr._time;
500 newCorr->_iod = clkCorr._iod;
501 newCorr->_acName = acName;
502 newCorr->_clkCorr = clkCorr;
503
504 // Check orbit correction
505 // ----------------------
506 if (!_orbCorrections.contains(acName)) {
507 delete newCorr;
508 continue;
509 }
510 else {
511 QMap<t_prn, t_orbCorr>& storage = _orbCorrections[acName];
512 if (!storage.contains(clkCorr._prn) ||
513 storage[clkCorr._prn]._iod != newCorr->_iod) {
514 delete newCorr;
515 continue;
516 }
517 else {
518 newCorr->_orbCorr = storage[clkCorr._prn];
519 }
520 }
521
522 // Check the Ephemeris
523 //--------------------
524 t_eph* ephLast = _ephUser.ephLast(prn);
525 t_eph* ephPrev = _ephUser.ephPrev(prn);
526 if (ephLast == 0) {
527#ifdef BNC_DEBUG_CMB
528 emit newMessage("bncComb: eph not found for " + prn.mid(0,3).toLatin1(), true);
529#endif
530 delete newCorr;
531 continue;
532 }
533 else if (ephLast->checkState() == t_eph::bad ||
534 ephLast->checkState() == t_eph::outdated ||
535 ephLast->checkState() == t_eph::unhealthy) {
536#ifdef BNC_DEBUG_CMB
537 emit newMessage("bncComb: ephLast not ok (checkState: " + ephLast->checkStateToString().toLatin1() + ") for " + prn.mid(0,3).toLatin1(), true);
538#endif
539 delete newCorr;
540 continue;
541 }
542 else {
543 if (ephLast->IOD() == newCorr->_iod) {
544 newCorr->_eph = ephLast;
545 }
546 else if (ephPrev && ephPrev->checkState() == t_eph::ok &&
547 ephPrev->IOD() == newCorr->_iod) {
548 newCorr->_eph = ephPrev;
549 switchToLastEph(ephLast, newCorr);
550 }
551 else {
552#ifdef BNC_DEBUG_CMB
553 emit newMessage("bncComb: eph not found for " + prn.mid(0,3).toLatin1() +
554 QString(" with IOD %1").arg(newCorr->_iod).toLatin1(), true);
555#endif
556 delete newCorr;
557 continue;
558 }
559 }
560
561 // Check satellite code biases
562 // ----------------------------
563 if (!_satCodeBiases.contains(acName)) {
564 delete newCorr;
565 continue;
566 }
567 else {
568 QMap<t_prn, t_satCodeBias>& storage = _satCodeBiases[acName];
569 if (!storage.contains(clkCorr._prn)) {
570 delete newCorr;
571 continue;
572 }
573 else {
574 newCorr->_satCodeBias = storage[clkCorr._prn];
575 QMap<t_frequency::type, double> codeBiasesRefSig;
576 for (unsigned ii = 1; ii < cmbRefSig::cIF; ii++) {
577 t_frequency::type frqType = cmbRefSig::toFreq(sys, static_cast<cmbRefSig::type>(ii));
578 char frqNum = t_frequency::toString(frqType)[1];
579 char attrib = cmbRefSig::toAttrib(sys, static_cast<cmbRefSig::type>(ii));
580 QString rnxType2ch = QString("%1%2").arg(frqNum).arg(attrib);
581 for (unsigned ii = 0; ii < newCorr->_satCodeBias._bias.size(); ii++) {
582 const t_frqCodeBias& bias = newCorr->_satCodeBias._bias[ii];
583 if (rnxType2ch.toStdString() == bias._rnxType2ch) {
584 codeBiasesRefSig[frqType] = bias._value;
585 }
586 }
587 }
588 if (codeBiasesRefSig.size() < 2) {
589 delete newCorr;
590 continue;
591 }
592 map<t_frequency::type, double> codeCoeff;
593 double channel = double(newCorr->_eph->slotNum());
594 cmbRefSig::coeff(sys, cmbRefSig::cIF, channel, codeCoeff);
595 map<t_frequency::type, double>::const_iterator it;
596 for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
597 t_frequency::type frqType = it->first;
598 newCorr->_codeBiasIF += it->second * codeBiasesRefSig[frqType];
599 }
600 newCorr->_satCodeBias._bias.clear();
601 }
602 }
603
604 // Store correction into the buffer
605 // --------------------------------
606 QVector<cmbCorr*>& corrs = _buffer[sys][newCorr->_time].corrs;
607 corrs.push_back(newCorr);
608 }
609
610 // Process previous Epoch(s)
611 // -------------------------
612 const double outWait = 1.0 * _cmbSampl;
613 QMapIterator<char, unsigned> itSys(_cmbSysPrn);
614 while (itSys.hasNext()) {
615 itSys.next();
616 char sys = itSys.key();
617 QListIterator<bncTime> itTime(_buffer[sys].keys());
618 while (itTime.hasNext()) {
619 bncTime epoTime = itTime.next();
620 if (epoTime < lastTime - outWait) {
621 _resTime = epoTime;
622 processEpoch(sys);
623 }
624 }
625 }
626}
627
628// Change the correction so that it refers to last received ephemeris
629////////////////////////////////////////////////////////////////////////////
630void bncComb::switchToLastEph(t_eph* lastEph, cmbCorr* corr) {
631
632 if (corr->_eph == lastEph) {
633 return;
634 }
635
636 ColumnVector oldXC(6);
637 ColumnVector oldVV(3);
638 if (corr->_eph->getCrd(corr->_time, oldXC, oldVV, false) != success) {
639 return;
640 }
641
642 ColumnVector newXC(6);
643 ColumnVector newVV(3);
644 if (lastEph->getCrd(corr->_time, newXC, newVV, false) != success) {
645 return;
646 }
647
648 ColumnVector dX = newXC.Rows(1,3) - oldXC.Rows(1,3);
649 ColumnVector dV = newVV - oldVV;
650 double dC = newXC(4) - oldXC(4);
651
652 ColumnVector dRAO(3);
653 XYZ_to_RSW(newXC.Rows(1,3), newVV, dX, dRAO);
654
655 ColumnVector dDotRAO(3);
656 XYZ_to_RSW(newXC.Rows(1,3), newVV, dV, dDotRAO);
657
658 QString msg = "switch corr " + corr->_prn.mid(0,3)
659 + QString(" %1 -> %2 %3").arg(corr->_iod,3).arg(lastEph->IOD(),3).arg(dC*t_CST::c, 8, 'f', 4);
660
661 emit newMessage(msg.toLatin1(), false);
662
663 corr->_iod = lastEph->IOD();
664 corr->_eph = lastEph;
665
666 corr->_orbCorr._xr += dRAO;
667 corr->_orbCorr._dotXr += dDotRAO;
668 corr->_clkCorr._dClk -= dC;
669}
670
671// Process Epoch
672////////////////////////////////////////////////////////////////////////////
673void bncComb::processEpoch(char sys) {
674
675 _log.clear();
676
677 QTextStream out(&_log, QIODevice::WriteOnly);
678
679 out << "\n" << "Combination: " << sys << "\n"
680 << "--------------------------------" << "\n";
681
682 // Observation Statistics
683 // ----------------------
684 bool masterPresent = false;
685 QListIterator<cmbAC*> icAC(_ACs);
686 while (icAC.hasNext()) {
687 cmbAC* AC = icAC.next();
688 AC->numObs[sys] = 0;
689 QVectorIterator<cmbCorr*> itCorr(corrs(sys));
690 while (itCorr.hasNext()) {
691 cmbCorr* corr = itCorr.next();
692 if (corr->_acName == AC->name) {
693 AC->numObs[sys] += 1;
694 if (AC->name == _masterOrbitAC[sys]) {
695 masterPresent = true;
696 }
697 }
698 }
699 out << AC->name.toLatin1().data() << ": " << AC->numObs[sys] << "\n";
700 }
701
702 // If Master not present, switch to another one
703 // --------------------------------------------
704 const unsigned switchMasterAfterGap = 1;
705 if (masterPresent) {
706 _masterMissingEpochs[sys] = 0;
707 }
708 else {
709 ++_masterMissingEpochs[sys];
710 if (_masterMissingEpochs[sys] < switchMasterAfterGap) {
711 out << "Missing Master, Epoch skipped" << "\n";
712 _buffer[sys].remove(_resTime);
713 emit newMessage(_log, false);
714 return;
715 }
716 else {
717 _masterMissingEpochs[sys] = 0;
718 QListIterator<cmbAC*> icAC(_ACs);
719 while (icAC.hasNext()) {
720 cmbAC* AC = icAC.next();
721 if (AC->numObs[sys] > 0) {
722 out << "Switching Master AC "
723 << _masterOrbitAC[sys].toLatin1().data() << " --> "
724 << AC->name.toLatin1().data() << " "
725 << _resTime.datestr().c_str() << " "
726 << _resTime.timestr().c_str() << "\n";
727 _masterOrbitAC[sys] = AC->name;
728 break;
729 }
730 }
731 }
732 }
733
734 QMap<QString, cmbCorr*> resCorr;
735
736 // Perform the actual Combination using selected Method
737 // ----------------------------------------------------
738 t_irc irc;
739 ColumnVector dx;
740 if (_method == filter) {
741 irc = processEpoch_filter(sys, out, resCorr, dx);
742 }
743 else {
744 irc = processEpoch_singleEpoch(sys, out, resCorr, dx);
745 }
746
747 // Update Parameter Values, Print Results
748 // --------------------------------------
749 if (irc == success) {
750 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
751 cmbParam* pp = _params[sys][iPar-1];
752 pp->xx += dx(iPar);
753 if (pp->type == cmbParam::clkSat) {
754 if (resCorr.find(pp->prn) != resCorr.end()) {
755 // set clock result
756 resCorr[pp->prn]->_dClkResult = pp->xx / t_CST::c;
757 // Add Code Biases from SINEX File
758 if (_bsx) {
759 map<t_frequency::type, double> codeCoeff;
760 double channel = double(resCorr[pp->prn]->_eph->slotNum());
761 cmbRefSig::coeff(sys, cmbRefSig::cIF, channel, codeCoeff);
762 t_frequency::type fType1 = cmbRefSig::toFreq(sys, cmbRefSig::c1);
763 t_frequency::type fType2 = cmbRefSig::toFreq(sys, cmbRefSig::c2);
764 _bsx->determineSsrSatCodeBiases(pp->prn.mid(0,3), codeCoeff[fType1], codeCoeff[fType2], resCorr[pp->prn]->_satCodeBias);
765 }
766 }
767 }
768 out << _resTime.datestr().c_str() << " "
769 << _resTime.timestr().c_str() << " ";
770 out.setRealNumberNotation(QTextStream::FixedNotation);
771 out.setFieldWidth(8);
772 out.setRealNumberPrecision(4);
773 out << pp->toString(sys) << " "
774 << pp->xx << " +- " << sqrt(_QQ[sys](pp->index,pp->index)) << "\n";
775 out.setFieldWidth(0);
776 }
777 printResults(out, resCorr);
778 dumpResults(resCorr);
779 }
780
781 // Delete Data, emit Message
782 // -------------------------
783 _buffer[sys].remove(_resTime);
784 emit newMessage(_log, false);
785}
786
787// Process Epoch - Filter Method
788////////////////////////////////////////////////////////////////////////////
789t_irc bncComb::processEpoch_filter(char sys, QTextStream& out,
790 QMap<QString, cmbCorr*>& resCorr,
791 ColumnVector& dx) {
792
793 // Prediction Step
794 // ---------------
795 int nPar = _params[sys].size();
796 ColumnVector x0(nPar);
797 for (int iPar = 1; iPar <= nPar; iPar++) {
798 cmbParam* pp = _params[sys][iPar-1];
799 if (pp->epoSpec) {
800 pp->xx = 0.0;
801 _QQ[sys].Row(iPar) = 0.0;
802 _QQ[sys].Column(iPar) = 0.0;
803 _QQ[sys](iPar,iPar) = pp->sig0 * pp->sig0;
804 }
805 else {
806 _QQ[sys](iPar,iPar) += pp->sigP * pp->sigP;
807 }
808 x0(iPar) = pp->xx;
809 }
810
811 // Check Satellite Positions for Outliers
812 // --------------------------------------
813 if (checkOrbits(sys, out) != success) {
814 return failure;
815 }
816
817 // Update and outlier detection loop
818 // ---------------------------------
819 SymmetricMatrix QQ_sav = _QQ[sys];
820 while (true) {
821
822 Matrix AA;
823 ColumnVector ll;
824 DiagonalMatrix PP;
825
826 if (createAmat(sys, AA, ll, PP, x0, resCorr) != success) {
827 return failure;
828 }
829
830 dx.ReSize(nPar); dx = 0.0;
831 kalman(AA, ll, PP, _QQ[sys], dx);
832
833 ColumnVector vv = ll - AA * dx;
834
835 int maxResIndex;
836 double maxRes = vv.maximum_absolute_value1(maxResIndex);
837 out.setRealNumberNotation(QTextStream::FixedNotation);
838 out.setRealNumberPrecision(3);
839 out << _resTime.datestr().c_str() << " " << _resTime.timestr().c_str()
840 << " Maximum Residuum " << maxRes << ' '
841 << corrs(sys)[maxResIndex-1]->_acName << ' ' << corrs(sys)[maxResIndex-1]->_prn.mid(0,3);
842 if (maxRes > _MAXRES) {
843 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
844 cmbParam* pp = _params[sys][iPar-1];
845 if (pp->type == cmbParam::offACSat &&
846 pp->AC == corrs(sys)[maxResIndex-1]->_acName &&
847 pp->prn == corrs(sys)[maxResIndex-1]->_prn.mid(0,3)) {
848 QQ_sav.Row(iPar) = 0.0;
849 QQ_sav.Column(iPar) = 0.0;
850 QQ_sav(iPar,iPar) = pp->sig0 * pp->sig0;
851 }
852 }
853
854 out << " Outlier" << "\n";
855 _QQ[sys] = QQ_sav;
856 corrs(sys).remove(maxResIndex-1);
857 }
858 else {
859 out << " OK" << "\n";
860 out.setRealNumberNotation(QTextStream::FixedNotation);
861 out.setRealNumberPrecision(4);
862 for (int ii = 0; ii < corrs(sys).size(); ii++) {
863 const cmbCorr* corr = corrs(sys)[ii];
864 out << _resTime.datestr().c_str() << ' '
865 << _resTime.timestr().c_str() << " "
866 << corr->_acName << ' ' << corr->_prn.mid(0,3);
867 out.setFieldWidth(10);
868 out << " res = " << vv[ii] << "\n";
869 out.setFieldWidth(0);
870 }
871 break;
872 }
873 }
874
875 return success;
876}
877
878// Print results
879////////////////////////////////////////////////////////////////////////////
880void bncComb::printResults(QTextStream& out,
881 const QMap<QString, cmbCorr*>& resCorr) {
882
883 QMapIterator<QString, cmbCorr*> it(resCorr);
884 while (it.hasNext()) {
885 it.next();
886 cmbCorr* corr = it.value();
887 const t_eph* eph = corr->_eph;
888 if (eph) {
889 ColumnVector xc(6);
890 ColumnVector vv(3);
891 if (eph->getCrd(_resTime, xc, vv, false) != success) {
892 continue;
893 }
894 out << _resTime.datestr().c_str() << " "
895 << _resTime.timestr().c_str() << " ";
896 out.setFieldWidth(3);
897 out << "Full Clock " << corr->_prn.mid(0,3) << " " << corr->_iod << " ";
898 out.setFieldWidth(14);
899 out << (xc(4) + corr->_dClkResult) * t_CST::c << "\n";
900 out.setFieldWidth(0);
901 }
902 else {
903 out << "bncComb::printResults bug" << "\n";
904 }
905 out.flush();
906 }
907}
908
909// Send results to RTNet Decoder and directly to PPP Client
910////////////////////////////////////////////////////////////////////////////
911void bncComb::dumpResults(const QMap<QString, cmbCorr*>& resCorr) {
912
913 QList<t_orbCorr> orbCorrections;
914 QList<t_clkCorr> clkCorrections;
915 QList<t_satCodeBias> satCodeBiasList;
916
917 unsigned year, month, day, hour, minute;
918 double sec;
919 _resTime.civil_date(year, month, day);
920 _resTime.civil_time(hour, minute, sec);
921
922 QString outLines = QString().asprintf("* %4d %2d %2d %d %d %12.8f\n",
923 year, month, day, hour, minute, sec);
924
925 QMapIterator<QString, cmbCorr*> it(resCorr);
926 while (it.hasNext()) {
927 it.next();
928 cmbCorr* corr = it.value();
929 // ORBIT
930 t_orbCorr orbCorr(corr->_orbCorr);
931 orbCorr._staID = "INTERNAL";
932 orbCorrections.push_back(orbCorr);
933 // CLOCK
934 t_clkCorr clkCorr(corr->_clkCorr);
935 clkCorr._staID = "INTERNAL";
936 clkCorr._dClk = corr->_dClkResult;
937 clkCorr._dotDClk = 0.0;
938 clkCorr._dotDotDClk = 0.0;
939 clkCorrections.push_back(clkCorr);
940 // CODE BIASES
941 t_satCodeBias satCodeBias(corr->_satCodeBias);
942 satCodeBias._staID = "INTERNAL";
943 satCodeBiasList.push_back(satCodeBias);
944
945 ColumnVector xc(6);
946 ColumnVector vv(3);
947 corr->_eph->setClkCorr(dynamic_cast<const t_clkCorr*>(&clkCorr));
948 corr->_eph->setOrbCorr(dynamic_cast<const t_orbCorr*>(&orbCorr));
949 if (corr->_eph->getCrd(_resTime, xc, vv, true) != success) {
950 delete corr;
951 continue;
952 }
953
954 // Correction Phase Center --> CoM
955 // -------------------------------
956 ColumnVector dx(3); dx = 0.0;
957 if (_antex) {
958 double Mjd = _resTime.mjd() + _resTime.daysec()/86400.0;
959 if (_antex->satCoMcorrection(corr->_prn, Mjd, xc.Rows(1,3), dx) != success) {
960 dx = 0;
961 _log += "antenna not found " + corr->_prn.mid(0,3).toLatin1() + '\n';
962 }
963 }
964
965 outLines += corr->_prn.mid(0,3);
966 QString hlp = QString().asprintf(" APC 3 %15.4f %15.4f %15.4f"
967 " Clk 1 %15.4f"
968 " Vel 3 %15.4f %15.4f %15.4f"
969 " CoM 3 %15.4f %15.4f %15.4f",
970 xc(1), xc(2), xc(3),
971 xc(4) * t_CST::c,
972 vv(1), vv(2), vv(3),
973 xc(1)-dx(1), xc(2)-dx(2), xc(3)-dx(3));
974 outLines += hlp;
975 hlp.clear();
976
977 if (satCodeBias._bias.size()) {
978 hlp = QString().asprintf(" CodeBias %2lu", satCodeBias._bias.size());
979 outLines += hlp;
980 hlp.clear();
981 for (unsigned ii = 0; ii < satCodeBias._bias.size(); ii++) {
982 const t_frqCodeBias& frqCodeBias = satCodeBias._bias[ii];
983 if (!frqCodeBias._rnxType2ch.empty()) {
984 hlp = QString().asprintf(" %s%10.6f", frqCodeBias._rnxType2ch.c_str(), frqCodeBias._value);
985 outLines += hlp;
986 hlp.clear();
987 }
988 }
989 }
990 outLines += "\n";
991 delete corr;
992 }
993
994 outLines += "EOE\n"; // End Of Epoch flag
995 //cout << outLines.toStdString();
996
997 if (!_rtnetDecoder) {
998 _rtnetDecoder = new bncRtnetDecoder();
999 }
1000
1001 vector<string> errmsg;
1002 _rtnetDecoder->Decode(outLines.toLatin1().data(), outLines.length(), errmsg);
1003
1004 // Send new Corrections to PPP etc.
1005 // --------------------------------
1006 if (orbCorrections.size() > 0 && clkCorrections.size() > 0) {
1007 emit newOrbCorrections(orbCorrections);
1008 emit newClkCorrections(clkCorrections);
1009 }
1010 if (satCodeBiasList.size()) {
1011 emit newCodeBiases(satCodeBiasList);
1012 }
1013}
1014
1015// Create First Design Matrix and Vector of Measurements
1016////////////////////////////////////////////////////////////////////////////
1017t_irc bncComb::createAmat(char sys, Matrix& AA, ColumnVector& ll, DiagonalMatrix& PP,
1018 const ColumnVector& x0, QMap<QString, cmbCorr*>& resCorr) {
1019
1020 unsigned nPar = _params[sys].size();
1021 unsigned nObs = corrs(sys).size();
1022
1023 if (nObs == 0) {
1024 return failure;
1025 }
1026
1027 int maxSat = _cmbSysPrn[sys];
1028
1029 const int nCon = (_method == filter) ? 1 + maxSat : 0;
1030
1031 AA.ReSize(nObs+nCon, nPar); AA = 0.0;
1032 ll.ReSize(nObs+nCon); ll = 0.0;
1033 PP.ReSize(nObs+nCon); PP = 1.0 / (sigObs * sigObs);
1034
1035 int iObs = 0;
1036 QVectorIterator<cmbCorr*> itCorr(corrs(sys));
1037 while (itCorr.hasNext()) {
1038 cmbCorr* corr = itCorr.next();
1039 QString prn = corr->_prn;
1040
1041 ++iObs;
1042
1043 if (corr->_acName == _masterOrbitAC[sys] && resCorr.find(prn) == resCorr.end()) {
1044 resCorr[prn] = new cmbCorr(*corr);
1045 }
1046
1047 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
1048 cmbParam* pp = _params[sys][iPar-1];
1049 AA(iObs, iPar) = pp->partial(sys, corr->_acName, prn);
1050 }
1051
1052 ll(iObs) = (corr->_clkCorr._dClk * t_CST::c - corr->_codeBiasIF) - DotProduct(AA.Row(iObs), x0);
1053 }
1054
1055 // Regularization
1056 // --------------
1057 if (_method == filter) {
1058 const double Ph = 1.e6;
1059 PP(nObs+1) = Ph;
1060 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
1061 cmbParam* pp = _params[sys][iPar-1];
1062 if ( AA.Column(iPar).maximum_absolute_value() > 0.0 &&
1063 pp->type == cmbParam::clkSat ) {
1064 AA(nObs+1, iPar) = 1.0;
1065 }
1066 }
1067 unsigned flag = 0;
1068 if (sys == 'E') {
1069 flag = 1;
1070 }
1071 if (sys == 'R') {
1072 return success;
1073 }
1074 int iCond = 1;
1075 // GNSS
1076 for (unsigned iGnss = 1; iGnss <= _cmbSysPrn[sys]; iGnss++) {
1077 QString prn = QString("%1%2_%3").arg(sys).arg(iGnss, 2, 10, QChar('0')).arg(flag);
1078 ++iCond;
1079 PP(nObs+iCond) = Ph;
1080 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
1081 cmbParam* pp = _params[sys][iPar-1];
1082 if ( pp &&
1083 AA.Column(iPar).maximum_absolute_value() > 0.0 &&
1084 pp->type == cmbParam::offACSat &&
1085 pp->prn == prn) {
1086 AA(nObs+iCond, iPar) = 1.0;
1087 }
1088 }
1089 }
1090 }
1091
1092 return success;
1093}
1094
1095// Process Epoch - Single-Epoch Method
1096////////////////////////////////////////////////////////////////////////////
1097t_irc bncComb::processEpoch_singleEpoch(char sys, QTextStream& out,
1098 QMap<QString, cmbCorr*>& resCorr,
1099 ColumnVector& dx) {
1100
1101 // Check Satellite Positions for Outliers
1102 // --------------------------------------
1103 if (checkOrbits(sys, out) != success) {
1104 return failure;
1105 }
1106
1107 // Outlier Detection Loop
1108 // ----------------------
1109 while (true) {
1110
1111 // Remove Satellites that are not in Master
1112 // ----------------------------------------
1113 QMutableVectorIterator<cmbCorr*> it(corrs(sys));
1114 while (it.hasNext()) {
1115 cmbCorr* corr = it.next();
1116 QString prn = corr->_prn;
1117 bool foundMaster = false;
1118 QVectorIterator<cmbCorr*> itHlp(corrs(sys));
1119 while (itHlp.hasNext()) {
1120 cmbCorr* corrHlp = itHlp.next();
1121 QString prnHlp = corrHlp->_prn;
1122 QString ACHlp = corrHlp->_acName;
1123 if (ACHlp == _masterOrbitAC[sys] && prn == prnHlp) {
1124 foundMaster = true;
1125 break;
1126 }
1127 }
1128 if (!foundMaster) {
1129 delete corr;
1130 it.remove();
1131 }
1132 }
1133
1134 // Count Number of Observations per Satellite and per AC
1135 // -----------------------------------------------------
1136 QMap<QString, int> numObsPrn;
1137 QMap<QString, int> numObsAC;
1138 QVectorIterator<cmbCorr*> itCorr(corrs(sys));
1139 while (itCorr.hasNext()) {
1140 cmbCorr* corr = itCorr.next();
1141 QString prn = corr->_prn;
1142 QString AC = corr->_acName;
1143 if (numObsPrn.find(prn) == numObsPrn.end()) {
1144 numObsPrn[prn] = 1;
1145 }
1146 else {
1147 numObsPrn[prn] += 1;
1148 }
1149 if (numObsAC.find(AC) == numObsAC.end()) {
1150 numObsAC[AC] = 1;
1151 }
1152 else {
1153 numObsAC[AC] += 1;
1154 }
1155 }
1156
1157 // Clean-Up the Parameters
1158 // -----------------------
1159 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
1160 delete _params[sys][iPar-1];
1161 }
1162 _params[sys].clear();
1163
1164 // Set new Parameters
1165 // ------------------
1166 int nextPar = 0;
1167
1168 QMapIterator<QString, int> itAC(numObsAC);
1169 while (itAC.hasNext()) {
1170 itAC.next();
1171 const QString& AC = itAC.key();
1172 int numObs = itAC.value();
1173 if (AC != _masterOrbitAC[sys] && numObs > 0) {
1174 _params[sys].push_back(new cmbParam(cmbParam::offACgnss, ++nextPar, AC, ""));
1175 }
1176 }
1177
1178 QMapIterator<QString, int> itPrn(numObsPrn);
1179 while (itPrn.hasNext()) {
1180 itPrn.next();
1181 const QString& prn = itPrn.key();
1182 int numObs = itPrn.value();
1183 if (numObs > 0) {
1184 _params[sys].push_back(new cmbParam(cmbParam::clkSat, ++nextPar, "", prn));
1185 }
1186 }
1187
1188 int nPar = _params[sys].size();
1189 ColumnVector x0(nPar);
1190 x0 = 0.0;
1191
1192 // Create First-Design Matrix
1193 // --------------------------
1194 Matrix AA;
1195 ColumnVector ll;
1196 DiagonalMatrix PP;
1197 if (createAmat(sys, AA, ll, PP, x0, resCorr) != success) {
1198 return failure;
1199 }
1200
1201 ColumnVector vv;
1202 try {
1203 Matrix ATP = AA.t() * PP;
1204 SymmetricMatrix NN; NN << ATP * AA;
1205 ColumnVector bb = ATP * ll;
1206 _QQ[sys] = NN.i();
1207 dx = _QQ[sys] * bb;
1208 vv = ll - AA * dx;
1209 }
1210 catch (Exception& exc) {
1211 out << exc.what() << "\n";
1212 return failure;
1213 }
1214
1215 int maxResIndex;
1216 double maxRes = vv.maximum_absolute_value1(maxResIndex);
1217 out.setRealNumberNotation(QTextStream::FixedNotation);
1218 out.setRealNumberPrecision(3);
1219 out << _resTime.datestr().c_str() << " " << _resTime.timestr().c_str()
1220 << " Maximum Residuum " << maxRes << ' '
1221 << corrs(sys)[maxResIndex-1]->_acName << ' ' << corrs(sys)[maxResIndex-1]->_prn.mid(0,3);
1222
1223 if (maxRes > _MAXRES) {
1224 out << " Outlier" << "\n";
1225 delete corrs(sys)[maxResIndex-1];
1226 corrs(sys).remove(maxResIndex-1);
1227 }
1228 else {
1229 out << " OK" << "\n";
1230 out.setRealNumberNotation(QTextStream::FixedNotation);
1231 out.setRealNumberPrecision(3);
1232 for (int ii = 0; ii < vv.Nrows(); ii++) {
1233 const cmbCorr* corr = corrs(sys)[ii];
1234 out << _resTime.datestr().c_str() << ' '
1235 << _resTime.timestr().c_str() << " "
1236 << corr->_acName << ' ' << corr->_prn.mid(0,3);
1237 out.setFieldWidth(6);
1238 out << " res = " << vv[ii] << "\n";
1239 out.setFieldWidth(0);
1240 }
1241 return success;
1242 }
1243
1244 }
1245
1246 return failure;
1247}
1248
1249// Check Satellite Positions for Outliers
1250////////////////////////////////////////////////////////////////////////////
1251t_irc bncComb::checkOrbits(char sys, QTextStream& out) {
1252
1253 const double MAX_DISPLACEMENT = 0.20;
1254
1255 // Switch to last ephemeris (if possible)
1256 // --------------------------------------
1257 QMutableVectorIterator<cmbCorr*> im(corrs(sys));
1258 while (im.hasNext()) {
1259 cmbCorr* corr = im.next();
1260 QString prn = corr->_prn;
1261
1262 t_eph* ephLast = _ephUser.ephLast(prn);
1263 t_eph* ephPrev = _ephUser.ephPrev(prn);
1264
1265 if (ephLast == 0) {
1266 out << "checkOrbit: missing eph (not found) " << corr->_prn.mid(0,3) << "\n";
1267 delete corr;
1268 im.remove();
1269 }
1270 else if (corr->_eph == 0) {
1271 out << "checkOrbit: missing eph (zero) " << corr->_prn.mid(0,3) << "\n";
1272 delete corr;
1273 im.remove();
1274 }
1275 else {
1276 if ( corr->_eph == ephLast || corr->_eph == ephPrev ) {
1277 switchToLastEph(ephLast, corr);
1278 }
1279 else {
1280 out << "checkOrbit: missing eph (deleted) " << corr->_prn.mid(0,3) << "\n";
1281 delete corr;
1282 im.remove();
1283 }
1284 }
1285 }
1286
1287 while (true) {
1288
1289 // Compute Mean Corrections for all Satellites
1290 // -------------------------------------------
1291 QMap<QString, int> numCorr;
1292 QMap<QString, ColumnVector> meanRao;
1293 QVectorIterator<cmbCorr*> it(corrs(sys));
1294 while (it.hasNext()) {
1295 cmbCorr* corr = it.next();
1296 QString prn = corr->_prn;
1297 if (meanRao.find(prn) == meanRao.end()) {
1298 meanRao[prn].ReSize(4);
1299 meanRao[prn].Rows(1,3) = corr->_orbCorr._xr;
1300 meanRao[prn](4) = 1;
1301 }
1302 else {
1303 meanRao[prn].Rows(1,3) += corr->_orbCorr._xr;
1304 meanRao[prn](4) += 1;
1305 }
1306 if (numCorr.find(prn) == numCorr.end()) {
1307 numCorr[prn] = 1;
1308 }
1309 else {
1310 numCorr[prn] += 1;
1311 }
1312 }
1313
1314 // Compute Differences wrt. Mean, find Maximum
1315 // -------------------------------------------
1316 QMap<QString, cmbCorr*> maxDiff;
1317 it.toFront();
1318 while (it.hasNext()) {
1319 cmbCorr* corr = it.next();
1320 QString prn = corr->_prn;
1321 if (meanRao[prn](4) != 0) {
1322 meanRao[prn] /= meanRao[prn](4);
1323 meanRao[prn](4) = 0;
1324 }
1325 corr->_diffRao = corr->_orbCorr._xr - meanRao[prn].Rows(1,3);
1326 if (maxDiff.find(prn) == maxDiff.end()) {
1327 maxDiff[prn] = corr;
1328 }
1329 else {
1330 double normMax = maxDiff[prn]->_diffRao.NormFrobenius();
1331 double norm = corr->_diffRao.NormFrobenius();
1332 if (norm > normMax) {
1333 maxDiff[prn] = corr;
1334 }
1335 }
1336 }
1337
1338 if (_ACs.size() == 1) {
1339 break;
1340 }
1341
1342 // Remove Outliers
1343 // ---------------
1344 bool removed = false;
1345 QMutableVectorIterator<cmbCorr*> im(corrs(sys));
1346 while (im.hasNext()) {
1347 cmbCorr* corr = im.next();
1348 QString prn = corr->_prn;
1349 if (numCorr[prn] < 2) {
1350 delete corr;
1351 im.remove();
1352 }
1353 else if (corr == maxDiff[prn]) {
1354 double norm = corr->_diffRao.NormFrobenius();
1355 if (norm > MAX_DISPLACEMENT) {
1356 out << _resTime.datestr().c_str() << " "
1357 << _resTime.timestr().c_str() << " "
1358 << "Orbit Outlier: "
1359 << corr->_acName.toLatin1().data() << " "
1360 << prn.mid(0,3).toLatin1().data() << " "
1361 << corr->_iod << " "
1362 << norm << "\n";
1363 delete corr;
1364 im.remove();
1365 removed = true;
1366 }
1367 }
1368 }
1369
1370 if (!removed) {
1371 break;
1372 }
1373 }
1374
1375 return success;
1376}
1377
1378//
1379////////////////////////////////////////////////////////////////////////////
1380void bncComb::slotProviderIDChanged(QString mountPoint) {
1381 QMutexLocker locker(&_mutex);
1382
1383 QTextStream out(&_log, QIODevice::WriteOnly);
1384
1385 // Find the AC Name
1386 // ----------------
1387 QString acName;
1388 QListIterator<cmbAC*> icAC(_ACs);
1389 while (icAC.hasNext()) {
1390 cmbAC* AC = icAC.next();
1391 if (AC->mountPoint == mountPoint) {
1392 acName = AC->name;
1393 out << "Provider ID changed: AC " << AC->name.toLatin1().data() << " "
1394 << _resTime.datestr().c_str() << " "
1395 << _resTime.timestr().c_str() << "\n";
1396 break;
1397 }
1398 }
1399 if (acName.isEmpty()) {
1400 return;
1401 }
1402 QMapIterator<char, unsigned> itSys(_cmbSysPrn);
1403 while (itSys.hasNext()) {
1404 itSys.next();
1405 char sys = itSys.key();
1406 // Remove all corrections of the corresponding AC
1407 // ----------------------------------------------
1408 QListIterator<bncTime> itTime(_buffer[sys].keys());
1409 while (itTime.hasNext()) {
1410 bncTime epoTime = itTime.next();
1411 QVector<cmbCorr*>& corrVec = _buffer[sys][epoTime].corrs;
1412 QMutableVectorIterator<cmbCorr*> it(corrVec);
1413 while (it.hasNext()) {
1414 cmbCorr* corr = it.next();
1415 if (acName == corr->_acName) {
1416 delete corr;
1417 it.remove();
1418 }
1419 }
1420 }
1421
1422 // Reset Satellite Offsets
1423 // -----------------------
1424 if (_method == filter) {
1425 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
1426 cmbParam* pp = _params[sys][iPar-1];
1427 if (pp->AC == acName && pp->type == cmbParam::offACSat) {
1428 pp->xx = 0.0;
1429 _QQ[sys].Row(iPar) = 0.0;
1430 _QQ[sys].Column(iPar) = 0.0;
1431 _QQ[sys](iPar,iPar) = pp->sig0 * pp->sig0;
1432 }
1433 }
1434 }
1435 }
1436}
Note: See TracBrowser for help on using the repository browser.