/* ------------------------------------------------------------------------- * BKG NTRIP Client * ------------------------------------------------------------------------- * * Class: t_pppClient * * Purpose: PPP Client processing starts here * * Author: L. Mervart * * Created: 29-Jul-2014 * * Changes: * * -----------------------------------------------------------------------*/ #include #include #include #include #include #include #include "pppClient.h" #include "pppEphPool.h" #include "pppObsPool.h" #include "bncconst.h" #include "bncutils.h" #include "pppStation.h" #include "bncantex.h" #include "pppFilter.h" using namespace BNC_PPP; using namespace std; // Global variable holding thread-specific pointers ////////////////////////////////////////////////////////////////////////////// QThreadStorage CLIENTS; // Static function returning thread-specific pointer ////////////////////////////////////////////////////////////////////////////// t_pppClient* t_pppClient::instance() { return CLIENTS.localData(); } // Constructor ////////////////////////////////////////////////////////////////////////////// t_pppClient::t_pppClient(const t_pppOptions* opt) { _output = 0; _opt = new t_pppOptions(*opt); _log = new ostringstream(); _ephPool = new t_pppEphPool(); _obsPool = new t_pppObsPool(); _staRover = new t_pppStation(); _filter = new t_pppFilter(_obsPool); _tides = new t_tides(); _antex = 0; if (!_opt->_antexFileName.empty()) { _antex = new bncAntex(_opt->_antexFileName.c_str()); } if (!_opt->_blqFileName.empty()) { if (_tides->readBlqFile(_opt->_blqFileName.c_str()) == success) { //_tides->printAllBlqSets(); } } if (_opt->_refSatRequired) { for (unsigned iSys = 0; iSys < _opt->systems().size(); iSys++) { char sys = _opt->systems()[iSys]; _refSatMap[sys] = new t_pppRefSat(); } } _offGR = 0.0; _offGE = 0.0; _offGC = 0.0; CLIENTS.setLocalData(this); // CLIENTS takes ownership over "this" } // Destructor ////////////////////////////////////////////////////////////////////////////// t_pppClient::~t_pppClient() { delete _log; delete _opt; delete _filter; delete _ephPool; delete _obsPool; delete _staRover; if (_antex) { delete _antex; } delete _tides; clearObs(); QMapIterator it(_refSatMap); while (it.hasNext()) { it.next(); delete it.value(); } _refSatMap.clear(); } // ////////////////////////////////////////////////////////////////////////////// void t_pppClient::putEphemeris(const t_eph* eph) { const t_ephGPS* ephGPS = dynamic_cast(eph); const t_ephGlo* ephGlo = dynamic_cast(eph); const t_ephGal* ephGal = dynamic_cast(eph); const t_ephBDS* ephBDS = dynamic_cast(eph); if (ephGPS) { _ephPool->putEphemeris(new t_ephGPS(*ephGPS)); } else if (ephGlo) { _ephPool->putEphemeris(new t_ephGlo(*ephGlo)); } else if (ephGal) { _ephPool->putEphemeris(new t_ephGal(*ephGal)); } else if (ephBDS) { _ephPool->putEphemeris(new t_ephBDS(*ephBDS)); } } // ////////////////////////////////////////////////////////////////////////////// void t_pppClient::putTec(const t_vTec* vTec) { _obsPool->putTec(new t_vTec(*vTec)); } // ////////////////////////////////////////////////////////////////////////////// void t_pppClient::putOrbCorrections(const vector& corr) { for (unsigned ii = 0; ii < corr.size(); ii++) { _ephPool->putOrbCorrection(new t_orbCorr(*corr[ii])); } } // ////////////////////////////////////////////////////////////////////////////// void t_pppClient::putClkCorrections(const vector& corr) { for (unsigned ii = 0; ii < corr.size(); ii++) { _ephPool->putClkCorrection(new t_clkCorr(*corr[ii])); } } // ////////////////////////////////////////////////////////////////////////////// void t_pppClient::putCodeBiases(const vector& biases) { for (unsigned ii = 0; ii < biases.size(); ii++) { _obsPool->putCodeBias(new t_satCodeBias(*biases[ii])); } } // ////////////////////////////////////////////////////////////////////////////// void t_pppClient::putPhaseBiases(const vector& biases) { for (unsigned ii = 0; ii < biases.size(); ii++) { _obsPool->putPhaseBias(new t_satPhaseBias(*biases[ii])); } } // ////////////////////////////////////////////////////////////////////////////// t_irc t_pppClient::prepareObs(const vector& satObs, vector& obsVector, bncTime& epoTime) { // Default // ------- epoTime.reset(); clearObs(); // Create vector of valid observations // ----------------------------------- for (unsigned ii = 0; ii < satObs.size(); ii++) { char sys = satObs[ii]->_prn.system(); if (_opt->useSystem(sys)) { t_pppSatObs* pppSatObs = new t_pppSatObs(*satObs[ii]); if (pppSatObs->isValid()) { obsVector.push_back(pppSatObs); } else { delete pppSatObs; } } } // Check whether data are synchronized, compute epoTime // ---------------------------------------------------- const double MAXSYNC = 0.05; // synchronization limit double meanDt = 0.0; for (unsigned ii = 0; ii < obsVector.size(); ii++) { const t_pppSatObs* satObs = obsVector.at(ii); if (epoTime.undef()) { epoTime = satObs->time(); } else { double dt = satObs->time() - epoTime; if (fabs(dt) > MAXSYNC) { LOG << "t_pppClient::prepareObs asynchronous observations" << endl; return failure; } meanDt += dt; } } if (obsVector.size() > 0) { epoTime += meanDt / obsVector.size(); } return success; } // ////////////////////////////////////////////////////////////////////////////// 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 sys = satObs->prn().system(); t_pppRefSat* refSat = _refSatMap[sys]; double stecRef = refSat->stecValue(); if (stecRef && !satObs->isReference()) { pseudoObsIono = true; satObs->setPseudoObsIono(t_frequency::G1, stecRef); } it++; } } /* vector::iterator it = obsVector.begin(); while (it != obsVector.end()) { t_pppSatObs* satObs = *it; satObs->printObsMinusComputed(); it++; } */ return pseudoObsIono; } // ////////////////////////////////////////////////////////////////////////////// void t_pppClient::useObsWithCodeBiasesOnly(std::vector& obsVector) { vector::iterator it = obsVector.begin(); while (it != obsVector.end()) { t_pppSatObs* pppSatObs = *it; char sys = pppSatObs->prn().system(); bool codeBiasesAvailable = false; t_frequency::type fType1 = t_lc::toFreq(sys,t_lc::c1); t_frequency::type fType2 = t_lc::toFreq(sys,t_lc::c2); if (pppSatObs->getCodeBias(fType1) && pppSatObs->getCodeBias(fType2)) { codeBiasesAvailable = true; } if (codeBiasesAvailable) { ++it; } else { it = obsVector.erase(it); delete pppSatObs; } } } // Compute the Bancroft position, check for blunders ////////////////////////////////////////////////////////////////////////////// t_irc t_pppClient::cmpBancroft(const bncTime& epoTime, vector& obsVector, ColumnVector& xyzc, bool print) { t_lc::type tLC = t_lc::dummy; while (true) { Matrix BB(obsVector.size(), 4); int iObs = -1; for (unsigned ii = 0; ii < obsVector.size(); ii++) { const t_pppSatObs* satObs = obsVector.at(ii); if (tLC == t_lc::dummy) { if (satObs->isValid(t_lc::cIF)) { tLC = t_lc::cIF; } else if (satObs->isValid(t_lc::c1)) { tLC = t_lc::c1; } else if (satObs->isValid(t_lc::c2)) { tLC = t_lc::c2; } } if ( satObs->isValid(tLC) && (!satObs->modelSet() || satObs->eleSat() >= _opt->_minEle) ) { ++iObs; BB[iObs][0] = satObs->xc()[0]; BB[iObs][1] = satObs->xc()[1]; BB[iObs][2] = satObs->xc()[2]; BB[iObs][3] = satObs->obsValue(tLC) - satObs->cmpValueForBanc(tLC); } } if (iObs + 1 < _opt->_minObs) { LOG << "t_pppClient::cmpBancroft not enough observations: " << iObs + 1; return failure; } BB = BB.Rows(1,iObs+1); if (bancroft(BB, xyzc) != success) { return failure; } xyzc[3] /= t_CST::c; // Check Blunders // -------------- const double BLUNDER = 100.0; double maxRes = 0.0; unsigned maxResIndex = 0; for (unsigned ii = 0; ii < obsVector.size(); ii++) { const t_pppSatObs* satObs = obsVector.at(ii); if (satObs->isValid() && (!satObs->modelSet() || satObs->eleSat() >= _opt->_minEle) ) { ColumnVector rr = satObs->xc().Rows(1,3) - xyzc.Rows(1,3); double res = rr.NormFrobenius() - satObs->obsValue(tLC) - (satObs->xc()[3] - xyzc[3]) * t_CST::c; if (fabs(res) > maxRes) { maxRes = fabs(res); maxResIndex = ii; } } } if (maxRes < BLUNDER) { if (print && _numEpoProcessing == 1) { LOG.setf(ios::fixed); LOG << "\nPPP of Epoch "; if (!_epoTimeRover.undef()) LOG << string(_epoTimeRover); LOG << "\n---------------------------------------------------------------\n"; LOG << string(epoTime) << " BANCROFT:" << ' ' << setw(14) << setprecision(3) << xyzc[0] << ' ' << setw(14) << setprecision(3) << xyzc[1] << ' ' << setw(14) << setprecision(3) << xyzc[2] << ' ' << setw(14) << setprecision(3) << xyzc[3] * t_CST::c << endl << endl; } break; } else { t_pppSatObs* satObs = obsVector.at(maxResIndex); LOG << "t_pppClient::cmpBancroft outlier " << satObs->prn().toString() << " " << maxRes << endl; delete satObs; obsVector.erase(obsVector.begin() + maxResIndex); } } return success; } // ////////////////////////////////////////////////////////////////////////////// void t_pppClient::initOutput(t_output* output) { _output = output; _output->_numSat = 0; _output->_hDop = 0.0; _output->_error = false; } // ////////////////////////////////////////////////////////////////////////////// void t_pppClient::clearObs() { for (unsigned ii = 0; ii < _obsRover.size(); ii++) { delete _obsRover.at(ii); } _obsRover.clear(); } // ////////////////////////////////////////////////////////////////////////////// void t_pppClient::finish(t_irc irc, int ind) { #ifdef BNC_DEBUG_PPP LOG << "t_pppClient::finish(" << ind << "): " << irc << endl; #endif clearObs(); _output->_epoTime = _epoTimeRover; if (irc == success) { _output->_xyzRover[0] = _staRover->xyzApr()[0] + _filter->x()[0]; _output->_xyzRover[1] = _staRover->xyzApr()[1] + _filter->x()[1]; _output->_xyzRover[2] = _staRover->xyzApr()[2] + _filter->x()[2]; xyz2neu(_staRover->ellApr().data(), _filter->x().data(), _output->_neu); copy(&_filter->Q().data()[0], &_filter->Q().data()[6], _output->_covMatrix); _output->_trp0 = t_tropo::delay_saast(_staRover->xyzApr(), M_PI/2.0); _output->_trp = _filter->trp(); _output->_trpStdev = _filter->trpStdev(); _output->_numSat = _filter->numSat(); _output->_hDop = _filter->HDOP(); _output->_error = false; } else { _output->_error = true; if (OPT->_obsModelType == OPT->DCMcodeBias || OPT->_obsModelType == OPT->DCMphaseBias) { if (ind == 1 || ind > 7) { reset(); } } } _output->_log = _log->str(); delete _log; _log = new ostringstream(); return; } // ////////////////////////////////////////////////////////////////////////////// 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); } else { station->setXyzApr(xyzc.Rows(1,3)); } station->setNeuEcc(_opt->_neuEccRover); // Receiver Clock // -------------- station->setDClk(xyzc[3]); // Tides // ----- station->setTideDsplEarth(_tides->earth(time, station->xyzApr())); station->setTideDsplOcean(_tides->ocean(time, station->xyzApr(), station->name())); // Observation model // ----------------- vector::iterator it = obsVector.begin(); while (it != obsVector.end()) { t_pppSatObs* satObs = *it; t_irc modelSetup; if (satObs->isValid()) { modelSetup = satObs->cmpModel(station); } if (satObs->isValid() && satObs->eleSat() >= _opt->_minEle && modelSetup == success) { ++it; } else { delete satObs; it = obsVector.erase(it); } } return success; } // ////////////////////////////////////////////////////////////////////////////// void t_pppClient::processEpoch(const vector& satObs, t_output* output) { try { initOutput(output); bool epochReProcessing = false; _numEpoProcessing = 0; _historicalRefSats.clear(); do { _numEpoProcessing++; #ifdef BNC_DEBUG_PPP LOG << "_numEpoProcessing " << _numEpoProcessing << endl; #endif if (_obsPool->refSatChanged()) { if(_filter->datumTransformation(_refSatMap) != success) { LOG << "t_pppFilter::datumTransformation() != success" << endl; return finish(failure,1); } else { LOG << "t_pppFilter::datumTransformation() == success" << endl; _filter->rememberState(1); } } // Prepare Observations of the Rover // --------------------------------- if (prepareObs(satObs, _obsRover, _epoTimeRover) != success) { return finish(failure,2); } 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,3); } if (cmpModel(_staRover, xyzc, _obsRover) != success) { return finish(failure,4); } } if (_opt->_refSatRequired) { if (handleRefSatellites(_obsRover) != success) { return finish(failure,6); } if (_obsPool->refSatChanged()) { LOG << "t_pppFilter: refSatChanged()" << endl; epochReProcessing = true; continue; } } // use observations only if satellite code biases are available // ------------------------------------------------------------ if (!_opt->_corrMount.empty() && (OPT->_obsModelType == OPT->DCMcodeBias || OPT->_obsModelType == OPT->DCMphaseBias)) { useObsWithCodeBiasesOnly(_obsRover); } if (int(_obsRover.size()) < _opt->_minObs) { LOG << "t_pppClient::processEpoch not enough observations" << endl; return finish(failure,5); } // Prepare Pseudo Observations of the Rover // ---------------------------------------- _pseudoObsIono = preparePseudoObs(_obsRover); // Store last epoch of data // ------------------------ _obsPool->putEpoch(_epoTimeRover, _obsRover, _pseudoObsIono, _refSatMap); // Process Epoch in Filter // ----------------------- if (_filter->processEpoch() != success) { LOG << "filter->processEpoch() != success" << endl; return finish(failure,7); } // Epoch re-processing required? // ----------------------------- if (_obsPool->refSatChangeRequired() && //SLIP _opt->_obsModelType != t_pppOptions::UncombPPP) { LOG << "pppClient: _obsPool->refSatChangeRequired() " << endl; epochReProcessing = true; _obsPool->deleteLastEpoch(); _filter->restoreState(0); setHistoricalRefSats(); } else { epochReProcessing = false; if (OPT->_obsModelType == OPT->DCMcodeBias || OPT->_obsModelType == OPT->DCMphaseBias) { _filter->rememberState(0); } } } while (epochReProcessing); } catch (Exception& exc) { LOG << exc.what() << endl; return finish(failure,8); } catch (t_except msg) { LOG << msg.what() << endl; return finish(failure,9); } catch (const char* msg) { LOG << msg << endl; return finish(failure,10); } catch (...) { LOG << "unknown exception" << endl; return finish(failure,11); } return finish(success,12); } // //////////////////////////////////////////////////////////////////////////// double lorentz(const ColumnVector& aa, const ColumnVector& bb) { return aa[0]*bb[0] + aa[1]*bb[1] + aa[2]*bb[2] - aa[3]*bb[3]; } // //////////////////////////////////////////////////////////////////////////// t_irc t_pppClient::bancroft(const Matrix& BBpass, ColumnVector& pos) { if (pos.Nrows() != 4) { pos.ReSize(4); } pos = 0.0; for (int iter = 1; iter <= 2; iter++) { Matrix BB = BBpass; int mm = BB.Nrows(); for (int ii = 1; ii <= mm; ii++) { double xx = BB(ii,1); double yy = BB(ii,2); double traveltime = 0.072; if (iter > 1) { double zz = BB(ii,3); double rho = sqrt( (xx-pos(1)) * (xx-pos(1)) + (yy-pos(2)) * (yy-pos(2)) + (zz-pos(3)) * (zz-pos(3)) ); traveltime = rho / t_CST::c; } double angle = traveltime * t_CST::omega; double cosa = cos(angle); double sina = sin(angle); BB(ii,1) = cosa * xx + sina * yy; BB(ii,2) = -sina * xx + cosa * yy; } Matrix BBB; if (mm > 4) { SymmetricMatrix hlp; hlp << BB.t() * BB; BBB = hlp.i() * BB.t(); } else { BBB = BB.i(); } ColumnVector ee(mm); ee = 1.0; ColumnVector alpha(mm); alpha = 0.0; for (int ii = 1; ii <= mm; ii++) { alpha(ii) = lorentz(BB.Row(ii).t(),BB.Row(ii).t())/2.0; } ColumnVector BBBe = BBB * ee; ColumnVector BBBalpha = BBB * alpha; double aa = lorentz(BBBe, BBBe); double bb = lorentz(BBBe, BBBalpha)-1; double cc = lorentz(BBBalpha, BBBalpha); double root = sqrt(bb*bb-aa*cc); Matrix hlpPos(4,2); hlpPos.Column(1) = (-bb-root)/aa * BBBe + BBBalpha; hlpPos.Column(2) = (-bb+root)/aa * BBBe + BBBalpha; ColumnVector omc(2); for (int pp = 1; pp <= 2; pp++) { hlpPos(4,pp) = -hlpPos(4,pp); omc(pp) = BB(1,4) - sqrt( (BB(1,1)-hlpPos(1,pp)) * (BB(1,1)-hlpPos(1,pp)) + (BB(1,2)-hlpPos(2,pp)) * (BB(1,2)-hlpPos(2,pp)) + (BB(1,3)-hlpPos(3,pp)) * (BB(1,3)-hlpPos(3,pp)) ) - hlpPos(4,pp); } if ( fabs(omc(1)) > fabs(omc(2)) ) { pos = hlpPos.Column(2); } else { pos = hlpPos.Column(1); } } if (pos.size() != 4 || std::isnan(pos(1)) || std::isnan(pos(2)) || std::isnan(pos(3)) || std::isnan(pos(4))) { return failure; } return success; } // ////////////////////////////////////////////////////////////////////////////// void t_pppClient::setRefSatellites(std::vector& obsVector) { // reference satellite definition per system for (unsigned iSys = 0; iSys < _opt->systems().size(); iSys++) { t_irc refSatReDefinition = success; char sys = _opt->systems()[iSys]; bool refSatDefined = false; t_pppRefSat* refSat = _refSatMap[sys]; for (unsigned ii = 0; ii < obsVector.size(); ii++) { t_pppSatObs* satObs = obsVector.at(ii); if (satObs->eleSat() < _opt->_minEle) { continue; } // reference satellite is unchanged // ================================ if ( !_obsPool->refSatChangeRequired(sys) && refSat->prn() == satObs->prn()) { refSatDefined = true; obsVector[ii]->setAsReference(); } // reference satellite has changed // =============================== else if ( _obsPool->refSatChangeRequired(sys) && refSat->prn() != satObs->prn()) { if (satObs->prn().system() == sys) { if (!_historicalRefSats[sys].contains(satObs->prn())) { refSatDefined = true; obsVector[ii]->setAsReference(); refSat->setPrn(satObs->prn()); } else if ( _historicalRefSats[sys].contains(satObs->prn())) { refSatReDefinition = failure; } } } if (refSatDefined) { if (OPT->_pseudoObsIono) { refSat->setStecValue(satObs->getIonoCodeDelay(t_frequency::G1)); } break; } } // all available satellites were already tried to use if ((!refSatDefined) && (refSatReDefinition == failure)) { refSat->setPrn(t_prn(sys, 99)); _obsPool->setRefSatChangeRequired(sys, false); return; } // reference satellite has to be initialized // ========================================= if (!refSatDefined) { for (unsigned ii = 0; ii < obsVector.size(); ii++) { t_pppSatObs* satObs = obsVector.at(ii); if (satObs->eleSat() < _opt->_minEle) { continue; } if (satObs->prn().system() == sys && !_historicalRefSats[sys].contains(satObs->prn())) { obsVector[ii]->setAsReference(); refSat->setPrn(satObs->prn()); if (OPT->_pseudoObsIono) { refSat->setStecValue(satObs->getIonoCodeDelay(t_frequency::G1)); } refSatDefined = true; break; } } } // no observations for that system // =============================== if (!refSatDefined) { refSat->setPrn(t_prn()); } _obsPool->setRefSatChangeRequired(sys, false); // done or impossible } return; } // ////////////////////////////////////////////////////////////////////////////// t_irc t_pppClient::handleRefSatellites(std::vector& obsVector) { // set refSats in current epoch // ============================ setRefSatellites(obsVector); // current epoch LOG.setf(ios::fixed); t_pppObsPool::t_epoch* epoch = _obsPool->lastEpoch(); const QMap& refSatMapLastEpoch = epoch->refSatMap(); QMapIterator it(_refSatMap); while (it.hasNext()) { it.next(); char sys = it.key(); t_prn prn = it.value()->prn(); t_prn refSatLastEpochPrn = t_prn(); if (epoch) { refSatLastEpochPrn = refSatMapLastEpoch[sys]->prn(); } if (prn.number() == 0) { // no obs for that system available continue; } else if (prn.number() == 99) { // refSat re-definition not possible LOG << " => refSat re-definition impossible: " << sys << endl; return failure; } QString str; if (prn == refSatLastEpochPrn || refSatLastEpochPrn == t_prn() ) { _obsPool->setRefSatChanged(sys, false); str = " SET "; } else { _obsPool->setRefSatChanged(sys, true); str = " RESET "; } LOG << string(_epoTimeRover) << str.toStdString() << " REFSAT " << sys << ": " << prn.toString() << endl; } return success; } void t_pppClient::setHistoricalRefSats() { QMapIterator it(_refSatMap); while (it.hasNext()) { it.next(); t_prn prn = it.value()->prn(); char sys = prn.system(); if (_obsPool->refSatChangeRequired(sys)) { _historicalRefSats[sys].append(prn); } } } // ////////////////////////////////////////////////////////////////////////////// void t_pppClient::reset() {LOG << "t_pppClient::reset()" << endl; // to delete old orbit and clock corrections delete _ephPool; _ephPool = new t_pppEphPool(); // to delete old code biases delete _obsPool; _obsPool = new t_pppObsPool(); // to delete all parameters delete _filter; _filter = new t_pppFilter(_obsPool); }