Index: /trunk/BNC/src/PPP/pppClient.cpp =================================================================== --- /trunk/BNC/src/PPP/pppClient.cpp (revision 5739) +++ /trunk/BNC/src/PPP/pppClient.cpp (revision 5740) @@ -1,2 +1,43 @@ + +// Part of BNC, a utility for retrieving decoding and +// converting GNSS data streams from NTRIP broadcasters. +// +// Copyright (C) 2007 +// German Federal Agency for Cartography and Geodesy (BKG) +// http://www.bkg.bund.de +// Czech Technical University Prague, Department of Geodesy +// http://www.fsv.cvut.cz +// +// Email: euref-ip@bkg.bund.de +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License +// as published by the Free Software Foundation, version 2. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +/* ------------------------------------------------------------------------- + * BKG NTRIP Client + * ------------------------------------------------------------------------- + * + * Class: t_pppClient + * + * Purpose: PPP Client processing starts here + * + * Author: L. Mervart + * + * Created: 29-Jul-2014 + * + * Changes: + * + * -----------------------------------------------------------------------*/ + #include @@ -10,9 +51,9 @@ #include "obspool.h" #include "satbias.h" -#include "genconst.h" -#include "utils.h" +#include "bncconst.h" +#include "bncutils.h" #include "station.h" -#include "tides.h" -#include "antex.h" +#include "bnctides.h" +#include "bncantex.h" #include "filter.h" @@ -33,6 +74,5 @@ _ephPool = new t_ephPool(); _obsPool = new t_obsPool(); - _staRover = new t_station(e_rover); - _staBase = new t_station(e_base); + _staRover = new t_station(); _filter = new t_filter(); } @@ -46,5 +86,4 @@ delete _obsPool; delete _staRover; - delete _staBase; delete _antex; delete _filter; @@ -54,14 +93,12 @@ // ////////////////////////////////////////////////////////////////////////////// -void t_pppClient::setOptions(const t_pppOpt* opt) { +void t_pppClient::setOptions(const t_options* opt) { delete _opt; - _opt = new t_options(opt); + _opt = new t_options(*opt); delete _antex; _antex = 0; - if (!_opt->antexFileName().empty()) { - _antex = new t_antex(_opt->antexFileName()); - } - - _opt->print(); + if (!_opt->_antexFile.empty()) { + _antex = new bncAntex(_opt->_antexFile.c_str()); + } } @@ -80,5 +117,5 @@ // ////////////////////////////////////////////////////////////////////////////// -void t_pppClient::putOrbCorrections(int numCorr, const t_pppOrbCorr* corr) { +void t_pppClient::putOrbCorrections(int numCorr, const t_orbCorr* corr) { for (int ii = 0; ii < numCorr; ii++) { _ephPool->putOrbCorrection(new t_orbCorr(corr[ii])); @@ -88,5 +125,5 @@ // ////////////////////////////////////////////////////////////////////////////// -void t_pppClient::putClkCorrections(int numCorr, const t_pppClkCorr* corr) { +void t_pppClient::putClkCorrections(int numCorr, const t_clkCorr* corr) { for (int ii = 0; ii < numCorr; ii++) { _ephPool->putClkCorrection(new t_clkCorr(corr[ii])); @@ -96,5 +133,5 @@ // ////////////////////////////////////////////////////////////////////////////// -void t_pppClient::putBiases(int numBiases, const t_pppSatBiases* biases) { +void t_pppClient::putBiases(int numBiases, const t_satBiases* biases) { for (int ii = 0; ii < numBiases; ii++) { _obsPool->putBiases(new t_satBias(biases[ii])); @@ -104,6 +141,6 @@ // ////////////////////////////////////////////////////////////////////////////// -t_irc::irc t_pppClient::prepareObs(int numSat, const t_pppSatObs* satObs, - vector& obsVector, t_time& epoTime) { +t_irc t_pppClient::prepareObs(int numSat, const t_pppSatObs* satObs, + vector& obsVector, bncTime& epoTime) { // Default // ------- @@ -114,5 +151,5 @@ int numValidGPS = 0; for (int ii = 0; ii < numSat; ii++) { - char system = satObs[ii]._satellite._system; + char system = satObs[ii]._prn.system(); if (system == 'G' || (system == 'R' && OPT->useGlonass())) { t_satObs* satObs = new t_satObs(satObs[ii]); @@ -142,5 +179,5 @@ if (fabs(dt) > MAXSYNC) { LOG << "t_pppClient::prepareObs asynchronous observations" << endl; - return t_irc::failure; + return failure; } meanDt += dt; @@ -149,10 +186,10 @@ epoTime += meanDt / obsVector.size(); - return t_irc::success; + return success; } // Compute the Bancroft position, check for blunders ////////////////////////////////////////////////////////////////////////////// -t_irc::irc t_pppClient::cmpBancroft(const t_time& epoTime, +t_irc t_pppClient::cmpBancroft(const bncTime& epoTime, vector& obsVector, ColumnVector& xyzc, bool print) { @@ -166,5 +203,5 @@ const t_satObs* satObs = obsVector.at(ii); if ( satObs->isValid() && satObs->prn().system() == 'G' && - (!satObs->modelSet() || satObs->eleSat() >= OPT->minEle()) ) { + (!satObs->modelSet() || satObs->eleSat() >= OPT->_minEle) ) { ++iObs; BB[iObs][0] = satObs->xc()[0]; @@ -174,12 +211,12 @@ } } - if (iObs + 1 < OPT->minobs()) { + if (iObs + 1 < OPT->_minObs) { LOG << "t_pppClient::cmpBancroft not enough observations" << endl; - return t_irc::failure; + return failure; } BB = BB.Rows(1,iObs+1); bancroft(BB, xyzc); - xyzc[3] /= t_genConst::c; + xyzc[3] /= t_CST::c; // Check Blunders @@ -191,8 +228,8 @@ const t_satObs* satObs = obsVector.at(ii); if ( satObs->isValid() && satObs->prn().system() == 'G' && - (!satObs->modelSet() || satObs->eleSat() >= OPT->minEle()) ) { + (!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) - - (satObs->xc()[3] - xyzc[3]) * t_genConst::c; + - (satObs->xc()[3] - xyzc[3]) * t_CST::c; if (fabs(res) > maxRes) { maxRes = fabs(res); @@ -208,5 +245,5 @@ << setw(14) << setprecision(3) << xyzc[1] << ' ' << setw(14) << setprecision(3) << xyzc[2] << ' ' - << setw(14) << setprecision(3) << xyzc[3] * t_genConst::c << endl << endl; + << setw(14) << setprecision(3) << xyzc[3] * t_CST::c << endl << endl; } break; @@ -221,5 +258,5 @@ } - return t_irc::success; + return success; } @@ -239,5 +276,5 @@ t_satObs* satObs = obsVector.at(ii); if ( !satObs->outlier() && satObs->isValid() && satObs->prn().system() == 'R' && - (!satObs->modelSet() || satObs->eleSat() >= OPT->minEle()) ) { + (!satObs->modelSet() || satObs->eleSat() >= OPT->_minEle) ) { ++nObs; double ll = satObs->obsValue(tLC) - satObs->cmpValue(tLC); @@ -270,11 +307,7 @@ void t_pppClient::initOutput(t_pppOutput* output) { _output = output; - _output->_epoTime._mjd = 0; - _output->_epoTime._sec = 0.0; - _output->_numSat = 0; - _output->_pDop = 0.0; - _output->_ambFixRate = 0.0; - _output->_log = 0; - _output->_error = false; + _output->_numSat = 0; + _output->_pDop = 0.0; + _output->_error = false; } @@ -286,20 +319,15 @@ } _obsRover.clear(); - for (unsigned ii = 0; ii < _obsBase.size(); ii++) { - delete _obsBase.at(ii); - } - _obsBase.clear(); -} - -// -////////////////////////////////////////////////////////////////////////////// -void t_pppClient::finish(t_irc::irc irc) { +} + +// +////////////////////////////////////////////////////////////////////////////// +void t_pppClient::finish(t_irc irc) { clearObs(); - _output->_epoTime._mjd = _epoTimeRover.mjd(); - _output->_epoTime._sec = _epoTimeRover.daysec(); - - if (irc == t_irc::success) { + _output->_epoTime = _epoTimeRover; + + if (irc == success) { _output->_xyzRover[0] = _staRover->xyzApr()[0] + _filter->x()[0]; _output->_xyzRover[1] = _staRover->xyzApr()[1] + _filter->x()[1]; @@ -308,5 +336,4 @@ _output->_numSat = _filter->numSat(); _output->_pDop = _filter->PDOP(); - _output->_ambFixRate = _filter->ambFixRate(); _output->_error = false; } @@ -314,7 +341,5 @@ _output->_error = true; } - if (OPT->logLevel() > 0) { - _output->_log = strdup(_log->str().c_str()); - } + _output->_log = _log->str(); delete _log; _log = new ostringstream(); } @@ -322,27 +347,18 @@ // ////////////////////////////////////////////////////////////////////////////// -t_irc::irc t_pppClient::cmpModel(t_station* station, const ColumnVector& xyzc, +t_irc t_pppClient::cmpModel(t_station* station, const ColumnVector& xyzc, vector& obsVector) { - t_time time; - if (station->roverBase() == e_rover) { - time = _epoTimeRover; - station->setName(OPT->roverName()); - station->setAntName(OPT->antNameRover()); - if (OPT->xyzAprRoverSet()) { - station->setXyzApr(OPT->xyzAprRover()); - } - else { - station->setXyzApr(xyzc.Rows(1,3)); - } - station->setNeuEcc(OPT->neuEccRover()); + bncTime time; + time = _epoTimeRover; + station->setName(OPT->_roverName); + station->setAntName(OPT->_antNameRover); + if (OPT->xyzAprRoverSet()) { + station->setXyzApr(OPT->_xyzAprRover); } else { - time = _epoTimeBase; - station->setName(OPT->baseName()); - station->setAntName(OPT->antNameBase()); - station->setXyzApr(OPT->xyzAprBase()); - station->setNeuEcc(OPT->neuEccBase()); - } + station->setXyzApr(xyzc.Rows(1,3)); + } + station->setNeuEcc(OPT->_neuEccRover); // Receiver Clock @@ -356,5 +372,7 @@ // Tides // ----- - station->setTideDspl(t_tides::displacement(time, station->xyzApr())); + ColumnVector hlp = station->xyzApr(); + tides(time, hlp); + station->setTideDspl(hlp - station->xyzApr()); // Observation model @@ -364,5 +382,5 @@ t_satObs* satObs = *it; satObs->cmpModel(station); - if (satObs->isValid() && satObs->eleSat() >= OPT->minEle()) { + if (satObs->isValid() && satObs->eleSat() >= OPT->_minEle) { ++it; } @@ -373,32 +391,5 @@ } - return t_irc::success; -} - -// -////////////////////////////////////////////////////////////////////////////// -t_irc::irc t_pppClient::createDifferences() { - - vector::iterator it = _obsRover.begin(); - while (it != _obsRover.end()) { - t_satObs* obsR = *it; - t_satObs* obsB = 0; - for (unsigned ii = 0; ii < _obsBase.size(); ii++) { - if (_obsBase[ii]->prn() == obsR->prn()) { - obsB = _obsBase[ii]; - break; - } - } - if (obsB) { - obsR->createDifference(obsB); - ++it; - } - else { - delete obsR; - it = _obsRover.erase(it); - } - } - - return t_irc::success; + return success; } @@ -414,6 +405,6 @@ // --------------------------------- if (prepareObs(numSatRover, satObsRover, _obsRover, - _epoTimeRover) != t_irc::success) { - return finish(t_irc::failure); + _epoTimeRover) != success) { + return finish(failure); } @@ -425,42 +416,13 @@ ColumnVector xyzc(4); xyzc = 0.0; bool print = (iter == 2); - if (cmpBancroft(_epoTimeRover, _obsRover, xyzc, print) != t_irc::success) { - return finish(t_irc::failure); - } - if (cmpModel(_staRover, xyzc, _obsRover) != t_irc::success) { - return finish(t_irc::failure); + if (cmpBancroft(_epoTimeRover, _obsRover, xyzc, print) != success) { + return finish(failure); + } + if (cmpModel(_staRover, xyzc, _obsRover) != success) { + return finish(failure); } } _offGG = cmpOffGG(_obsRover); - - // Prepare Observations of the Base - // -------------------------------- - if (OPT->isDifferential()) { - if (prepareObs(numSatBase, satObsBase, _obsBase, - _epoTimeBase) != t_irc::success) { - return finish(t_irc::failure); - } - ColumnVector xyzc(4); xyzc = 0.0; - if (cmpModel(_staBase, xyzc, _obsBase) != t_irc::success) { - return finish(t_irc::failure); - } - if (cmpBancroft(_epoTimeBase, _obsBase, xyzc, true) != t_irc::success) { - return finish(t_irc::failure); - } - _staBase->setDClk(xyzc[3]); - - _offGG -= cmpOffGG(_obsBase); - - // Create single-differences of observations - // ----------------------------------------- - if (createDifferences() != t_irc::success) { - return finish(t_irc::failure); - } - for (unsigned ii = 0; ii < _obsBase.size(); ii++) { - delete _obsBase[ii]; - } - _obsBase.clear(); - } // Store last epoch of data @@ -470,29 +432,108 @@ // Process Epoch in Filter // ----------------------- - if (_filter->processEpoch(_obsPool) != t_irc::success) { - return finish(t_irc::failure); + if (_filter->processEpoch(_obsPool) != success) { + return finish(failure); } } catch (Exception& exc) { LOG << exc.what() << endl; - return finish(t_irc::failure); + return finish(failure); } catch (string& msg) { LOG << "Exception: " << msg << endl; - return finish(t_irc::failure); + return finish(failure); } catch (logic_error exc) { LOG << exc.what() << endl; - return finish(t_irc::failure); + return finish(failure); } catch (const char* msg) { LOG << msg << endl; - return finish(t_irc::failure); + return finish(failure); } catch (...) { LOG << "unknown exception" << endl; - return finish(t_irc::failure); - } - - return finish(t_irc::success); -} + return finish(failure); + } + + return finish(success); +} + +// +//////////////////////////////////////////////////////////////////////////// +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]; +} + +// +//////////////////////////////////////////////////////////////////////////// +void 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); + } + } +} + Index: /trunk/BNC/src/PPP/pppClient.h =================================================================== --- /trunk/BNC/src/PPP/pppClient.h (revision 5739) +++ /trunk/BNC/src/PPP/pppClient.h (revision 5740) @@ -4,5 +4,9 @@ #include #include -#include "RTCM3/ephemeris.h" +#include "ppp.h" +#include "ephemeris.h" +#include "options.h" + +class bncAntex; namespace BNC { @@ -12,5 +16,4 @@ class t_satObs; class t_station; -class t_antex; class t_filter; @@ -19,17 +22,18 @@ t_pppClient(); ~t_pppClient(); - void setOptions(const t_pppOpt* opt); + void setOptions(const t_options* opt); void putGPSEphemeris(const t_ephGPS* eph); void putGloEphemeris(const t_ephGlo* eph); - void putOrbCorrections(int numCorr, const t_pppOrbCorr* corr); - void putClkCorrections(int numCorr, const t_pppClkCorr* corr); - void putBiases(int numBiases, const t_pppSatBiases* biases); + void putOrbCorrections(int numCorr, const t_orbCorr* corr); + void putClkCorrections(int numCorr, const t_clkCorr* corr); + void putBiases(int numBiases, const t_satBiases* biases); void processEpoch(int numSatRover, const t_pppSatObs* satObsRover, t_pppOutput* output); const t_ephPool* ephPool() const {return _ephPool;} const t_obsPool* obsPool() const {return _obsPool;} - const t_antex* antex() const {return _antex;} + const bncAntex* antex() const {return _antex;} const t_station* staRover() const {return _staRover;} double offGG() const {return _offGG;} + static void bancroft(const Matrix& BBpass, ColumnVector& pos); std::ostringstream* _log; @@ -38,14 +42,13 @@ private: void initOutput(t_pppOutput* output); - void finish(t_irc::irc irc); + void finish(t_irc irc); void clearObs(); - t_irc::irc prepareObs(int numSat, const t_pppSatObs* satObs, - std::vector& obsVector, t_time& epoTime); - t_irc::irc cmpModel(t_station* station, const ColumnVector& xyzc, + t_irc prepareObs(int numSat, const t_pppSatObs* satObs, + std::vector& obsVector, bncTime& epoTime); + t_irc cmpModel(t_station* station, const ColumnVector& xyzc, std::vector& obsVector); - t_irc::irc cmpBancroft(const t_time& epoTime, + t_irc cmpBancroft(const bncTime& epoTime, std::vector& obsVector, ColumnVector& xyzc, bool print); - t_irc::irc createDifferences(); double cmpOffGG(std::vector& obsVector); @@ -53,13 +56,10 @@ t_ephPool* _ephPool; t_obsPool* _obsPool; - t_time _epoTimeRover; - t_time _epoTimeBase; + bncTime _epoTimeRover; t_station* _staRover; - t_station* _staBase; - t_antex* _antex; + bncAntex* _antex; t_filter* _filter; double _offGG; std::vector _obsRover; - std::vector _obsBase; };