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

Last change on this file since 10495 was 10495, checked in by stuerze, 6 months ago

changes regarding Combi

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