Index: trunk/BNC/src/PPP/pppClient.cpp =================================================================== --- trunk/BNC/src/PPP/pppClient.cpp (revision 8904) +++ trunk/BNC/src/PPP/pppClient.cpp (revision 8905) @@ -55,21 +55,23 @@ _obsPool = new t_pppObsPool(); _staRover = new t_pppStation(); - _filter = new t_pppFilter(); + _filter = new t_pppFilter(_obsPool); _tides = new t_tides(); - + _antex = 0; if (!_opt->_antexFileName.empty()) { _antex = new bncAntex(_opt->_antexFileName.c_str()); } - else { - _antex = 0; - } - if (!_opt->_blqFileName.empty()) { - _loading = new t_loading(_opt->_blqFileName.c_str()); - } - else { - _loading = 0; - } - + if (_tides->readBlqFile(_opt->_blqFileName.c_str()) == success) { + //_tides->printAllBlqSets(); + } + } + + if (OPT->_refSatRequired) { + for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) { + char system = OPT->systems()[iSys]; + _obsPool->initRefSatMapElement(system); + } + } + _offGG = 0.0; CLIENTS.setLocalData(this); // CLIENTS takes ownership over "this" } @@ -80,4 +82,5 @@ delete _log; delete _opt; + delete _filter; delete _ephPool; delete _obsPool; @@ -86,9 +89,5 @@ delete _antex; } - delete _filter; delete _tides; - if (_loading) { - delete _loading; - } clearObs(); } @@ -157,4 +156,5 @@ t_irc t_pppClient::prepareObs(const vector& satObs, vector& obsVector, bncTime& epoTime) { + // Default // ------- @@ -176,4 +176,48 @@ } + // (re)set reference satellites per system if required + // --------------------------------------------------- + if (_opt->_refSatRequired) { + // reference satellite definition per system + for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) { + char system = OPT->systems()[iSys]; + bool refSatDefined = false; + t_pppRefSat* refSat = _obsPool->getRefSatMapElement(system); + for (unsigned ii = 0; ii < obsVector.size(); ii++) { + const t_pppSatObs* satObs = obsVector.at(ii); + // reference satellite is unchanged + if (!_obsPool->epoReProcessing() && refSat->prn() == satObs->prn()) { + refSatDefined = true; + obsVector[ii]->setAsReference(); + refSat->setStatus(t_pppRefSat::unchanged); + } + // reference satellite has changed + else if (_obsPool->epoReProcessing() && refSat->prn() != satObs->prn()) { + if (satObs->prn().system() == system) { + refSatDefined = true; + obsVector[ii]->setAsReference(); + refSat->setStatus(t_pppRefSat::changed); + refSat->setPrn(satObs->prn()); + } + } + if (refSatDefined) { + break; + } + } + // reference satellite has to be initialized + if (!refSatDefined) { + for (unsigned ii = 0; ii < obsVector.size(); ii++) { + const t_pppSatObs* satObs = obsVector.at(ii); + if (satObs->prn().system() == system) { + obsVector[ii]->setAsReference(); + refSat->setStatus(t_pppRefSat::initialized); + refSat->setPrn(satObs->prn()); + } + } + } + } + _obsPool->setEpoReProcessing(false); //TODO: später erst nach Trafo false setzen + } + // Check whether data are synchronized, compute epoTime // ---------------------------------------------------- @@ -202,4 +246,29 @@ } +// +////////////////////////////////////////////////////////////////////////////// +bool t_pppClient::preparePseudoObs(std::vector& obsVector) { + + bool pseudoObsIono = false; + + if (OPT->_pseudoObsIono) { + vector::iterator it = obsVector.begin(); + while (it != obsVector.end()) { + t_pppSatObs* satObs = *it; + char system = satObs->prn().system(); + t_pppRefSat* refSat = _obsPool->getRefSatMapElement(system); + double stecRef = refSat->stecValue(); + if (stecRef && !satObs->isReference()) { + pseudoObsIono = true; + satObs->setPseudoObsIono(t_frequency::G1, stecRef); + } + satObs->printObsMinusComputed(); + it++; + } + } + + return pseudoObsIono; +} + // Compute the Bancroft position, check for blunders ////////////////////////////////////////////////////////////////////////////// @@ -250,8 +319,9 @@ for (unsigned ii = 0; ii < obsVector.size(); ii++) { const t_pppSatObs* satObs = obsVector.at(ii); - if ( satObs->isValid() && satObs->prn().system() == 'G' && - (!satObs->modelSet() || satObs->eleSat() >= OPT->_minEle) ) { + if (satObs->isValid() && + satObs->prn().system() == 'G' && + (!satObs->modelSet() || satObs->eleSat() >= OPT->_minEle) ) { ColumnVector rr = satObs->xc().Rows(1,3) - xyzc.Rows(1,3); - double res = rr.norm_Frobenius() - satObs->obsValue(tLC) + double res = rr.NormFrobenius() - satObs->obsValue(tLC) - (satObs->xc()[3] - xyzc[3]) * t_CST::c; if (std::isnan(res) || fabs(res) > maxRes) { @@ -393,9 +463,9 @@ t_irc t_pppClient::cmpModel(t_pppStation* station, const ColumnVector& xyzc, vector& obsVector) { - bncTime time; time = _epoTimeRover; station->setName(OPT->_roverName); station->setAntName(OPT->_antNameRover); + station->setEpochTime(time); if (OPT->xyzAprRoverSet()) { station->setXyzApr(OPT->_xyzAprRover); @@ -412,9 +482,6 @@ // Tides // ----- - station->setTideDspl( _tides->displacement(time, station->xyzApr()) ); - - // Ionosphere - // ---------- - station->setIonoEpochTime(time); + station->setTideDsplEarth(_tides->earth(time, station->xyzApr())); + station->setTideDsplOcean(_tides->ocean(time, station->xyzApr(), station->name())); // Observation model @@ -430,4 +497,9 @@ satObs->eleSat() >= OPT->_minEle && modelSetup == success) { + if (satObs->isReference() && OPT->_pseudoObsIono) { + char system = satObs->prn().system(); + t_pppRefSat* refSat = _obsPool->getRefSatMapElement(system); + refSat->setStecValue(satObs->getIonoCodeDelay(t_frequency::G1)); + } ++it; } @@ -447,42 +519,65 @@ try { initOutput(output); - - // Prepare Observations of the Rover - // --------------------------------- - if (prepareObs(satObs, _obsRover, _epoTimeRover) != success) { - return finish(failure); - } - - LOG << "\nPPP of Epoch "; - if (!_epoTimeRover.undef()) LOG << string(_epoTimeRover); - LOG << "\n---------------------------------------------------------------\n"; - - for (int iter = 1; iter <= 2; iter++) { - ColumnVector xyzc(4); xyzc = 0.0; - bool print = (iter == 2); - if (cmpBancroft(_epoTimeRover, _obsRover, xyzc, print) != success) { + _num = 0; + _obsPool->setEpoReProcessing(false); // initialize for epoch + + do { + _num++; + + // Prepare Observations of the Rover + // --------------------------------- + if (prepareObs(satObs, _obsRover, _epoTimeRover) != success) { return finish(failure); } - if (cmpModel(_staRover, xyzc, _obsRover) != success) { + + LOG << "\nPPP of Epoch "; + if (!_epoTimeRover.undef()) LOG << string(_epoTimeRover); + LOG << "\n---------------------------------------------------------------\n"; + + for (int iter = 1; iter <= 2; iter++) { + ColumnVector xyzc(4); xyzc = 0.0; + bool print = (iter == 2); + if (cmpBancroft(_epoTimeRover, _obsRover, xyzc, print) != success) { + return finish(failure); + } + if (cmpModel(_staRover, xyzc, _obsRover) != success) { + return finish(failure); + } + } + + _offGG = cmpOffGG(_obsRover); + + // Prepare Pseudo Observations of the Rover + // ---------------------------------------- + _pseudoObsIono = preparePseudoObs(_obsRover);//qDebug() << "_pseudoObsIonoAvailable: " << _pseudoObsIono; + + if (int(_obsRover.size()) < OPT->_minObs) { + LOG << "t_pppClient::processEpoch not enough observations" << endl; return finish(failure); } - } - - _offGG = cmpOffGG(_obsRover); - - if (int(_obsRover.size()) < OPT->_minObs) { - LOG << "t_pppClient::processEpoch not enough observations" << endl; - return finish(failure); - } - - // Store last epoch of data - // ------------------------ - _obsPool->putEpoch(_epoTimeRover, _obsRover); - - // Process Epoch in Filter - // ----------------------- - if (_filter->processEpoch(_obsPool) != success) { - return finish(failure); - } + + if (OPT->_refSatRequired) { + LOG.setf(ios::fixed); + QMapIterator it(_obsPool->getRefSatMap()); + while (it.hasNext()) { + it.next(); + char sys = it.key(); + string prn = it.value()->prn().toString(); + LOG << string(_epoTimeRover) << " REFSAT " << sys << ": " << prn << endl; + } + } + + // Store last epoch of data + // ------------------------ + //_obsRover.resize(2); + _obsPool->putEpoch(_epoTimeRover, _obsRover, _pseudoObsIono); + + // Process Epoch in Filter + // ----------------------- + if (_filter->processEpoch(_num) != success) { + return finish(failure); + } + // if num > 1 und !obsPool->epoReProcessing() => filter ->datumTransformation + } while (_obsPool->epoReProcessing()); } catch (Exception& exc) { @@ -588,8 +683,4 @@ void t_pppClient::reset() { - // to delete all parameters - delete _filter; - _filter = new t_pppFilter(); - // to delete old orbit and clock corrections delete _ephPool; @@ -599,3 +690,8 @@ delete _obsPool; _obsPool = new t_pppObsPool(); -} + + // to delete all parameters + delete _filter; + _filter = new t_pppFilter(_obsPool); + +} Index: trunk/BNC/src/PPP/pppClient.h =================================================================== --- trunk/BNC/src/PPP/pppClient.h (revision 8904) +++ trunk/BNC/src/PPP/pppClient.h (revision 8905) @@ -10,4 +10,5 @@ class bncAntex; +class t_pppRefSat; namespace BNC_PPP { @@ -34,7 +35,7 @@ const t_pppEphPool* ephPool() const {return _ephPool;} const t_pppObsPool* obsPool() const {return _obsPool;} - const bncAntex* antex() const {return _antex;} + const bncAntex* antex() const {return _antex;} const t_pppStation* staRover() const {return _staRover;} - double offGG() const {return _offGG;} + double offGG() const {return _offGG;} std::ostringstream& log() {return *_log;} @@ -53,4 +54,5 @@ t_irc prepareObs(const std::vector& satObs, std::vector& obsVector, bncTime& epoTime); + bool preparePseudoObs(std::vector& obsVector); t_irc cmpModel(t_pppStation* station, const ColumnVector& xyzc, std::vector& obsVector); @@ -71,5 +73,6 @@ t_pppOptions* _opt; t_tides* _tides; - t_loading* _loading; + bool _pseudoObsIono; + int _num; }; Index: trunk/BNC/src/PPP/pppFilter.cpp =================================================================== --- trunk/BNC/src/PPP/pppFilter.cpp (revision 8904) +++ trunk/BNC/src/PPP/pppFilter.cpp (revision 8905) @@ -34,6 +34,8 @@ // Constructor //////////////////////////////////////////////////////////////////////////// -t_pppFilter::t_pppFilter() { +t_pppFilter::t_pppFilter(t_pppObsPool* obsPool) { _parlist = 0; + _numSat = 0; + _obsPool = obsPool; } @@ -46,6 +48,5 @@ // Process Single Epoch //////////////////////////////////////////////////////////////////////////// -t_irc t_pppFilter::processEpoch(t_pppObsPool* obsPool) { - +t_irc t_pppFilter::processEpoch(int num) { _numSat = 0; const double maxSolGap = 60.0; @@ -57,5 +58,5 @@ // Vector of all Observations // -------------------------- - t_pppObsPool::t_epoch* epoch = obsPool->lastEpoch(); + t_pppObsPool::t_epoch* epoch = _obsPool->lastEpoch(); if (!epoch) { return failure; @@ -75,4 +76,5 @@ string epoTimeStr = string(_epoTime); + //-- // Set Parameters // -------------- @@ -91,7 +93,5 @@ for (unsigned ii = 0; ii < params.size(); ii++) { const t_pppParam* par1 = params[ii]; - _x0[ii] = par1->x0(); - int iOld = par1->indexOld(); if (iOld < 0) { @@ -113,8 +113,32 @@ predictCovCrdPart(QFltOld); + + // Pre-Process Satellite Systems separately + // ---------------------------------------- + bool preProcessing = false; + if (OPT->_obsModelType == OPT->DCMcodeBias || + OPT->_obsModelType == OPT->DCMphaseBias) { + preProcessing = true; + for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) { + char system = OPT->systems()[iSys]; + t_prn refPrn = (_obsPool->getRefSatMapElement(system))->prn(); + vector obsVector; + for (unsigned jj = 0; jj < allObs.size(); jj++) { + if (allObs[jj]->prn().system() == system) { + obsVector.push_back(allObs[jj]); + } + } + if (processSystem(OPT->LCs(system), obsVector, refPrn, + epoch->pseudoObsIono(), preProcessing) != success) { + return failure; + } + } + } + // Process Satellite Systems separately // ------------------------------------ for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) { char system = OPT->systems()[iSys]; + t_prn refPrn = (_obsPool->getRefSatMapElement(system))->prn(); unsigned int num = 0; vector obsVector; @@ -126,13 +150,20 @@ } LOG << epoTimeStr << " SATNUM " << system << ' ' << right << setw(2) << num << endl; - if ( processSystem(OPT->LCs(system), obsVector) != success ) { + if (processSystem(OPT->LCs(system), obsVector, refPrn, + epoch->pseudoObsIono(), preProcessing) != success) { return failure; } } - - cmpDOP(allObs); - - _parlist->printResult(_epoTime, _QFlt, _xFlt); - _lastEpoTimeOK = _epoTime; // remember time of last successful epoch processing + if (_obsPool->epoReProcessing()) { + // set A1 und A2 abhängig von num + // if num == 1 => A1 + // if num > 1 => A2 + } + else { + cmpDOP(allObs); + _parlist->printResult(_epoTime, _QFlt, _xFlt); + _lastEpoTimeOK = _epoTime; // remember time of last successful epoch processing + } + return success; } @@ -141,23 +172,38 @@ //////////////////////////////////////////////////////////////////////////// t_irc t_pppFilter::processSystem(const vector& LCs, - const vector& obsVector) { - + const vector& obsVector, + const t_prn& refPrn, + bool pseudoObsIonoAvailable, + bool preProcessing) { + qDebug() << "======t_pppFilter::processSystem======="; LOG.setf(ios::fixed); // Detect Cycle Slips // ------------------ - if (detectCycleSlips(LCs, obsVector) != success) { + if (detectCycleSlips(LCs, obsVector, refPrn, preProcessing) != success) { return failure; } - ColumnVector xSav = _xFlt; - SymmetricMatrix QSav = _QFlt; - string epoTimeStr = string(_epoTime); + unsigned usedLCs = LCs.size(); //qDebug() << "usedLCs: " << usedLCs; + if (OPT->_pseudoObsIono && !pseudoObsIonoAvailable) { + usedLCs -= 1; // GIM not used + } + ColumnVector xSav = _xFlt; + SymmetricMatrix QSav = _QFlt; + string epoTimeStr = string(_epoTime); const vector& params = _parlist->params(); - unsigned maxObs = obsVector.size() * LCs.size(); + unsigned maxObs = obsVector.size() * usedLCs; + //unsigned maxObs = 2 * usedLCs; + + if (OPT->_pseudoObsIono && pseudoObsIonoAvailable) { // stecDiff w.r.t refSat + maxObs -= 1; + } + qDebug() << "par.size() : " << _parlist->nPar() << " LCs.size(): " << usedLCs; + qDebug() << "obsVector.size(): " << obsVector.size() << " maxObs : " << maxObs; // Outlier Detection Loop // ---------------------- for (unsigned iOutlier = 0; iOutlier < maxObs; iOutlier++) { + qDebug() << "iOutlier: " << iOutlier; if (iOutlier > 0) { @@ -168,7 +214,14 @@ // First-Design Matrix, Terms Observed-Computed, Weight Matrix // ----------------------------------------------------------- + qDebug() << "A(" << maxObs << "," << _parlist->nPar() << ")"; Matrix AA(maxObs, _parlist->nPar()); ColumnVector ll(maxObs); DiagonalMatrix PP(maxObs); PP = 0.0; + //TETSPLOT + for (unsigned iPar = 0; iPar < params.size(); iPar++) { + const t_pppParam* par = params[iPar]; + cout << " " << par->toString() < usedTypes; for (unsigned ii = 0; ii < obsVector.size(); ii++) { - t_pppSatObs* obs = obsVector[ii]; + t_pppSatObs* obs = obsVector[ii]; qDebug() << "SATELLITE: " << obs->prn().toString().c_str() << "isRef: " << obs->isReference(); if (!obs->outlier()) { - for (unsigned jj = 0; jj < LCs.size(); jj++) { + for (unsigned jj = 0; jj < usedLCs; jj++) { const t_lc::type tLC = LCs[jj]; + if (tLC == t_lc::GIM && obs->isReference()) {continue;} ++iObs; usedObs.push_back(obs); @@ -186,7 +240,11 @@ const t_pppParam* par = params[iPar]; AA[iObs][iPar] = par->partial(_epoTime, obs, tLC); + cout << setw(10) << par->partial(_epoTime, obs, tLC); } + cout << endl; ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1)); PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC)); + //qDebug() << "ll[iObs]: " << ll[iObs]; + //qDebug() << "PP[iObs]: " << PP[iObs]; } } @@ -229,15 +287,23 @@ t_pppSatObs* obs = usedObs[maxOutlierIndex]; t_pppParam* par = 0; - LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' ' - << obs->prn().toString() << ' ' - << setw(8) << setprecision(4) << maxOutlier << endl; + if (!preProcessing) { + LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' ' + << obs->prn().toString() << ' ' + << setw(8) << setprecision(4) << maxOutlier << endl; + } for (unsigned iPar = 0; iPar < params.size(); iPar++) { t_pppParam* hlp = params[iPar]; - if (hlp->type() == t_pppParam::amb && hlp->prn() == obs->prn() && + if (hlp->type() == t_pppParam::amb && + hlp->prn() == obs->prn() && hlp->tLC() == usedTypes[maxOutlierIndex]) { par = hlp; } } - if (par) { + if (par && preProcessing) { + if (par->prn() == refPrn) { + _obsPool->setEpoReProcessing(true); + } + } + else if (par && !preProcessing) { if (par->ambResetCandidate()) { resetAmb(par->prn(), obsVector, &QSav, &xSav); @@ -248,5 +314,5 @@ } } - else { + else if (!par && !preProcessing){ obs->setOutlier(); } @@ -256,14 +322,16 @@ // --------------- else { - for (unsigned jj = 0; jj < LCs.size(); jj++) { - for (unsigned ii = 0; ii < usedObs.size(); ii++) { - const t_lc::type tLC = usedTypes[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 << ' ' - << obs->prn().toString().substr(0,3) << ' ' - << setw(8) << setprecision(4) << vv[ii] << endl; + if (!preProcessing) { + for (unsigned jj = 0; jj < LCs.size(); jj++) { + for (unsigned ii = 0; ii < usedObs.size(); ii++) { + const t_lc::type tLC = usedTypes[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 << ' ' + << obs->prn().toString().substr(0,3) << ' ' + << setw(8) << setprecision(4) << vv[ii] << endl; + } } } @@ -279,9 +347,11 @@ //////////////////////////////////////////////////////////////////////////// t_irc t_pppFilter::detectCycleSlips(const vector& LCs, - const vector& obsVector) { + const vector& obsVector, + const t_prn& refPrn, + bool preProcessing) { const double SLIP = 20.0; // slip threshold string epoTimeStr = string(_epoTime); - const vector& params = _parlist->params(); + const vector& params = _parlist->params(); for (unsigned ii = 0; ii < LCs.size(); ii++) { @@ -294,5 +364,4 @@ // --------------------------------- bool slip = false; - if (obs->slip()) { LOG << epoTimeStr << "cycle slip set (obs) " << obs->prn().toString() << endl; @@ -312,12 +381,17 @@ slip = true; } - _slips[obs->prn()]._biasJumpCounter = obs->biasJumpCounter(); // Slip Set // -------- if (slip) { - resetAmb(obs->prn(), obsVector); - } - + if (preProcessing) { + if (obs->prn() == refPrn) { + _obsPool->setEpoReProcessing(true); + } + } + else { + resetAmb(obs->prn(), obsVector); + } + } // Check Pre-Fit Residuals // ----------------------- @@ -328,12 +402,17 @@ AA[iPar] = par->partial(_epoTime, obs, tLC); } - 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; - resetAmb(obs->prn(), obsVector); + if (preProcessing) { + if (obs->prn() == refPrn) { + _obsPool->setEpoReProcessing(true); + } + } + else { + LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' ' + << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl; + resetAmb(obs->prn(), obsVector); + } } } Index: trunk/BNC/src/PPP/pppFilter.h =================================================================== --- trunk/BNC/src/PPP/pppFilter.h (revision 8904) +++ trunk/BNC/src/PPP/pppFilter.h (revision 8905) @@ -17,8 +17,8 @@ class t_pppFilter { public: - t_pppFilter(); + t_pppFilter(t_pppObsPool* obsPool); ~t_pppFilter(); - t_irc processEpoch(t_pppObsPool* obsPool); + t_irc processEpoch(int num); const ColumnVector& x() const {return _xFlt;} @@ -74,8 +74,13 @@ t_irc processSystem(const std::vector& LCs, - const std::vector& obsVector); + const std::vector& obsVector, + const t_prn& refPrn, + bool pseudoObsIonoAvailable, + bool preProcessing); t_irc detectCycleSlips(const std::vector& LCs, - const std::vector& obsVector); + const std::vector& obsVector, + const t_prn& refPrn, + bool preProcessing); t_irc resetAmb(t_prn prn, const std::vector& obsVector, @@ -88,4 +93,5 @@ bncTime _epoTime; t_pppParlist* _parlist; + t_pppObsPool* _obsPool; SymmetricMatrix _QFlt; ColumnVector _xFlt; Index: trunk/BNC/src/PPP/pppObsPool.cpp =================================================================== --- trunk/BNC/src/PPP/pppObsPool.cpp (revision 8904) +++ trunk/BNC/src/PPP/pppObsPool.cpp (revision 8905) @@ -22,6 +22,8 @@ // Constructor ///////////////////////////////////////////////////////////////////////////// -t_pppObsPool::t_epoch::t_epoch(const bncTime& epoTime, vector& obsVector) { - _epoTime = epoTime; +t_pppObsPool::t_epoch::t_epoch(const bncTime& epoTime, vector& obsVector, + bool pseudoObsIono) { + _epoTime = epoTime; + _pseudoObsIono = pseudoObsIono; for (unsigned ii = 0; ii < obsVector.size(); ii++) { _obsVector.push_back(obsVector[ii]); @@ -48,4 +50,5 @@ } _vTec = 0; + _epoReProcessing = false; } @@ -64,4 +67,5 @@ _epochs.pop_front(); } + clearRefSatMap(); } @@ -92,7 +96,13 @@ // ///////////////////////////////////////////////////////////////////////////// -void t_pppObsPool::putEpoch(const bncTime& epoTime, vector& obsVector) { +void t_pppObsPool::putEpoch(const bncTime& epoTime, vector& obsVector, + bool pseudoObsIono) { const unsigned MAXSIZE = 2; - _epochs.push_back(new t_epoch(epoTime, obsVector)); + + if (epoReProcessing()) { + delete _epochs.back(); + _epochs.pop_back(); + } + _epochs.push_back(new t_epoch(epoTime, obsVector,pseudoObsIono)); if (_epochs.size() > MAXSIZE) { delete _epochs.front(); Index: trunk/BNC/src/PPP/pppObsPool.h =================================================================== --- trunk/BNC/src/PPP/pppObsPool.h (revision 8904) +++ trunk/BNC/src/PPP/pppObsPool.h (revision 8905) @@ -4,6 +4,8 @@ #include #include +#include #include "pppSatObs.h" #include "bnctime.h" +#include "pppRefSat.h" namespace BNC_PPP { @@ -14,12 +16,15 @@ class t_epoch { public: - t_epoch(const bncTime& epoTime, std::vector& obsVector); + t_epoch(const bncTime& epoTime, std::vector& obsVector, + bool pseudoObsIono); ~t_epoch(); - std::vector& obsVector() {return _obsVector;} + std::vector& obsVector() {return _obsVector;} const std::vector& obsVector() const {return _obsVector;} const bncTime& epoTime() const {return _epoTime;} + bool pseudoObsIono() const {return _pseudoObsIono;} private: - bncTime _epoTime; + bncTime _epoTime; std::vector _obsVector; + bool _pseudoObsIono; }; @@ -30,5 +35,6 @@ void putTec(t_vTec* _vTec); - void putEpoch(const bncTime& epoTime, std::vector& obsVector); + void putEpoch(const bncTime& epoTime, std::vector& obsVector, + bool pseudoObs); const t_satCodeBias* satCodeBias(const t_prn& prn) const { @@ -49,9 +55,32 @@ } + void initRefSatMapElement(char system) { + _refSatMap[system] = new t_pppRefSat(); + } + void clearRefSatMap() { + QMapIterator it(_refSatMap); + while (it.hasNext()) { + it.next(); + delete it.value(); + } + _refSatMap.clear(); + } + t_pppRefSat* getRefSatMapElement(char system) { + return _refSatMap[system]; + } + QMap getRefSatMap() {return _refSatMap;} + + void setEpoReProcessing(bool epoReProcessing) { + _epoReProcessing = epoReProcessing; + } + bool epoReProcessing() {return _epoReProcessing;} + private: - t_satCodeBias* _satCodeBiases[t_prn::MAXPRN+1]; - t_satPhaseBias* _satPhaseBiases[t_prn::MAXPRN+1]; - t_vTec* _vTec; - std::deque _epochs; + t_satCodeBias* _satCodeBiases[t_prn::MAXPRN+1]; + t_satPhaseBias* _satPhaseBiases[t_prn::MAXPRN+1]; + t_vTec* _vTec; + std::deque _epochs; + QMap _refSatMap; + bool _epoReProcessing; }; Index: trunk/BNC/src/PPP/pppParlist.cpp =================================================================== --- trunk/BNC/src/PPP/pppParlist.cpp (revision 8904) +++ trunk/BNC/src/PPP/pppParlist.cpp (revision 8905) @@ -94,4 +94,21 @@ _noise = OPT->_noiseTrp; break; + case ion: + _epoSpec = false; + _sigma0 = OPT->_aprSigIon; + _noise = OPT->_noiseIon; + break; + case cBias1: + case cBias2: + _epoSpec = false; + _sigma0 = OPT->_aprSigCodeBias; + _noise = OPT->_noiseCodeBias; + break; + case pBias1: + case pBias2: + _epoSpec = false; + _sigma0 = OPT->_aprSigPhaseBias; + _noise = OPT->_noisePhaseBias; + break; } } @@ -106,5 +123,5 @@ //////////////////////////////////////////////////////////////////////////// double t_pppParam::partial(const bncTime& /* epoTime */, const t_pppSatObs* obs, - const t_lc::type& tLC) const { + const t_lc::type& tLC) const {//qDebug() << "t_pppParam::partial: " << tLC; // Special Case - Melbourne-Wuebbena @@ -115,42 +132,110 @@ const t_pppStation* sta = PPP_CLIENT->staRover(); - ColumnVector rhoV = sta->xyzApr() - obs->xc().Rows(1,3); + ColumnVector rhoV = sta->xyzApr() - obs->xc().Rows(1,3); + + map codeCoeff; + map phaseCoeff; + map ionoCoeff; + obs->lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff); switch (_type) { case crdX: - return (sta->xyzApr()[0] - obs->xc()[0]) / rhoV.norm_Frobenius(); + if (tLC == t_lc::GIM) {return 0.0;} + return (sta->xyzApr()[0] - obs->xc()[0]) / rhoV.NormFrobenius(); case crdY: - return (sta->xyzApr()[1] - obs->xc()[1]) / rhoV.norm_Frobenius(); + if (tLC == t_lc::GIM) {return 0.0;} + return (sta->xyzApr()[1] - obs->xc()[1]) / rhoV.NormFrobenius(); case crdZ: - return (sta->xyzApr()[2] - obs->xc()[2]) / rhoV.norm_Frobenius(); + if (tLC == t_lc::GIM) {return 0.0;} + return (sta->xyzApr()[2] - obs->xc()[2]) / rhoV.NormFrobenius(); case clkR: + if (tLC == t_lc::GIM) {return 0.0;} return 1.0; case offGG: + if (tLC == t_lc::GIM) {return 0.0;} return (obs->prn().system() == 'R') ? 1.0 : 0.0; case amb: + if (tLC == t_lc::GIM) {return 0.0;} + else if ((OPT->_obsModelType == OPT->IF) || + (OPT->_obsModelType == OPT->PPPRTK) || + (OPT->_obsModelType == OPT->UncombPPP) || + (OPT->_obsModelType == OPT->DCMcodeBias && !obs->isReference()) || + (OPT->_obsModelType == OPT->DCMphaseBias && !obs->isReference()) ) { + if (obs->prn() == _prn) { + if (tLC == _tLC) { + return (obs->lambda(tLC)); + } + else if (tLC == t_lc::lIF && _tLC == t_lc::MW) { + return obs->lambda(t_lc::lIF) * obs->lambda(t_lc::MW) / obs->lambda(t_lc::l2); + } + else { + if (_tLC == t_lc::l1) { + return obs->lambda(t_lc::l1) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)]; + } + else if (_tLC == t_lc::l2) { + return obs->lambda(t_lc::l2) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)]; + } + } + } + } + break; + case trp: + if (tLC == t_lc::GIM) {return 0.0;} + return 1.0 / sin(obs->eleSat()); + case ion: +// qDebug() << "refPrn: " << _refPrn.toString().c_str(); if (obs->prn() == _prn) { - if (tLC == _tLC) { - return (obs->lambda(tLC)); - } - else if (tLC == t_lc::lIF && _tLC == t_lc::MW) { - return obs->lambda(t_lc::lIF) * obs->lambda(t_lc::MW) / obs->lambda(t_lc::l2); - } - else { - map codeCoeff; - map phaseCoeff; - obs->lcCoeff(tLC, codeCoeff, phaseCoeff); - if (_tLC == t_lc::l1) { - return obs->lambda(t_lc::l1) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)]; - } - else if (_tLC == t_lc::l2) { - return obs->lambda(t_lc::l2) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)]; - } - } - } - return 0.0; - case trp: - return 1.0 / sin(obs->eleSat()); - } - + if (tLC == t_lc::c1) { + return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::c1)]; + } + else if (tLC == t_lc::c2) { + return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::c2)]; + } + else if (tLC == t_lc::l1) { + return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)]; + } + else if (tLC == t_lc::l2) { + return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)]; + } + else if (tLC == t_lc::GIM) { + return -1.0; + } + } + if (tLC == t_lc::GIM && _prn == _refPrn) { + return 1.0; + } + break; + case cBias1: + if (tLC == t_lc::c1) { + return 1.0; + } + else { + return 0.0; + } + break; + case cBias2: + if (tLC == t_lc::c2) { + return 1.0; + } + else { + return 0.0; + } + break; + case pBias1: + if (tLC == t_lc::l1) { + return 1.0; + } + else { + return 0.0; + } + break; + case pBias2: + if (tLC == t_lc::l2) { + return 1.0; + } + else { + return 0.0; + } + } return 0.0; } @@ -171,8 +256,5 @@ break; case clkR: - ss << "CLK "; - break; - case amb: - ss << "AMB " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString(); + ss << "REC_CLK "; break; case offGG: @@ -181,4 +263,16 @@ case trp: ss << "TRP "; + break; + case amb: + ss << "AMB " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString(); + break; + case ion: + ss << "ION " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString(); + break; + case cBias1: + case cBias2: + case pBias1: + case pBias2: + ss << "BIAS " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << "REC"; break; } @@ -201,5 +295,6 @@ // //////////////////////////////////////////////////////////////////////////// -t_irc t_pppParlist::set(const bncTime& epoTime, const std::vector& obsVector) { +t_irc t_pppParlist::set(const bncTime& epoTime, + const std::vector& obsVector) { // Remove some Parameters @@ -271,5 +366,5 @@ required.push_back(new t_pppParam(t_pppParam::clkR, t_prn(), t_lc::dummy)); - // GPS-Glonass Clock Offset + // GPS-GLONASS Clock Offset // ------------------------ if (OPT->useSystem('R')) { @@ -283,12 +378,44 @@ } + // Ionosphere + // ---------- + if (OPT->_obsModelType == OPT->UncombPPP || + OPT->_obsModelType == OPT->DCMcodeBias || + OPT->_obsModelType == OPT->DCMphaseBias ) { + for (unsigned jj = 0; jj < obsVector.size(); jj++) { + const t_pppSatObs* satObs = obsVector[jj]; + required.push_back(new t_pppParam(t_pppParam::ion, satObs->prn(), t_lc::dummy)); + } + } + // Ambiguities // ----------- for (unsigned jj = 0; jj < obsVector.size(); jj++) { const t_pppSatObs* satObs = obsVector[jj]; - const vector& ambLCs = OPT->ambLCs(satObs->prn().system()); - for (unsigned ii = 0; ii < ambLCs.size(); ii++) { - required.push_back(new t_pppParam(t_pppParam::amb, satObs->prn(), ambLCs[ii], &obsVector)); - } + if ((OPT->_obsModelType == OPT->IF) || + (OPT->_obsModelType == OPT->PPPRTK) || + (OPT->_obsModelType == OPT->UncombPPP) || + (OPT->_obsModelType == OPT->DCMcodeBias && !satObs->isReference()) || + (OPT->_obsModelType == OPT->DCMphaseBias && !satObs->isReference()) ) { + const vector& ambLCs = OPT->ambLCs(satObs->prn().system()); + for (unsigned ii = 0; ii < ambLCs.size(); ii++) { + required.push_back(new t_pppParam(t_pppParam::amb, satObs->prn(), ambLCs[ii], &obsVector)); + } + } + } + + // Receiver Code Biases + // -------------------- + if (OPT->_obsModelType == OPT->DCMcodeBias) { + required.push_back(new t_pppParam(t_pppParam::cBias1, t_prn(), t_lc::c1)); + required.push_back(new t_pppParam(t_pppParam::cBias2, t_prn(), t_lc::c2)); + } + + // Receiver Phase Biases + // --------------------- + if ((OPT->_obsModelType == OPT->DCMphaseBias) || + (OPT->_obsModelType == OPT->PPPRTK) ) { + required.push_back(new t_pppParam(t_pppParam::pBias1, t_prn(), t_lc::l1)); + required.push_back(new t_pppParam(t_pppParam::pBias2, t_prn(), t_lc::l2)); } Index: trunk/BNC/src/PPP/pppParlist.h =================================================================== --- trunk/BNC/src/PPP/pppParlist.h (revision 8904) +++ trunk/BNC/src/PPP/pppParlist.h (revision 8905) @@ -14,5 +14,5 @@ class t_pppParam { public: - enum e_type {crdX, crdY, crdZ, clkR, amb, offGG, trp}; + enum e_type {crdX, crdY, crdZ, clkR, amb, offGG, trp, ion, cBias1, cBias2, pBias1, pBias2}; t_pppParam(e_type type, const t_prn& prn, t_lc::type tLC, @@ -20,8 +20,8 @@ ~t_pppParam(); + e_type type() const {return _type;} double x0() const {return _x0;} - double partial(const bncTime& epoTime, const t_pppSatObs* obs, - const t_lc::type& tLC) const; + double partial(const bncTime& epoTime, const t_pppSatObs* obs, const t_lc::type& tLC) const; bool epoSpec() const {return _epoSpec;} bool isEqual(const t_pppParam* par2) const { @@ -51,4 +51,5 @@ unsigned ambNumEpo() const {return _ambInfo ? _ambInfo->_numEpo : 0;} void stepAmbNumEpo() {if (_ambInfo) _ambInfo->_numEpo += 1;} + void setRefPrn(t_prn prn) {_refPrn = prn;} static bool sortFunction(const t_pppParam* p1, const t_pppParam* p2) { @@ -78,16 +79,17 @@ unsigned _numEpo; }; - e_type _type; - t_prn _prn; - t_lc::type _tLC; - double _x0; - bool _epoSpec; - int _indexOld; - int _indexNew; - double _sigma0; - double _noise; - t_ambInfo* _ambInfo; - bncTime _lastObsTime; - bncTime _firstObsTime; + e_type _type; + t_prn _prn; + t_prn _refPrn; + t_lc::type _tLC; + double _x0; + bool _epoSpec; + int _indexOld; + int _indexNew; + double _sigma0; + double _noise; + t_ambInfo* _ambInfo; + bncTime _lastObsTime; + bncTime _firstObsTime; }; @@ -100,7 +102,8 @@ unsigned nPar() const {return _params.size();} const std::vector& params() const {return _params;} - std::vector& params() {return _params;} - void printResult(const bncTime& epoTime, const SymmetricMatrix& QQ, + std::vector& params() {return _params;} + void printResult(const bncTime& epoTime, const SymmetricMatrix& QQ, const ColumnVector& xx) const; + private: std::vector _params; Index: trunk/BNC/src/PPP/pppRefSat.cpp =================================================================== --- trunk/BNC/src/PPP/pppRefSat.cpp (revision 8905) +++ trunk/BNC/src/PPP/pppRefSat.cpp (revision 8905) @@ -0,0 +1,23 @@ +/* + * pppRefSat.cpp + * + * Created on: Feb 27, 2020 + * Author: A. Stuerze + */ + +#include "pppRefSat.h" + +using namespace BNC_PPP; + + +// Constructor +//////////////////////////////////////////////////////////////////////////// +t_pppRefSat::t_pppRefSat() { + _status = undefined; + _prn = t_prn(); + _stecValue = 0.0; +} + +t_pppRefSat::~t_pppRefSat(){ + +} Index: trunk/BNC/src/PPP/pppRefSat.h =================================================================== --- trunk/BNC/src/PPP/pppRefSat.h (revision 8905) +++ trunk/BNC/src/PPP/pppRefSat.h (revision 8905) @@ -0,0 +1,36 @@ +/* + * pppRefSat.h + * + * Created on: Feb 27, 2020 + * Author: stuerze + */ + +#ifndef PPPREFSAT_H_ +#define PPPREFSAT_H_ + +#include "t_prn.h" + +namespace BNC_PPP { + +class t_pppRefSat { +public: + t_pppRefSat (); + ~t_pppRefSat(); + enum e_status {undefined, initialized, changed, unchanged};// ev. wieder löschen + e_status status() const {return _status;} + void setStatus(e_status status) {_status = status;} + t_prn prn() {return _prn;} + void setPrn(t_prn prn) {_prn = prn;} + double stecValue() {return _stecValue;} + void setStecValue(double stecValue) {_stecValue = stecValue;} +private: + e_status _status; + t_prn _prn; + double _stecValue; +}; + + +} + + +#endif /* PPPREFSAT_H_ */ Index: trunk/BNC/src/PPP/pppSatObs.cpp =================================================================== --- trunk/BNC/src/PPP/pppSatObs.cpp (revision 8904) +++ trunk/BNC/src/PPP/pppSatObs.cpp (revision 8905) @@ -35,8 +35,11 @@ //////////////////////////////////////////////////////////////////////////// t_pppSatObs::t_pppSatObs(const t_satObs& pppSatObs) { - _prn = pppSatObs._prn; - _time = pppSatObs._time; - _outlier = false; - _valid = true; + _prn = pppSatObs._prn; + _time = pppSatObs._time; + _outlier = false; + _valid = true; + _reference = false; + _stecRefSat = 0.0; + _stecSat = 0.0; for (unsigned ii = 0; ii < t_frequency::max; ii++) { _obs[ii] = 0; @@ -61,6 +64,5 @@ // Select pseudoranges and phase observations // ------------------------------------------ - const string preferredAttrib = "CWPXI_"; - //const string preferredAttrib = "G:12&PWCSLXYN G:5&IQX R:12&PC R:3&IQX E:16&BCX E:578&IQX J:1&SLXCZ J:26&SLX J:5&IQX C:IQX I:ABCX S:1&C S:5&IQX"; + const string preferredAttrib = "G:12&PWCSLXYN G:5&IQX R:12&PC R:3&IQX E:16&BCX E:578&IQX J:1&SLXCZ J:26&SLX J:5&IQX C:IQX I:ABCX S:1&C S:5&IQX"; for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) { @@ -90,4 +92,5 @@ for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) { t_lc::type tLC = OPT->LCs(_prn.system())[ii]; + if (tLC == t_lc::GIM) {continue;} if (!isValid(tLC)) { _valid = false; @@ -121,5 +124,5 @@ ColumnVector dx = _xcSat - satPosOld; dx[3] *= t_CST::c; - if (dx.norm_Frobenius() < 1.e-4) { + if (dx.NormFrobenius() < 1.e-4) { totOK = true; break; @@ -140,18 +143,23 @@ void t_pppSatObs::lcCoeff(t_lc::type tLC, map& codeCoeff, - map& phaseCoeff) const { + map& phaseCoeff, + map& ionoCoeff) const { codeCoeff.clear(); phaseCoeff.clear(); + ionoCoeff.clear(); double f1 = t_CST::freq(_fType1, _channel); double f2 = t_CST::freq(_fType2, _channel); + double f1GPS = t_CST::freq(t_frequency::G1, 0); switch (tLC) { case t_lc::l1: - phaseCoeff[_fType1] = 1.0; + phaseCoeff[_fType1] = 1.0; + ionoCoeff [_fType1] = -1.0 * pow(f1GPS, 2) / pow(f1, 2); return; case t_lc::l2: - phaseCoeff[_fType2] = 1.0; + phaseCoeff[_fType2] = 1.0; + ionoCoeff [_fType2] = -1.0 * pow(f1GPS, 2) / pow(f2, 2); return; case t_lc::lIF: @@ -167,11 +175,13 @@ case t_lc::CL: phaseCoeff[_fType1] = 0.5; - codeCoeff[_fType1] = 0.5; + codeCoeff [_fType1] = 0.5; return; case t_lc::c1: codeCoeff[_fType1] = 1.0; + ionoCoeff[_fType1] = pow(f1GPS, 2) / pow(f1, 2); return; case t_lc::c2: codeCoeff[_fType2] = 1.0; + ionoCoeff[_fType2] = pow(f1GPS, 2) / pow(f2, 2); return; case t_lc::cIF: @@ -179,4 +189,5 @@ codeCoeff[_fType2] = -f2 * f2 / (f1 * f1 - f2 * f2); return; + case t_lc::GIM: case t_lc::dummy: case t_lc::maxLc: @@ -190,4 +201,5 @@ bool valid = true; obsValue(tLC, &valid); + //qDebug() << "tLC: " << tLC << " valid: " << valid; return valid; } @@ -195,13 +207,27 @@ //////////////////////////////////////////////////////////////////////////// double t_pppSatObs::obsValue(t_lc::type tLC, bool* valid) const { + + double retVal = 0.0; + if (valid) *valid = true; + + // Pseudo observations + if (tLC == t_lc::GIM) { + if (_stecRefSat == 0.0 || _stecSat == 0.0) { + if (valid) *valid = false; + return 0.0; + } + else { + return _stecRefSat - _stecSat; + } + } map codeCoeff; map phaseCoeff; - lcCoeff(tLC, codeCoeff, phaseCoeff); - - double retVal = 0.0; - if (valid) *valid = true; + map ionoCoeff; + lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff); map::const_iterator it; + + // Code observations for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) { t_frequency::type tFreq = it->first; @@ -214,4 +240,5 @@ } } + // Phase observations for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) { t_frequency::type tFreq = it->first; @@ -224,5 +251,4 @@ } } - return retVal; } @@ -258,15 +284,20 @@ double t_pppSatObs::sigma(t_lc::type tLC) const { + double retVal = 0.0; map codeCoeff; map phaseCoeff; - lcCoeff(tLC, codeCoeff, phaseCoeff); - - double retVal = 0.0; + map ionoCoeff; + lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff); + + if (tLC == t_lc::GIM) { + retVal = OPT->_sigmaGIMdiff * OPT->_sigmaGIMdiff; + } map::const_iterator it; - for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) { + for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {//qDebug() << "codeCoeff : " << t_frequency::toString(it->first).c_str() << ": " << it->second; retVal += it->second * it->second * OPT->_sigmaC1 * OPT->_sigmaC1; } - for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) { + + for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {//qDebug() << "phaseCoeff: " << t_frequency::toString(it->first).c_str() << ": " << it->second; retVal += it->second * it->second * OPT->_sigmaL1 * OPT->_sigmaL1; } @@ -295,20 +326,22 @@ // //////////////////////////////////////////////////////////////////////////// -double t_pppSatObs::maxRes(t_lc::type tLC) const { +double t_pppSatObs::maxRes(t_lc::type tLC) const {//qDebug() << "t_pppSatObs::maxRes(t_lc::type tLC)"; + double retVal = 0.0; map codeCoeff; map phaseCoeff; - lcCoeff(tLC, codeCoeff, phaseCoeff); - - double retVal = 0.0; + map ionoCoeff; + lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff); map::const_iterator it; - for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) { + for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {//qDebug() << "codeCoeff: " << it->first << ": " << it->second; retVal += it->second * it->second * OPT->_maxResC1 * OPT->_maxResC1; } - for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) { + for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {//qDebug() << "phaseCoeff: " << it->first << ": " << it->second; retVal += it->second * it->second * OPT->_maxResL1 * OPT->_maxResL1; } - + if (tLC == t_lc::GIM) { + retVal = 3 * (OPT->_sigmaGIMdiff * OPT->_sigmaGIMdiff); + } return sqrt(retVal); } @@ -326,6 +359,7 @@ // ------------------------------ ColumnVector rSat = _xcSat.Rows(1,3); - ColumnVector rhoV = rSat - station->xyzApr(); - _model._rho = rhoV.norm_Frobenius(); + ColumnVector rRec = station->xyzApr(); + ColumnVector rhoV = rSat - rRec; + _model._rho = rhoV.NormFrobenius(); ColumnVector vSat = _vvSat; @@ -334,5 +368,5 @@ xyz2neu(station->ellApr().data(), rhoV.data(), neu.data()); - _model._eleSat = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / _model._rho ); + _model._eleSat = acos(sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / _model._rho); if (neu[2] < 0) { _model._eleSat *= -1.0; @@ -354,5 +388,5 @@ Omega[1] = 0.0; Omega[2] = t_CST::omega / t_CST::c; - _model._sagnac = DotProduct(Omega, crossproduct(rSat, station->xyzApr())); + _model._sagnac = DotProduct(Omega, crossproduct(rSat, rRec)); // Antenna Eccentricity @@ -373,5 +407,5 @@ // Tropospheric Delay // ------------------ - _model._tropo = t_tropo::delay_saast(station->xyzApr(), _model._eleSat); + _model._tropo = t_tropo::delay_saast(rRec, _model._eleSat); // Code Biases @@ -392,10 +426,13 @@ // Phase Biases // ----------- - // TODO: consideration of fix indicators and jump counter const t_satPhaseBias* satPhaseBias = PPP_CLIENT->obsPool()->satPhaseBias(_prn); double yaw = 0.0; bool ssr = false; if (satPhaseBias) { - yaw = satPhaseBias->_yaw; + double dt = station->epochTime() - satPhaseBias->_time; + if (satPhaseBias->_updateInt) { + dt -= (0.5 * ssrUpdateInt[satPhaseBias->_updateInt]); + } + yaw = satPhaseBias->_yaw + satPhaseBias->_yawRate * dt; ssr = true; for (unsigned ii = 0; ii < satPhaseBias->_bias.size(); ii++) { @@ -414,8 +451,15 @@ _model._windUp = station->windUp(_time, _prn, rSat, ssr, yaw, vSat) ; + // Relativistic effect due to earth gravity + // ---------------------------------------- + double a = rSat.NormFrobenius() + rRec.NormFrobenius(); + double b = (rSat - rRec).NormFrobenius(); + double gm = 3.986004418e14; // m3/s2 + _model._rel = 2 * gm / t_CST::c / t_CST::c * log((a + b) / (a - b)); // Tidal Correction // ---------------- - _model._tide = -DotProduct(station->tideDspl(), rhoV) / _model._rho; + _model._tideEarth = -DotProduct(station->tideDsplEarth(), rhoV) / _model._rho; + _model._tideOcean = -DotProduct(station->tideDsplOcean(), rhoV) / _model._rho; // Ionospheric Delay @@ -429,19 +473,21 @@ } } + if (vTecUsage && vTec) { - double stec = station->stec(vTec, _signalPropagationTime, rSat); + double stec = station->stec(vTec, _signalPropagationTime, rSat); + double f1GPS = t_CST::freq(t_frequency::G1, 0); for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) { - t_frequency::type frqType = static_cast(iFreq); - _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(t_CST::freq(frqType, _channel), 2) * stec; - } - } - - // Relativistic effect due to earth gravity - // ---------------------------------------- - // TODO - - // Ocean Loading - // ------------- - // TODO + if (OPT->_pseudoObsIono) { // DCMcodeBias, DCMphaseBias + // For scaling the slant ionospheric delays the trick is to be consistent with units! + // The conversion of TECU into meters requires the frequency of the signal. + // Hence, GPS L1 frequency is used for all systems. The same is true for mu_i in lcCoeff(). + _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(f1GPS, 2) * stec; + } + else { // PPP-RTK + t_frequency::type frqType = static_cast(iFreq); + _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(t_CST::freq(frqType, _channel), 2) * stec; + } + } + } // Set Model Set Flag @@ -449,5 +495,5 @@ _model._set = true; - //printModel(); + printModel(); return success; @@ -457,41 +503,54 @@ //////////////////////////////////////////////////////////////////////////// void t_pppSatObs::printModel() const { - - LOG.setf(ios::fixed); - LOG << "MODEL for Satellite " << _prn.toString() << endl - << "RHO: " << setw(12) << setprecision(3) << _model._rho << endl - << "ELE: " << setw(12) << setprecision(3) << _model._eleSat * 180.0 / M_PI << endl - << "AZI: " << setw(12) << setprecision(3) << _model._azSat * 180.0 / M_PI << endl - << "SATCLK: " << setw(12) << setprecision(3) << _model._satClkM << endl - << "RECCLK: " << setw(12) << setprecision(3) << _model._recClkM << endl - << "SAGNAC: " << setw(12) << setprecision(3) << _model._sagnac << endl - << "ANTECC: " << setw(12) << setprecision(3) << _model._antEcc << endl - << "TROPO: " << setw(12) << setprecision(3) << _model._tropo << endl - << "WINDUP: " << setw(12) << setprecision(3) << _model._windUp << endl - << "TIDES: " << setw(12) << setprecision(3) << _model._tide << endl; +// TODO: cout should be LOG + cout.setf(ios::fixed); + cout << "\nMODEL for Satellite " << _prn.toString() << (isReference() ? " (Reference Satellite)" : "") << endl + << "======================= " << endl + << "PPP STRATEGY : " << OPT->_obsmodelTypeStr.at((int)OPT->_obsModelType).toLocal8Bit().constData() + << ((OPT->_pseudoObsIono) ? " with pseudo-observations for STEC" : "") << endl + << "RHO : " << setw(12) << setprecision(3) << _model._rho << endl + << "ELE : " << setw(12) << setprecision(3) << _model._eleSat * RHO_DEG << endl + << "AZI : " << setw(12) << setprecision(3) << _model._azSat * RHO_DEG << endl + << "SATCLK : " << setw(12) << setprecision(3) << _model._satClkM << endl + << "RECCLK : " << setw(12) << setprecision(3) << _model._recClkM << endl + << "SAGNAC : " << setw(12) << setprecision(3) << _model._sagnac << endl + << "ANTECC : " << setw(12) << setprecision(3) << _model._antEcc << endl + << "TROPO : " << setw(12) << setprecision(3) << _model._tropo << endl + << "WINDUP : " << setw(12) << setprecision(3) << _model._windUp << endl + << "REL : " << setw(12) << setprecision(3) << _model._rel << endl + << "EARTH TIDES : " << setw(12) << setprecision(3) << _model._tideEarth << endl + << "OCEAN TIDES : " << setw(12) << setprecision(3) << _model._tideOcean << endl + << endl + << "FREQUENCY DEPENDENT CORRECTIONS:" << endl + << "-------------------------------" << endl; for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) { if (_obs[iFreq]) { string frqStr = t_frequency::toString(t_frequency::type(iFreq)); if (_prn.system() == frqStr[0]) { - LOG << "PCO : " << frqStr << setw(12) << setprecision(3) << _model._antPCO[iFreq] << endl - << "BIAS CODE : " << frqStr << setw(12) << setprecision(3) << _model._codeBias[iFreq] << endl - << "BIAS PHASE : " << frqStr << setw(12) << setprecision(3) << _model._phaseBias[iFreq] << endl - << "IONO CODEDELAY: " << frqStr << setw(12) << setprecision(3) << _model._ionoCodeDelay[iFreq] << endl; - } - } - } + cout << "PCO : " << frqStr << setw(12) << setprecision(3) << _model._antPCO[iFreq] << endl + << "BIAS CODE : " << frqStr << setw(12) << setprecision(3) << _model._codeBias[iFreq] << endl + << "BIAS PHASE : " << frqStr << setw(12) << setprecision(3) << _model._phaseBias[iFreq] << endl + << "IONO CODEDELAY: " << frqStr << setw(12) << setprecision(3) << _model._ionoCodeDelay[iFreq]<< endl; + } + } + } +} + +// +//////////////////////////////////////////////////////////////////////////// +void t_pppSatObs::printObsMinusComputed() const { +// TODO: cout should be LOG + cout.setf(ios::fixed); + cout << "\nOBS-COMP for Satellite " << _prn.toString() << (isReference() ? " (Reference Satellite)" : "") << endl + << "========================== " << endl; for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) { t_lc::type tLC = OPT->LCs(_prn.system())[ii]; - LOG << "OBS-CMP " << t_lc::toString(tLC) << ": " << _prn.toString() << " " + cout << "OBS-CMP " << setw(4) << t_lc::toString(tLC) << ": " << _prn.toString() << " " << setw(12) << setprecision(3) << obsValue(tLC) << " " << setw(12) << setprecision(3) << cmpValue(tLC) << " " << setw(12) << setprecision(3) << obsValue(tLC) - cmpValue(tLC) << endl; - - } - LOG << "OBS-CMP MW: " << _prn.toString() << " " - << setw(12) << setprecision(3) << obsValue(t_lc::MW) << " " - << setw(12) << setprecision(3) << cmpValue(t_lc::MW) << " " - << setw(12) << setprecision(3) << obsValue(t_lc::MW) - cmpValue(t_lc::MW) << endl; -} + } +} + // @@ -504,37 +563,47 @@ //////////////////////////////////////////////////////////////////////////// double t_pppSatObs::cmpValue(t_lc::type tLC) const { - - if (!isValid(tLC)) { - return 0.0; - } - - // Non-Dispersive Part - // ------------------- - double nonDisp = _model._rho + _model._recClkM - _model._satClkM - + _model._sagnac + _model._antEcc + _model._tropo - + _model._tide; - - // Add Dispersive Part - // ------------------- - map codeCoeff; - map phaseCoeff; - lcCoeff(tLC, codeCoeff, phaseCoeff); - - double dispPart = 0.0; - - map::const_iterator it; - for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) { - t_frequency::type tFreq = it->first; - dispPart += it->second * (_model._antPCO[tFreq] - _model._codeBias[tFreq] + - _model._ionoCodeDelay[tFreq]); - } - for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) { - t_frequency::type tFreq = it->first; - dispPart += it->second * (_model._antPCO[tFreq] - _model._phaseBias[tFreq] + - _model._windUp * t_CST::lambda(tFreq, _channel) - - _model._ionoCodeDelay[tFreq]); - } - - return nonDisp + dispPart; + double cmpValue; + + if (!isValid(tLC)) { + cmpValue = 0.0; + } + else if (tLC == t_lc::GIM) { + cmpValue = 0.0; + } + else { + // Non-Dispersive Part + // ------------------- + double nonDisp = _model._rho + + _model._recClkM - _model._satClkM + + _model._sagnac + _model._antEcc + _model._tropo + + _model._tideEarth + _model._tideOcean + _model._rel; + + // Add Dispersive Part + // ------------------- + double dispPart = 0.0; + map codeCoeff; + map phaseCoeff; + map ionoCoeff; + lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff); + map::const_iterator it; + for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) { + t_frequency::type tFreq = it->first; + dispPart += it->second * (_model._antPCO[tFreq] - _model._codeBias[tFreq]); + if (OPT->PPPRTK) { + dispPart += it->second * (_model._ionoCodeDelay[tFreq]); + } + } + for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) { + t_frequency::type tFreq = it->first; + dispPart += it->second * (_model._antPCO[tFreq] - _model._phaseBias[tFreq] + + _model._windUp * t_CST::lambda(tFreq, _channel)); + if (OPT->PPPRTK) { + dispPart += it->second * (- _model._ionoCodeDelay[tFreq]); + } + } + cmpValue = nonDisp + dispPart; + } + + return cmpValue; } @@ -556,2 +625,9 @@ } } + +// +//////////////////////////////////////////////////////////////////////////// +void t_pppSatObs::setPseudoObsIono(t_frequency::type freq, double stecRefSat) { + _stecSat = _model._ionoCodeDelay[freq]; + _stecRefSat = stecRefSat; +} Index: trunk/BNC/src/PPP/pppSatObs.h =================================================================== --- trunk/BNC/src/PPP/pppSatObs.h (revision 8904) +++ trunk/BNC/src/PPP/pppSatObs.h (revision 8905) @@ -19,4 +19,6 @@ bool isValid() const {return _valid;}; bool isValid(t_lc::type tLC) const; + bool isReference() const {return _reference;}; + void setAsReference() {_reference = true;}; const t_prn& prn() const {return _prn;} const ColumnVector& xc() const {return _xcSat;} @@ -31,7 +33,9 @@ bool modelSet() const {return _model._set;} void printModel() const; + void printObsMinusComputed() const; void lcCoeff(t_lc::type tLC, std::map& codeCoeff, - std::map& phaseCoeff) const; + std::map& phaseCoeff, + std::map& ionoCoeff) const; double lambda(t_lc::type tLC) const; double sigma(t_lc::type tLC) const; @@ -41,4 +45,6 @@ void setRes(t_lc::type tLC, double res); double getRes(t_lc::type tLC) const; + void setPseudoObsIono(t_frequency::type freq, double stecRefSat); + double getIonoCodeDelay(t_frequency::type freq) {return _model._ionoCodeDelay[freq];} // RINEX @@ -79,16 +85,17 @@ ~t_model() {} void reset() { - _set = false; - _rho = 0.0; - _eleSat = 0.0; - _azSat = 0.0; - _recClkM = 0.0; - _satClkM = 0.0; - _sagnac = 0.0; - _antEcc = 0.0; - _tropo = 0.0; - _tide = 0.0; - _windUp = 0.0; - _rel = 0.0; + _set = false; + _rho = 0.0; + _eleSat = 0.0; + _azSat = 0.0; + _recClkM = 0.0; + _satClkM = 0.0; + _sagnac = 0.0; + _antEcc = 0.0; + _tropo = 0.0; + _tideEarth = 0.0; + _tideOcean = 0.0; + _windUp = 0.0; + _rel = 0.0; for (unsigned ii = 0; ii < t_frequency::max; ii++) { _antPCO[ii] = 0.0; @@ -107,5 +114,6 @@ double _antEcc; double _tropo; - double _tide; + double _tideEarth; + double _tideOcean; double _windUp; double _rel; @@ -119,4 +127,5 @@ bool _valid; + bool _reference; t_frequency::type _fType1; t_frequency::type _fType2; @@ -131,4 +140,6 @@ std::map _res; double _signalPropagationTime; + double _stecRefSat; + double _stecSat; }; Index: trunk/BNC/src/PPP/pppStation.h =================================================================== --- trunk/BNC/src/PPP/pppStation.h (revision 8904) +++ trunk/BNC/src/PPP/pppStation.h (revision 8905) @@ -21,6 +21,7 @@ void setNeuEcc(const ColumnVector& neuEcc); void setDClk(double dClk) {_dClk = dClk;} - void setTideDspl(const ColumnVector& tideDspl) {_tideDspl = tideDspl;} - void setIonoEpochTime(const bncTime& epochTime) {_epochTime = epochTime;} + void setTideDsplEarth(const ColumnVector& tideDsplEarth) {_tideDsplEarth = tideDsplEarth;} + void setTideDsplOcean(const ColumnVector& tideDsplOcean) {_tideDsplOcean = tideDsplOcean;} + void setEpochTime(const bncTime& epochTime) {_epochTime = epochTime;} const std::string& name() const {return _name;} const std::string& antName() const {return _antName;} @@ -29,5 +30,7 @@ const ColumnVector& neuEcc() const {return _neuEcc;} const ColumnVector& xyzEcc() const {return _xyzEcc;} - const ColumnVector& tideDspl() const {return _tideDspl;} + const ColumnVector& tideDsplEarth() const {return _tideDsplEarth;} + const ColumnVector& tideDsplOcean() const {return _tideDsplOcean;} + const bncTime epochTime() const {return _epochTime;} double dClk() const {return _dClk;} double windUp(const bncTime& time, t_prn prn, const ColumnVector& rSat, bool ssr, @@ -42,5 +45,6 @@ ColumnVector _neuEcc; ColumnVector _xyzEcc; - ColumnVector _tideDspl; + ColumnVector _tideDsplEarth; + ColumnVector _tideDsplOcean; double _dClk; mutable t_windUp* _windUp; Index: trunk/BNC/src/bncwindow.cpp =================================================================== --- trunk/BNC/src/bncwindow.cpp (revision 8904) +++ trunk/BNC/src/bncwindow.cpp (revision 8905) @@ -1023,13 +1023,18 @@ pppLayout3->addWidget(_pppWidgets._minEle, ir, 7); ++ir; - pppLayout3->addWidget(new QLabel("Wait for clock corr."), ir, 0, Qt::AlignLeft); - pppLayout3->addWidget(_pppWidgets._corrWaitTime, ir, 1); - pppLayout3->addWidget(new QLabel("Seeding (sec)"), ir, 3, Qt::AlignLeft); - pppLayout3->addWidget(_pppWidgets._seedingTime, ir, 4); + pppLayout3->addWidget(new QLabel("Model Obs"), ir, 0, Qt::AlignLeft); + pppLayout3->addWidget(_pppWidgets._modelObs, ir, 1); + pppLayout3->addWidget(new QLabel("Wait for clock corr."), ir, 3, Qt::AlignLeft); + pppLayout3->addWidget(_pppWidgets._corrWaitTime, ir, 4); + pppLayout3->addWidget(new QLabel("Seeding (sec)"), ir, 6, Qt::AlignLeft); + pppLayout3->addWidget(_pppWidgets._seedingTime, ir, 7); + ++ir; + pppLayout3->addWidget(new QLabel("Pseudo Obs"), ir, 0, Qt::AlignLeft); + pppLayout3->addWidget(_pppWidgets._pseudoObs, ir, 1); ++ir; pppLayout3->addWidget(new QLabel(""), ir, 8); pppLayout3->setColumnStretch(8, 999); ++ir; - pppLayout3->addWidget(new QLabel(""), ir, 1); + pppLayout3->addWidget(new QLabel(""), ir, 1); pppLayout3->setRowStretch(ir, 999); @@ -1402,5 +1407,6 @@ _pppWidgets._lcGLONASS->setWhatsThis(tr("

Specify which kind of GLONASS observations you want to use and on which kind of ionosphere-free linear combination of GLONASS observations you want to base ambiguity resolutions.

  • Specifying 'P3' means that you request BNC to use code data and so-called P3 ionosphere-free linear combination of code observations.
  • 'L3' means that you request BNC to use phase data and so-called L3 ionosphere-free linear combination of phase observations.
  • 'P3&L3' means that you request BNC to use both, code and phase data and so-called P3 and L3 ionosphere-free linear combination of code and phase observations.

Specifying 'no' means that you don't want BNC to use GLONASS data. [key: PPP/lcGLONASS]

")); _pppWidgets._lcGalileo->setWhatsThis(tr("

Specify which kind of Galileo observations you want to use and on which kind of of ionosphere-free linear combination of Galileo observations you want to base ambiguity resolutions.

  • Specifying 'P3' means that you request BNC to use code data and so-called P3 ionosphere-free linear combination of code observations.
  • 'L3' means that you request BNC to use phase data and so-called L3 ionosphere-free linear combination of phase observations.
  • 'P3&L3' means that you request BNC to use both, code and phase data and so-called P3 and L3 ionosphere-free linear combination of code and phase observations.

Specifying on of these options makes only sense if Galileo data are part of the processed observation stream.

Specifying 'no' means that you don't want BNC to use Galileo data. [key: PPP/lcGalileo]

")); - _pppWidgets._lcBDS->setWhatsThis(tr("

Specify which kind of BDS observations you want to use and on which kind of ionosphere-free linear combination of BDS observations you want to base ambiguity resolutions.

  • Specifying 'P3' means that you request BNC to use code data and so-called P3 ionosphere-free linear combination of code observations.
  • 'L3' means that you request BNC to use phase data and so-called L3 ionosphere-free linear combination of phase observations.
  • 'P3&L3' means that you request BNC to use both, code and phase data and so-called P3 and L3 ionosphere-free linear combination of code and phase observations.

Specifying on of these options makes only sense if BDS data are part of the processed observation stream.

Specifying 'no' means that you don't want to use BDS data. [key: PPP/lcBDS]

")); + _pppWidgets._lcBDS->setWhatsThis(tr("

Specify which kind of model for observations you want to use and on which kind of ionosphere-free linear combination of BDS observations you want to base ambiguity resolutions.

  • Specifying 'P3' means that you request BNC to use code data and so-called P3 ionosphere-free linear combination of code observations.
  • 'L3' means that you request BNC to use phase data and so-called L3 ionosphere-free linear combination of phase observations.
  • 'P3&L3' means that you request BNC to use both, code and phase data and so-called P3 and L3 ionosphere-free linear combination of code and phase observations.

Specifying on of these options makes only sense if BDS data are part of the processed observation stream.

Specifying 'no' means that you don't want to use BDS data. [key: PPP/lcBDS]

")); + _pppWidgets._modelObs->setWhatsThis(tr("

Specify which kind of PPP model you want to use. [key: PPP/modelObs]

")); _pppWidgets._sigmaC1->setWhatsThis(tr("

Enter a Sigma for GNSS C1 code observations in meters.

The higher the sigma you enter, the less the contribution of C1 code observations to a PPP solution from combined code and phase data. 2.0 is likely to be an appropriate choice.

Default is an empty option field, meaning
'Sigma C1 = 2.0' [key: PPP/sigmaC1]

")); _pppWidgets._sigmaL1->setWhatsThis(tr("

Enter a Sigma for GNSS L1 phase observations in meters.

The higher the sigma you enter, the less the contribution of L1 phase observations to a PPP solutions from combined code and phase data. 0.01 is likely to be an appropriate choice.

Default is an empty option field, meaning
'Sigma L1 = 0.01' [key: PPP/sigmaL1]

")); @@ -2855,5 +2861,5 @@ QComboBox* system = new QComboBox(); system->setEditable(false); - system->addItems(QString("ALL,GPS,GLONASS,Galileo,BDS,QZSS,SBAS").split(",")); + system->addItems(QString("ALL,GPS,GLONASS,Galileo,BDS,QZSS,SBAS,IRNSS").split(",")); system->setFrame(false); _uploadEphTable->setCellWidget(iRow, iCol, system); @@ -2941,5 +2947,5 @@ QComboBox* system = new QComboBox(); system->setEditable(false); - system->addItems(QString("ALL,GPS,GLONASS,Galileo,BDS,QZSS,SBAS").split(",")); + system->addItems(QString("ALL,GPS,GLONASS,Galileo,BDS,QZSS,SBAS,IRNSS").split(",")); system->setFrame(false); system->setCurrentIndex(system->findText(hlp[iCol])); Index: trunk/BNC/src/pppInclude.h =================================================================== --- trunk/BNC/src/pppInclude.h (revision 8904) +++ trunk/BNC/src/pppInclude.h (revision 8905) @@ -42,5 +42,5 @@ class t_lc { public: - enum type {dummy = 0, l1, l2, c1, c2, lIF, cIF, MW, CL, maxLc}; + enum type {dummy = 0, l1, l2, c1, c2, lIF, cIF, MW, CL, GIM, maxLc}; static bool includesPhase(type tt) { @@ -56,5 +56,8 @@ case cIF: return false; - case dummy: case maxLc: return false; + case dummy: + case maxLc: + case GIM: + return false; } return false; @@ -73,5 +76,8 @@ case lIF: return false; - case dummy: case maxLc: return false; + case dummy: + case maxLc: + case GIM: + return false; } return false; @@ -94,5 +100,8 @@ case lIF: case cIF: case MW: case CL: return t_frequency::dummy; - case dummy: case maxLc: return t_frequency::dummy; + case dummy: + case maxLc: + case GIM: + return t_frequency::dummy; } return t_frequency::dummy; @@ -109,5 +118,8 @@ case c2: return "c2"; case cIF: return "cIF"; - case dummy: case maxLc: return ""; + case GIM: return "GIM"; + case dummy: + case maxLc: + return ""; } return ""; Index: trunk/BNC/src/pppMain.cpp =================================================================== --- trunk/BNC/src/pppMain.cpp (revision 8904) +++ trunk/BNC/src/pppMain.cpp (revision 8905) @@ -109,9 +109,9 @@ } delete pppThread; - } + } #endif } _pppThreads.clear(); - } + } _running = false; @@ -183,103 +183,198 @@ opt->_corrWaitTime = 0; } + opt->_obsModelType = t_pppOptions::IF; + opt->_pseudoObsIono = false; + opt->_refSatRequired = false; +#ifdef USE_PPP + // Pseudo Observations + if (settings.value("PPP/pseudoObs").toString() == "Ionosphere") { + opt->_pseudoObsIono = true; + } + else if (settings.value("PPP/pseudoObs").toString() == "no") { + opt->_pseudoObsIono = false; + } + // Observation Model + if (settings.value("PPP/modelObs").toString() == "Ionosphere-free PPP") { + opt->_obsModelType = t_pppOptions::IF; + opt->_pseudoObsIono = false; + } + else if (settings.value("PPP/modelObs").toString() == "PPP-RTK") { + opt->_obsModelType = t_pppOptions::PPPRTK; + opt->_pseudoObsIono = false; + } + else if (settings.value("PPP/modelObs").toString() == "Uncombined PPP") { + opt->_obsModelType = t_pppOptions::UncombPPP; + if (opt->_pseudoObsIono) { + opt->_refSatRequired = true; + } + } + else if (settings.value("PPP/modelObs").toString() == "DCM with Code Biases") { + opt->_obsModelType = t_pppOptions::DCMcodeBias; + opt->_refSatRequired = true; + } + else if (settings.value("PPP/modelObs").toString() == "DCM with Phase Biases") { + opt->_obsModelType = t_pppOptions::DCMphaseBias; + opt->_refSatRequired = true; + } +#endif // GPS - if (settings.value("PPP/lcGPS").toString() == "Pi") { - opt->_LCsGPS.push_back(t_lc::c1); - opt->_LCsGPS.push_back(t_lc::c2); + if (settings.value("PPP/lcGPS").toString() == "Pi") { + if (opt->_obsModelType == t_pppOptions::IF) { + opt->_LCsGPS.push_back(t_lc::cIF); + } + else { + opt->_LCsGPS.push_back(t_lc::c1); + opt->_LCsGPS.push_back(t_lc::c2); + if (opt->_pseudoObsIono) { + opt->_LCsGPS.push_back(t_lc::GIM); + } + } } else if (settings.value("PPP/lcGPS").toString() == "Li") { - opt->_LCsGPS.push_back(t_lc::l1); - opt->_LCsGPS.push_back(t_lc::l2); + if (opt->_obsModelType == t_pppOptions::IF) { + opt->_LCsGPS.push_back(t_lc::lIF); + } + else { + opt->_LCsGPS.push_back(t_lc::l1); + opt->_LCsGPS.push_back(t_lc::l2); + if (opt->_pseudoObsIono) { + opt->_LCsGPS.push_back(t_lc::GIM); + } + } } else if (settings.value("PPP/lcGPS").toString() == "Pi&Li") { - opt->_LCsGPS.push_back(t_lc::c1); - opt->_LCsGPS.push_back(t_lc::c2); - opt->_LCsGPS.push_back(t_lc::l1); - opt->_LCsGPS.push_back(t_lc::l2); - } - if (settings.value("PPP/lcGPS").toString() == "P3") { - opt->_LCsGPS.push_back(t_lc::cIF); - } - else if (settings.value("PPP/lcGPS").toString() == "L3") { - opt->_LCsGPS.push_back(t_lc::lIF); - } - else if (settings.value("PPP/lcGPS").toString() == "P3&L3") { - opt->_LCsGPS.push_back(t_lc::cIF); - opt->_LCsGPS.push_back(t_lc::lIF); + if (opt->_obsModelType == t_pppOptions::IF) { + opt->_LCsGPS.push_back(t_lc::cIF); + opt->_LCsGPS.push_back(t_lc::lIF); + } + else { + opt->_LCsGPS.push_back(t_lc::c1); + opt->_LCsGPS.push_back(t_lc::c2); + opt->_LCsGPS.push_back(t_lc::l1); + opt->_LCsGPS.push_back(t_lc::l2); + if (opt->_pseudoObsIono) { + opt->_LCsGPS.push_back(t_lc::GIM); + } + } } // GLONASS - if (settings.value("PPP/lcGLONASS").toString() == "Pi") { - opt->_LCsGLONASS.push_back(t_lc::c1); - opt->_LCsGLONASS.push_back(t_lc::c2); + if (settings.value("PPP/lcGLONASS").toString() == "Pi") { + if (opt->_obsModelType == t_pppOptions::IF) { + opt->_LCsGLONASS.push_back(t_lc::cIF); + } + else { + opt->_LCsGLONASS.push_back(t_lc::c1); + opt->_LCsGLONASS.push_back(t_lc::c2); + if (opt->_pseudoObsIono) { + opt->_LCsGLONASS.push_back(t_lc::GIM); + } + } } else if (settings.value("PPP/lcGLONASS").toString() == "Li") { - opt->_LCsGLONASS.push_back(t_lc::l1); - opt->_LCsGLONASS.push_back(t_lc::l2); + if (opt->_obsModelType == t_pppOptions::IF) { + opt->_LCsGLONASS.push_back(t_lc::lIF); + } + else { + opt->_LCsGLONASS.push_back(t_lc::l1); + opt->_LCsGLONASS.push_back(t_lc::l2); + if (opt->_obsModelType == t_pppOptions::IF) { + opt->_LCsGLONASS.push_back(t_lc::GIM); + } + } } else if (settings.value("PPP/lcGLONASS").toString() == "Pi&Li") { - opt->_LCsGLONASS.push_back(t_lc::c1); - opt->_LCsGLONASS.push_back(t_lc::c2); - opt->_LCsGLONASS.push_back(t_lc::l1); - opt->_LCsGLONASS.push_back(t_lc::l2); - } - if (settings.value("PPP/lcGLONASS").toString() == "P3") { - opt->_LCsGLONASS.push_back(t_lc::cIF); - } - else if (settings.value("PPP/lcGLONASS").toString() == "L3") { - opt->_LCsGLONASS.push_back(t_lc::lIF); - } - else if (settings.value("PPP/lcGLONASS").toString() == "P3&L3") { - opt->_LCsGLONASS.push_back(t_lc::cIF); - opt->_LCsGLONASS.push_back(t_lc::lIF); + if (opt->_obsModelType == t_pppOptions::IF) { + opt->_LCsGLONASS.push_back(t_lc::cIF); + opt->_LCsGLONASS.push_back(t_lc::lIF); + } + else { + opt->_LCsGLONASS.push_back(t_lc::c1); + opt->_LCsGLONASS.push_back(t_lc::c2); + opt->_LCsGLONASS.push_back(t_lc::l1); + opt->_LCsGLONASS.push_back(t_lc::l2); + if (opt->_pseudoObsIono) { + opt->_LCsGLONASS.push_back(t_lc::GIM); + } + } } // Galileo - if (settings.value("PPP/lcGalileo").toString() == "Pi") { - opt->_LCsGalileo.push_back(t_lc::c1); - opt->_LCsGalileo.push_back(t_lc::c2); + if (settings.value("PPP/lcGalileo").toString() == "Pi") { + if (opt->_obsModelType == t_pppOptions::IF) { + opt->_LCsGalileo.push_back(t_lc::cIF); + } + else { + opt->_LCsGalileo.push_back(t_lc::c1); + opt->_LCsGalileo.push_back(t_lc::c2); + if (opt->_pseudoObsIono) { + opt->_LCsGalileo.push_back(t_lc::GIM); + } + } } else if (settings.value("PPP/lcGalileo").toString() == "Li") { - opt->_LCsGalileo.push_back(t_lc::l1); - opt->_LCsGalileo.push_back(t_lc::l2); + if (opt->_obsModelType == t_pppOptions::IF) { + opt->_LCsGalileo.push_back(t_lc::lIF); + } + else { + opt->_LCsGalileo.push_back(t_lc::l1); + opt->_LCsGalileo.push_back(t_lc::l2); + if (opt->_pseudoObsIono) { + opt->_LCsGalileo.push_back(t_lc::GIM); + } + } } else if (settings.value("PPP/lcGalileo").toString() == "Pi&Li") { - opt->_LCsGalileo.push_back(t_lc::c1); - opt->_LCsGalileo.push_back(t_lc::c2); - opt->_LCsGalileo.push_back(t_lc::l1); - opt->_LCsGalileo.push_back(t_lc::l2); - } - if (settings.value("PPP/lcGalileo").toString() == "P3") { - opt->_LCsGalileo.push_back(t_lc::cIF); - } - else if (settings.value("PPP/lcGalileo").toString() == "L3") { - opt->_LCsGalileo.push_back(t_lc::lIF); - } - else if (settings.value("PPP/lcGalileo").toString() == "P3&L3") { - opt->_LCsGalileo.push_back(t_lc::cIF); - opt->_LCsGalileo.push_back(t_lc::lIF); + if (opt->_obsModelType == t_pppOptions::IF) { + opt->_LCsGalileo.push_back(t_lc::cIF); + opt->_LCsGalileo.push_back(t_lc::lIF); + } + else { + opt->_LCsGalileo.push_back(t_lc::c1); + opt->_LCsGalileo.push_back(t_lc::c2); + opt->_LCsGalileo.push_back(t_lc::l1); + opt->_LCsGalileo.push_back(t_lc::l2); + if (opt->_pseudoObsIono) { + opt->_LCsGalileo.push_back(t_lc::GIM); + } + } } // BDS - if (settings.value("PPP/lcBDS").toString() == "Pi") { - opt->_LCsBDS.push_back(t_lc::c1); - opt->_LCsBDS.push_back(t_lc::c2); + if (settings.value("PPP/lcBDS").toString() == "Pi") { + if (opt->_obsModelType == t_pppOptions::IF) { + opt->_LCsBDS.push_back(t_lc::cIF); + } + else { + opt->_LCsBDS.push_back(t_lc::c1); + opt->_LCsBDS.push_back(t_lc::c2); + if (opt->_pseudoObsIono) { + opt->_LCsBDS.push_back(t_lc::GIM); + } + } } else if (settings.value("PPP/lcBDS").toString() == "Li") { - opt->_LCsBDS.push_back(t_lc::l1); - opt->_LCsBDS.push_back(t_lc::l2); + if (opt->_obsModelType == t_pppOptions::IF) { + opt->_LCsBDS.push_back(t_lc::lIF); + } + else { + opt->_LCsBDS.push_back(t_lc::l1); + opt->_LCsBDS.push_back(t_lc::l2); + if (opt->_pseudoObsIono) { + opt->_LCsBDS.push_back(t_lc::GIM); + } + } } else if (settings.value("PPP/lcBDS").toString() == "Pi&Li") { - opt->_LCsBDS.push_back(t_lc::c1); - opt->_LCsBDS.push_back(t_lc::c2); - opt->_LCsBDS.push_back(t_lc::l1); - opt->_LCsBDS.push_back(t_lc::l2); - } - if (settings.value("PPP/lcBDS").toString() == "P3") { - opt->_LCsBDS.push_back(t_lc::cIF); - } - else if (settings.value("PPP/lcBDS").toString() == "L3") { - opt->_LCsBDS.push_back(t_lc::lIF); - } - else if (settings.value("PPP/lcBDS").toString() == "P3&L3") { - opt->_LCsBDS.push_back(t_lc::cIF); - opt->_LCsBDS.push_back(t_lc::lIF); + if (opt->_obsModelType == t_pppOptions::IF) { + opt->_LCsBDS.push_back(t_lc::cIF); + opt->_LCsBDS.push_back(t_lc::lIF); + } + else { + opt->_LCsBDS.push_back(t_lc::c1); + opt->_LCsBDS.push_back(t_lc::c2); + opt->_LCsBDS.push_back(t_lc::l1); + opt->_LCsBDS.push_back(t_lc::l2); + if (opt->_pseudoObsIono) { + opt->_LCsBDS.push_back(t_lc::GIM); + } + } } @@ -316,6 +411,13 @@ // Some default values // ------------------- - opt->_aprSigAmb = 1000.0; - opt->_noiseClk = 1000.0; + opt->_aprSigAmb = 1000.0; + opt->_aprSigIon = 1000.0; + opt->_noiseClk = 1000.0; + opt->_aprSigCodeBias = 1000.0; + opt->_aprSigPhaseBias = 1000.0; + opt->_noiseIon = 1.0; + opt->_noiseCodeBias = 1.0; + opt->_noisePhaseBias = 0.1; + opt->_sigmaGIMdiff = 0.05; // pseudo observation GIM: STEC(ref_sat) - STEC(sat) _options << opt; Index: trunk/BNC/src/pppModel.cpp =================================================================== --- trunk/BNC/src/pppModel.cpp (revision 8904) +++ trunk/BNC/src/pppModel.cpp (revision 8905) @@ -1,3 +1,2 @@ - // Part of BNC, a utility for retrieving decoding and // converting GNSS data streams from NTRIP broadcasters. @@ -40,5 +39,4 @@ * -----------------------------------------------------------------------*/ - #include @@ -48,15 +46,19 @@ using namespace std; -const double t_astro::RHO_DEG = 180.0 / M_PI; -const double t_astro::RHO_SEC = 3600.0 * 180.0 / M_PI; -const double t_astro::MJD_J2000 = 51544.5; + Matrix t_astro::rotX(double Angle) { const double C = cos(Angle); const double S = sin(Angle); - Matrix UU(3,3); - UU[0][0] = 1.0; UU[0][1] = 0.0; UU[0][2] = 0.0; - UU[1][0] = 0.0; UU[1][1] = +C; UU[1][2] = +S; - UU[2][0] = 0.0; UU[2][1] = -S; UU[2][2] = +C; + Matrix UU(3, 3); + UU[0][0] = 1.0; + UU[0][1] = 0.0; + UU[0][2] = 0.0; + UU[1][0] = 0.0; + UU[1][1] = +C; + UU[1][2] = +S; + UU[2][0] = 0.0; + UU[2][1] = -S; + UU[2][2] = +C; return UU; } @@ -65,8 +67,14 @@ const double C = cos(Angle); const double S = sin(Angle); - Matrix UU(3,3); - UU[0][0] = +C; UU[0][1] = 0.0; UU[0][2] = -S; - UU[1][0] = 0.0; UU[1][1] = 1.0; UU[1][2] = 0.0; - UU[2][0] = +S; UU[2][1] = 0.0; UU[2][2] = +C; + Matrix UU(3, 3); + UU[0][0] = +C; + UU[0][1] = 0.0; + UU[0][2] = -S; + UU[1][0] = 0.0; + UU[1][1] = 1.0; + UU[1][2] = 0.0; + UU[2][0] = +S; + UU[2][1] = 0.0; + UU[2][2] = +C; return UU; } @@ -75,8 +83,14 @@ const double C = cos(Angle); const double S = sin(Angle); - Matrix UU(3,3); - UU[0][0] = +C; UU[0][1] = +S; UU[0][2] = 0.0; - UU[1][0] = -S; UU[1][1] = +C; UU[1][2] = 0.0; - UU[2][0] = 0.0; UU[2][1] = 0.0; UU[2][2] = 1.0; + Matrix UU(3, 3); + UU[0][0] = +C; + UU[0][1] = +S; + UU[0][2] = 0.0; + UU[1][0] = -S; + UU[1][1] = +C; + UU[1][2] = 0.0; + UU[2][0] = 0.0; + UU[2][1] = 0.0; + UU[2][2] = 1.0; return UU; } @@ -89,12 +103,12 @@ double Mjd_0 = floor(Mjd_UT1); - double UT1 = Secs*(Mjd_UT1-Mjd_0); - double T_0 = (Mjd_0 -MJD_J2000)/36525.0; - double T = (Mjd_UT1-MJD_J2000)/36525.0; - - double gmst = 24110.54841 + 8640184.812866*T_0 + 1.002737909350795*UT1 - + (0.093104-6.2e-6*T)*T*T; - - return 2.0*M_PI*Frac(gmst/Secs); + double UT1 = Secs * (Mjd_UT1 - Mjd_0); + double T_0 = (Mjd_0 - MJD_J2000) / 36525.0; + double T = (Mjd_UT1 - MJD_J2000) / 36525.0; + + double gmst = 24110.54841 + 8640184.812866 * T_0 + 1.002737909350795 * UT1 + + (0.093104 - 6.2e-6 * T) * T * T; + + return 2.0 * M_PI * Frac(gmst / Secs); } @@ -103,19 +117,21 @@ Matrix t_astro::NutMatrix(double Mjd_TT) { - const double T = (Mjd_TT-MJD_J2000)/36525.0; - - double ls = 2.0*M_PI*Frac(0.993133+ 99.997306*T); - double D = 2.0*M_PI*Frac(0.827362+1236.853087*T); - double F = 2.0*M_PI*Frac(0.259089+1342.227826*T); - double N = 2.0*M_PI*Frac(0.347346- 5.372447*T); - - double dpsi = ( -17.200*sin(N) - 1.319*sin(2*(F-D+N)) - 0.227*sin(2*(F+N)) - + 0.206*sin(2*N) + 0.143*sin(ls) ) / RHO_SEC; - double deps = ( + 9.203*cos(N) + 0.574*cos(2*(F-D+N)) + 0.098*cos(2*(F+N)) - - 0.090*cos(2*N) ) / RHO_SEC; - - double eps = 0.4090928-2.2696E-4*T; - - return rotX(-eps-deps)*rotZ(-dpsi)*rotX(+eps); + const double T = (Mjd_TT - MJD_J2000) / 36525.0; + + double ls = 2.0 * M_PI * Frac(0.993133 + 99.997306 * T); + double D = 2.0 * M_PI * Frac(0.827362 + 1236.853087 * T); + double F = 2.0 * M_PI * Frac(0.259089 + 1342.227826 * T); + double N = 2.0 * M_PI * Frac(0.347346 - 5.372447 * T); + + double dpsi = (-17.200 * sin(N) - 1.319 * sin(2 * (F - D + N)) + - 0.227 * sin(2 * (F + N)) + + 0.206 * sin(2 * N) + 0.143 * sin(ls)) / RHO_SEC; + double deps = (+9.203 * cos(N) + 0.574 * cos(2 * (F - D + N)) + + 0.098 * cos(2 * (F + N)) + - 0.090 * cos(2 * N)) / RHO_SEC; + + double eps = 0.4090928 - 2.2696E-4 * T; + + return rotX(-eps - deps) * rotZ(-dpsi) * rotX(+eps); } @@ -124,12 +140,13 @@ Matrix t_astro::PrecMatrix(double Mjd_1, double Mjd_2) { - const double T = (Mjd_1-MJD_J2000)/36525.0; - const double dT = (Mjd_2-Mjd_1)/36525.0; - - double zeta = ( (2306.2181+(1.39656-0.000139*T)*T)+ - ((0.30188-0.000344*T)+0.017998*dT)*dT )*dT/RHO_SEC; - double z = zeta + ( (0.79280+0.000411*T)+0.000205*dT)*dT*dT/RHO_SEC; - double theta = ( (2004.3109-(0.85330+0.000217*T)*T)- - ((0.42665+0.000217*T)+0.041833*dT)*dT )*dT/RHO_SEC; + const double T = (Mjd_1 - MJD_J2000) / 36525.0; + const double dT = (Mjd_2 - Mjd_1) / 36525.0; + + double zeta = ((2306.2181 + (1.39656 - 0.000139 * T) * T) + + ((0.30188 - 0.000344 * T) + 0.017998 * dT) * dT) * dT / RHO_SEC; + double z = zeta + + ((0.79280 + 0.000411 * T) + 0.000205 * dT) * dT * dT / RHO_SEC; + double theta = ((2004.3109 - (0.85330 + 0.000217 * T) * T) - + ((0.42665 + 0.000217 * T) + 0.041833 * dT) * dT) * dT / RHO_SEC; return rotZ(-z) * rotY(theta) * rotZ(-zeta); @@ -140,19 +157,20 @@ ColumnVector t_astro::Sun(double Mjd_TT) { - const double eps = 23.43929111/RHO_DEG; - const double T = (Mjd_TT-MJD_J2000)/36525.0; - - double M = 2.0*M_PI * Frac ( 0.9931267 + 99.9973583*T); - double L = 2.0*M_PI * Frac ( 0.7859444 + M/2.0/M_PI + - (6892.0*sin(M)+72.0*sin(2.0*M)) / 1296.0e3); - double r = 149.619e9 - 2.499e9*cos(M) - 0.021e9*cos(2*M); + const double eps = 23.43929111 / RHO_DEG; + const double T = (Mjd_TT - MJD_J2000) / 36525.0; + + double M = 2.0 * M_PI * Frac(0.9931267 + 99.9973583 * T); + double L = 2.0 * M_PI * Frac(0.7859444 + M / 2.0 / M_PI + + (6892.0 * sin(M) + 72.0 * sin(2.0 * M)) / 1296.0e3); + double r = 149.619e9 - 2.499e9 * cos(M) - 0.021e9 * cos(2 * M); ColumnVector r_Sun(3); - r_Sun << r*cos(L) << r*sin(L) << 0.0; r_Sun = rotX(-eps) * r_Sun; - - return rotZ(GMST(Mjd_TT)) - * NutMatrix(Mjd_TT) - * PrecMatrix(MJD_J2000, Mjd_TT) - * r_Sun; + r_Sun << r * cos(L) << r * sin(L) << 0.0; + r_Sun = rotX(-eps) * r_Sun; + + return rotZ(GMST(Mjd_TT)) + * NutMatrix(Mjd_TT) + * PrecMatrix(MJD_J2000, Mjd_TT) + * r_Sun; } @@ -161,49 +179,56 @@ ColumnVector t_astro::Moon(double Mjd_TT) { - const double eps = 23.43929111/RHO_DEG; - const double T = (Mjd_TT-MJD_J2000)/36525.0; - - double L_0 = Frac ( 0.606433 + 1336.851344*T ); - double l = 2.0*M_PI*Frac ( 0.374897 + 1325.552410*T ); - double lp = 2.0*M_PI*Frac ( 0.993133 + 99.997361*T ); - double D = 2.0*M_PI*Frac ( 0.827361 + 1236.853086*T ); - double F = 2.0*M_PI*Frac ( 0.259086 + 1342.227825*T ); - - double dL = +22640*sin(l) - 4586*sin(l-2*D) + 2370*sin(2*D) + 769*sin(2*l) - -668*sin(lp) - 412*sin(2*F) - 212*sin(2*l-2*D)- 206*sin(l+lp-2*D) - +192*sin(l+2*D) - 165*sin(lp-2*D) - 125*sin(D) - 110*sin(l+lp) - +148*sin(l-lp) - 55*sin(2*F-2*D); - - double L = 2.0*M_PI * Frac( L_0 + dL/1296.0e3 ); - - double S = F + (dL+412*sin(2*F)+541*sin(lp)) / RHO_SEC; - double h = F-2*D; - double N = -526*sin(h) + 44*sin(l+h) - 31*sin(-l+h) - 23*sin(lp+h) - +11*sin(-lp+h) - 25*sin(-2*l+F) + 21*sin(-l+F); - - double B = ( 18520.0*sin(S) + N ) / RHO_SEC; + const double eps = 23.43929111 / RHO_DEG; + const double T = (Mjd_TT - MJD_J2000) / 36525.0; + + double L_0 = Frac(0.606433 + 1336.851344 * T); + double l = 2.0 * M_PI * Frac(0.374897 + 1325.552410 * T); + double lp = 2.0 * M_PI * Frac(0.993133 + 99.997361 * T); + double D = 2.0 * M_PI * Frac(0.827361 + 1236.853086 * T); + double F = 2.0 * M_PI * Frac(0.259086 + 1342.227825 * T); + + double dL = +22640 * sin(l) - 4586 * sin(l - 2 * D) + 2370 * sin(2 * D) + + 769 * sin(2 * l) + - 668 * sin(lp) - 412 * sin(2 * F) - 212 * sin(2 * l - 2 * D) + - 206 * sin(l + lp - 2 * D) + + 192 * sin(l + 2 * D) - 165 * sin(lp - 2 * D) - 125 * sin(D) + - 110 * sin(l + lp) + + 148 * sin(l - lp) - 55 * sin(2 * F - 2 * D); + + double L = 2.0 * M_PI * Frac(L_0 + dL / 1296.0e3); + + double S = F + (dL + 412 * sin(2 * F) + 541 * sin(lp)) / RHO_SEC; + double h = F - 2 * D; + double N = -526 * sin(h) + 44 * sin(l + h) - 31 * sin(-l + h) + - 23 * sin(lp + h) + + 11 * sin(-lp + h) - 25 * sin(-2 * l + F) + 21 * sin(-l + F); + + double B = (18520.0 * sin(S) + N) / RHO_SEC; double cosB = cos(B); - double R = 385000e3 - 20905e3*cos(l) - 3699e3*cos(2*D-l) - 2956e3*cos(2*D) - -570e3*cos(2*l) + 246e3*cos(2*l-2*D) - 205e3*cos(lp-2*D) - -171e3*cos(l+2*D) - 152e3*cos(l+lp-2*D); + double R = 385000e3 - 20905e3 * cos(l) - 3699e3 * cos(2 * D - l) + - 2956e3 * cos(2 * D) + - 570e3 * cos(2 * l) + 246e3 * cos(2 * l - 2 * D) + - 205e3 * cos(lp - 2 * D) + - 171e3 * cos(l + 2 * D) - 152e3 * cos(l + lp - 2 * D); ColumnVector r_Moon(3); - r_Moon << R*cos(L)*cosB << R*sin(L)*cosB << R*sin(B); + r_Moon << R * cos(L) * cosB << R * sin(L) * cosB << R * sin(B); r_Moon = rotX(-eps) * r_Moon; - return rotZ(GMST(Mjd_TT)) - * NutMatrix(Mjd_TT) - * PrecMatrix(MJD_J2000, Mjd_TT) - * r_Moon; + return rotZ(GMST(Mjd_TT)) + * NutMatrix(Mjd_TT) + * PrecMatrix(MJD_J2000, Mjd_TT) + * r_Moon; } // Tidal Correction //////////////////////////////////////////////////////////////////////////// -ColumnVector t_tides::displacement(const bncTime& time, const ColumnVector& xyz) { +ColumnVector t_tides::earth(const bncTime& time, const ColumnVector& xyz) { if (time.undef()) { - ColumnVector dX(3); dX = 0.0; + ColumnVector dX(3); + dX = 0.0; return dX; } @@ -214,12 +239,12 @@ _lastMjd = Mjd; _xSun = t_astro::Sun(Mjd); - _rSun = sqrt(DotProduct(_xSun,_xSun)); + _rSun = sqrt(DotProduct(_xSun, _xSun)); _xSun /= _rSun; _xMoon = t_astro::Moon(Mjd); - _rMoon = sqrt(DotProduct(_xMoon,_xMoon)); + _rMoon = sqrt(DotProduct(_xMoon, _xMoon)); _xMoon /= _rMoon; } - double rRec = sqrt(DotProduct(xyz, xyz)); + double rRec = sqrt(DotProduct(xyz, xyz)); ColumnVector xyzUnit = xyz / rRec; @@ -231,25 +256,25 @@ // Tidal Displacement // ------------------ - double scSun = DotProduct(xyzUnit, _xSun); + double scSun = DotProduct(xyzUnit, _xSun); double scMoon = DotProduct(xyzUnit, _xMoon); - double p2Sun = 3.0 * (H2/2.0-L2) * scSun * scSun - H2/2.0; - double p2Moon = 3.0 * (H2/2.0-L2) * scMoon * scMoon - H2/2.0; - - double x2Sun = 3.0 * L2 * scSun; + double p2Sun = 3.0 * (H2 / 2.0 - L2) * scSun * scSun - H2 / 2.0; + double p2Moon = 3.0 * (H2 / 2.0 - L2) * scMoon * scMoon - H2 / 2.0; + + double x2Sun = 3.0 * L2 * scSun; double x2Moon = 3.0 * L2 * scMoon; const double gmWGS = 398.6005e12; - const double gms = 1.3271250e20; - const double gmm = 4.9027890e12; - - double facSun = gms / gmWGS * - (rRec * rRec * rRec * rRec) / (_rSun * _rSun * _rSun); + const double gms = 1.3271250e20; + const double gmm = 4.9027890e12; + + double facSun = gms / gmWGS * + (rRec * rRec * rRec * rRec) / (_rSun * _rSun * _rSun); double facMoon = gmm / gmWGS * - (rRec * rRec * rRec * rRec) / (_rMoon * _rMoon * _rMoon); - - ColumnVector dX = facSun * (x2Sun * _xSun + p2Sun * xyzUnit) + - facMoon * (x2Moon * _xMoon + p2Moon * xyzUnit); + (rRec * rRec * rRec * rRec) / (_rMoon * _rMoon * _rMoon); + + ColumnVector dX = facSun * (x2Sun * _xSun + p2Sun * xyzUnit) + + facMoon * (x2Moon * _xMoon + p2Moon * xyzUnit); return dX; @@ -258,12 +283,12 @@ // Constructor /////////////////////////////////////////////////////////////////////////// -t_loading::t_loading(const QString& fileName) { - +t_tides::t_tides() { + _lastMjd = 0.0; + _rSun = 0.0; + _rMoon = 0.0; newBlqData = 0; - readFile(fileName); - printAll(); -} - -t_loading::~t_loading() { +} + +t_tides::~t_tides() { if (newBlqData) { @@ -278,5 +303,5 @@ } -t_irc t_loading::readFile(const QString& fileName) { +t_irc t_tides::readBlqFile(const char* fileName) { QFile inFile(fileName); inFile.open(QIODevice::ReadOnly | QIODevice::Text); @@ -285,43 +310,45 @@ QString site = QString(); - while ( !in.atEnd() ) { + while (!in.atEnd()) { QString line = in.readLine(); // skip empty lines and comments - if (line.indexOf("$$") != -1 || line.size() == 0) { + if (line.indexOf("$$") != -1) { continue; } - line = line.trimmed(); QTextStream inLine(line.toLatin1(), QIODevice::ReadOnly); - switch (row) { case 0: site = line; site = site.toUpper(); - if (!newBlqData) { - newBlqData = new t_blqData; - newBlqData->amplitudes.ReSize(3,11); - newBlqData->phases.ReSize(3,11); + newBlqData = new t_blqData; + newBlqData->amplitudes.ReSize(3, 11); + newBlqData->phases.ReSize(3, 11); + break; + case 1: + case 2: + case 3: + for (int ii = 0; ii < 11; ii++) { + inLine >> newBlqData->amplitudes[row - 1][ii]; } break; - case 1: case 2: case 3: + case 4: + case 5: for (int ii = 0; ii < 11; ii++) { - inLine >> newBlqData->amplitudes[row-1][ii]; + inLine >> newBlqData->phases[row - 4][ii]; } break; - case 4: case 5: case 6: + case 6: for (int ii = 0; ii < 11; ii++) { - inLine >> newBlqData->phases[row-4][ii]; + inLine >> newBlqData->phases[row - 4][ii]; } - break; - case 7: if (newBlqData && !site.isEmpty()) { blqMap[site] = newBlqData; site = QString(); - row = -1; newBlqData = 0; } + row = -1; break; } @@ -332,7 +359,134 @@ } +ColumnVector t_tides::ocean(const bncTime& time, const ColumnVector& xyz, + const std::string& station) { + ColumnVector dX(3); dX = 0.0; + if (time.undef()) { + return dX; + } + QString stationQ = station.c_str(); + if (blqMap.find(stationQ) == blqMap.end()) { + return dX; + } + t_blqData* blqSet = blqMap[stationQ]; //printBlqSet(station, blqSet); + + // angular argument: see arg2.f from IERS Conventions software collection + double speed[11] = {1.40519e-4, 1.45444e-4, 1.3788e-4, 1.45842e-4, 7.2921e-5, + 6.7598e-5, 7.2523e-5, 6.4959e-5, 5.3234e-6, 2.6392e-6, 3.982e-7}; + + double angfac[4][11]; + angfac[0][0] = 2.0; + angfac[1][0] =-2.0; + angfac[2][0] = 0.0; + angfac[3][0] = 0.0; + + angfac[0][1] = 0.0; + angfac[1][1] = 0.0; + angfac[2][1] = 0.0; + angfac[3][1] = 0.0; + + angfac[0][2] = 2.0; + angfac[1][2] =-3.0; + angfac[2][2] = 1.0; + angfac[3][2] = 0.0; + + angfac[0][3] = 2.0; + angfac[1][3] = 0.0; + angfac[2][3] = 0.0; + angfac[3][3] = 0.0; + + angfac[0][4] = 1.0; + angfac[1][4] = 0.0; + angfac[2][4] = 0.0; + angfac[3][4] = .25; + + angfac[0][5] = 1.0; + angfac[1][5] =-2.0; + angfac[2][5] = 0.0; + angfac[3][5] =-.25; + + angfac[0][6] =-1.0; + angfac[1][6] = 0.0; + angfac[2][6] = 0.0; + angfac[3][6] =-.25; + + angfac[0][7] = 1.0; + angfac[1][7] =-3.0; + angfac[2][7] = 1.0; + angfac[3][7] =-.25; + + angfac[0][8] = 0.0; + angfac[1][8] = 2.0; + angfac[2][8] = 0.0; + angfac[3][8] = 0.0; + + angfac[0][9] = 0.0; + angfac[1][9] = 1.0; + angfac[2][9] =-1.0; + angfac[3][9] = 0.0; + + angfac[0][10] = 2.0; + angfac[1][10] = 0.0; + angfac[2][10] = 0.0; + angfac[3][10] = 0.0; + + double twopi = 6.283185307179586476925287e0; + double dtr = 0.0174532925199; + + // fractional part of the day in seconds + unsigned int year, month, day; + time.civil_date(year, month, day); + int iyear = year - 2000; + QDateTime datTim = QDateTime::fromString(QString::fromStdString(time.datestr()), Qt::ISODate); + int doy = datTim.date().dayOfYear(); + double fday = time.daysec(); + int icapd = doy + 365 * (iyear - 75) + ((iyear - 73) / 4); + double capt = (icapd * 1.000000035 + 27392.500528) / 36525.0; + + // mean longitude of the sun at the beginning of the day + double h0 = (279.69668e0 + (36000.768930485e0 + 3.03e-4 * capt) * capt) * dtr; + + // mean longitude of moon at the beginning of the day + double s0 = (((1.9e-6 * capt - .001133e0) * capt + 481267.88314137e0) * capt + 270.434358e0) * dtr; + + // mean longitude of lunar perigee at the beginning of the day + double p0 = (((-1.2e-5 * capt - .010325e0) * capt + 4069.0340329577e0) * capt + 334.329653e0) * dtr; + + // tidal angle arguments + double angle[11]; + for (int k = 0; k < 11; ++k) { + angle[k] = speed[k] * fday + + angfac[0][k] * h0 + + angfac[1][k] * s0 + + angfac[2][k] * p0 + + angfac[3][k] * twopi; + angle[k] = fmod(angle[k], twopi); + if (angle[k] < 0.0) { + angle[k] += twopi; + } + } + + // displacement by 11 constituents + ColumnVector rwsSta(3); rwsSta = 0.0; // radial, west, south + for (int rr = 0; rr < 3; rr++) { + for (int cc = 0; cc < 11; cc++) { + rwsSta[rr] += blqSet->amplitudes[rr][cc] * cos((angle[cc] - (blqSet->phases[rr][cc]/RHO_DEG))); + } + } + + // neu2xyz + ColumnVector dneu(3); // neu + dneu[0] = -rwsSta[2]; + dneu[1] = -rwsSta[1]; + dneu[2] = rwsSta[0]; + double recEll[3]; xyz2ell(xyz.data(), recEll) ; + neu2xyz(recEll, dneu.data(), dX.data()); + + return dX; +} + // Print //////////////////////////////////////////////////////////////////////////// -void t_loading::printAll() const { +void t_tides::printAllBlqSets() const { QMapIterator it(blqMap); @@ -341,5 +495,5 @@ t_blqData* blq = it.value(); QString site = it.key(); - cout << site.toStdString().c_str() << "\n"; + cout << site.toStdString().c_str() << "\n===============\n"; for (int rr = 0; rr < 3; rr++) { for (int cc = 0; cc < 11; cc++) { @@ -357,9 +511,27 @@ } +// Print +//////////////////////////////////////////////////////////////////////////// +void t_tides::printBlqSet(const std::string& station, t_blqData* blq) { + cout << station << endl; + for (int rr = 0; rr < 3; rr++) { + for (int cc = 0; cc < 11; cc++) { + cout << blq->amplitudes[rr][cc] << " "; + } + cout << endl; + } + for (int rr = 0; rr < 3; rr++) { + for (int cc = 0; cc < 11; cc++) { + cout << blq->phases[rr][cc] << " "; + } + cout << endl; + } +} + // Constructor /////////////////////////////////////////////////////////////////////////// t_windUp::t_windUp() { for (unsigned ii = 0; ii <= t_prn::MAXPRN; ii++) { - sumWind[ii] = 0.0; + sumWind[ii] = 0.0; lastEtime[ii] = 0.0; } @@ -369,6 +541,6 @@ /////////////////////////////////////////////////////////////////////////// double t_windUp::value(const bncTime& etime, const ColumnVector& rRec, - t_prn prn, const ColumnVector& rSat, bool ssr, - double yaw, const ColumnVector& vSat) { + t_prn prn, const ColumnVector& rSat, bool ssr, + double yaw, const ColumnVector& vSat) { if (etime.mjddec() != lastEtime[prn.toInt()]) { @@ -377,5 +549,5 @@ // -------------------------------------- ColumnVector rho = rRec - rSat; - rho /= rho.norm_Frobenius(); + rho /= rho.NormFrobenius(); // GPS Satellite unit Vectors sz, sy, sx @@ -392,16 +564,25 @@ sHlp = vSat + crossproduct(Omega, rSat); } - sHlp /= sHlp.norm_Frobenius(); - - ColumnVector sz = -rSat / rSat.norm_Frobenius(); + sHlp /= sHlp.NormFrobenius(); + + ColumnVector sz = -rSat / rSat.NormFrobenius(); ColumnVector sy = crossproduct(sz, sHlp); ColumnVector sx = crossproduct(sy, sz); - // Yaw consideration - sx = t_astro::rotZ(yaw) * sx; - + if (ssr) { + // Yaw angle consideration + Matrix SXYZ(3, 3); + SXYZ.Column(1) = sx; + SXYZ.Column(2) = sy; + SXYZ.Column(3) = sz; + SXYZ = DotProduct(t_astro::rotZ(yaw), SXYZ); + sx = SXYZ.Column(1); + sy = SXYZ.Column(2); + sz = SXYZ.Column(3); + } // Effective Dipole of the GPS Satellite Antenna // --------------------------------------------- - ColumnVector dipSat = sx - rho * DotProduct(rho,sx) - crossproduct(rho, sy); + ColumnVector dipSat = sx - rho * DotProduct(rho, sx) + - crossproduct(rho, sy); // Receiver unit Vectors rx, ry @@ -409,5 +590,6 @@ ColumnVector rx(3); ColumnVector ry(3); - double recEll[3]; xyz2ell(rRec.data(), recEll) ; + double recEll[3]; + xyz2ell(rRec.data(), recEll); double neu[3]; @@ -417,24 +599,27 @@ neu2xyz(recEll, neu, rx.data()); - neu[0] = 0.0; + neu[0] = 0.0; neu[1] = -1.0; - neu[2] = 0.0; + neu[2] = 0.0; neu2xyz(recEll, neu, ry.data()); // Effective Dipole of the Receiver Antenna // ---------------------------------------- - ColumnVector dipRec = rx - rho * DotProduct(rho,rx) + crossproduct(rho, ry); + ColumnVector dipRec = rx - rho * DotProduct(rho, rx) + + crossproduct(rho, ry); // Resulting Effect // ---------------- - double alpha = DotProduct(dipSat,dipRec) / - (dipSat.norm_Frobenius() * dipRec.norm_Frobenius()); -/* - if (alpha > 1.0) alpha = 1.0; - if (alpha < -1.0) alpha = -1.0; -*/ + double alpha = DotProduct(dipSat, dipRec) + / (dipSat.NormFrobenius() * dipRec.NormFrobenius()); + + if (alpha > 1.0) + alpha = 1.0; + if (alpha < -1.0) + alpha = -1.0; + double dphi = acos(alpha) / 2.0 / M_PI; // in cycles - if ( DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0 ) { + if (DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0) { dphi = -dphi; } @@ -467,14 +652,19 @@ double height = ell[2]; - double pp = 1013.25 * pow(1.0 - 2.26e-5 * height, 5.225); - double TT = 18.0 - height * 0.0065 + 273.15; - double hh = 50.0 * exp(-6.396e-4 * height); - double ee = hh / 100.0 * exp(-37.2465 + 0.213166*TT - 0.000256908*TT*TT); + double pp = 1013.25 * pow(1.0 - 2.26e-5 * height, 5.225); + double TT = 18.0 - height * 0.0065 + 273.15; + double hh = 50.0 * exp(-6.396e-4 * height); + double ee = hh / 100.0 + * exp(-37.2465 + 0.213166 * TT - 0.000256908 * TT * TT); double h_km = height / 1000.0; - if (h_km < 0.0) h_km = 0.0; - if (h_km > 5.0) h_km = 5.0; - int ii = int(h_km + 1); if (ii > 5) ii = 5; + if (h_km < 0.0) + h_km = 0.0; + if (h_km > 5.0) + h_km = 5.0; + int ii = int(h_km + 1); + if (ii > 5) + ii = 5; double href = ii - 1; @@ -487,9 +677,10 @@ bCor[5] = 0.563; - double BB = bCor[ii-1] + (bCor[ii]-bCor[ii-1]) * (h_km - href); - - double zen = M_PI/2.0 - Ele; - - return (0.002277/cos(zen)) * (pp + ((1255.0/TT)+0.05)*ee - BB*(tan(zen)*tan(zen))); + double BB = bCor[ii - 1] + (bCor[ii] - bCor[ii - 1]) * (h_km - href); + + double zen = M_PI / 2.0 - Ele; + + return (0.002277 / cos(zen)) + * (pp + ((1255.0 / TT) + 0.05) * ee - BB * (tan(zen) * tan(zen))); } @@ -500,9 +691,10 @@ } -t_iono::~t_iono() {} +t_iono::~t_iono() { +} double t_iono::stec(const t_vTec* vTec, double signalPropagationTime, - const ColumnVector& rSat, const bncTime& epochTime, - const ColumnVector& xyzSta) { + const ColumnVector& rSat, const bncTime& epochTime, + const ColumnVector& xyzSta) { // Latitude, longitude, height are defined with respect to a spherical earth model @@ -524,8 +716,8 @@ // ------------------------------------------------------------- ColumnVector rhoV = xyzSat - xyzSta; - double rho = rhoV.norm_Frobenius(); + double rho = rhoV.NormFrobenius(); ColumnVector neu(3); xyz2neu(geocSta.data(), rhoV.data(), neu.data()); - double sphEle = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / rho ); + double sphEle = acos(sqrt(neu[0] * neu[0] + neu[1] * neu[1]) / rho); if (neu[2] < 0) { sphEle *= -1.0; @@ -537,5 +729,6 @@ double stec = 0.0; for (unsigned ii = 0; ii < vTec->_layers.size(); ii++) { - piercePoint(vTec->_layers[ii]._height, epoch, geocSta.data(), sphEle, sphAzi); + piercePoint(vTec->_layers[ii]._height, epoch, geocSta.data(), sphEle, + sphAzi); double vtec = vtecSingleLayerContribution(vTec->_layers[ii]); stec += vtec * sin(sphEle + _psiPP); @@ -547,6 +740,6 @@ double vtec = 0.0; - int N = vTecLayer._C.Nrows()-1; - int M = vTecLayer._C.Ncols()-1; + int N = vTecLayer._C.Nrows() - 1; + int M = vTecLayer._C.Ncols() - 1; double fac; @@ -576,21 +769,29 @@ } -void t_iono::piercePoint(double layerHeight, double epoch, const double* geocSta, +void t_iono::piercePoint(double layerHeight, double epoch, + const double* geocSta, double sphEle, double sphAzi) { double q = (t_CST::rgeoc + geocSta[2]) / (t_CST::rgeoc + layerHeight); - _psiPP = M_PI/2 - sphEle - asin(q * cos(sphEle)); - - _phiPP = asin(sin(geocSta[0]) * cos(_psiPP) + cos(geocSta[0]) * sin(_psiPP) * cos(sphAzi)); - - if (( (geocSta[0]*180.0/M_PI > 0) && ( tan(_psiPP) * cos(sphAzi) > tan(M_PI/2 - geocSta[0])) ) || - ( (geocSta[0]*180.0/M_PI < 0) && (-(tan(_psiPP) * cos(sphAzi)) > tan(M_PI/2 + geocSta[0])) )) { - _lambdaPP = geocSta[1] + M_PI - asin((sin(_psiPP) * sin(sphAzi) / cos(_phiPP))); - } else { - _lambdaPP = geocSta[1] + asin((sin(_psiPP) * sin(sphAzi) / cos(_phiPP))); - } - - _lonS = fmod((_lambdaPP + (epoch - 50400) * M_PI / 43200), 2*M_PI); + _psiPP = M_PI / 2 - sphEle - asin(q * cos(sphEle)); + + _phiPP = asin( + sin(geocSta[0]) * cos(_psiPP) + + cos(geocSta[0]) * sin(_psiPP) * cos(sphAzi)); + + if (((geocSta[0] * 180.0 / M_PI > 0) + && (tan(_psiPP) * cos(sphAzi) > tan(M_PI / 2 - geocSta[0]))) + || + ((geocSta[0] * 180.0 / M_PI < 0) + && (-(tan(_psiPP) * cos(sphAzi)) > tan(M_PI / 2 + geocSta[0])))) { + _lambdaPP = geocSta[1] + M_PI + - asin((sin(_psiPP) * sin(sphAzi) / cos(_phiPP))); + } + else { + _lambdaPP = geocSta[1] + asin((sin(_psiPP) * sin(sphAzi) / cos(_phiPP))); + } + + _lonS = fmod((_lambdaPP + (epoch - 50400) * M_PI / 43200), 2 * M_PI); return; Index: trunk/BNC/src/pppModel.h =================================================================== --- trunk/BNC/src/pppModel.h (revision 8904) +++ trunk/BNC/src/pppModel.h (revision 8905) @@ -5,4 +5,5 @@ #include #include +#include #include "bnctime.h" #include "t_prn.h" @@ -21,8 +22,4 @@ private: - static const double RHO_DEG; - static const double RHO_SEC; - static const double MJD_J2000; - static double GMST(double Mjd_UT1); static Matrix NutMatrix(double Mjd_TT); @@ -32,11 +29,10 @@ class t_tides { public: - t_tides() { - _lastMjd = 0.0; - _rSun = 0.0; - _rMoon = 0.0; - } - ~t_tides() {} - ColumnVector displacement(const bncTime& time, const ColumnVector& xyz); + t_tides(); + ~t_tides(); + ColumnVector earth(const bncTime& time, const ColumnVector& xyz); + ColumnVector ocean(const bncTime& time, const ColumnVector& xyz, const std::string& station); + t_irc readBlqFile(const char* fileName); + void printAllBlqSets() const; private: double _lastMjd; @@ -45,14 +41,5 @@ double _rSun; double _rMoon; -}; -class t_loading { - public: - t_loading(const QString& fileName); - ~t_loading(); - t_irc readFile(const QString& fileName); - void printAll() const; - - private: class t_blqData { public: @@ -63,5 +50,5 @@ t_blqData* newBlqData; QMap blqMap; - + void printBlqSet(const std::string& station, t_blqData* blq); }; @@ -98,6 +85,4 @@ double _lambdaPP; double _lonS; - - }; Index: trunk/BNC/src/pppOptions.cpp =================================================================== --- trunk/BNC/src/pppOptions.cpp (revision 8904) +++ trunk/BNC/src/pppOptions.cpp (revision 8905) @@ -35,5 +35,5 @@ * Created: 29-Jul-2014 * - * Changes: + * Changes: * * -----------------------------------------------------------------------*/ @@ -59,5 +59,5 @@ } -// +// ////////////////////////////////////////////////////////////////////////////// const std::vector& t_pppOptions::LCs(char system) const { @@ -77,5 +77,5 @@ } -// +// ////////////////////////////////////////////////////////////////////////////// bool t_pppOptions::useOrbClkCorr() const { @@ -99,8 +99,8 @@ } -// +// ///////////////////////////////////////////////////////////////////////////// vector t_pppOptions::ambLCs(char system) const { - + const vector& allLCs = LCs(system); vector phaseLCs; @@ -146,2 +146,4 @@ return answ; } + + Index: trunk/BNC/src/pppOptions.h =================================================================== --- trunk/BNC/src/pppOptions.h (revision 8904) +++ trunk/BNC/src/pppOptions.h (revision 8905) @@ -11,4 +11,5 @@ class t_pppOptions { public: + enum e_type {IF, UncombPPP, PPPRTK, DCMcodeBias, DCMphaseBias}; t_pppOptions(); ~t_pppOptions(); @@ -18,4 +19,5 @@ std::vector ambLCs(char system) const; std::vector codeLCs(char system) const; + std::vector ionoLCs(char system) const; bool useSystem(char system) const {return LCs(system).size() > 0;} bool useOrbClkCorr() const; @@ -25,4 +27,11 @@ } + e_type _obsModelType; + QStringList _obsmodelTypeStr = QStringList() + << "IF PPP" + << "Uncombined PPP" + << "PPP-RTK" + << "DCM with Code Biases" + << "DCM with Phase Biases"; bool _realTime; std::string _crdFile; @@ -43,4 +52,6 @@ double _maxResC1; double _maxResL1; + double _sigmaGIMdiff; + double _maxResGIMdiff; bool _eleWgtCode; bool _eleWgtPhase; @@ -52,4 +63,10 @@ double _aprSigTrp; double _noiseTrp; + double _aprSigIon; + double _noiseIon; + double _aprSigCodeBias; + double _noiseCodeBias; + double _aprSigPhaseBias; + double _noisePhaseBias; int _nmeaPort; double _aprSigAmb; @@ -59,4 +76,6 @@ std::vector _LCsGalileo; std::vector _LCsBDS; + bool _pseudoObsIono; + bool _refSatRequired; }; Index: trunk/BNC/src/pppWidgets.cpp =================================================================== --- trunk/BNC/src/pppWidgets.cpp (revision 8904) +++ trunk/BNC/src/pppWidgets.cpp (revision 8905) @@ -82,4 +82,6 @@ _lcGalileo = new QComboBox(); _lcGalileo ->setObjectName("PPP/lcGalileo"); _widgets << _lcGalileo; _lcBDS = new QComboBox(); _lcBDS ->setObjectName("PPP/lcBDS"); _widgets << _lcBDS; + _modelObs = new QComboBox(); _modelObs ->setObjectName("PPP/modelObs"); _widgets << _modelObs; + _pseudoObs = new QComboBox(); _pseudoObs ->setObjectName("PPP/pseudoObs"); _widgets << _pseudoObs; _sigmaC1 = new QLineEdit(); _sigmaC1 ->setObjectName("PPP/sigmaC1"); _widgets << _sigmaC1; _sigmaL1 = new QLineEdit(); _sigmaL1 ->setObjectName("PPP/sigmaL1"); _widgets << _sigmaL1; @@ -111,13 +113,8 @@ _dataSource->addItems(QString(",Real-Time Streams,RINEX Files").split(",")); connect(_dataSource, SIGNAL(currentIndexChanged(const QString&)), this, SLOT(slotEnableWidgets())); - - connect(_snxtroPath, SIGNAL(textChanged(const QString &)), - this, SLOT(slotPPPTextChanged())); - - connect(_snxtroAc, SIGNAL(textChanged(const QString &)), - this, SLOT(slotPPPTextChanged())); - - connect(_snxtroSol, SIGNAL(textChanged(const QString &)), - this, SLOT(slotPPPTextChanged())); + connect(_modelObs, SIGNAL(currentIndexChanged(const QString&)), this, SLOT(slotEnableWidgets())); + connect(_snxtroPath, SIGNAL(textChanged(const QString &)), this, SLOT(slotPPPTextChanged())); + connect(_snxtroAc, SIGNAL(textChanged(const QString &)), this, SLOT(slotPPPTextChanged())); + connect(_snxtroSol, SIGNAL(textChanged(const QString &)), this, SLOT(slotPPPTextChanged())); slotEnableWidgets(); @@ -127,5 +124,5 @@ _lcGPS->addItems(QString("P3,P3&L3").split(",")); #else - _lcGPS->addItems(QString("no,Pi,Li,Pi&Li,P3,L3,P3&L3").split(",")); + _lcGPS->addItems(QString("no,Pi,Li,Pi&Li").split(",")); #endif @@ -134,5 +131,5 @@ _lcGLONASS->addItems(QString("no,P3,L3,P3&L3").split(",")); #else - _lcGLONASS->addItems(QString("no,Pi,Li,Pi&Li,P3,L3,P3&L3").split(",")); + _lcGLONASS->addItems(QString("no,Pi,Li,Pi&Li").split(",")); #endif @@ -141,5 +138,5 @@ _lcGalileo->addItems(QString("no,P3,L3,P3&L3").split(",")); #else - _lcGalileo->addItems(QString("no,Pi,Li,Pi&Li,P3,L3,P3&L3").split(",")); + _lcGalileo->addItems(QString("no,Pi,Li,Pi&Li").split(",")); #endif @@ -148,5 +145,15 @@ _lcBDS->addItems(QString("no,P3,L3,P3&L3").split(",")); #else - _lcBDS->addItems(QString("no,Pi,Li,Pi&Li,P3,L3,P3&L3").split(",")); + _lcBDS->addItems(QString("no,Pi,Li,Pi&Li").split(",")); +#endif + + _modelObs->setEditable(false); + _pseudoObs->setEditable(false); +#ifdef USE_PPP_SSR_I + _modelObs->addItems(QString("Ionosphere-free PPP").split(",")); + _pseudoObs->addItems(QString("no").split(",")); +#else + _modelObs->addItems(QString("Ionosphere-free PPP,Uncombined PPP,PPP-RTK,DCM with Code Biases,DCM with Phase Biases").split(",")); + _pseudoObs->addItems(QString("no,Ionosphere").split(",")); #endif @@ -247,4 +254,6 @@ delete _lcGalileo; delete _lcBDS; + delete _modelObs; + delete _pseudoObs; delete _sigmaC1; delete _sigmaL1; @@ -297,4 +306,12 @@ _lcBDS->setCurrentIndex(ii); } + ii = _modelObs->findText(settings.value(_modelObs->objectName()).toString()); + if (ii != -1) { + _modelObs->setCurrentIndex(ii); + } + ii = _pseudoObs->findText(settings.value(_pseudoObs->objectName()).toString()); + if (ii != -1) { + _pseudoObs->setCurrentIndex(ii); + } ii = _snxtroIntr->findText(settings.value(_snxtroIntr->objectName()).toString()); if (ii != -1) { @@ -433,4 +450,6 @@ settings.setValue(_lcGalileo ->objectName(), _lcGalileo ->currentText()); settings.setValue(_lcBDS ->objectName(), _lcBDS ->currentText()); + settings.setValue(_modelObs ->objectName(), _modelObs ->currentText()); + settings.setValue(_pseudoObs ->objectName(), _pseudoObs ->currentText()); settings.setValue(_sigmaC1 ->objectName(), _sigmaC1 ->text()); settings.setValue(_sigmaL1 ->objectName(), _sigmaL1 ->text()); @@ -477,4 +496,12 @@ bool realTime = _dataSource->currentText() == "Real-Time Streams"; bool rinexFiles = _dataSource->currentText() == "RINEX Files"; + bool enablePseudoObs; + if (_modelObs->currentText() == "PPP-RTK" || + _modelObs->currentText() == "Ionosphere-free PPP") { + enablePseudoObs = false; + } + else { + enablePseudoObs = true; + } QListIterator it(_widgets); @@ -490,7 +517,6 @@ } else if (rinexFiles) { - _corrMount->setEnabled(false); -// _plotCoordinates->setEnabled(false); -// _audioResponse->setEnabled(false); + _corrMount ->setEnabled(false); + _audioResponse->setEnabled(false); } @@ -506,4 +532,10 @@ _snxtroAc ->setEnabled(false); _snxtroSol ->setEnabled(false); + } + + if (enablePseudoObs) { + _pseudoObs->setEnabled(true); + } else { + _pseudoObs->setEnabled(false); } @@ -592,3 +624,5 @@ } } -} + + +} Index: trunk/BNC/src/pppWidgets.h =================================================================== --- trunk/BNC/src/pppWidgets.h (revision 8904) +++ trunk/BNC/src/pppWidgets.h (revision 8905) @@ -66,4 +66,6 @@ QComboBox* _lcGalileo; QComboBox* _lcBDS; + QComboBox* _modelObs; + QComboBox* _pseudoObs; QLineEdit* _sigmaC1; QLineEdit* _sigmaL1; Index: trunk/BNC/src/src.pri =================================================================== --- trunk/BNC/src/src.pri (revision 8904) +++ trunk/BNC/src/src.pri (revision 8905) @@ -122,8 +122,8 @@ HEADERS += PPP/pppClient.h PPP/pppObsPool.h PPP/pppEphPool.h \ PPP/pppStation.h PPP/pppFilter.h PPP/pppParlist.h \ - PPP/pppSatObs.h + PPP/pppSatObs.h PPP/pppRefSat.h SOURCES += PPP/pppClient.cpp PPP/pppObsPool.cpp PPP/pppEphPool.cpp \ PPP/pppStation.cpp PPP/pppFilter.cpp PPP/pppParlist.cpp \ - PPP/pppSatObs.cpp + PPP/pppSatObs.cpp PPP/pppRefSat.cpp } else {