source: ntrip/trunk/BNC/combination/bnccomb.cpp@ 3471

Last change on this file since 3471 was 3471, checked in by mervart, 13 years ago
File size: 21.5 KB
RevLine 
[2898]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
[2933]17#include <newmatio.h>
[2921]18#include <iomanip>
[3214]19#include <sstream>
[2898]20
21#include "bnccomb.h"
22#include "bncapp.h"
[3201]23#include "upload/bncrtnetdecoder.h"
[2910]24#include "bncsettings.h"
[2933]25#include "bncmodel.h"
[2988]26#include "bncutils.h"
[3027]27#include "bncpppclient.h"
[3181]28#include "bncsp3.h"
[3052]29#include "bncantex.h"
[3055]30#include "bnctides.h"
[2898]31
[3455]32const int moduloTime = 10;
[2898]33
[3455]34const double sig0_offAC = 1000.0;
35const double sig0_offACSat = 100.0;
[3468]36const double sigP_offACSat = 0.01;
[3455]37const double sig0_clkSat = 100.0;
38
39const double sigObs = 0.05;
40
[3134]41const int MAXPRN_GPS = 32;
42
[3455]43using namespace std;
44
[2898]45// Constructor
46////////////////////////////////////////////////////////////////////////////
[3429]47cmbParam::cmbParam(parType type_, int index_,
48 const QString& ac_, const QString& prn_) {
[3032]49
[3433]50 type = type_;
51 index = index_;
52 AC = ac_;
53 prn = prn_;
54 xx = 0.0;
[3455]55 eph = 0;
[3429]56
57 if (type == offAC) {
[3455]58 epoSpec = true;
59 sig0 = sig0_offAC;
60 sigP = sig0;
[3429]61 }
62 else if (type == offACSat) {
[3455]63 epoSpec = false;
64 sig0 = sig0_offACSat;
65 sigP = sigP_offACSat;
[3429]66 }
67 else if (type == clkSat) {
[3455]68 epoSpec = true;
69 sig0 = sig0_clkSat;
70 sigP = sig0;
[3429]71 }
[2933]72}
73
74// Destructor
75////////////////////////////////////////////////////////////////////////////
76cmbParam::~cmbParam() {
77}
78
79// Partial
80////////////////////////////////////////////////////////////////////////////
[3429]81double cmbParam::partial(const QString& AC_, const QString& prn_) {
[2933]82
[3429]83 if (type == offAC) {
[3032]84 if (AC == AC_) {
[2934]85 return 1.0;
86 }
87 }
[3429]88 else if (type == offACSat) {
89 if (AC == AC_ && prn == prn_) {
[2933]90 return 1.0;
91 }
92 }
[3429]93 else if (type == clkSat) {
94 if (prn == prn_) {
[2933]95 return 1.0;
96 }
97 }
98
99 return 0.0;
100}
101
[2990]102//
103////////////////////////////////////////////////////////////////////////////
104QString cmbParam::toString() const {
[2933]105
[2990]106 QString outStr;
107
[3429]108 if (type == offAC) {
[2992]109 outStr = "AC offset " + AC;
[2990]110 }
[3429]111 else if (type == offACSat) {
[2992]112 outStr = "Sat Offset " + AC + " " + prn;
[2990]113 }
[3429]114 else if (type == clkSat) {
[2992]115 outStr = "Clk Corr " + prn;
[2990]116 }
[2933]117
[2990]118 return outStr;
119}
120
[2933]121// Constructor
122////////////////////////////////////////////////////////////////////////////
[2898]123bncComb::bncComb() {
[2910]124
125 bncSettings settings;
126
127 QStringList combineStreams = settings.value("combineStreams").toStringList();
[2918]128
[3455]129 _masterMissingEpochs = 0;
130
[3143]131 if (combineStreams.size() >= 1 && !combineStreams[0].isEmpty()) {
[2910]132 QListIterator<QString> it(combineStreams);
133 while (it.hasNext()) {
134 QStringList hlp = it.next().split(" ");
[2918]135 cmbAC* newAC = new cmbAC();
136 newAC->mountPoint = hlp[0];
137 newAC->name = hlp[1];
138 newAC->weight = hlp[2].toDouble();
[3453]139 if (_masterOrbitAC.isEmpty()) {
140 _masterOrbitAC = newAC->name;
141 }
[3429]142 _ACs.append(newAC);
[2910]143 }
144 }
[2927]145
[3211]146 _rtnetDecoder = 0;
[3201]147
[2990]148 connect(this, SIGNAL(newMessage(QByteArray,bool)),
149 ((bncApp*)qApp), SLOT(slotMessage(const QByteArray,bool)));
[2933]150
[3470]151 // Combination Method
152 // ------------------
153 if (settings.value("cmbMethod").toString() == "Filter") {
154 _method = filter;
155 }
156 else {
157 _method = singleEpoch;
158 }
[3032]159
160 // Initialize Parameters (model: Clk_Corr = AC_Offset + Sat_Offset + Clk)
161 // ----------------------------------------------------------------------
[3470]162 if (_method == filter) {
163 int nextPar = 0;
164 QListIterator<cmbAC*> it(_ACs);
165 while (it.hasNext()) {
166 cmbAC* AC = it.next();
167 _params.push_back(new cmbParam(cmbParam::offAC, ++nextPar, AC->name, ""));
168 for (int iGps = 1; iGps <= MAXPRN_GPS; iGps++) {
169 QString prn = QString("G%1").arg(iGps, 2, 10, QChar('0'));
170 _params.push_back(new cmbParam(cmbParam::offACSat, ++nextPar,
171 AC->name, prn));
172 }
173 }
[3134]174 for (int iGps = 1; iGps <= MAXPRN_GPS; iGps++) {
175 QString prn = QString("G%1").arg(iGps, 2, 10, QChar('0'));
[3470]176 _params.push_back(new cmbParam(cmbParam::clkSat, ++nextPar, "", prn));
[2991]177 }
[3470]178
179 // Initialize Variance-Covariance Matrix
180 // -------------------------------------
181 _QQ.ReSize(_params.size());
182 _QQ = 0.0;
183 for (int iPar = 1; iPar <= _params.size(); iPar++) {
184 cmbParam* pp = _params[iPar-1];
185 _QQ(iPar,iPar) = pp->sig0 * pp->sig0;
186 }
[2991]187 }
[2933]188
[3052]189 // ANTEX File
190 // ----------
191 _antex = 0;
192 QString antexFileName = settings.value("pppAntex").toString();
193 if (!antexFileName.isEmpty()) {
194 _antex = new bncAntex();
195 if (_antex->readFile(antexFileName) != success) {
196 emit newMessage("wrong ANTEX file", true);
197 delete _antex;
198 _antex = 0;
199 }
200 }
[3138]201
[3329]202 // Maximal Residuum
203 // ----------------
204 _MAXRES = settings.value("cmbMaxres").toDouble();
205 if (_MAXRES <= 0.0) {
206 _MAXRES = 999.0;
207 }
[2898]208}
209
210// Destructor
211////////////////////////////////////////////////////////////////////////////
212bncComb::~bncComb() {
[3429]213 QListIterator<cmbAC*> icAC(_ACs);
214 while (icAC.hasNext()) {
215 delete icAC.next();
[2918]216 }
[3201]217 delete _rtnetDecoder;
[3052]218 delete _antex;
[3296]219 for (int iPar = 1; iPar <= _params.size(); iPar++) {
220 delete _params[iPar-1];
221 }
[3434]222 QVectorIterator<cmbCorr*> itCorr(corrs());
[3429]223 while (itCorr.hasNext()) {
224 delete itCorr.next();
225 }
[2898]226}
227
[2928]228// Read and store one correction line
[2898]229////////////////////////////////////////////////////////////////////////////
230void bncComb::processCorrLine(const QString& staID, const QString& line) {
[2906]231 QMutexLocker locker(&_mutex);
[2913]232
[3431]233 // Find the AC Name
234 // ----------------
235 QString acName;
236 QListIterator<cmbAC*> icAC(_ACs);
237 while (icAC.hasNext()) {
238 cmbAC* AC = icAC.next();
239 if (AC->mountPoint == staID) {
240 acName = AC->name;
241 break;
242 }
243 }
244 if (acName.isEmpty()) {
245 return;
246 }
247
[2918]248 // Read the Correction
249 // -------------------
[3429]250 cmbCorr* newCorr = new cmbCorr();
[3431]251 newCorr->acName = acName;
[2913]252 if (!newCorr->readLine(line) == success) {
253 delete newCorr;
254 return;
255 }
256
[3274]257 // Check Modulo Time
258 // -----------------
259 if (int(newCorr->tt.gpssec()) % moduloTime != 0.0) {
260 delete newCorr;
261 return;
262 }
263
[3299]264 // Delete old corrections
265 // ----------------------
[3437]266 if (_resTime.valid() && newCorr->tt <= _resTime) {
[2924]267 delete newCorr;
268 return;
269 }
270
[3029]271 // Check the Ephemeris
272 //--------------------
[2986]273 if (_eph.find(newCorr->prn) == _eph.end()) {
274 delete newCorr;
275 return;
276 }
277 else {
278 t_eph* lastEph = _eph[newCorr->prn]->last;
279 t_eph* prevEph = _eph[newCorr->prn]->prev;
[3029]280 if (lastEph && lastEph->IOD() == newCorr->iod) {
281 newCorr->eph = lastEph;
[2986]282 }
[3029]283 else if (prevEph && prevEph->IOD() == newCorr->iod) {
284 newCorr->eph = prevEph;
[3429]285 switchToLastEph(lastEph, newCorr);
[3029]286 }
287 else {
[2986]288 delete newCorr;
289 return;
290 }
291 }
292
[3435]293 // Process previous Epoch(s)
294 // -------------------------
295 QListIterator<bncTime> itTime(_buffer.keys());
296 while (itTime.hasNext()) {
297 bncTime epoTime = itTime.next();
[3438]298 if (epoTime < newCorr->tt - moduloTime) {
[3435]299 _resTime = epoTime;
300 processEpoch();
301 }
[2927]302 }
303
[3429]304 // Merge or add the correction
305 // ---------------------------
[3435]306 QVector<cmbCorr*>& corrs = _buffer[newCorr->tt].corrs;
[3429]307 cmbCorr* existingCorr = 0;
[3435]308 QVectorIterator<cmbCorr*> itCorr(corrs);
[3429]309 while (itCorr.hasNext()) {
310 cmbCorr* hlp = itCorr.next();
311 if (hlp->prn == newCorr->prn && hlp->acName == newCorr->prn) {
312 existingCorr = hlp;
[2918]313 break;
314 }
315 }
[3429]316 if (existingCorr) {
[3298]317 delete newCorr;
[3429]318 existingCorr->readLine(line); // merge (multiple messages)
[2918]319 }
[2919]320 else {
[3435]321 corrs.append(newCorr);
[2919]322 }
[2898]323}
324
[2986]325// Change the correction so that it refers to last received ephemeris
326////////////////////////////////////////////////////////////////////////////
[3029]327void bncComb::switchToLastEph(const t_eph* lastEph, t_corr* corr) {
[3028]328
[3429]329 if (corr->eph == lastEph) {
330 return;
331 }
332
[2987]333 ColumnVector oldXC(4);
334 ColumnVector oldVV(3);
[3029]335 corr->eph->position(corr->tt.gpsw(), corr->tt.gpssec(),
336 oldXC.data(), oldVV.data());
[2988]337
[2987]338 ColumnVector newXC(4);
339 ColumnVector newVV(3);
[3029]340 lastEph->position(corr->tt.gpsw(), corr->tt.gpssec(),
[2987]341 newXC.data(), newVV.data());
[2988]342
343 ColumnVector dX = newXC.Rows(1,3) - oldXC.Rows(1,3);
[2989]344 ColumnVector dV = newVV - oldVV;
345 double dC = newXC(4) - oldXC(4);
346
[2988]347 ColumnVector dRAO(3);
348 XYZ_to_RSW(newXC.Rows(1,3), newVV, dX, dRAO);
[2989]349
350 ColumnVector dDotRAO(3);
351 XYZ_to_RSW(newXC.Rows(1,3), newVV, dV, dDotRAO);
352
[3455]353 QString msg = "switch corr " + corr->prn
[3029]354 + QString(" %1 -> %2 %3").arg(corr->iod,3)
[3013]355 .arg(lastEph->IOD(),3).arg(dC*t_CST::c, 8, 'f', 4);
356
[3029]357 emit newMessage(msg.toAscii(), false);
[3028]358
[3029]359 corr->iod = lastEph->IOD();
360 corr->eph = lastEph;
[3031]361 corr->rao += dRAO;
362 corr->dotRao += dDotRAO;
[3029]363 corr->dClk -= dC;
[2986]364}
365
[3429]366// Process Epoch
[2928]367////////////////////////////////////////////////////////////////////////////
[3429]368void bncComb::processEpoch() {
[3470]369 if (_method == filter) {
370 return processEpoch_filter();
371 }
372 else {
373 return processEpoch_singleEpoch();
374 }
375}
[2918]376
[3470]377// Process Epoch - Filter Method
378////////////////////////////////////////////////////////////////////////////
379void bncComb::processEpoch_filter() {
[3429]380 int nPar = _params.size();
381
[2990]382 _log.clear();
383
384 QTextStream out(&_log, QIODevice::WriteOnly);
385
[3471]386 out << endl << "Combination (Filter):" << endl
387 << "---------------------------------------" << endl;
[2990]388
[3433]389 // Observation Statistics
390 // ----------------------
[3455]391 bool masterPresent = false;
[3433]392 QListIterator<cmbAC*> icAC(_ACs);
393 while (icAC.hasNext()) {
394 cmbAC* AC = icAC.next();
395 AC->numObs = 0;
[3434]396 QVectorIterator<cmbCorr*> itCorr(corrs());
[3433]397 while (itCorr.hasNext()) {
398 cmbCorr* corr = itCorr.next();
399 if (corr->acName == AC->name) {
400 AC->numObs += 1;
[3453]401 if (AC->name == _masterOrbitAC) {
402 masterPresent = true;
403 }
[3433]404 }
405 }
406 out << AC->name.toAscii().data() << ": " << AC->numObs << endl;
407 }
408
[3453]409 // If Master not present, switch to another one
410 // --------------------------------------------
[3456]411 if (masterPresent) {
412 _masterMissingEpochs = 0;
413 }
414 else {
[3455]415 ++_masterMissingEpochs;
[3459]416 if (_masterMissingEpochs < 10) {
[3457]417 out << "Missing Master, Epoch skipped" << endl;
[3455]418 _buffer.remove(_resTime);
419 emit newMessage(_log, false);
420 return;
[3453]421 }
[3455]422 else {
423 _masterMissingEpochs = 0;
424 QListIterator<cmbAC*> icAC(_ACs);
425 while (icAC.hasNext()) {
426 cmbAC* AC = icAC.next();
427 if (AC->numObs > 0) {
428 out << "Switching Master AC "
429 << _masterOrbitAC.toAscii().data() << " --> "
430 << AC->name.toAscii().data() << " "
431 << _resTime.datestr().c_str() << " "
432 << _resTime.timestr().c_str() << endl;
433 _masterOrbitAC = AC->name;
434 break;
435 }
[3452]436 }
437 }
438 }
439
[3429]440 // Prediction Step
441 // ---------------
[2933]442 ColumnVector x0(nPar);
443 for (int iPar = 1; iPar <= _params.size(); iPar++) {
[3455]444 cmbParam* pp = _params[iPar-1];
445 if (pp->epoSpec) {
[3451]446 pp->xx = 0.0;
[3455]447 _QQ.Row(iPar) = 0.0;
448 _QQ.Column(iPar) = 0.0;
449 _QQ(iPar,iPar) = pp->sig0 * pp->sig0;
[3451]450 }
[3455]451 else {
452 _QQ(iPar,iPar) += pp->sigP * pp->sigP;
453 }
[2933]454 x0(iPar) = pp->xx;
455 }
456
[3455]457 SymmetricMatrix QQ_sav = _QQ;
[3441]458
[3455]459 ColumnVector dx;
[3429]460 QMap<QString, t_corr*> resCorr;
[3455]461
462 // Update and outlier detection loop
463 // ---------------------------------
464 while (true) {
[3135]465
[3455]466 Matrix AA;
467 ColumnVector ll;
468 DiagonalMatrix PP;
[3429]469
[3455]470 if (createAmat(AA, ll, PP, x0, resCorr) != success) {
471 _buffer.remove(_resTime);
472 emit newMessage(_log, false);
473 return;
[3441]474 }
[2933]475
[3429]476 bncModel::kalman(AA, ll, PP, _QQ, dx);
477 ColumnVector vv = ll - AA * dx;
[3080]478
[3429]479 int maxResIndex;
480 double maxRes = vv.maximum_absolute_value1(maxResIndex);
481 out.setRealNumberNotation(QTextStream::FixedNotation);
482 out.setRealNumberPrecision(3);
483 out << _resTime.datestr().c_str() << " " << _resTime.timestr().c_str()
[3432]484 << " Maximum Residuum " << maxRes << ' '
[3434]485 << corrs()[maxResIndex-1]->acName << ' ' << corrs()[maxResIndex-1]->prn;
[3080]486
[3432]487 if (maxRes > _MAXRES) {
488 for (int iPar = 1; iPar <= _params.size(); iPar++) {
489 cmbParam* pp = _params[iPar-1];
490 if (pp->type == cmbParam::offACSat &&
[3434]491 pp->AC == corrs()[maxResIndex-1]->acName &&
492 pp->prn == corrs()[maxResIndex-1]->prn) {
[3432]493 QQ_sav.Row(iPar) = 0.0;
494 QQ_sav.Column(iPar) = 0.0;
[3455]495 QQ_sav(iPar,iPar) = pp->sig0 * pp->sig0;
[3432]496 }
497 }
498
499 out << " Outlier" << endl;
500 _QQ = QQ_sav;
[3455]501 corrs().remove(maxResIndex-1);
[3432]502 }
503 else {
504 out << " OK" << endl;
505 break;
506 }
507 }
508
509 // Update Parameter Values
510 // -----------------------
[3429]511 for (int iPar = 1; iPar <= _params.size(); iPar++) {
512 cmbParam* pp = _params[iPar-1];
513 pp->xx += dx(iPar);
514 if (pp->type == cmbParam::clkSat) {
515 if (resCorr.find(pp->prn) != resCorr.end()) {
516 resCorr[pp->prn]->dClk = pp->xx / t_CST::c;
[3080]517 }
518 }
[3429]519 out << _resTime.datestr().c_str() << " "
520 << _resTime.timestr().c_str() << " ";
521 out.setRealNumberNotation(QTextStream::FixedNotation);
522 out.setFieldWidth(8);
523 out.setRealNumberPrecision(4);
524 out << pp->toString() << " "
525 << pp->xx << " +- " << sqrt(_QQ(pp->index,pp->index)) << endl;
526 out.setFieldWidth(0);
[2918]527 }
[2919]528
[3432]529 // Print Results
530 // -------------
[3429]531 printResults(out, resCorr);
532 dumpResults(resCorr);
[2928]533
[2990]534 emit newMessage(_log, false);
[3296]535
[3432]536 // Delete Data
537 // -----------
[3434]538 _buffer.remove(_resTime);
[2918]539}
[3011]540
[3201]541// Print results
[3011]542////////////////////////////////////////////////////////////////////////////
[3429]543void bncComb::printResults(QTextStream& out,
[3011]544 const QMap<QString, t_corr*>& resCorr) {
545
546 QMapIterator<QString, t_corr*> it(resCorr);
547 while (it.hasNext()) {
548 it.next();
549 t_corr* corr = it.value();
[3029]550 const t_eph* eph = corr->eph;
[3015]551 if (eph) {
552 double xx, yy, zz, cc;
[3429]553 eph->position(_resTime.gpsw(), _resTime.gpssec(), xx, yy, zz, cc);
[3015]554
[3429]555 out << _resTime.datestr().c_str() << " "
556 << _resTime.timestr().c_str() << " ";
[3015]557 out.setFieldWidth(3);
558 out << "Full Clock " << corr->prn << " " << corr->iod << " ";
559 out.setFieldWidth(14);
560 out << (cc + corr->dClk) * t_CST::c << endl;
[3370]561 out.setFieldWidth(0);
[3015]562 }
563 else {
564 out << "bncComb::printResuls bug" << endl;
565 }
[3011]566 }
567}
[3202]568
569// Send results to RTNet Decoder and directly to PPP Client
570////////////////////////////////////////////////////////////////////////////
[3429]571void bncComb::dumpResults(const QMap<QString, t_corr*>& resCorr) {
[3202]572
[3214]573 ostringstream out; out.setf(std::ios::fixed);
[3216]574 QStringList corrLines;
[3214]575
576 unsigned year, month, day, hour, minute;
577 double sec;
[3429]578 _resTime.civil_date(year, month, day);
579 _resTime.civil_time(hour, minute, sec);
[3214]580
581 out << "* "
582 << setw(4) << year << " "
583 << setw(2) << month << " "
584 << setw(2) << day << " "
585 << setw(2) << hour << " "
586 << setw(2) << minute << " "
587 << setw(12) << setprecision(8) << sec << " "
588 << endl;
589
[3202]590 QMapIterator<QString, t_corr*> it(resCorr);
591 while (it.hasNext()) {
592 it.next();
593 t_corr* corr = it.value();
594
[3214]595 double dT = 60.0;
[3202]596
[3214]597 for (int iTime = 1; iTime <= 2; iTime++) {
598
[3429]599 bncTime time12 = (iTime == 1) ? _resTime : _resTime + dT;
[3214]600
601 ColumnVector xc(4);
602 ColumnVector vv(3);
603 corr->eph->position(time12.gpsw(), time12.gpssec(),
604 xc.data(), vv.data());
605 bncPPPclient::applyCorr(time12, corr, xc, vv);
606
607 // Relativistic Correction
608 // -----------------------
609 double relCorr = - 2.0 * DotProduct(xc.Rows(1,3),vv) / t_CST::c / t_CST::c;
610 xc(4) -= relCorr;
611
612 // Code Biases
613 // -----------
614 double dcbP1C1 = 0.0;
615 double dcbP1P2 = 0.0;
616
617 // Correction Phase Center --> CoM
618 // -------------------------------
619 ColumnVector dx(3); dx = 0.0;
620 if (_antex) {
621 double Mjd = time12.mjd() + time12.daysec()/86400.0;
622 if (_antex->satCoMcorrection(corr->prn, Mjd, xc.Rows(1,3), dx) == success) {
623 xc(1) -= dx(1);
624 xc(2) -= dx(2);
625 xc(3) -= dx(3);
626 }
627 else {
628 cout << "antenna not found" << endl;
629 }
630 }
631
632 if (iTime == 1) {
633 out << 'P' << corr->prn.toAscii().data()
634 << setw(14) << setprecision(6) << xc(1) / 1000.0
635 << setw(14) << setprecision(6) << xc(2) / 1000.0
636 << setw(14) << setprecision(6) << xc(3) / 1000.0
637 << setw(14) << setprecision(6) << xc(4) * 1e6
638 << setw(14) << setprecision(6) << relCorr * 1e6
639 << setw(8) << setprecision(3) << dx(1)
640 << setw(8) << setprecision(3) << dx(2)
641 << setw(8) << setprecision(3) << dx(3)
642 << setw(8) << setprecision(3) << dcbP1C1
643 << setw(8) << setprecision(3) << dcbP1P2
644 << setw(6) << setprecision(1) << dT;
[3216]645
[3218]646 QString line;
[3217]647 int messageType = COTYPE_GPSCOMBINED;
[3218]648 int updateInt = 0;
649 line.sprintf("%d %d %d %.1f %s"
650 " %3d"
651 " %8.3f %8.3f %8.3f %8.3f"
652 " %10.5f %10.5f %10.5f %10.5f"
[3262]653 " %10.5f %10.5f %10.5f %10.5f INTERNAL",
[3218]654 messageType, updateInt, time12.gpsw(), time12.gpssec(),
655 corr->prn.toAscii().data(),
656 corr->iod,
[3219]657 corr->dClk * t_CST::c,
[3218]658 corr->rao[0],
659 corr->rao[1],
660 corr->rao[2],
[3219]661 corr->dotDClk * t_CST::c,
[3218]662 corr->dotRao[0],
663 corr->dotRao[1],
664 corr->dotRao[2],
[3262]665 corr->dotDotDClk * t_CST::c,
666 corr->dotDotRao[0],
667 corr->dotDotRao[1],
668 corr->dotDotRao[2]);
[3220]669 corrLines << line;
[3214]670 }
671 else {
672 out << setw(14) << setprecision(6) << xc(1) / 1000.0
673 << setw(14) << setprecision(6) << xc(2) / 1000.0
674 << setw(14) << setprecision(6) << xc(3) / 1000.0 << endl;
675 }
676 }
677
678 delete corr;
679 }
[3228]680 out << "EOE" << endl; // End Of Epoch flag
[3214]681
[3215]682 if (!_rtnetDecoder) {
683 _rtnetDecoder = new bncRtnetDecoder();
684 }
[3221]685
[3215]686 vector<string> errmsg;
[3221]687 _rtnetDecoder->Decode((char*) out.str().data(), out.str().size(), errmsg);
[3215]688
[3216]689 // Optionally send new Corrections to PPP
690 // --------------------------------------
691 bncApp* app = (bncApp*) qApp;
692 if (app->_bncPPPclient) {
693 app->_bncPPPclient->slotNewCorrections(corrLines);
694 }
[3202]695}
[3455]696
697// Create First Design Matrix and Vector of Measurements
698////////////////////////////////////////////////////////////////////////////
699t_irc bncComb::createAmat(Matrix& AA, ColumnVector& ll, DiagonalMatrix& PP,
700 const ColumnVector& x0,
701 QMap<QString, t_corr*>& resCorr) {
702
703 unsigned nPar = _params.size();
704 unsigned nObs = corrs().size();
705
706 if (nObs == 0) {
707 return failure;
708 }
709
[3461]710 const int nCon = 1 + MAXPRN_GPS;
[3455]711
712 AA.ReSize(nObs+nCon, nPar); AA = 0.0;
713 ll.ReSize(nObs+nCon); ll = 0.0;
714 PP.ReSize(nObs+nCon); PP = 1.0 / (sigObs * sigObs);
715
716 int iObs = 0;
717
718 QVectorIterator<cmbCorr*> itCorr(corrs());
719 while (itCorr.hasNext()) {
720 cmbCorr* corr = itCorr.next();
721 QString prn = corr->prn;
722 switchToLastEph(_eph[prn]->last, corr);
723 ++iObs;
724
725 if (corr->acName == _masterOrbitAC && resCorr.find(prn) == resCorr.end()) {
726 resCorr[prn] = new t_corr(*corr);
727 }
728
729 for (int iPar = 1; iPar <= _params.size(); iPar++) {
730 cmbParam* pp = _params[iPar-1];
731 AA(iObs, iPar) = pp->partial(corr->acName, prn);
732 }
733
734 ll(iObs) = corr->dClk * t_CST::c - DotProduct(AA.Row(iObs), x0);
735 }
736
737 // Regularization
738 // --------------
739 const double Ph = 1.e6;
[3461]740 PP(nObs+1) = Ph;
[3455]741 for (int iPar = 1; iPar <= _params.size(); iPar++) {
742 cmbParam* pp = _params[iPar-1];
[3465]743 if ( AA.Column(iPar).maximum_absolute_value() > 0.0 &&
744 pp->type == cmbParam::clkSat ) {
745 AA(nObs+1, iPar) = 1.0;
[3455]746 }
747 }
[3461]748 int iCond = 1;
749 for (int iGps = 1; iGps <= MAXPRN_GPS; iGps++) {
750 QString prn = QString("G%1").arg(iGps, 2, 10, QChar('0'));
751 ++iCond;
752 PP(nObs+iCond) = Ph;
753 for (int iPar = 1; iPar <= _params.size(); iPar++) {
754 cmbParam* pp = _params[iPar-1];
[3465]755 if ( AA.Column(iPar).maximum_absolute_value() > 0.0 &&
756 pp->type == cmbParam::offACSat &&
757 pp->prn == prn) {
[3461]758 AA(nObs+iCond, iPar) = 1.0;
[3455]759 }
760 }
761 }
762
763 return success;
764}
[3470]765
766// Process Epoch - Single-Epoch Method
767////////////////////////////////////////////////////////////////////////////
768void bncComb::processEpoch_singleEpoch() {
769
[3471]770 _log.clear();
[3470]771
[3471]772 QTextStream out(&_log, QIODevice::WriteOnly);
773
774 out << endl << "Combination (Single-Epoch):" << endl
775 << "---------------------------------------------" << endl;
776
777 // Observation Statistics
778 // ----------------------
779 QListIterator<cmbAC*> icAC(_ACs);
780 while (icAC.hasNext()) {
781 cmbAC* AC = icAC.next();
782 AC->numObs = 0;
783 QVectorIterator<cmbCorr*> itCorr(corrs());
784 while (itCorr.hasNext()) {
785 cmbCorr* corr = itCorr.next();
786 if (corr->acName == AC->name) {
787 AC->numObs += 1;
788 }
789 }
790 out << AC->name.toAscii().data() << ": " << AC->numObs << endl;
791 }
792
793 QMap<QString, t_corr*> resCorr;
794
795
796 // Print Results
797 // -------------
798 printResults(out, resCorr);
799 dumpResults(resCorr);
800
801 emit newMessage(_log, false);
802
803 // Delete Data
804 // -----------
805 _buffer.remove(_resTime);
[3470]806}
Note: See TracBrowser for help on using the repository browser.