Index: /trunk/BNC/src/PPP/pppFilter.cpp =================================================================== --- /trunk/BNC/src/PPP/pppFilter.cpp (revision 9641) +++ /trunk/BNC/src/PPP/pppFilter.cpp (revision 9642) @@ -33,8 +33,8 @@ // Constructor //////////////////////////////////////////////////////////////////////////// -t_pppFilter::t_pppFilter(t_pppObsPool* obsPool) { +t_pppFilter::t_pppFilter(t_pppObsPool *obsPool) { _numSat = 0; _obsPool = obsPool; - _refPrn = t_prn(); + _refPrn = t_prn(); _datumTrafo = new t_datumTrafo(); } @@ -49,14 +49,14 @@ //////////////////////////////////////////////////////////////////////////// t_irc t_pppFilter::processEpoch() { - _numSat = 0; + _numSat = 0; const double maxSolGap = 60.0; // Vector of all Observations // -------------------------- - t_pppObsPool::t_epoch* epoch = _obsPool->lastEpoch(); + t_pppObsPool::t_epoch *epoch = _obsPool->lastEpoch(); if (!epoch) { return failure; } - vector& allObs = epoch->obsVector(); + vector &allObs = epoch->obsVector(); // Time of the Epoch @@ -64,7 +64,6 @@ _epoTime = epoch->epoTime(); - if (!_firstEpoTime.valid() || - !_lastEpoTimeOK.valid() || - (maxSolGap > 0.0 && _epoTime - _lastEpoTimeOK > maxSolGap)) { + if (!_firstEpoTime.valid() || !_lastEpoTimeOK.valid() + || (maxSolGap > 0.0 && _epoTime - _lastEpoTimeOK > maxSolGap)) { _firstEpoTime = _epoTime; } @@ -72,5 +71,5 @@ string epoTimeStr = string(_epoTime); - const QMap& refSatMap = epoch->refSatMap(); + const QMap &refSatMap = epoch->refSatMap(); //LOG << "processEpoch: printParams before set" << endl; @@ -84,5 +83,5 @@ if (OPT->_obsModelType == OPT->DCMcodeBias || - OPT->_obsModelType == OPT->DCMphaseBias) { + OPT->_obsModelType == OPT->DCMphaseBias) { #ifdef BNC_DEBUG_PPP _parlist.printParams(_epoTime); @@ -91,5 +90,5 @@ // Status Vector, Variance-Covariance Matrix // ----------------------------------------- - ColumnVector xFltOld = _xFlt; + ColumnVector xFltOld = _xFlt; SymmetricMatrix QFltOld = _QFlt; setStateVectorAndVarCovMatrix(xFltOld, QFltOld); @@ -99,7 +98,7 @@ bool preProcessing = false; if (OPT->_obsModelType == OPT->DCMcodeBias || - OPT->_obsModelType == OPT->DCMphaseBias) { + OPT->_obsModelType == OPT->DCMphaseBias) { preProcessing = true; - const QList& usedSystems = _parlist.usedSystems(); + const QList &usedSystems = _parlist.usedSystems(); for (int iSys = 0; iSys < usedSystems.size(); iSys++) { char sys = usedSystems[iSys]; @@ -108,5 +107,5 @@ _refPrn = refSatMap[sys]->prn(); } - vector obsVector; + vector obsVector; for (unsigned jj = 0; jj < allObs.size(); jj++) { if (allObs[jj]->prn().system() == sys) { @@ -118,5 +117,5 @@ } if (processSystem(OPT->LCs(sys), obsVector, _refPrn, - epoch->pseudoObsIono(), preProcessing) != success) { + epoch->pseudoObsIono(), preProcessing) != success) { LOG << sys << ": processSystem != success (pre-processing)" << endl; _xFlt = xFltOld; @@ -133,6 +132,5 @@ _QFlt = QFltOld; return success; - } - else if (!_obsPool->refSatChangeRequired()) { + } else if (!_obsPool->refSatChangeRequired()) { initDatumTransformation(allObs, epoch->pseudoObsIono()); } @@ -142,5 +140,5 @@ // ------------------------------------ preProcessing = false; - const QList& usedSystems = _parlist. usedSystems(); + const QList &usedSystems = _parlist.usedSystems(); for (int iSys = 0; iSys < usedSystems.size(); iSys++) { char sys = usedSystems[iSys]; @@ -160,10 +158,11 @@ _datumTrafo->setFirstSystem(sys); } - LOG << epoTimeStr << " SATNUM " << sys << ' ' << right << setw(2) << num << endl; - if (processSystem(OPT->LCs(sys), obsVector, _refPrn, - epoch->pseudoObsIono(), preProcessing) != success) { + LOG << epoTimeStr << " SATNUM " << sys << ' ' << right << setw(2) << num + << endl; + if (processSystem(OPT->LCs(sys), obsVector, _refPrn, epoch->pseudoObsIono(), + preProcessing) != success) { LOG << "processSystem != success (fin-processing)" << endl; if (OPT->_obsModelType == OPT->DCMcodeBias || - OPT->_obsModelType == OPT->DCMphaseBias) { + OPT->_obsModelType == OPT->DCMphaseBias) { _xFlt = xFltOld; _QFlt = QFltOld; @@ -179,5 +178,5 @@ cmpDOP(allObs, refSatMap); _parlist.printResult(_epoTime, _QFlt, _xFlt); - _lastEpoTimeOK = _epoTime; // remember time of last successful epoch processing + _lastEpoTimeOK = _epoTime; // remember time of last successful epoch processing return success; } @@ -185,9 +184,7 @@ // Process Selected LCs //////////////////////////////////////////////////////////////////////////// -t_irc t_pppFilter::processSystem(const vector& LCs, - const vector& obsVector, - const t_prn& refPrn, - bool pseudoObsIonoAvailable, - bool preProcessing) { +t_irc t_pppFilter::processSystem(const vector &LCs, + const vector &obsVector, const t_prn &refPrn, + bool pseudoObsIonoAvailable, bool preProcessing) { LOG.setf(ios::fixed); char sys = refPrn.system(); @@ -202,13 +199,13 @@ } - ColumnVector xSav = _xFlt; - SymmetricMatrix QSav = _QFlt; - string epoTimeStr = string(_epoTime); - const vector& params = _parlist.params(); + ColumnVector xSav = _xFlt; + SymmetricMatrix QSav = _QFlt; + string epoTimeStr = string(_epoTime); + const vector ¶ms = _parlist.params(); unsigned nPar = _parlist.nPar(); - unsigned usedLCs = LCs.size(); + unsigned usedLCs = LCs.size(); if (OPT->_pseudoObsIono && !pseudoObsIonoAvailable) { - usedLCs -= 1; // GIM not used + usedLCs -= 1; // GIM not used } // max Obs @@ -230,11 +227,12 @@ // First-Design Matrix, Terms Observed-Computed, Weight Matrix // ----------------------------------------------------------- - Matrix AA(maxObs, nPar); - ColumnVector ll(maxObs); - DiagonalMatrix PP(maxObs); PP = 0.0; + Matrix AA(maxObs, nPar); + ColumnVector ll(maxObs); + DiagonalMatrix PP(maxObs); + PP = 0.0; int iObs = -1; vector usedObs; - vector usedTypes; + vector usedTypes; // Real Observations @@ -242,5 +240,5 @@ double nSat = 0; for (unsigned ii = 0; ii < obsVector.size(); ii++) { - t_pppSatObs* obs = obsVector[ii]; + t_pppSatObs *obs = obsVector[ii]; if (iOutlier == 0 && !preProcessing) { obs->resetOutlier(); @@ -250,13 +248,16 @@ for (unsigned jj = 0; jj < usedLCs; jj++) { const t_lc::type tLC = LCs[jj]; - if (tLC == t_lc::GIM) {continue;} + if (tLC == t_lc::GIM) { + continue; + } ++iObs; usedObs.push_back(obs); usedTypes.push_back(tLC); for (unsigned iPar = 0; iPar < nPar; iPar++) { - const t_pppParam* par = params[iPar]; + const t_pppParam *par = params[iPar]; AA[iObs][iPar] = par->partial(_epoTime, obs, tLC, refPrn); } - ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1)); + ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) + - DotProduct(_x0, AA.Row(iObs + 1)); PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC)); } @@ -271,10 +272,8 @@ } - if ((!iOutlier) && - (OPT->_obsModelType == OPT->DCMcodeBias || - OPT->_obsModelType == OPT->DCMphaseBias) && - (!preProcessing)) { - _datumTrafo->updateIndices(sys, iObs+1); - _datumTrafo->prepareAA(AA.SubMatrix(1, iObs+1 , 1, _parlist.nPar()), 1); + if ((!iOutlier) && (OPT->_obsModelType == OPT->DCMcodeBias || + OPT->_obsModelType == OPT->DCMphaseBias) && (!preProcessing)) { + _datumTrafo->updateIndices(sys, iObs + 1); + _datumTrafo->prepareAA(AA.SubMatrix(1, iObs + 1, 1, _parlist.nPar()), 1); } @@ -283,5 +282,5 @@ if (OPT->_pseudoObsIono && pseudoObsIonoAvailable) { for (unsigned ii = 0; ii < obsVector.size(); ii++) { - t_pppSatObs* obs = obsVector[ii]; + t_pppSatObs *obs = obsVector[ii]; if (!obs->outlier()) { for (unsigned jj = 0; jj < usedLCs; jj++) { @@ -289,12 +288,15 @@ if (tLC == t_lc::GIM && !obs->isReference()) { ++iObs; - } else {continue;} + } else { + continue; + } usedObs.push_back(obs); usedTypes.push_back(tLC); for (unsigned iPar = 0; iPar < nPar; iPar++) { - const t_pppParam* par = params[iPar]; + const t_pppParam *par = params[iPar]; AA[iObs][iPar] = par->partial(_epoTime, obs, tLC, refPrn); } - ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1)); + ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) + - DotProduct(_x0, AA.Row(iObs + 1)); PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC)); } @@ -305,7 +307,7 @@ // Truncate matrices // ----------------- - AA = AA.Rows(1, iObs+1); - ll = ll.Rows(1, iObs+1); - PP = PP.SymSubMatrix(1, iObs+1); + AA = AA.Rows(1, iObs + 1); + ll = ll.Rows(1, iObs + 1); + PP = PP.SymSubMatrix(1, iObs + 1); // Kalman update step @@ -316,7 +318,7 @@ // --------------- ColumnVector vv = AA * _xFlt - ll; - double maxOutlier = 0.0; - int maxOutlierIndex = -1; - t_lc::type maxOutlierLC = t_lc::dummy; + double maxOutlier = 0.0; + int maxOutlierIndex = -1; + t_lc::type maxOutlierLC = t_lc::dummy; for (unsigned ii = 0; ii < usedObs.size(); ii++) { const t_lc::type tLC = usedTypes[ii]; @@ -324,7 +326,7 @@ if (res > usedObs[ii]->maxRes(tLC)) { if (res > fabs(maxOutlier)) { - maxOutlier = vv[ii]; + maxOutlier = vv[ii]; maxOutlierIndex = ii; - maxOutlierLC = tLC; + maxOutlierLC = tLC; } } @@ -334,11 +336,10 @@ // -------------------------------------------- if (maxOutlierIndex > -1) { - t_pppSatObs* obs = usedObs[maxOutlierIndex]; - t_pppParam* par = 0; + t_pppSatObs *obs = usedObs[maxOutlierIndex]; + t_pppParam *par = 0; for (unsigned iPar = 0; iPar < nPar; iPar++) { - t_pppParam* hlp = params[iPar]; - if (hlp->type() == t_pppParam::amb && - hlp->prn() == obs->prn() && - hlp->tLC() == usedTypes[maxOutlierIndex]) { + t_pppParam *hlp = params[iPar]; + if (hlp->type() == t_pppParam::amb && hlp->prn() == obs->prn() + && hlp->tLC() == usedTypes[maxOutlierIndex]) { par = hlp; } @@ -346,24 +347,23 @@ if (preProcessing) { // for refSats no ambiguity parameter exists - if ((obs->prn() == refPrn) && - (t_lc::toString(maxOutlierLC) == "l1" || - t_lc::toString(maxOutlierLC) == "l2" )) { - _obsPool->setRefSatChangeRequired(sys, true); - LOG << epoTimeStr << " Outlier (" - << ((preProcessing) ? "pre-processing) " : "fin-processing) ") << t_lc::toString(maxOutlierLC) << ' ' - << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << maxOutlier << endl; - break; - } - else { - obs->setOutlier(); - } - } - else {// fin-processing + if ((obs->prn() == refPrn) + && (t_lc::toString(maxOutlierLC) == "l1" + || t_lc::toString(maxOutlierLC) == "l2")) { + _obsPool->setRefSatChangeRequired(sys, true); + LOG << epoTimeStr << " Outlier (" + << ((preProcessing) ? "pre-processing) " : "fin-processing) ") + << t_lc::toString(maxOutlierLC) << ' ' << obs->prn().toString() + << ' ' << setw(8) << setprecision(4) << maxOutlier << endl; + break; + } else { + obs->setOutlier(); + } + } else { // fin-processing LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' ' - << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << maxOutlier << endl; + << obs->prn().toString() << ' ' << setw(8) << setprecision(4) + << maxOutlier << endl; if (par) { resetAmb(par->prn(), obsVector, &QSav, &xSav); - } - else { + } else { obs->setOutlier(); } @@ -377,11 +377,14 @@ for (unsigned ii = 0; ii < usedObs.size(); ii++) { const t_lc::type tLC = usedTypes[ii]; - t_pppSatObs* obs = usedObs[ii]; + t_pppSatObs *obs = usedObs[ii]; if (tLC == LCs[jj]) { obs->setRes(tLC, vv[ii]); - LOG << epoTimeStr << " RES " - << left << setw(3) << t_lc::toString(tLC) << right << ' '; - if (t_lc::toString(tLC) == "Tz0") {LOG << sys << " ";} - else {LOG << obs->prn().toString();} + LOG << epoTimeStr << " RES " << left << setw(3) + << t_lc::toString(tLC) << right << ' '; + if (t_lc::toString(tLC) == "Tz0") { + LOG << sys << " "; + } else { + LOG << obs->prn().toString(); + } LOG << setw(9) << setprecision(4) << vv[ii] << endl; } @@ -397,19 +400,18 @@ // Cycle-Slip Detection //////////////////////////////////////////////////////////////////////////// -t_irc t_pppFilter::detectCycleSlips(const vector& LCs, - const vector& obsVector, - const t_prn& refPrn, - bool preProcessing) { +t_irc t_pppFilter::detectCycleSlips(const vector &LCs, + const vector &obsVector, const t_prn &refPrn, + bool preProcessing) { const double SLIP = 20.0; char sys = refPrn.system(); string epoTimeStr = string(_epoTime); - const vector& params = _parlist.params(); + const vector ¶ms = _parlist.params(); unsigned nPar = _parlist.nPar(); for (unsigned ii = 0; ii < LCs.size(); ii++) { - const t_lc::type& tLC = LCs[ii]; + const t_lc::type &tLC = LCs[ii]; if (t_lc::includesPhase(tLC)) { for (unsigned iObs = 0; iObs < obsVector.size(); iObs++) { - const t_pppSatObs* obs = obsVector[iObs]; + const t_pppSatObs *obs = obsVector[iObs]; // Check set Slips and Jump Counters @@ -423,11 +425,13 @@ if (obs->slip()) { - LOG << epoTimeStr << "cycle slip set (obs) " << obs->prn().toString() << endl; + LOG << epoTimeStr << "cycle slip set (obs) " << obs->prn().toString() + << endl; slip = true; } - if (_slips[obs->prn()]._obsSlipCounter != -1 && - _slips[obs->prn()]._obsSlipCounter != obs->slipCounter()) { - LOG << epoTimeStr << "cycle slip set (obsSlipCounter) " << obs->prn().toString() << endl; + if (_slips[obs->prn()]._obsSlipCounter != -1 + && _slips[obs->prn()]._obsSlipCounter != obs->slipCounter()) { + LOG << epoTimeStr << "cycle slip set (obsSlipCounter) " + << obs->prn().toString() << endl; slip = true; } @@ -435,7 +439,8 @@ _slips[obs->prn()]._obsSlipCounter = obs->slipCounter(); } - if (_slips[obs->prn()]._biasJumpCounter != -1 && - _slips[obs->prn()]._biasJumpCounter != obs->biasJumpCounter()) { - LOG << epoTimeStr << "cycle slip set (biasJumpCounter) " << obs->prn().toString() << endl; + if (_slips[obs->prn()]._biasJumpCounter != -1 + && _slips[obs->prn()]._biasJumpCounter != obs->biasJumpCounter()) { + LOG << epoTimeStr << "cycle slip set (biasJumpCounter) " + << obs->prn().toString() << endl; slip = true; } @@ -448,6 +453,5 @@ if (preProcessing) { _obsPool->setRefSatChangeRequired(sys, true); - } - else { + } else { resetAmb(obs->prn(), obsVector); } @@ -456,19 +460,22 @@ // ----------------------- else { - if (refPrn != t_prn()) {return success;} + if (refPrn != t_prn()) { + return success; + } ColumnVector AA(nPar); for (unsigned iPar = 0; iPar < nPar; iPar++) { - const t_pppParam* par = params[iPar]; + const t_pppParam *par = params[iPar]; AA[iPar] = par->partial(_epoTime, obs, tLC, refPrn); } - double ll = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA); + double ll = obs->obsValue(tLC) - obs->cmpValue(tLC) + - DotProduct(_x0, AA); double vv = DotProduct(AA, _xFlt) - ll; if (fabs(vv) > SLIP) { - LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' ' - << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl; + LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) + << ' ' << obs->prn().toString() << ' ' << setw(8) + << setprecision(4) << vv << endl; if (preProcessing) { _obsPool->setRefSatChangeRequired(sys, true); - } - else { + } else { resetAmb(obs->prn(), obsVector); } @@ -483,27 +490,28 @@ // Reset Ambiguity Parameter (cycle slip) //////////////////////////////////////////////////////////////////////////// -t_irc t_pppFilter::resetAmb(t_prn prn, const vector& obsVector, - SymmetricMatrix* QSav, ColumnVector* xSav) { +t_irc t_pppFilter::resetAmb(t_prn prn, const vector &obsVector, + SymmetricMatrix *QSav, ColumnVector *xSav) { t_irc irc = failure; - vector& params = _parlist.params(); + vector ¶ms = _parlist.params(); for (unsigned iPar = 0; iPar < params.size(); iPar++) { - t_pppParam* par = params[iPar]; + t_pppParam *par = params[iPar]; if (par->type() == t_pppParam::amb && par->prn() == prn) { int ind = par->indexNew(); t_lc::type tLC = par->tLC(); LOG << string(_epoTime) << " RESET " << par->toString() << endl; - delete par; par = new t_pppParam(t_pppParam::amb, prn, tLC, &obsVector); + delete par; + par = new t_pppParam(t_pppParam::amb, prn, tLC, &obsVector); par->setIndex(ind); params[iPar] = par; for (unsigned ii = 1; ii <= params.size(); ii++) { - _QFlt(ii, ind+1) = 0.0; + _QFlt(ii, ind + 1) = 0.0; if (QSav) { - (*QSav)(ii, ind+1) = 0.0; - } - } - _QFlt(ind+1,ind+1) = par->sigma0() * par->sigma0(); + (*QSav)(ii, ind + 1) = 0.0; + } + } + _QFlt(ind + 1, ind + 1) = par->sigma0() * par->sigma0(); if (QSav) { - (*QSav)(ind+1,ind+1) = _QFlt(ind+1,ind+1); + (*QSav)(ind + 1, ind + 1) = _QFlt(ind + 1, ind + 1); } _xFlt[ind] = 0.0; @@ -521,15 +529,16 @@ // Add infinite noise to iono //////////////////////////////////////////////////////////////////////////// -t_irc t_pppFilter::addNoiseToIono(char sys) { +t_irc t_pppFilter::addNoiseToPar(t_pppParam::e_type parType, char sys) { t_irc irc = failure; - vector& params = _parlist.params(); + vector ¶ms = _parlist.params(); for (unsigned iPar = 0; iPar < params.size(); iPar++) { - t_pppParam* par = params[iPar]; - if (par->type() == t_pppParam::ion && par->prn().system() == sys) { + t_pppParam *par = params[iPar]; + if (par->type() == parType && par->prn().system() == sys) { int ind = par->indexNew(); - LOG << string(_epoTime) << " ADD IONO_NOISE TO " << par->prn().toString() << endl; + LOG << string(_epoTime) << " ADD IONO_NOISE TO " << par->prn().toString() + << endl; par->setIndex(ind); - _QFlt(ind+1,ind+1) = par->sigma0() * par->sigma0(); + _QFlt(ind + 1, ind + 1) = par->sigma0() * par->sigma0(); irc = success; } @@ -541,6 +550,6 @@ // Compute various DOP Values //////////////////////////////////////////////////////////////////////////// -void t_pppFilter::cmpDOP(const vector& obsVector, - const QMap& refSatMap) { +void t_pppFilter::cmpDOP(const vector &obsVector, + const QMap &refSatMap) { _dop.reset(); @@ -551,5 +560,5 @@ _numSat = 0; for (unsigned ii = 0; ii < obsVector.size(); ii++) { - t_pppSatObs* obs = obsVector[ii]; + t_pppSatObs *obs = obsVector[ii]; char sys = obs->prn().system(); t_prn refPrn = t_prn(); @@ -560,6 +569,6 @@ ++_numSat; for (unsigned iPar = 0; iPar < numPar; iPar++) { - const t_pppParam* par = _parlist.params()[iPar]; - AA[_numSat-1][iPar] = par->partial(_epoTime, obs, t_lc::c1, refPrn); + const t_pppParam *par = _parlist.params()[iPar]; + AA[_numSat - 1][iPar] = par->partial(_epoTime, obs, t_lc::c1, refPrn); } } @@ -569,14 +578,14 @@ } AA = AA.Rows(1, _numSat); - SymmetricMatrix NN; NN << AA.t() * AA; + SymmetricMatrix NN; + NN << AA.t() * AA; SymmetricMatrix QQ = NN.i(); - _dop.H = sqrt(QQ(1,1) + QQ(2,2)); - _dop.V = sqrt(QQ(3,3)); - _dop.P = sqrt(QQ(1,1) + QQ(2,2) + QQ(3,3)); - _dop.T = sqrt(QQ(4,4)); - _dop.G = sqrt(QQ(1,1) + QQ(2,2) + QQ(3,3) + QQ(4,4)); - } - catch (...) { + _dop.H = sqrt(QQ(1, 1) + QQ(2, 2)); + _dop.V = sqrt(QQ(3, 3)); + _dop.P = sqrt(QQ(1, 1) + QQ(2, 2) + QQ(3, 3)); + _dop.T = sqrt(QQ(4, 4)); + _dop.G = sqrt(QQ(1, 1) + QQ(2, 2) + QQ(3, 3) + QQ(4, 4)); + } catch (...) { } } @@ -584,7 +593,7 @@ // //////////////////////////////////////////////////////////////////////////// -void t_pppFilter::predictCovCrdPart(const SymmetricMatrix& QFltOld) { - - const vector& params = _parlist.params(); +void t_pppFilter::predictCovCrdPart(const SymmetricMatrix &QFltOld) { + + const vector ¶ms = _parlist.params(); if (params.size() < 3) { return; @@ -593,29 +602,27 @@ bool first = (params[0]->indexOld() < 0); - SymmetricMatrix Qneu(3); Qneu = 0.0; + SymmetricMatrix Qneu(3); + Qneu = 0.0; for (unsigned ii = 0; ii < 3; ii++) { - const t_pppParam* par = params[ii]; + const t_pppParam *par = params[ii]; if (first) { Qneu[ii][ii] = par->sigma0() * par->sigma0(); - } - else { + } else { Qneu[ii][ii] = par->noise() * par->noise(); } } - const t_pppStation* sta = PPP_CLIENT->staRover(); + const t_pppStation *sta = PPP_CLIENT->staRover(); SymmetricMatrix Qxyz(3); covariNEU_XYZ(Qneu, sta->ellApr().data(), Qxyz); if (first) { - _QFlt.SymSubMatrix(1,3) = Qxyz; - } - else { + _QFlt.SymSubMatrix(1, 3) = Qxyz; + } else { double dt = _epoTime - _firstEpoTime; if (dt < OPT->_seedingTime) { - _QFlt.SymSubMatrix(1,3) = QFltOld.SymSubMatrix(1,3); - } - else { - _QFlt.SymSubMatrix(1,3) = QFltOld.SymSubMatrix(1,3) + Qxyz; + _QFlt.SymSubMatrix(1, 3) = QFltOld.SymSubMatrix(1, 3); + } else { + _QFlt.SymSubMatrix(1, 3) = QFltOld.SymSubMatrix(1, 3) + Qxyz; } } @@ -624,16 +631,19 @@ // //////////////////////////////////////////////////////////////////////////// -void t_pppFilter::setStateVectorAndVarCovMatrix(const ColumnVector& xFltOld, - const SymmetricMatrix& QFltOld) { - - const vector& params = _parlist.params(); +void t_pppFilter::setStateVectorAndVarCovMatrix(const ColumnVector &xFltOld, + const SymmetricMatrix &QFltOld) { + + const vector ¶ms = _parlist.params(); unsigned nPar = params.size(); - _QFlt.ReSize(nPar); _QFlt = 0.0; - _xFlt.ReSize(nPar); _xFlt = 0.0; - _x0.ReSize(nPar); _x0 = 0.0; + _QFlt.ReSize(nPar); + _QFlt = 0.0; + _xFlt.ReSize(nPar); + _xFlt = 0.0; + _x0.ReSize(nPar); + _x0 = 0.0; for (unsigned ii = 0; ii < nPar; ii++) { - t_pppParam* par1 = params[ii]; + t_pppParam *par1 = params[ii]; if (QFltOld.size() == 0) { par1->resetIndex(); @@ -643,13 +653,12 @@ if (iOld < 0) { _QFlt[ii][ii] = par1->sigma0() * par1->sigma0(); // new parameter - } - else { + } else { _QFlt[ii][ii] = QFltOld[iOld][iOld] + par1->noise() * par1->noise(); - _xFlt[ii] = xFltOld[iOld]; + _xFlt[ii] = xFltOld[iOld]; for (unsigned jj = 0; jj < ii; jj++) { - t_pppParam* par2 = params[jj]; - int jOld = par2->indexOld(); + t_pppParam *par2 = params[jj]; + int jOld = par2->indexOld(); if (jOld >= 0) { - _QFlt[ii][jj] = QFltOld(iOld+1,jOld+1); + _QFlt[ii][jj] = QFltOld(iOld + 1, jOld + 1); } } @@ -661,8 +670,9 @@ // Compute datum transformation //////////////////////////////////////////////////////////////////////////// -t_irc t_pppFilter::datumTransformation(const QMap& refSatMap) { +t_irc t_pppFilter::datumTransformation( + const QMap &refSatMap) { // get last epoch - t_pppObsPool::t_epoch* epoch = _obsPool->lastEpoch(); + t_pppObsPool::t_epoch *epoch = _obsPool->lastEpoch(); if (!epoch) { LOG << "t_pppFilter::datumTransformation: !lastEpoch" << endl; @@ -673,7 +683,7 @@ LOG << string(_epoTime) << " DATUM TRANSFORMATION " << endl; - vector& allObs = epoch->obsVector(); - - const QMap& refSatMapLastEpoch = epoch->refSatMap(); + vector &allObs = epoch->obsVector(); + + const QMap &refSatMapLastEpoch = epoch->refSatMap(); bool pseudoObsIono = epoch->pseudoObsIono(); @@ -681,6 +691,8 @@ // reset old and set new refSats in last epoch (ambiguities/GIM) // ============================================================= - if (resetRefSatellitesLastEpoch(allObs, refSatMap, refSatMapLastEpoch) != true) { - LOG << "datumTransformation: resetRefSatellitesLastEpoch != success" << endl; + if (resetRefSatellitesLastEpoch(allObs, refSatMap, refSatMapLastEpoch) + != true) { + LOG << "datumTransformation: resetRefSatellitesLastEpoch != success" + << endl; return failure; } @@ -703,5 +715,5 @@ #endif - const QList& usedSystems = _parlist.usedSystems(); + const QList &usedSystems = _parlist.usedSystems(); for (int iSys = 0; iSys < usedSystems.size(); iSys++) { char sys = usedSystems[iSys]; @@ -720,5 +732,5 @@ unsigned usedLCs = LCs.size(); if (OPT->_pseudoObsIono && !pseudoObsIono) { - usedLCs -= 1; // GIM not used + usedLCs -= 1; // GIM not used } // max Obs @@ -729,8 +741,8 @@ } - const vector& params = _parlist.params(); + const vector ¶ms = _parlist.params(); unsigned nPar = _parlist.nPar(); - Matrix AA(maxObs, nPar); + Matrix AA(maxObs, nPar); // Real Observations @@ -738,11 +750,13 @@ int iObs = -1; for (unsigned ii = 0; ii < obsVector.size(); ii++) { - t_pppSatObs* obs = obsVector[ii]; + t_pppSatObs *obs = obsVector[ii]; for (unsigned jj = 0; jj < usedLCs; jj++) { const t_lc::type tLC = LCs[jj]; - if (tLC == t_lc::GIM) {continue;} + if (tLC == t_lc::GIM) { + continue; + } ++iObs; for (unsigned iPar = 0; iPar < nPar; iPar++) { - const t_pppParam* par = params[iPar]; + const t_pppParam *par = params[iPar]; AA[iObs][iPar] = par->partial(_epoTime, obs, tLC, refPrn); } @@ -753,6 +767,7 @@ continue; } - _datumTrafo->updateIndices(sys, iObs+1); //LOG << "AA Ncols/Nrows: " << AA.Ncols() << "/" << AA.Nrows() << " nPar: " << nPar << endl; //LOG << "AA.SubMatrix(1 .. " << iObs+1 << " , 1 .. " << nPar << ")" << endl; - if(_datumTrafo->prepareAA(AA.SubMatrix(1, iObs+1 , 1, nPar), 2) != success) { + _datumTrafo->updateIndices(sys, iObs + 1); //LOG << "AA Ncols/Nrows: " << AA.Ncols() << "/" << AA.Nrows() << " nPar: " << nPar << endl; //LOG << "AA.SubMatrix(1 .. " << iObs+1 << " , 1 .. " << nPar << ")" << endl; + if (_datumTrafo->prepareAA(AA.SubMatrix(1, iObs + 1, 1, nPar), 2) + != success) { return failure; } @@ -766,5 +781,5 @@ //LOG << "AA2\n"; _datumTrafo->printMatrix(_datumTrafo->AA2(), _datumTrafo->numObs(), _datumTrafo->numPar()); #endif - if(_datumTrafo->computeTrafoMatrix() != success) { + if (_datumTrafo->computeTrafoMatrix() != success) { return failure; } @@ -772,9 +787,9 @@ //LOG << "D21" << endl; _datumTrafo->printMatrix(_datumTrafo->D21(), _datumTrafo->numObs(), _datumTrafo->numPar()); #endif - ColumnVector xFltOld = _xFlt; + ColumnVector xFltOld = _xFlt; SymmetricMatrix QFltOld = _QFlt; _QFlt << _datumTrafo->D21() * QFltOld * _datumTrafo->D21().t(); - _xFlt = _datumTrafo->D21() * xFltOld; + _xFlt = _datumTrafo->D21() * xFltOld; #ifdef BNC_DEBUG_PPP @@ -790,6 +805,23 @@ t_prn refPrnNew = refSatMap[sys]->prn(); if (refPrnNew != refPrnOld) { - t_irc irc = resetAmb(refPrnOld, allObs); - if (OPT->_obsModelType == OPT->DCMcodeBias) {if (irc == success) {addNoiseToIono(sys);}} + if (resetAmb(refPrnOld, allObs) == success) { + if (OPT->_obsModelType == OPT->DCMcodeBias) { + addNoiseToPar(t_pppParam::ion, sys); + } else if (OPT->_obsModelType == OPT->DCMphaseBias) { + if (sys == 'G') { + addNoiseToPar(t_pppParam::pBiasG1, sys); + addNoiseToPar(t_pppParam::pBiasG2, sys); + } else if (sys == 'R') { + addNoiseToPar(t_pppParam::pBiasR1, sys); + addNoiseToPar(t_pppParam::pBiasR2, sys); + } else if (sys == 'E') { + addNoiseToPar(t_pppParam::pBiasE1, sys); + addNoiseToPar(t_pppParam::pBiasE2, sys); + } else if (sys == 'C') { + addNoiseToPar(t_pppParam::pBiasC1, sys); + addNoiseToPar(t_pppParam::pBiasC2, sys); + } + } + } } } @@ -806,8 +838,8 @@ // Init datum transformation //////////////////////////////////////////////////////////////////////////// -void t_pppFilter::initDatumTransformation(const std::vector& allObs, - bool pseudoObsIono) { +void t_pppFilter::initDatumTransformation( + const std::vector &allObs, bool pseudoObsIono) { unsigned trafoObs = 0; - const QList& usedSystems = _parlist. usedSystems(); + const QList &usedSystems = _parlist.usedSystems(); for (int iSys = 0; iSys < usedSystems.size(); iSys++) { char sys = usedSystems[iSys]; @@ -838,10 +870,11 @@ // ////////////////////////////////////////////////////////////////////////////// -bool t_pppFilter::resetRefSatellitesLastEpoch(std::vector& obsVector, - const QMap& refSatMap, - const QMap& refSatMapLastEpoch) { +bool t_pppFilter::resetRefSatellitesLastEpoch( + std::vector &obsVector, + const QMap &refSatMap, + const QMap &refSatMapLastEpoch) { bool resetRefSat; // reference satellite definition per system - const QList& usedSystems = _parlist.usedSystems(); + const QList &usedSystems = _parlist.usedSystems(); for (int iSys = 0; iSys < usedSystems.size(); iSys++) { resetRefSat = false; @@ -856,6 +889,6 @@ vector::iterator it = obsVector.begin(); while (it != obsVector.end()) { - t_pppSatObs* satObs = *it; - if (satObs->prn() == newPrn) { + t_pppSatObs *satObs = *it; + if (satObs->prn() == newPrn) { resetRefSat = true; satObs->setAsReference(); @@ -863,5 +896,5 @@ satObs->resetReference(); } - it++; + it++; } if (!resetRefSat) { Index: /trunk/BNC/src/PPP/pppFilter.h =================================================================== --- /trunk/BNC/src/PPP/pppFilter.h (revision 9641) +++ /trunk/BNC/src/PPP/pppFilter.h (revision 9642) @@ -195,5 +195,5 @@ void predictCovCrdPart(const SymmetricMatrix& QFltOld); - t_irc addNoiseToIono(char sys); + t_irc addNoiseToPar(t_pppParam::e_type parType, char sys); bool resetRefSatellitesLastEpoch(std::vector& obsVector, Index: /trunk/BNC/src/PPP/pppParlist.cpp =================================================================== --- /trunk/BNC/src/PPP/pppParlist.cpp (revision 9641) +++ /trunk/BNC/src/PPP/pppParlist.cpp (revision 9642) @@ -102,15 +102,27 @@ _noise = OPT->_noiseIon; break; - case cBiasG1: case cBiasR1: case cBiasE1: case cBiasC1: - case cBiasG2: case cBiasR2: case cBiasE2: case cBiasC2: + case cBiasG1: case cBiasE1: case cBiasC1: + case cBiasG2: case cBiasE2: case cBiasC2: _epoSpec = false; _sigma0 = OPT->_aprSigCodeBias; _noise = OPT->_noiseCodeBias; break; - case pBiasG1: case pBiasR1: case pBiasE1: case pBiasC1: - case pBiasG2: case pBiasR2: case pBiasE2: case pBiasC2: + case cBiasR1: + case cBiasR2: + _epoSpec = false; + _sigma0 = OPT->_aprSigCodeBias; + _noise = 10.0 * OPT->_noiseCodeBias; + break; + case pBiasG1: case pBiasE1: case pBiasC1: + case pBiasG2: case pBiasE2: case pBiasC2: _epoSpec = false; _sigma0 = OPT->_aprSigPhaseBias; _noise = OPT->_noisePhaseBias; + break; + case pBiasR1: + case pBiasR2: + _epoSpec = false; + _sigma0 = OPT->_aprSigPhaseBias; + _noise = 10.0 * OPT->_noisePhaseBias; break; }