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

Last change on this file since 9690 was 9690, checked in by stuerze, 2 years ago

minor changes to read the bsx file every 6 hours
File size: 41.4 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("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("wrong Bias SINEX file", true);
298 delete _bsx;
299 _bsx = 0;
300 }
301 }
302 if (_bsx) {
303 int 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 if (!bsxFileName.isEmpty()) {
349 _bsx = new bncBiasSnx();
350 }
351
352 if (_bsx->readFile(bsxFileName) != success) {
353 emit newMessage("wrong Bias SINEX file", true);
354 delete _bsx;
355 _bsx = 0;
356 }
357//#ifdef BNC_DEBUG_CMB
358 else {
359 emit newMessage("Successfully read Bias SINEX file", true);
360 }
361//#endif
362}
363
364
365// Remember orbit corrections
366////////////////////////////////////////////////////////////////////////////
367void bncComb::slotNewOrbCorrections(QList<t_orbCorr> orbCorrections) {
368 QMutexLocker locker(&_mutex);
369
370 for (int ii = 0; ii < orbCorrections.size(); ii++) {
371 t_orbCorr& orbCorr = orbCorrections[ii];
372 QString staID(orbCorr._staID.c_str());
373 char sys = orbCorr._prn.system();
374
375 if (!_cmbSysPrn.contains(sys)){
376 continue;
377 }
378
379 // Find/Check the AC Name
380 // ----------------------
381 QString acName;
382 QListIterator<cmbAC*> icAC(_ACs);
383 while (icAC.hasNext()) {
384 cmbAC* AC = icAC.next();
385 if (AC->mountPoint == staID) {
386 acName = AC->name;
387 break;
388 }
389 }
390 if (acName.isEmpty()) {
391 continue;
392 }
393
394 // Store the correction
395 // --------------------
396 QMap<t_prn, t_orbCorr>& storage = _orbCorrections[acName];
397 storage[orbCorr._prn] = orbCorr;
398 }
399}
400
401// Remember Satellite Code Biases
402////////////////////////////////////////////////////////////////////////////
403void bncComb::slotNewCodeBiases(QList<t_satCodeBias> satCodeBiases) {
404 QMutexLocker locker(&_mutex);
405
406 for (int ii = 0; ii < satCodeBiases.size(); ii++) {
407 t_satCodeBias& satCodeBias = satCodeBiases[ii];
408 QString staID(satCodeBias._staID.c_str());
409 char sys = satCodeBias._prn.system();
410
411 if (!_cmbSysPrn.contains(sys)){
412 continue;
413 }
414
415 // Find/Check the AC Name
416 // ----------------------
417 QString acName;
418 QListIterator<cmbAC*> icAC(_ACs);
419 while (icAC.hasNext()) {
420 cmbAC* AC = icAC.next();
421 if (AC->mountPoint == staID) {
422 acName = AC->name;
423 break;
424 }
425 }
426 if (acName.isEmpty()) {
427 continue;
428 }
429
430 // Store the correction
431 // --------------------
432 QMap<t_prn, t_satCodeBias>& storage = _satCodeBiases[acName];
433 storage[satCodeBias._prn] = satCodeBias;
434 }
435}
436
437// Process clock corrections
438////////////////////////////////////////////////////////////////////////////
439void bncComb::slotNewClkCorrections(QList<t_clkCorr> clkCorrections) {
440 QMutexLocker locker(&_mutex);
441
442 bncTime lastTime;
443
444 for (int ii = 0; ii < clkCorrections.size(); ii++) {
445 t_clkCorr& clkCorr = clkCorrections[ii];
446 QString staID(clkCorr._staID.c_str());
447 QString prn(clkCorr._prn.toInternalString().c_str());
448 char sys = clkCorr._prn.system();
449
450 if (!_cmbSysPrn.contains(sys)){
451 continue;
452 }
453
454 // Set the last time
455 // -----------------
456 if (lastTime.undef() || clkCorr._time > lastTime) {
457 lastTime = clkCorr._time;
458 }
459
460 // Find/Check the AC Name
461 // ----------------------
462 QString acName;
463 QListIterator<cmbAC*> icAC(_ACs);
464 while (icAC.hasNext()) {
465 cmbAC* AC = icAC.next();
466 if (AC->mountPoint == staID) {
467 acName = AC->name;
468 break;
469 }
470 }
471 if (acName.isEmpty()) {
472 continue;
473 }
474
475 // Check Modulo Time
476 // -----------------
477 int sec = int(nint(clkCorr._time.gpssec()*10));
478 if (sec % (_cmbSampl*10) != 0.0) {
479 continue;
480 }
481
482 // Check Correction Age
483 // --------------------
484 if (_resTime.valid() && clkCorr._time <= _resTime) {
485#ifdef BNC_DEBUG_CMB
486 emit newMessage("bncComb: old correction: " + acName.toLatin1() + " " + prn.mid(0,3).toLatin1(), true);
487#endif
488 continue;
489 }
490
491 // Create new correction
492 // ---------------------
493 cmbCorr* newCorr = new cmbCorr();
494 newCorr->_prn = prn;
495 newCorr->_time = clkCorr._time;
496 newCorr->_iod = clkCorr._iod;
497 newCorr->_acName = acName;
498 newCorr->_clkCorr = clkCorr;
499
500 // Check orbit correction
501 // ----------------------
502 if (!_orbCorrections.contains(acName)) {
503 delete newCorr;
504 continue;
505 }
506 else {
507 QMap<t_prn, t_orbCorr>& storage = _orbCorrections[acName];
508 if (!storage.contains(clkCorr._prn) ||
509 storage[clkCorr._prn]._iod != newCorr->_iod) {
510 delete newCorr;
511 continue;
512 }
513 else {
514 newCorr->_orbCorr = storage[clkCorr._prn];
515 }
516 }
517
518 // Check the Ephemeris
519 //--------------------
520 t_eph* ephLast = _ephUser.ephLast(prn);
521 t_eph* ephPrev = _ephUser.ephPrev(prn);
522 if (ephLast == 0) {
523#ifdef BNC_DEBUG_CMB
524 emit newMessage("bncComb: eph not found for " + prn.mid(0,3).toLatin1(), true);
525#endif
526 delete newCorr;
527 continue;
528 }
529 else if (ephLast->checkState() == t_eph::bad ||
530 ephLast->checkState() == t_eph::outdated ||
531 ephLast->checkState() == t_eph::unhealthy) {
532#ifdef BNC_DEBUG_CMB
533 emit newMessage("bncComb: ephLast not ok (checkState: " + ephLast->checkStateToString().toLatin1() + ") for " + prn.mid(0,3).toLatin1(), true);
534#endif
535 delete newCorr;
536 continue;
537 }
538 else {
539 if (ephLast->IOD() == newCorr->_iod) {
540 newCorr->_eph = ephLast;
541 }
542 else if (ephPrev && ephPrev->checkState() == t_eph::ok &&
543 ephPrev->IOD() == newCorr->_iod) {
544 newCorr->_eph = ephPrev;
545 switchToLastEph(ephLast, newCorr);
546 }
547 else {
548#ifdef BNC_DEBUG_CMB
549 emit newMessage("bncComb: eph not found for " + prn.mid(0,3).toLatin1() +
550 QString(" with IOD %1").arg(newCorr->_iod).toLatin1(), true);
551#endif
552 delete newCorr;
553 continue;
554 }
555 }
556
557 // Check satellite code biases
558 // ----------------------------
559 if (!_satCodeBiases.contains(acName)) {
560 delete newCorr;
561 continue;
562 }
563 else {
564 QMap<t_prn, t_satCodeBias>& storage = _satCodeBiases[acName];
565 if (!storage.contains(clkCorr._prn)) {
566 delete newCorr;
567 continue;
568 }
569 else {
570 newCorr->_satCodeBias = storage[clkCorr._prn];
571 QMap<t_frequency::type, double> codeBiasesRefSig;
572 for (unsigned ii = 1; ii < cmbRefSig::cIF; ii++) {
573 t_frequency::type frqType = cmbRefSig::toFreq(sys, static_cast<cmbRefSig::type>(ii));
574 char frqNum = t_frequency::toString(frqType)[1];
575 char attrib = cmbRefSig::toAttrib(sys, static_cast<cmbRefSig::type>(ii));
576 QString rnxType2ch = QString("%1%2").arg(frqNum).arg(attrib);
577 for (unsigned ii = 0; ii < newCorr->_satCodeBias._bias.size(); ii++) {
578 const t_frqCodeBias& bias = newCorr->_satCodeBias._bias[ii];
579 if (rnxType2ch.toStdString() == bias._rnxType2ch) {
580 codeBiasesRefSig[frqType] = bias._value;
581 }
582 }
583 }
584 if (codeBiasesRefSig.size() < 2) {
585 delete newCorr;
586 continue;
587 }
588 map<t_frequency::type, double> codeCoeff;
589 double channel = double(newCorr->_eph->slotNum());
590 cmbRefSig::coeff(sys, cmbRefSig::cIF, channel, codeCoeff);
591 map<t_frequency::type, double>::const_iterator it;
592 for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
593 t_frequency::type frqType = it->first;
594 newCorr->_codeBiasIF += it->second * codeBiasesRefSig[frqType];
595 }
596 newCorr->_satCodeBias._bias.clear();
597 }
598 //cout << acName.toStdString() << " " << clkCorr._prn.toString().c_str() << " _codeBiasIF: " << newCorr->_codeBiasIF << endl;
599 }
600
601 // Store correction into the buffer
602 // --------------------------------
603 QVector<cmbCorr*>& corrs = _buffer[sys][newCorr->_time].corrs;
604 corrs.push_back(newCorr);
605 }
606
607 // Process previous Epoch(s)
608 // -------------------------
609 const double outWait = 1.0 * _cmbSampl;
610 QMapIterator<char, unsigned> itSys(_cmbSysPrn);
611 while (itSys.hasNext()) {
612 itSys.next();
613 char sys = itSys.key();
614 QListIterator<bncTime> itTime(_buffer[sys].keys());
615 while (itTime.hasNext()) {
616 bncTime epoTime = itTime.next();
617 if (epoTime < lastTime - outWait) {
618 _resTime = epoTime;
619 processEpoch(sys);
620 }
621 }
622 }
623}
624
625// Change the correction so that it refers to last received ephemeris
626////////////////////////////////////////////////////////////////////////////
627void bncComb::switchToLastEph(t_eph* lastEph, cmbCorr* corr) {
628
629 if (corr->_eph == lastEph) {
630 return;
631 }
632
633 ColumnVector oldXC(6);
634 ColumnVector oldVV(3);
635 if (corr->_eph->getCrd(corr->_time, oldXC, oldVV, false) != success) {
636 return;
637 }
638
639 ColumnVector newXC(6);
640 ColumnVector newVV(3);
641 if (lastEph->getCrd(corr->_time, newXC, newVV, false) != success) {
642 return;
643 }
644
645 ColumnVector dX = newXC.Rows(1,3) - oldXC.Rows(1,3);
646 ColumnVector dV = newVV - oldVV;
647 double dC = newXC(4) - oldXC(4);
648
649 ColumnVector dRAO(3);
650 XYZ_to_RSW(newXC.Rows(1,3), newVV, dX, dRAO);
651
652 ColumnVector dDotRAO(3);
653 XYZ_to_RSW(newXC.Rows(1,3), newVV, dV, dDotRAO);
654
655 QString msg = "switch corr " + corr->_prn.mid(0,3)
656 + QString(" %1 -> %2 %3").arg(corr->_iod,3).arg(lastEph->IOD(),3).arg(dC*t_CST::c, 8, 'f', 4);
657
658 emit newMessage(msg.toLatin1(), false);
659
660 corr->_iod = lastEph->IOD();
661 corr->_eph = lastEph;
662
663 corr->_orbCorr._xr += dRAO;
664 corr->_orbCorr._dotXr += dDotRAO;
665 corr->_clkCorr._dClk -= dC;
666}
667
668// Process Epoch
669////////////////////////////////////////////////////////////////////////////
670void bncComb::processEpoch(char sys) {
671
672 _log.clear();
673
674 QTextStream out(&_log, QIODevice::WriteOnly);
675
676 out << "\n" << "Combination: " << sys << "\n"
677 << "--------------------------------" << "\n";
678
679 // Observation Statistics
680 // ----------------------
681 bool masterPresent = false;
682 QListIterator<cmbAC*> icAC(_ACs);
683 while (icAC.hasNext()) {
684 cmbAC* AC = icAC.next();
685 AC->numObs[sys] = 0;
686 QVectorIterator<cmbCorr*> itCorr(corrs(sys));
687 while (itCorr.hasNext()) {
688 cmbCorr* corr = itCorr.next();
689 if (corr->_acName == AC->name) {
690 AC->numObs[sys] += 1;
691 if (AC->name == _masterOrbitAC[sys]) {
692 masterPresent = true;
693 }
694 }
695 }
696 out << AC->name.toLatin1().data() << ": " << AC->numObs[sys] << "\n";
697 }
698
699 // If Master not present, switch to another one
700 // --------------------------------------------
701 const unsigned switchMasterAfterGap = 1;
702 if (masterPresent) {
703 _masterMissingEpochs[sys] = 0;
704 }
705 else {
706 ++_masterMissingEpochs[sys];
707 if (_masterMissingEpochs[sys] < switchMasterAfterGap) {
708 out << "Missing Master, Epoch skipped" << "\n";
709 _buffer[sys].remove(_resTime);
710 emit newMessage(_log, false);
711 return;
712 }
713 else {
714 _masterMissingEpochs[sys] = 0;
715 QListIterator<cmbAC*> icAC(_ACs);
716 while (icAC.hasNext()) {
717 cmbAC* AC = icAC.next();
718 if (AC->numObs[sys] > 0) {
719 out << "Switching Master AC "
720 << _masterOrbitAC[sys].toLatin1().data() << " --> "
721 << AC->name.toLatin1().data() << " "
722 << _resTime.datestr().c_str() << " "
723 << _resTime.timestr().c_str() << "\n";
724 _masterOrbitAC[sys] = AC->name;
725 break;
726 }
727 }
728 }
729 }
730
731 QMap<QString, cmbCorr*> resCorr;
732
733 // Perform the actual Combination using selected Method
734 // ----------------------------------------------------
735 t_irc irc;
736 ColumnVector dx;
737 if (_method == filter) {
738 irc = processEpoch_filter(sys, out, resCorr, dx);
739 }
740 else {
741 irc = processEpoch_singleEpoch(sys, out, resCorr, dx);
742 }
743
744 // Update Parameter Values, Print Results
745 // --------------------------------------
746 if (irc == success) {
747 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
748 cmbParam* pp = _params[sys][iPar-1];
749 pp->xx += dx(iPar);
750 if (pp->type == cmbParam::clkSat) {
751 if (resCorr.find(pp->prn) != resCorr.end()) {
752 // set clock result
753 resCorr[pp->prn]->_dClkResult = pp->xx / t_CST::c;
754 // Add Code Biases from SINEX File
755 if (_bsx) {
756 map<t_frequency::type, double> codeCoeff;
757 double channel = double(resCorr[pp->prn]->_eph->slotNum());
758 cmbRefSig::coeff(sys, cmbRefSig::cIF, channel, codeCoeff);
759 t_frequency::type fType1 = cmbRefSig::toFreq(sys, cmbRefSig::c1);
760 t_frequency::type fType2 = cmbRefSig::toFreq(sys, cmbRefSig::c2);
761 _bsx->determineSsrSatCodeBiases(pp->prn.mid(0,3), codeCoeff[fType1], codeCoeff[fType2], resCorr[pp->prn]->_satCodeBias);
762 }
763 }
764 }
765 out << _resTime.datestr().c_str() << " "
766 << _resTime.timestr().c_str() << " ";
767 out.setRealNumberNotation(QTextStream::FixedNotation);
768 out.setFieldWidth(8);
769 out.setRealNumberPrecision(4);
770 out << pp->toString(sys) << " "
771 << pp->xx << " +- " << sqrt(_QQ[sys](pp->index,pp->index)) << "\n";
772 out.setFieldWidth(0);
773 }
774 printResults(out, resCorr);
775 dumpResults(resCorr);
776 }
777
778 // Delete Data, emit Message
779 // -------------------------
780 _buffer[sys].remove(_resTime);
781 emit newMessage(_log, false);
782}
783
784// Process Epoch - Filter Method
785////////////////////////////////////////////////////////////////////////////
786t_irc bncComb::processEpoch_filter(char sys, QTextStream& out,
787 QMap<QString, cmbCorr*>& resCorr,
788 ColumnVector& dx) {
789
790 // Prediction Step
791 // ---------------
792 int nPar = _params[sys].size();
793 ColumnVector x0(nPar);
794 for (int iPar = 1; iPar <= nPar; iPar++) {
795 cmbParam* pp = _params[sys][iPar-1];
796 if (pp->epoSpec) {
797 pp->xx = 0.0;
798 _QQ[sys].Row(iPar) = 0.0;
799 _QQ[sys].Column(iPar) = 0.0;
800 _QQ[sys](iPar,iPar) = pp->sig0 * pp->sig0;
801 }
802 else {
803 _QQ[sys](iPar,iPar) += pp->sigP * pp->sigP;
804 }
805 x0(iPar) = pp->xx;
806 }
807
808 // Check Satellite Positions for Outliers
809 // --------------------------------------
810 if (checkOrbits(sys, out) != success) {
811 return failure;
812 }
813
814 // Update and outlier detection loop
815 // ---------------------------------
816 SymmetricMatrix QQ_sav = _QQ[sys];
817 while (true) {
818
819 Matrix AA;
820 ColumnVector ll;
821 DiagonalMatrix PP;
822
823 if (createAmat(sys, AA, ll, PP, x0, resCorr) != success) {
824 return failure;
825 }
826
827 dx.ReSize(nPar); dx = 0.0;
828 kalman(AA, ll, PP, _QQ[sys], dx);
829
830 ColumnVector vv = ll - AA * dx;
831
832 int maxResIndex;
833 double maxRes = vv.maximum_absolute_value1(maxResIndex);
834 out.setRealNumberNotation(QTextStream::FixedNotation);
835 out.setRealNumberPrecision(3);
836 out << _resTime.datestr().c_str() << " " << _resTime.timestr().c_str()
837 << " Maximum Residuum " << maxRes << ' '
838 << corrs(sys)[maxResIndex-1]->_acName << ' ' << corrs(sys)[maxResIndex-1]->_prn.mid(0,3);
839 if (maxRes > _MAXRES) {
840 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
841 cmbParam* pp = _params[sys][iPar-1];
842 if (pp->type == cmbParam::offACSat &&
843 pp->AC == corrs(sys)[maxResIndex-1]->_acName &&
844 pp->prn == corrs(sys)[maxResIndex-1]->_prn.mid(0,3)) {
845 QQ_sav.Row(iPar) = 0.0;
846 QQ_sav.Column(iPar) = 0.0;
847 QQ_sav(iPar,iPar) = pp->sig0 * pp->sig0;
848 }
849 }
850
851 out << " Outlier" << "\n";
852 _QQ[sys] = QQ_sav;
853 corrs(sys).remove(maxResIndex-1);
854 }
855 else {
856 out << " OK" << "\n";
857 out.setRealNumberNotation(QTextStream::FixedNotation);
858 out.setRealNumberPrecision(4);
859 for (int ii = 0; ii < corrs(sys).size(); ii++) {
860 const cmbCorr* corr = corrs(sys)[ii];
861 out << _resTime.datestr().c_str() << ' '
862 << _resTime.timestr().c_str() << " "
863 << corr->_acName << ' ' << corr->_prn.mid(0,3);
864 out.setFieldWidth(10);
865 out << " res = " << vv[ii] << "\n";
866 out.setFieldWidth(0);
867 }
868 break;
869 }
870 }
871
872 return success;
873}
874
875// Print results
876////////////////////////////////////////////////////////////////////////////
877void bncComb::printResults(QTextStream& out,
878 const QMap<QString, cmbCorr*>& resCorr) {
879
880 QMapIterator<QString, cmbCorr*> it(resCorr);
881 while (it.hasNext()) {
882 it.next();
883 cmbCorr* corr = it.value();
884 const t_eph* eph = corr->_eph;
885 if (eph) {
886 ColumnVector xc(6);
887 ColumnVector vv(3);
888 if (eph->getCrd(_resTime, xc, vv, false) != success) {
889 continue;
890 }
891 out << _resTime.datestr().c_str() << " "
892 << _resTime.timestr().c_str() << " ";
893 out.setFieldWidth(3);
894 out << "Full Clock " << corr->_prn.mid(0,3) << " " << corr->_iod << " ";
895 out.setFieldWidth(14);
896 out << (xc(4) + corr->_dClkResult) * t_CST::c << "\n";
897 out.setFieldWidth(0);
898 }
899 else {
900 out << "bncComb::printResuls bug" << "\n";
901 }
902 out.flush();
903 }
904}
905
906// Send results to RTNet Decoder and directly to PPP Client
907////////////////////////////////////////////////////////////////////////////
908void bncComb::dumpResults(const QMap<QString, cmbCorr*>& resCorr) {
909
910 QList<t_orbCorr> orbCorrections;
911 QList<t_clkCorr> clkCorrections;
912 QList<t_satCodeBias> satCodeBiasList;
913
914 unsigned year, month, day, hour, minute;
915 double sec;
916 _resTime.civil_date(year, month, day);
917 _resTime.civil_time(hour, minute, sec);
918
919 QString outLines = QString().asprintf("* %4d %2d %2d %d %d %12.8f\n",
920 year, month, day, hour, minute, sec);
921
922 QMapIterator<QString, cmbCorr*> it(resCorr);
923 while (it.hasNext()) {
924 it.next();
925 cmbCorr* corr = it.value();
926 // ORBIT
927 t_orbCorr orbCorr(corr->_orbCorr);
928 orbCorr._staID = "INTERNAL";
929 orbCorrections.push_back(orbCorr);
930 // CLOCK
931 t_clkCorr clkCorr(corr->_clkCorr);
932 clkCorr._staID = "INTERNAL";
933 clkCorr._dClk = corr->_dClkResult;
934 clkCorr._dotDClk = 0.0;
935 clkCorr._dotDotDClk = 0.0;
936 clkCorrections.push_back(clkCorr);
937
938 ColumnVector xc(6);
939 ColumnVector vv(3);
940 corr->_eph->setClkCorr(dynamic_cast<const t_clkCorr*>(&clkCorr));
941 corr->_eph->setOrbCorr(dynamic_cast<const t_orbCorr*>(&orbCorr));
942 if (corr->_eph->getCrd(_resTime, xc, vv, true) != success) {
943 delete corr;
944 continue;
945 }
946
947 // Correction Phase Center --> CoM
948 // -------------------------------
949 ColumnVector dx(3); dx = 0.0;
950 if (_antex) {
951 double Mjd = _resTime.mjd() + _resTime.daysec()/86400.0;
952 if (_antex->satCoMcorrection(corr->_prn, Mjd, xc.Rows(1,3), dx) != success) {
953 dx = 0;
954 _log += "antenna not found " + corr->_prn.mid(0,3).toLatin1() + '\n';
955 }
956 }
957
958 outLines += corr->_prn.mid(0,3);
959 QString hlp = QString().asprintf(" APC 3 %15.4f %15.4f %15.4f"
960 " Clk 1 %15.4f"
961 " Vel 3 %15.4f %15.4f %15.4f"
962 " CoM 3 %15.4f %15.4f %15.4f",
963 xc(1), xc(2), xc(3),
964 xc(4) * t_CST::c,
965 vv(1), vv(2), vv(3),
966 xc(1)-dx(1), xc(2)-dx(2), xc(3)-dx(3));
967 outLines += hlp;
968 hlp.clear();
969
970 // CODE BIASES
971 t_satCodeBias satCodeBias(corr->_satCodeBias);
972 satCodeBias._staID = "INTERNAL";
973 satCodeBiasList.push_back(satCodeBias);
974
975 if (satCodeBias._bias.size()) {
976 hlp = QString().asprintf(" CodeBias %2lu", satCodeBias._bias.size());
977 outLines += hlp;
978 hlp.clear();
979 for (unsigned ii = 0; ii < satCodeBias._bias.size(); ii++) {
980 const t_frqCodeBias& frqCodeBias = satCodeBias._bias[ii];
981 if (!frqCodeBias._rnxType2ch.empty()) {
982 hlp = QString().asprintf(" %s%10.6f", frqCodeBias._rnxType2ch.c_str(), frqCodeBias._value);
983 outLines += hlp;
984 hlp.clear();
985 }
986 }
987 }
988 outLines += "\n";
989 //cout << outLines.toStdString();
990 delete corr;
991 }
992
993 outLines += "EOE\n"; // End Of Epoch flag
994
995 if (!_rtnetDecoder) {
996 _rtnetDecoder = new bncRtnetDecoder();
997 }
998
999 vector<string> errmsg;
1000 _rtnetDecoder->Decode(outLines.toLatin1().data(), outLines.length(), errmsg);
1001
1002 // Send new Corrections to PPP etc.
1003 // --------------------------------
1004 if (orbCorrections.size() > 0 && clkCorrections.size() > 0) {
1005 emit newOrbCorrections(orbCorrections);
1006 emit newClkCorrections(clkCorrections);
1007 }
1008 if (satCodeBiasList.size()) {
1009 emit newCodeBiases(satCodeBiasList);
1010 }
1011}
1012
1013// Create First Design Matrix and Vector of Measurements
1014////////////////////////////////////////////////////////////////////////////
1015t_irc bncComb::createAmat(char sys, Matrix& AA, ColumnVector& ll, DiagonalMatrix& PP,
1016 const ColumnVector& x0, QMap<QString, cmbCorr*>& resCorr) {
1017
1018 unsigned nPar = _params[sys].size();
1019 unsigned nObs = corrs(sys).size();
1020
1021 if (nObs == 0) {
1022 return failure;
1023 }
1024
1025 int maxSat = _cmbSysPrn[sys];
1026
1027 const int nCon = (_method == filter) ? 1 + maxSat : 0;
1028
1029 AA.ReSize(nObs+nCon, nPar); AA = 0.0;
1030 ll.ReSize(nObs+nCon); ll = 0.0;
1031 PP.ReSize(nObs+nCon); PP = 1.0 / (sigObs * sigObs);
1032
1033 int iObs = 0;
1034 QVectorIterator<cmbCorr*> itCorr(corrs(sys));
1035 while (itCorr.hasNext()) {
1036 cmbCorr* corr = itCorr.next();
1037 QString prn = corr->_prn;
1038
1039 ++iObs;
1040
1041 if (corr->_acName == _masterOrbitAC[sys] && resCorr.find(prn) == resCorr.end()) {
1042 resCorr[prn] = new cmbCorr(*corr);
1043 }
1044
1045 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
1046 cmbParam* pp = _params[sys][iPar-1];
1047 AA(iObs, iPar) = pp->partial(sys, corr->_acName, prn);
1048 }
1049
1050 ll(iObs) = (corr->_clkCorr._dClk * t_CST::c - corr->_codeBiasIF) - DotProduct(AA.Row(iObs), x0);
1051 }
1052
1053 // Regularization
1054 // --------------
1055 if (_method == filter) {
1056 const double Ph = 1.e6;
1057 PP(nObs+1) = Ph;
1058 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
1059 cmbParam* pp = _params[sys][iPar-1];
1060 if ( AA.Column(iPar).maximum_absolute_value() > 0.0 &&
1061 pp->type == cmbParam::clkSat ) {
1062 AA(nObs+1, iPar) = 1.0;
1063 }
1064 }
1065 unsigned flag = 0;
1066 if (sys == 'E') {
1067 flag = 1;
1068 }
1069 if (sys == 'R') {
1070 return success;
1071 }
1072 int iCond = 1;
1073 // GNSS
1074 for (unsigned iGnss = 1; iGnss <= _cmbSysPrn[sys]; iGnss++) {
1075 QString prn = QString("%1%2_%3").arg(sys).arg(iGnss, 2, 10, QChar('0')).arg(flag);
1076 ++iCond;
1077 PP(nObs+iCond) = Ph;
1078 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
1079 cmbParam* pp = _params[sys][iPar-1];
1080 if ( pp &&
1081 AA.Column(iPar).maximum_absolute_value() > 0.0 &&
1082 pp->type == cmbParam::offACSat &&
1083 pp->prn == prn) {
1084 AA(nObs+iCond, iPar) = 1.0;
1085 }
1086 }
1087 }
1088 }
1089
1090 return success;
1091}
1092
1093// Process Epoch - Single-Epoch Method
1094////////////////////////////////////////////////////////////////////////////
1095t_irc bncComb::processEpoch_singleEpoch(char sys, QTextStream& out,
1096 QMap<QString, cmbCorr*>& resCorr,
1097 ColumnVector& dx) {
1098
1099 // Check Satellite Positions for Outliers
1100 // --------------------------------------
1101 if (checkOrbits(sys, out) != success) {
1102 return failure;
1103 }
1104
1105 // Outlier Detection Loop
1106 // ----------------------
1107 while (true) {
1108
1109 // Remove Satellites that are not in Master
1110 // ----------------------------------------
1111 QMutableVectorIterator<cmbCorr*> it(corrs(sys));
1112 while (it.hasNext()) {
1113 cmbCorr* corr = it.next();
1114 QString prn = corr->_prn;
1115 bool foundMaster = false;
1116 QVectorIterator<cmbCorr*> itHlp(corrs(sys));
1117 while (itHlp.hasNext()) {
1118 cmbCorr* corrHlp = itHlp.next();
1119 QString prnHlp = corrHlp->_prn;
1120 QString ACHlp = corrHlp->_acName;
1121 if (ACHlp == _masterOrbitAC[sys] && prn == prnHlp) {
1122 foundMaster = true;
1123 break;
1124 }
1125 }
1126 if (!foundMaster) {
1127 delete corr;
1128 it.remove();
1129 }
1130 }
1131
1132 // Count Number of Observations per Satellite and per AC
1133 // -----------------------------------------------------
1134 QMap<QString, int> numObsPrn;
1135 QMap<QString, int> numObsAC;
1136 QVectorIterator<cmbCorr*> itCorr(corrs(sys));
1137 while (itCorr.hasNext()) {
1138 cmbCorr* corr = itCorr.next();
1139 QString prn = corr->_prn;
1140 QString AC = corr->_acName;
1141 if (numObsPrn.find(prn) == numObsPrn.end()) {
1142 numObsPrn[prn] = 1;
1143 }
1144 else {
1145 numObsPrn[prn] += 1;
1146 }
1147 if (numObsAC.find(AC) == numObsAC.end()) {
1148 numObsAC[AC] = 1;
1149 }
1150 else {
1151 numObsAC[AC] += 1;
1152 }
1153 }
1154
1155 // Clean-Up the Parameters
1156 // -----------------------
1157 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
1158 delete _params[sys][iPar-1];
1159 }
1160 _params[sys].clear();
1161
1162 // Set new Parameters
1163 // ------------------
1164 int nextPar = 0;
1165
1166 QMapIterator<QString, int> itAC(numObsAC);
1167 while (itAC.hasNext()) {
1168 itAC.next();
1169 const QString& AC = itAC.key();
1170 int numObs = itAC.value();
1171 if (AC != _masterOrbitAC[sys] && numObs > 0) {
1172 _params[sys].push_back(new cmbParam(cmbParam::offACgnss, ++nextPar, AC, ""));
1173 }
1174 }
1175
1176 QMapIterator<QString, int> itPrn(numObsPrn);
1177 while (itPrn.hasNext()) {
1178 itPrn.next();
1179 const QString& prn = itPrn.key();
1180 int numObs = itPrn.value();
1181 if (numObs > 0) {
1182 _params[sys].push_back(new cmbParam(cmbParam::clkSat, ++nextPar, "", prn));
1183 }
1184 }
1185
1186 int nPar = _params[sys].size();
1187 ColumnVector x0(nPar);
1188 x0 = 0.0;
1189
1190 // Create First-Design Matrix
1191 // --------------------------
1192 Matrix AA;
1193 ColumnVector ll;
1194 DiagonalMatrix PP;
1195 if (createAmat(sys, AA, ll, PP, x0, resCorr) != success) {
1196 return failure;
1197 }
1198
1199 ColumnVector vv;
1200 try {
1201 Matrix ATP = AA.t() * PP;
1202 SymmetricMatrix NN; NN << ATP * AA;
1203 ColumnVector bb = ATP * ll;
1204 _QQ[sys] = NN.i();
1205 dx = _QQ[sys] * bb;
1206 vv = ll - AA * dx;
1207 }
1208 catch (Exception& exc) {
1209 out << exc.what() << "\n";
1210 return failure;
1211 }
1212
1213 int maxResIndex;
1214 double maxRes = vv.maximum_absolute_value1(maxResIndex);
1215 out.setRealNumberNotation(QTextStream::FixedNotation);
1216 out.setRealNumberPrecision(3);
1217 out << _resTime.datestr().c_str() << " " << _resTime.timestr().c_str()
1218 << " Maximum Residuum " << maxRes << ' '
1219 << corrs(sys)[maxResIndex-1]->_acName << ' ' << corrs(sys)[maxResIndex-1]->_prn.mid(0,3);
1220
1221 if (maxRes > _MAXRES) {
1222 out << " Outlier" << "\n";
1223 delete corrs(sys)[maxResIndex-1];
1224 corrs(sys).remove(maxResIndex-1);
1225 }
1226 else {
1227 out << " OK" << "\n";
1228 out.setRealNumberNotation(QTextStream::FixedNotation);
1229 out.setRealNumberPrecision(3);
1230 for (int ii = 0; ii < vv.Nrows(); ii++) {
1231 const cmbCorr* corr = corrs(sys)[ii];
1232 out << _resTime.datestr().c_str() << ' '
1233 << _resTime.timestr().c_str() << " "
1234 << corr->_acName << ' ' << corr->_prn.mid(0,3);
1235 out.setFieldWidth(6);
1236 out << " res = " << vv[ii] << "\n";
1237 out.setFieldWidth(0);
1238 }
1239 return success;
1240 }
1241
1242 }
1243
1244 return failure;
1245}
1246
1247// Check Satellite Positions for Outliers
1248////////////////////////////////////////////////////////////////////////////
1249t_irc bncComb::checkOrbits(char sys, QTextStream& out) {
1250
1251 const double MAX_DISPLACEMENT = 0.20;
1252
1253 // Switch to last ephemeris (if possible)
1254 // --------------------------------------
1255 QMutableVectorIterator<cmbCorr*> im(corrs(sys));
1256 while (im.hasNext()) {
1257 cmbCorr* corr = im.next();
1258 QString prn = corr->_prn;
1259
1260 t_eph* ephLast = _ephUser.ephLast(prn);
1261 t_eph* ephPrev = _ephUser.ephPrev(prn);
1262
1263 if (ephLast == 0) {
1264 out << "checkOrbit: missing eph (not found) " << corr->_prn.mid(0,3) << "\n";
1265 delete corr;
1266 im.remove();
1267 }
1268 else if (corr->_eph == 0) {
1269 out << "checkOrbit: missing eph (zero) " << corr->_prn.mid(0,3) << "\n";
1270 delete corr;
1271 im.remove();
1272 }
1273 else {
1274 if ( corr->_eph == ephLast || corr->_eph == ephPrev ) {
1275 switchToLastEph(ephLast, corr);
1276 }
1277 else {
1278 out << "checkOrbit: missing eph (deleted) " << corr->_prn.mid(0,3) << "\n";
1279 delete corr;
1280 im.remove();
1281 }
1282 }
1283 }
1284
1285 while (true) {
1286
1287 // Compute Mean Corrections for all Satellites
1288 // -------------------------------------------
1289 QMap<QString, int> numCorr;
1290 QMap<QString, ColumnVector> meanRao;
1291 QVectorIterator<cmbCorr*> it(corrs(sys));
1292 while (it.hasNext()) {
1293 cmbCorr* corr = it.next();
1294 QString prn = corr->_prn;
1295 if (meanRao.find(prn) == meanRao.end()) {
1296 meanRao[prn].ReSize(4);
1297 meanRao[prn].Rows(1,3) = corr->_orbCorr._xr;
1298 meanRao[prn](4) = 1;
1299 }
1300 else {
1301 meanRao[prn].Rows(1,3) += corr->_orbCorr._xr;
1302 meanRao[prn](4) += 1;
1303 }
1304 if (numCorr.find(prn) == numCorr.end()) {
1305 numCorr[prn] = 1;
1306 }
1307 else {
1308 numCorr[prn] += 1;
1309 }
1310 }
1311
1312 // Compute Differences wrt. Mean, find Maximum
1313 // -------------------------------------------
1314 QMap<QString, cmbCorr*> maxDiff;
1315 it.toFront();
1316 while (it.hasNext()) {
1317 cmbCorr* corr = it.next();
1318 QString prn = corr->_prn;
1319 if (meanRao[prn](4) != 0) {
1320 meanRao[prn] /= meanRao[prn](4);
1321 meanRao[prn](4) = 0;
1322 }
1323 corr->_diffRao = corr->_orbCorr._xr - meanRao[prn].Rows(1,3);
1324 if (maxDiff.find(prn) == maxDiff.end()) {
1325 maxDiff[prn] = corr;
1326 }
1327 else {
1328 double normMax = maxDiff[prn]->_diffRao.NormFrobenius();
1329 double norm = corr->_diffRao.NormFrobenius();
1330 if (norm > normMax) {
1331 maxDiff[prn] = corr;
1332 }
1333 }
1334 }
1335
1336 if (_ACs.size() == 1) {
1337 break;
1338 }
1339
1340 // Remove Outliers
1341 // ---------------
1342 bool removed = false;
1343 QMutableVectorIterator<cmbCorr*> im(corrs(sys));
1344 while (im.hasNext()) {
1345 cmbCorr* corr = im.next();
1346 QString prn = corr->_prn;
1347 if (numCorr[prn] < 2) {
1348 delete corr;
1349 im.remove();
1350 }
1351 else if (corr == maxDiff[prn]) {
1352 double norm = corr->_diffRao.NormFrobenius();
1353 if (norm > MAX_DISPLACEMENT) {
1354 out << _resTime.datestr().c_str() << " "
1355 << _resTime.timestr().c_str() << " "
1356 << "Orbit Outlier: "
1357 << corr->_acName.toLatin1().data() << " "
1358 << prn.mid(0,3).toLatin1().data() << " "
1359 << corr->_iod << " "
1360 << norm << "\n";
1361 delete corr;
1362 im.remove();
1363 removed = true;
1364 }
1365 }
1366 }
1367
1368 if (!removed) {
1369 break;
1370 }
1371 }
1372
1373 return success;
1374}
1375
1376//
1377////////////////////////////////////////////////////////////////////////////
1378void bncComb::slotProviderIDChanged(QString mountPoint) {
1379 QMutexLocker locker(&_mutex);
1380
1381 QTextStream out(&_log, QIODevice::WriteOnly);
1382
1383 // Find the AC Name
1384 // ----------------
1385 QString acName;
1386 QListIterator<cmbAC*> icAC(_ACs);
1387 while (icAC.hasNext()) {
1388 cmbAC* AC = icAC.next();
1389 if (AC->mountPoint == mountPoint) {
1390 acName = AC->name;
1391 out << "Provider ID changed: AC " << AC->name.toLatin1().data() << " "
1392 << _resTime.datestr().c_str() << " "
1393 << _resTime.timestr().c_str() << "\n";
1394 break;
1395 }
1396 }
1397 if (acName.isEmpty()) {
1398 return;
1399 }
1400 QMapIterator<char, unsigned> itSys(_cmbSysPrn);
1401 while (itSys.hasNext()) {
1402 itSys.next();
1403 char sys = itSys.key();
1404 // Remove all corrections of the corresponding AC
1405 // ----------------------------------------------
1406 QListIterator<bncTime> itTime(_buffer[sys].keys());
1407 while (itTime.hasNext()) {
1408 bncTime epoTime = itTime.next();
1409 QVector<cmbCorr*>& corrVec = _buffer[sys][epoTime].corrs;
1410 QMutableVectorIterator<cmbCorr*> it(corrVec);
1411 while (it.hasNext()) {
1412 cmbCorr* corr = it.next();
1413 if (acName == corr->_acName) {
1414 delete corr;
1415 it.remove();
1416 }
1417 }
1418 }
1419
1420 // Reset Satellite Offsets
1421 // -----------------------
1422 if (_method == filter) {
1423 for (int iPar = 1; iPar <= _params[sys].size(); iPar++) {
1424 cmbParam* pp = _params[sys][iPar-1];
1425 if (pp->AC == acName && pp->type == cmbParam::offACSat) {
1426 pp->xx = 0.0;
1427 _QQ[sys].Row(iPar) = 0.0;
1428 _QQ[sys].Column(iPar) = 0.0;
1429 _QQ[sys](iPar,iPar) = pp->sig0 * pp->sig0;
1430 }
1431 }
1432 }
1433 }
1434}
Note: See TracBrowser for help on using the repository browser.