Changeset 7203 in ntrip

Timestamp:

Aug 17, 2015, 12:30:54 PM (9 years ago)

Author:

stuerze

Message:

some renaming regarding PPP

Location:

trunk/BNC/src

Files:

• PPP/pppClient.cpp (modified) (4 diffs)
• PPP/pppClient.h (modified) (1 diff)
• PPP/pppFilter.cpp (modified) (2 diffs)
• PPP/pppFilter.h (modified) (1 diff)
• PPP_AR (added)
• PPP_AR/pppClient.cpp (added)
• PPP_AR/pppClient.h (added)
• PPP_AR/pppEphPool.cpp (added)
• PPP_AR/pppEphPool.h (added)
• PPP_AR/pppFilter.cpp (added)
• PPP_AR/pppFilter.h (added)
• PPP_AR/pppObsPool.cpp (added)
• PPP_AR/pppObsPool.h (added)
• PPP_AR/pppParlist.cpp (added)
• PPP_AR/pppParlist.h (added)
• PPP_AR/pppSatObs.cpp (added)
• PPP_AR/pppSatObs.h (added)
• PPP_AR/pppStation.cpp (added)
• PPP_AR/pppStation.h (added)
• src.pri (modified) (1 diff)

trunk/BNC/src/PPP/pppClient.cpp

@@ +1 @@
+// 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
@@ -5 +29 @@
  * Class:      t_pppClient
  *
- * Purpose:    PPP Client processing starts here
+ * Purpose:    Precise Point Positioning
  *
  * Author:     L. Mervart
  *
- * Created:    29-Jul-2014
+ * Created:    21-Nov-2009
  *
  * Changes:    
@@ -15 +39 @@
  * -----------------------------------------------------------------------*/
 
-#include <QThreadStorage>
-
-#include <iostream>
+#include <newmatio.h>
 #include <iomanip>
-#include <stdlib.h>
-#include <string.h>
-#include <stdexcept>
+#include <sstream>
 
 #include "pppClient.h"
-#include "pppEphPool.h"
-#include "pppObsPool.h"
-#include "bncconst.h"
+#include "bncephuser.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
+// Constructor
+////////////////////////////////////////////////////////////////////////////
+t_pppClient::t_pppClient(const t_pppOptions* opt) {
+
+  _opt     = new t_pppOptions(*opt);
+  _filter  = new t_pppFilter(this);
+  _epoData = new t_epoData();
+  _log     = new ostringstream();
+  _ephUser = new bncEphUser(false);
+
+  for (unsigned ii = 0; ii <= t_prn::MAXPRN; ii++) {
+    _satCodeBiases[ii] = 0;
+  }
+
+}
+
+// Destructor
+////////////////////////////////////////////////////////////////////////////
+t_pppClient::~t_pppClient() {
+
+  for (unsigned ii = 0; ii <= t_prn::MAXPRN; ii++) {
+    delete _satCodeBiases[ii];
+  }
+
+  delete _filter;
+  delete _epoData;
+  delete _opt;
+  delete _ephUser;
+  delete _log;
+}
+
+//
+////////////////////////////////////////////////////////////////////////////
+void t_pppClient::processEpoch(const vector<t_satObs*>& satObs, t_output* output) {
+  
+  // Convert and store observations
+  // ------------------------------
+  _epoData->clear();
+  for (unsigned ii = 0; ii < satObs.size(); ii++) {
+    const t_satObs* obs     = satObs[ii]; 
+    t_prn prn = obs->_prn;
+    if (prn.system() == 'E') {prn.setFlags(1);} // force I/NAV usage
+    t_satData*   satData = new t_satData();
+
+    if (_epoData->tt.undef()) {
+      _epoData->tt = obs->_time;
+    }
+
+    satData->tt       = obs->_time;
+    satData->prn      = QString(prn.toInternalString().c_str());
+    satData->slipFlag = false;
+    satData->P1       = 0.0;
+    satData->P2       = 0.0;
+    satData->P5       = 0.0;
+    satData->P7       = 0.0;
+    satData->L1       = 0.0;
+    satData->L2       = 0.0;
+    satData->L5       = 0.0;
+    satData->L7       = 0.0;
+    for (unsigned ifrq = 0; ifrq < obs->_obs.size(); ifrq++) {
+      t_frqObs* frqObs = obs->_obs[ifrq];
+      double cb = 0.0;
+      const t_satCodeBias* satCB = satCodeBias(prn);
+      if (satCB && satCB->_bias.size()) {
+        for (unsigned ii = 0; ii < satCB->_bias.size(); ii++) {
+
+          const t_frqCodeBias& bias = satCB->_bias[ii];
+          if (frqObs && frqObs->_rnxType2ch == bias._rnxType2ch) {
+            cb  = bias._value;
+          }
+        }
+      }
+      if      (frqObs->_rnxType2ch[0] == '1') {
+        if (frqObs->_codeValid)  satData->P1       = frqObs->_code + cb;
+        if (frqObs->_phaseValid) satData->L1       = frqObs->_phase;
+        if (frqObs->_slip)       satData->slipFlag = true;
+      }
+      else if (frqObs->_rnxType2ch[0] == '2') {
+        if (frqObs->_codeValid)  satData->P2       = frqObs->_code + cb;
+        if (frqObs->_phaseValid) satData->L2       = frqObs->_phase;
+        if (frqObs->_slip)       satData->slipFlag = true;
+      }
+      else if (frqObs->_rnxType2ch[0] == '5') {
+        if (frqObs->_codeValid)  satData->P5       = frqObs->_code + cb;
+        if (frqObs->_phaseValid) satData->L5       = frqObs->_phase;
+        if (frqObs->_slip)       satData->slipFlag = true;
+      }
+      else if (frqObs->_rnxType2ch[0] == '7') {
+        if (frqObs->_codeValid)  satData->P7       = frqObs->_code + cb;
+        if (frqObs->_phaseValid) satData->L7       = frqObs->_phase;
+        if (frqObs->_slip)       satData->slipFlag = true;
+      }
+    }
+    putNewObs(satData);
+  }
+
+  // Data Pre-Processing
+  // -------------------
+  QMutableMapIterator<QString, t_satData*> it(_epoData->satData);
+  while (it.hasNext()) {
+    it.next();
+    QString    prn     = it.key();
+    t_satData* satData = it.value();
+    
+    if (cmpToT(satData) != success) {
+      delete satData;
+      it.remove();
+      continue;
+    }
+
+  }
+
+  // Filter Solution
+  // ---------------
+  if (_filter->update(_epoData) == success) {
+    output->_error = false;
+    output->_epoTime     = _filter->time();
+    output->_xyzRover[0] = _filter->x();
+    output->_xyzRover[1] = _filter->y();
+    output->_xyzRover[2] = _filter->z();
+    copy(&_filter->Q().data()[0], &_filter->Q().data()[6], output->_covMatrix);
+    output->_neu[0]      = _filter->neu()[0];
+    output->_neu[1]      = _filter->neu()[1];
+    output->_neu[2]      = _filter->neu()[2];
+    output->_numSat      = _filter->numSat();
+    output->_pDop        = _filter->PDOP();
+    output->_trp0        = _filter->trp0();
+    output->_trp         = _filter->trp();
+  }
+  else {
+    output->_error = true;
+  }
+
+  output->_log = _log->str();
+  delete _log; _log = new ostringstream();
+}
+
+//
+////////////////////////////////////////////////////////////////////////////
+void t_pppClient::putNewObs(t_satData* satData) {
+
+  // Set Observations GPS
+  // --------------------
+  if      (satData->system() == 'G') {
+    if (satData->P1 != 0.0 && satData->P2 != 0.0 &&
+        satData->L1 != 0.0 && satData->L2 != 0.0 ) {
+      t_frequency::type fType1 = t_lc::toFreq(satData->system(), t_lc::l1);
+      t_frequency::type fType2 = t_lc::toFreq(satData->system(), t_lc::l2);
+      double f1 = t_CST::freq(fType1, 0);
+      double f2 = t_CST::freq(fType2, 0);
+      double a1 =   f1 * f1 / (f1 * f1 - f2 * f2);
+      double a2 = - f2 * f2 / (f1 * f1 - f2 * f2);
+      satData->L1      = satData->L1 * t_CST::c / f1;
+      satData->L2      = satData->L2 * t_CST::c / f2;
+      satData->P3      = a1 * satData->P1 + a2 * satData->P2;
+      satData->L3      = a1 * satData->L1 + a2 * satData->L2;
+      satData->lambda3 = a1 * t_CST::c / f1 + a2 * t_CST::c / f2;
+      satData->lkA     = a1;
+      satData->lkB     = a2;
+      _epoData->satData[satData->prn] = satData;
+    }
+    else {
+      delete satData;
+    }
+  }
+
+  // Set Observations Glonass
+  // ------------------------
+  else if (satData->system() == 'R' && _opt->useSystem('R')) {
+    if (satData->P1 != 0.0 && satData->P2 != 0.0 && 
+        satData->L1 != 0.0 && satData->L2 != 0.0 ) {
+
+      int channel = 0;
+      if (satData->system() == 'R') {
+        const t_eph* eph = _ephUser->ephLast(satData->prn);
+        if (eph) {
+          channel = eph->slotNum();
+        }
+        else {
+          delete satData;
+          return;
+        }
+      }
+
+      t_frequency::type fType1 = t_lc::toFreq(satData->system(), t_lc::l1);
+      t_frequency::type fType2 = t_lc::toFreq(satData->system(), t_lc::l2);
+      double f1 = t_CST::freq(fType1, channel);
+      double f2 = t_CST::freq(fType2, channel);
+      double a1 =   f1 * f1 / (f1 * f1 - f2 * f2);
+      double a2 = - f2 * f2 / (f1 * f1 - f2 * f2);
+      satData->L1      = satData->L1 * t_CST::c / f1;
+      satData->L2      = satData->L2 * t_CST::c / f2;
+      satData->P3      = a1 * satData->P1 + a2 * satData->P2;
+      satData->L3      = a1 * satData->L1 + a2 * satData->L2;
+      satData->lambda3 = a1 * t_CST::c / f1 + a2 * t_CST::c / f2;
+      satData->lkA     = a1;
+      satData->lkB     = a2;
+      _epoData->satData[satData->prn] = satData;
+    }
+    else {
+      delete satData;
+    }
+  }
+
+  // Set Observations Galileo
+  // ------------------------
+  else if (satData->system() == 'E' && _opt->useSystem('E')) {
+    if (satData->P1 != 0.0 && satData->P5 != 0.0 && 
+        satData->L1 != 0.0 && satData->L5 != 0.0 ) {
+      double f1 = t_CST::freq(t_frequency::E1, 0);
+      double f5 = t_CST::freq(t_frequency::E5, 0);
+      double a1 =   f1 * f1 / (f1 * f1 - f5 * f5);
+      double a5 = - f5 * f5 / (f1 * f1 - f5 * f5);
+      satData->L1      = satData->L1 * t_CST::c / f1;
+      satData->L5      = satData->L5 * t_CST::c / f5;
+      satData->P3      = a1 * satData->P1 + a5 * satData->P5;
+      satData->L3      = a1 * satData->L1 + a5 * satData->L5;
+      satData->lambda3 = a1 * t_CST::c / f1 + a5 * t_CST::c / f5;
+      satData->lkA     = a1;
+      satData->lkB     = a5;
+      _epoData->satData[satData->prn] = satData;
+    }
+    else {
+      delete satData;
+    }
+  }
+
+  // Set Observations BDS
+  // ---------------------
+  else if (satData->system() == 'C' && _opt->useSystem('C')) {
+    if (satData->P2 != 0.0 && satData->P7 != 0.0 &&
+        satData->L2 != 0.0 && satData->L7 != 0.0 ) {
+      double f2 = t_CST::freq(t_frequency::C2, 0);
+      double f7 = t_CST::freq(t_frequency::C7, 0);
+      double a2 =   f2 * f2 / (f2 * f2 - f7 * f7);
+      double a7 = - f7 * f7 / (f2 * f2 - f7 * f7);
+      satData->L2      = satData->L2 * t_CST::c / f2;
+      satData->L7      = satData->L7 * t_CST::c / f7;
+      satData->P3      = a2 * satData->P2 + a7 * satData->P7;
+      satData->L3      = a2 * satData->L2 + a7 * satData->L7;
+      satData->lambda3 = a2 * t_CST::c / f2 + a7 * t_CST::c / f7;
+      satData->lkA     = a2;
+      satData->lkB     = a7;
+      _epoData->satData[satData->prn] = satData;
+    }
+    else {
+      delete satData;
+    }
+  }
+}
+
+// 
+////////////////////////////////////////////////////////////////////////////
+void t_pppClient::putOrbCorrections(const std::vector<t_orbCorr*>& corr) {
+  for (unsigned ii = 0; ii < corr.size(); ii++) {
+    QString prn = QString(corr[ii]->_prn.toInternalString().c_str());
+    t_eph* eLast = _ephUser->ephLast(prn);
+    t_eph* ePrev = _ephUser->ephPrev(prn);
+    if      (eLast && eLast->IOD() == corr[ii]->_iod) {
+      eLast->setOrbCorr(corr[ii]);
+    }
+    else if (ePrev && ePrev->IOD() == corr[ii]->_iod) {
+      ePrev->setOrbCorr(corr[ii]);
+    }
+  }
+}
+
+// 
+////////////////////////////////////////////////////////////////////////////
+void t_pppClient::putClkCorrections(const std::vector<t_clkCorr*>& corr) {
+  for (unsigned ii = 0; ii < corr.size(); ii++) {
+    QString prn = QString(corr[ii]->_prn.toInternalString().c_str());
+    t_eph* eLast = _ephUser->ephLast(prn);
+    t_eph* ePrev = _ephUser->ephPrev(prn);
+    if      (eLast && eLast->IOD() == corr[ii]->_iod) {
+      eLast->setClkCorr(corr[ii]);
+    }
+    else if (ePrev && ePrev->IOD() == corr[ii]->_iod) {
+      ePrev->setClkCorr(corr[ii]);
+    }
+  }
+}
+
+// 
 //////////////////////////////////////////////////////////////////////////////
-QThreadStorage<t_pppClient*> 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();
-  _tides    = new t_tides();
-
-  if (!_opt->_antexFileName.empty()) {
-    _antex = new bncAntex(_opt->_antexFileName.c_str());
-  }
-  else {
-    _antex = 0;
-  }
-
-  CLIENTS.setLocalData(this);  // CLIENTS takes ownership over "this"
-}
-
-// Destructor
-//////////////////////////////////////////////////////////////////////////////
-t_pppClient::~t_pppClient() {
-  delete _log;
-  delete _opt;
-  delete _ephPool;
-  delete _obsPool;
-  delete _staRover;
-  delete _antex;
-  delete _filter;
-  delete _tides;
-  clearObs();
+void t_pppClient::putCodeBiases(const std::vector<t_satCodeBias*>& satCodeBias) {
+  for (unsigned ii = 0; ii < satCodeBias.size(); ii++) {
+    putCodeBias(new t_satCodeBias(*satCodeBias[ii]));
+  }
 }
 
@@ -84 +344 @@
 //////////////////////////////////////////////////////////////////////////////
 void t_pppClient::putEphemeris(const t_eph* eph) {
+  bool check = _opt->_realTime;
   const t_ephGPS* ephGPS = dynamic_cast<const t_ephGPS*>(eph);
   const t_ephGlo* ephGlo = dynamic_cast<const t_ephGlo*>(eph);
   const t_ephGal* ephGal = dynamic_cast<const t_ephGal*>(eph);
+  const t_ephBDS* ephBDS = dynamic_cast<const t_ephBDS*>(eph);
   if      (ephGPS) {
-    _ephPool->putEphemeris(new t_ephGPS(*ephGPS));
+    _ephUser->putNewEph(new t_ephGPS(*ephGPS), check);
   }
   else if (ephGlo) {
-    _ephPool->putEphemeris(new t_ephGlo(*ephGlo));
+    _ephUser->putNewEph(new t_ephGlo(*ephGlo), check);
   }
   else if (ephGal) {
-    _ephPool->putEphemeris(new t_ephGal(*ephGal));
-  }
-}
-
-// 
-//////////////////////////////////////////////////////////////////////////////
-void t_pppClient::putOrbCorrections(const vector<t_orbCorr*>& corr) {
-  for (unsigned ii = 0; ii < corr.size(); ii++) {
-    _ephPool->putOrbCorrection(new t_orbCorr(*corr[ii]));
-  }
-}
-
-// 
-//////////////////////////////////////////////////////////////////////////////
-void t_pppClient::putClkCorrections(const vector<t_clkCorr*>& corr) {
-  for (unsigned ii = 0; ii < corr.size(); ii++) {
-    _ephPool->putClkCorrection(new t_clkCorr(*corr[ii]));
-  }
-}
-
-// 
-//////////////////////////////////////////////////////////////////////////////
-void t_pppClient::putCodeBiases(const vector<t_satCodeBias*>& biases) {
-  for (unsigned ii = 0; ii < biases.size(); ii++) {
-    _obsPool->putCodeBias(new t_satCodeBias(*biases[ii]));
-  }
-}
-
-// 
-//////////////////////////////////////////////////////////////////////////////
-t_irc t_pppClient::prepareObs(const vector<t_satObs*>& satObs,
-                              vector<t_pppSatObs*>& obsVector, bncTime& epoTime) {
-  // Default 
-  // -------
-  epoTime.reset();
-
-  // Create vector of valid observations
-  // -----------------------------------
-  for (unsigned ii = 0; ii < satObs.size(); ii++) {
-    char system = satObs[ii]->_prn.system();
-    if (OPT->useSystem(system)) {
-      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;
-}
-
-// Compute the Bancroft position, check for blunders
-//////////////////////////////////////////////////////////////////////////////
-t_irc t_pppClient::cmpBancroft(const bncTime& epoTime, 
-                                  vector<t_pppSatObs*>& 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 (satObs->prn().system() == 'G') {
-        if (tLC == t_lc::dummy) {
-          tLC = satObs->isValid(t_lc::cIF) ? t_lc::cIF : t_lc::c1;
-        }
-        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" << endl;
+    _ephUser->putNewEph(new t_ephGal(*ephGal), check);
+  }
+  else if (ephBDS) {
+      _ephUser->putNewEph(new t_ephBDS(*ephBDS), check);
+    }
+}
+
+// Satellite Position
+////////////////////////////////////////////////////////////////////////////
+t_irc t_pppClient::getSatPos(const bncTime& tt, const QString& prn, 
+                              ColumnVector& xc, ColumnVector& vv) {
+
+  t_eph* eLast = _ephUser->ephLast(prn);
+  t_eph* ePrev = _ephUser->ephPrev(prn);
+  if      (eLast && eLast->getCrd(tt, xc, vv, _opt->useOrbClkCorr()) == success) {
+    return success;
+  }
+  else if (ePrev && ePrev->getCrd(tt, xc, vv, _opt->useOrbClkCorr()) == success) {
+    return success;
+  }
+  return failure;
+}
+
+// Correct Time of Transmission
+////////////////////////////////////////////////////////////////////////////
+t_irc t_pppClient::cmpToT(t_satData* satData) {
+
+  double prange = satData->P3;
+  if (prange == 0.0) {
+    return failure;
+  }
+
+  double clkSat = 0.0;
+  for (int ii = 1; ii <= 10; ii++) {
+
+    bncTime ToT = satData->tt - prange / t_CST::c - clkSat;
+
+    ColumnVector xc(4);
+    ColumnVector vv(3);
+    if (getSatPos(ToT, satData->prn, xc, vv) != success) {
       return failure;
     }
-    BB = BB.Rows(1,iObs+1);
-    bancroft(BB, xyzc);
-
-    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->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) 
-          - (satObs->xc()[3] - xyzc[3]) * t_CST::c;
-        if (fabs(res) > maxRes) {
-          maxRes      = fabs(res);
-          maxResIndex = ii;
-        }
-      }
-    }
-    if (maxRes < BLUNDER) {
-      if (print) {
-        LOG.setf(ios::fixed);
-        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;
-}
-
-// Compute A Priori GPS-Glonass Offset
-//////////////////////////////////////////////////////////////////////////////
-double t_pppClient::cmpOffGG(vector<t_pppSatObs*>& obsVector) {
-
-  t_lc::type tLC   = t_lc::dummy;
-  double     offGG = 0.0;
-
-  if (OPT->useSystem('R')) {
-
-    while (obsVector.size() > 0) {
-      offGG = 0.0;
-      double   maxRes      = 0.0;
-      int      maxResIndex = -1;
-      t_prn    maxResPrn;
-      unsigned nObs        = 0;
-      for (unsigned ii = 0; ii < obsVector.size(); ii++) {
-        t_pppSatObs* satObs = obsVector.at(ii);
-        if (satObs->prn().system() == 'R') {
-          if (tLC == t_lc::dummy) {
-            tLC = satObs->isValid(t_lc::cIF) ? t_lc::cIF : t_lc::c1;
-          }
-          if (satObs->isValid(tLC) && (!satObs->modelSet() || satObs->eleSat() >= OPT->_minEle)) {
-            double ll = satObs->obsValue(tLC) - satObs->cmpValue(tLC);
-            ++nObs;
-            offGG += ll;
-            if (fabs(ll) > fabs(maxRes)) {
-              maxRes      = ll;
-              maxResIndex = ii;
-              maxResPrn   = satObs->prn();
-            }
-          }
-        }
-      }
-
-      if (nObs > 0) {
-        offGG = offGG / nObs;
-      }
-      else {
-        offGG = 0.0;
-      }
-
-      if (fabs(maxRes) > 1000.0) {
-        LOG << "t_pppClient::cmpOffGG outlier " << maxResPrn.toString() << " " << maxRes << endl;
-        obsVector.erase(obsVector.begin() + maxResIndex);
-      }
-      else {
-        break;
-      }
-    }
-  }
-
-  return offGG;
-}
-
-// 
-//////////////////////////////////////////////////////////////////////////////
-void t_pppClient::initOutput(t_output* output) {
-  _output = output;
-  _output->_numSat = 0;
-  _output->_pDop   = 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) {
-
-  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->_pDop       = _filter->PDOP();
-    _output->_error = false;
-  }
-  else {
-    _output->_error = true;
-  }
-  _output->_log = _log->str();
-  delete _log; _log = new ostringstream();
-}
-
-// 
-//////////////////////////////////////////////////////////////////////////////
-t_irc t_pppClient::cmpModel(t_pppStation* station, const ColumnVector& xyzc,
-                               vector<t_pppSatObs*>& obsVector) {
-
-  bncTime time;
-  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);
-
-  // Receiver Clock
-  // --------------
-  station->setDClk(xyzc[3]);
-
-  // Tides
-  // -----
-  station->setTideDspl( _tides->displacement(time, station->xyzApr()) );
-  
-  // Observation model
-  // -----------------
-  vector<t_pppSatObs*>::iterator it = obsVector.begin();
-  while (it != obsVector.end()) {
-    t_pppSatObs* satObs = *it;
-    if (satObs->isValid()) {
-      satObs->cmpModel(station);
-    }
-    if (satObs->isValid() && satObs->eleSat() >= OPT->_minEle) {
-      ++it;
-    }
-    else {
-      delete satObs;
-      it = obsVector.erase(it);
-    }
-  }
-
-  return success;
-}
-
-// 
-//////////////////////////////////////////////////////////////////////////////
-void t_pppClient::processEpoch(const vector<t_satObs*>& satObs, t_output* output) {
-
-  try {
-    initOutput(output);
-
-    // Prepare Observations of the Rover
-    // ---------------------------------    
-    if (prepareObs(satObs, _obsRover, _epoTimeRover) != success) {
-      return finish(failure);
-    }
-
-    LOG << "\nResults 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);
-
-    // Store last epoch of data
-    // ------------------------    
-    _obsPool->putEpoch(_epoTimeRover, _obsRover);
-
-    // Process Epoch in Filter
-    // -----------------------
-    if (_filter->processEpoch(_obsPool) != success) {
-      return finish(failure);
-    }
-  }
-  catch (Exception& exc) {
-    LOG << exc.what() << endl;
-    return finish(failure);
-  }
-  catch (t_except msg) {
-    LOG << msg.what() << endl;
-    return finish(failure);
-  }
-  catch (const char* msg) {
-    LOG << msg << endl;
-    return finish(failure);
-  }
-  catch (...) {
-    LOG << "unknown exception" << endl;
-    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);
-    }
-  }
-}
-
+
+    double clkSatOld = clkSat;
+    clkSat = xc(4);
+
+    if ( fabs(clkSat-clkSatOld) * t_CST::c < 1.e-4 ) {
+      satData->xx      = xc.Rows(1,3);
+      satData->vv      = vv;
+      satData->clk     = clkSat * t_CST::c;
+      return success;
+    } 
+  }
+
+  return failure;
+}
+
+void t_pppClient::putCodeBias(t_satCodeBias* satCodeBias) {
+  int iPrn = satCodeBias->_prn.toInt();
+  delete _satCodeBiases[iPrn];
+  _satCodeBiases[iPrn] = satCodeBias;
+}

trunk/BNC/src/PPP/pppClient.h

@@ +1 @@
+// 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.
+
 #ifndef PPPCLIENT_H
 #define PPPCLIENT_H
 
-#include <sstream>
 #include <vector>
+#include <QtCore>
+
 #include "pppInclude.h"
-#include "ephemeris.h"
 #include "pppOptions.h"
-#include "pppModel.h"
-#include "satObs.h"
+#include "pppFilter.h"
 
-class bncAntex;
+class bncEphUser;
+class t_eph;
 
 namespace BNC_PPP {
-
-class t_pppEphPool;
-class t_pppObsPool;
-class t_pppSatObs;
-class t_pppStation;
-class t_pppFilter;
-
+  
 class t_pppClient : public interface_pppClient {
  public:
-  t_pppClient(const t_pppOptions* opt);                                                      
-  ~t_pppClient();                                                     
+  t_pppClient(const t_pppOptions* opt);
+  ~t_pppClient();
+  void                processEpoch(const std::vector<t_satObs*>& satObs, t_output* output);
+  void                putEphemeris(const t_eph* eph);                  
+  void                putOrbCorrections(const std::vector<t_orbCorr*>& corr); 
+  void                putClkCorrections(const std::vector<t_clkCorr*>& corr); 
+  void                putCodeBiases(const std::vector<t_satCodeBias*>& satCodeBias);   
+  std::ostringstream& log() {return *_log;}
+  const t_pppOptions* opt() const {return _opt;}
+  void putCodeBias(t_satCodeBias* satCodeBias);
+  const t_satCodeBias* satCodeBias(const t_prn& prn) const {
+    return _satCodeBiases[prn.toInt()];
+  }
+ private:
+  t_irc getSatPos(const bncTime& tt, const QString& prn, ColumnVector& xc, ColumnVector& vv);
+  void  putNewObs(t_satData* satData);
+  t_irc cmpToT(t_satData* satData);
 
-  void putEphemeris(const t_eph* eph);                  
-  void putOrbCorrections(const std::vector<t_orbCorr*>& corr); 
-  void putClkCorrections(const std::vector<t_clkCorr*>& corr); 
-  void putCodeBiases(const std::vector<t_satCodeBias*>& satBias);   
-  void processEpoch(const std::vector<t_satObs*>& satObs, t_output* output);
-
-  const t_pppEphPool* ephPool() const {return _ephPool;}
-  const t_pppObsPool* obsPool() const {return _obsPool;}
-  const bncAntex*  antex() const {return _antex;}
-  const t_pppStation* staRover() const {return _staRover;}
-  double           offGG() const {return _offGG;}
-
-  std::ostringstream& log() {return *_log;}
-  const t_pppOptions*    opt() const {return _opt;}
-
-  static void bancroft(const Matrix& BBpass, ColumnVector& pos);
-
-  static t_pppClient* instance();
-
- private:
-  void initOutput(t_output* output);
-  void finish(t_irc irc);
-  void clearObs();
-  t_irc prepareObs(const std::vector<t_satObs*>& satObs,
-                   std::vector<t_pppSatObs*>& obsVector, bncTime& epoTime);
-  t_irc cmpModel(t_pppStation* station, const ColumnVector& xyzc,
-                 std::vector<t_pppSatObs*>& obsVector);
-  t_irc cmpBancroft(const bncTime& epoTime, std::vector<t_pppSatObs*>& obsVector, 
-                    ColumnVector& xyzc, bool print);
-  double cmpOffGG(std::vector<t_pppSatObs*>& obsVector);
-
-  t_output*                 _output;
-  t_pppEphPool*             _ephPool;
-  t_pppObsPool*             _obsPool;
-  bncTime                   _epoTimeRover;
-  t_pppStation*             _staRover;
-  bncAntex*                 _antex;
-  t_pppFilter*              _filter;
-  double                    _offGG;
-  std::vector<t_pppSatObs*> _obsRover;
-  std::ostringstream*       _log; 
-  t_pppOptions*             _opt; 
-  t_tides*                  _tides;
+  t_satCodeBias*       _satCodeBiases[t_prn::MAXPRN+1];
+  bncEphUser*         _ephUser;
+  t_pppOptions*       _opt;
+  t_epoData*          _epoData;
+  t_pppFilter*        _filter;
+  std::ostringstream* _log; 
 };
 
-}; // namespace BNC_PPP
-
-#define PPP_CLIENT (BNC_PPP::t_pppClient::instance())
-#define LOG        (BNC_PPP::t_pppClient::instance()->log())
-#define OPT        (BNC_PPP::t_pppClient::instance()->opt())
+} // namespace
 
 #endif

trunk/BNC/src/PPP/pppFilter.cpp

@@ +1 @@
+// 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_pppFilter
+ * Class:      t_pppParam, t_pppFilter
  *
- * Purpose:    Filter Adjustment
+ * Purpose:    Model for PPP
  *
  * Author:     L. Mervart
  *
- * Created:    29-Jul-2014
+ * Created:    01-Dec-2009
  *
  * Changes:    
@@ -15 +39 @@
  * -----------------------------------------------------------------------*/
 
-#include <iostream>
 #include <iomanip>
 #include <cmath>
-#include <newmat.h>
+#include <sstream>
 #include <newmatio.h>
 #include <newmatap.h>
 
 #include "pppFilter.h"
+#include "pppClient.h"
 #include "bncutils.h"
-#include "pppParlist.h"
-#include "pppObsPool.h"
-#include "pppStation.h"
-#include "pppClient.h"
+#include "bncantex.h"
+#include "pppOptions.h"
+#include "pppModel.h"
 
 using namespace BNC_PPP;
 using namespace std;
 
+const double   MAXRES_CODE           = 2.98 * 3.0;
+const double   MAXRES_PHASE_GPS      = 2.98 * 0.03;
+const double   MAXRES_PHASE_GLONASS  = 2.98 * 0.03;
+const double   GLONASS_WEIGHT_FACTOR = 1.0;
+const double   BDS_WEIGHT_FACTOR     = 2.0;
+
+#define LOG (_pppClient->log())
+#define OPT (_pppClient->opt())
+
 // Constructor
 ////////////////////////////////////////////////////////////////////////////
-t_pppFilter::t_pppFilter() {
-  _parlist = 0;
+t_pppParam::t_pppParam(t_pppParam::parType typeIn, int indexIn, 
+                   const QString& prnIn) {
+  type      = typeIn;
+  index     = indexIn;
+  prn       = prnIn;
+  index_old = 0;
+  xx        = 0.0;
+  numEpo    = 0;
 }
 
 // Destructor
 ////////////////////////////////////////////////////////////////////////////
+t_pppParam::~t_pppParam() {
+}
+
+// Partial
+////////////////////////////////////////////////////////////////////////////
+double t_pppParam::partial(t_satData* satData, bool phase) {
+
+  Tracer tracer("t_pppParam::partial");
+
+  // Coordinates
+  // -----------
+  if      (type == CRD_X) {
+    return (xx - satData->xx(1)) / satData->rho; 
+  }
+  else if (type == CRD_Y) {
+    return (xx - satData->xx(2)) / satData->rho; 
+  }
+  else if (type == CRD_Z) {
+    return (xx - satData->xx(3)) / satData->rho; 
+  }
+
+  // Receiver Clocks
+  // ---------------
+  else if (type == RECCLK) {
+    return 1.0;
+  }
+
+  // Troposphere
+  // -----------
+  else if (type == TROPO) {
+    return 1.0 / sin(satData->eleSat); 
+  }
+
+  // Glonass Offset
+  // --------------
+  else if (type == GLONASS_OFFSET) {
+    if (satData->prn[0] == 'R') {
+      return 1.0;
+    }
+    else {
+      return 0.0;
+    }
+  }
+
+  // Galileo Offset
+  // --------------
+  else if (type == GALILEO_OFFSET) {
+    if (satData->prn[0] == 'E') {
+      return 1.0;
+    }
+    else {
+      return 0.0;
+    }
+  }
+
+  // BDS Offset
+  // ----------
+  else if (type == BDS_OFFSET) {
+    if (satData->prn[0] == 'C') {
+      return 1.0;
+    }
+    else {
+      return 0.0;
+    }
+  }
+
+  // Ambiguities
+  // -----------
+  else if (type == AMB_L3) {
+    if (phase && satData->prn == prn) {
+      return 1.0;
+    }
+    else {
+      return 0.0;
+    }
+  }
+
+  // Default return
+  // --------------
+  return 0.0;
+}
+
+// Constructor
+////////////////////////////////////////////////////////////////////////////
+t_pppFilter::t_pppFilter(t_pppClient* pppClient) {
+
+  _pppClient = pppClient;
+  _tides     = new t_tides();
+
+  // Antenna Name, ANTEX File
+  // ------------------------
+  _antex = 0;
+  if (!OPT->_antexFileName.empty()) {
+    _antex = new bncAntex(OPT->_antexFileName.c_str());
+  }
+
+  // Bancroft Coordinates
+  // --------------------
+  _xcBanc.ReSize(4);  _xcBanc  = 0.0;
+  _ellBanc.ReSize(3); _ellBanc = 0.0;
+
+  // Save copy of data (used in outlier detection)
+  // ---------------------------------------------
+  _epoData_sav = new t_epoData();
+
+  // Some statistics
+  // ---------------
+  _neu.ReSize(3); _neu = 0.0;
+  _numSat = 0;
+  _pDop   = 0.0;
+}
+
+// Destructor
+////////////////////////////////////////////////////////////////////////////
 t_pppFilter::~t_pppFilter() {
-  delete _parlist;
-}
-
-// Process Single Epoch
-////////////////////////////////////////////////////////////////////////////
-t_irc t_pppFilter::processEpoch(t_pppObsPool* obsPool) {
-
-  _numSat     = 0;
-
-  if (!_parlist) {
-    _parlist = new t_pppParlist();
-  }
-
-  // Vector of all Observations
-  // --------------------------
-  t_pppObsPool::t_epoch* epoch = obsPool->lastEpoch();
-  if (!epoch) {
+  delete _tides;
+  delete _antex;
+  for (int iPar = 1; iPar <= _params.size(); iPar++) {
+    delete _params[iPar-1];
+  }
+  for (int iPar = 1; iPar <= _params_sav.size(); iPar++) {
+    delete _params_sav[iPar-1];
+  }
+  delete _epoData_sav;
+}
+
+// Reset Parameters and Variance-Covariance Matrix
+////////////////////////////////////////////////////////////////////////////
+void t_pppFilter::reset() {
+
+  Tracer tracer("t_pppFilter::reset");
+
+  double lastTrp = 0.0;
+  for (int ii = 0; ii < _params.size(); ii++) {
+    t_pppParam* pp = _params[ii];
+    if (pp->type == t_pppParam::TROPO) {
+      lastTrp = pp->xx;
+    }
+    delete pp;
+  }
+  _params.clear();
+
+  int nextPar = 0;
+  _params.push_back(new t_pppParam(t_pppParam::CRD_X,  ++nextPar, ""));
+  _params.push_back(new t_pppParam(t_pppParam::CRD_Y,  ++nextPar, ""));
+  _params.push_back(new t_pppParam(t_pppParam::CRD_Z,  ++nextPar, ""));
+  _params.push_back(new t_pppParam(t_pppParam::RECCLK, ++nextPar, ""));
+  if (OPT->estTrp()) {
+    _params.push_back(new t_pppParam(t_pppParam::TROPO, ++nextPar, ""));
+  }
+  if (OPT->useSystem('R')) {
+    _params.push_back(new t_pppParam(t_pppParam::GLONASS_OFFSET, ++nextPar, ""));
+  }
+  if (OPT->useSystem('E')) {
+    _params.push_back(new t_pppParam(t_pppParam::GALILEO_OFFSET, ++nextPar, ""));
+  }
+  if (OPT->useSystem('C')) {
+    _params.push_back(new t_pppParam(t_pppParam::BDS_OFFSET, ++nextPar, ""));
+  }
+
+  _QQ.ReSize(_params.size()); 
+  _QQ = 0.0;
+  for (int iPar = 1; iPar <= _params.size(); iPar++) {
+    t_pppParam* pp = _params[iPar-1];
+    pp->xx = 0.0;
+    if      (pp->isCrd()) {
+      _QQ(iPar,iPar) = OPT->_aprSigCrd(1) * OPT->_aprSigCrd(1); 
+    }
+    else if (pp->type == t_pppParam::RECCLK) {
+      _QQ(iPar,iPar) = OPT->_noiseClk * OPT->_noiseClk; 
+    }
+    else if (pp->type == t_pppParam::TROPO) {
+      _QQ(iPar,iPar) = OPT->_aprSigTrp * OPT->_aprSigTrp; 
+      pp->xx = lastTrp;
+    }
+    else if (pp->type == t_pppParam::GLONASS_OFFSET) {
+      _QQ(iPar,iPar) = 1000.0 * 1000.0;
+    }
+    else if (pp->type == t_pppParam::GALILEO_OFFSET) {
+      _QQ(iPar,iPar) = 1000.0 * 1000.0;
+    }
+    else if (pp->type == t_pppParam::BDS_OFFSET) {
+      _QQ(iPar,iPar) = 1000.0 * 1000.0;
+    }
+  }
+}
+
+// Bancroft Solution
+////////////////////////////////////////////////////////////////////////////
+t_irc t_pppFilter::cmpBancroft(t_epoData* epoData) {
+
+  Tracer tracer("t_pppFilter::cmpBancroft");
+
+  if (int(epoData->sizeSys('G')) < OPT->_minObs) {
+    LOG << "t_pppFilter::cmpBancroft: not enough data\n";
     return failure;
   }
-  vector<t_pppSatObs*>& allObs = epoch->obsVector();
-
-  // Time of the Epoch
-  // -----------------
-  _epoTime = epoch->epoTime();
-
-  if (!_firstEpoTime.valid()) {
-    _firstEpoTime = _epoTime;
-  }
-
-  // Set Parameters
-  // --------------
-  _parlist->set(_epoTime, allObs);
-  const vector<t_pppParam*>& params = _parlist->params();
-
-  // Status Vector, Variance-Covariance Matrix
-  // -----------------------------------------
-  ColumnVector    xFltOld = _xFlt;
-  SymmetricMatrix QFltOld = _QFlt;
-
-  _QFlt.ReSize(_parlist->nPar()); _QFlt = 0.0;
-  _xFlt.ReSize(_parlist->nPar()); _xFlt = 0.0;
-  _x0.ReSize(_parlist->nPar());   _x0   = 0.0;
+
+  Matrix BB(epoData->sizeSys('G'), 4);
+
+  QMapIterator<QString, t_satData*> it(epoData->satData);
+  int iObsBanc = 0;
+  while (it.hasNext()) {
+    it.next();
+    t_satData* satData = it.value();
+    if (satData->system() == 'G') {
+      ++iObsBanc;
+      QString    prn     = it.key();
+      BB(iObsBanc, 1) = satData->xx(1);
+      BB(iObsBanc, 2) = satData->xx(2);
+      BB(iObsBanc, 3) = satData->xx(3);
+      BB(iObsBanc, 4) = satData->P3 + satData->clk;
+    }
+  }
+
+  bancroft(BB, _xcBanc);
+
+  // Ellipsoidal Coordinates
+  // ------------------------
+  xyz2ell(_xcBanc.data(), _ellBanc.data());
+
+  // Compute Satellite Elevations
+  // ----------------------------
+  QMutableMapIterator<QString, t_satData*> im(epoData->satData);
+  while (im.hasNext()) {
+    im.next();
+    t_satData* satData = im.value();
+    cmpEle(satData);
+    if (satData->eleSat < OPT->_minEle) {
+      delete satData;
+      im.remove();
+    }
+  }
+
+  return success;
+}
+
+// Computed Value
+////////////////////////////////////////////////////////////////////////////
+double t_pppFilter::cmpValue(t_satData* satData, bool phase) {
+
+  Tracer tracer("t_pppFilter::cmpValue");
+
+  ColumnVector xRec(3);
+  xRec(1) = x();
+  xRec(2) = y();
+  xRec(3) = z();
+
+  double rho0 = (satData->xx - xRec).norm_Frobenius();
+  double dPhi = t_CST::omega * rho0 / t_CST::c; 
+
+  xRec(1) = x() * cos(dPhi) - y() * sin(dPhi); 
+  xRec(2) = y() * cos(dPhi) + x() * sin(dPhi); 
+  xRec(3) = z();
+
+  xRec += _tides->displacement(_time, xRec);
+
+  satData->rho = (satData->xx - xRec).norm_Frobenius();
+
+  double tropDelay = delay_saast(satData->eleSat) + 
+                     trp() / sin(satData->eleSat);
+
+  double wind = 0.0;
+  if (phase) {
+    wind = windUp(satData->prn, satData->xx, xRec) * satData->lambda3;
+  }
+
+  double offset = 0.0;
+  t_frequency::type frqA = t_frequency::G1;
+  t_frequency::type frqB = t_frequency::G2;
+  if      (satData->prn[0] == 'R') {
+    offset = Glonass_offset();
+    frqA = t_frequency::R1;
+    frqB = t_frequency::R2;
+  }
+  else if (satData->prn[0] == 'E') {
+    offset = Galileo_offset();
+    //frqA = t_frequency::E1; as soon as available
+    //frqB = t_frequency::E5; -"-
+  }
+  else if (satData->prn[0] == 'C') {
+    offset = Bds_offset();
+    //frqA = t_frequency::C2; as soon as available
+    //frqB = t_frequency::C7; -"-
+  }
+  double phaseCenter = 0.0;
+  if (_antex) { 
+    bool found;
+    phaseCenter = satData->lkA * _antex->rcvCorr(OPT->_antNameRover, frqA,
+                                                 satData->eleSat, satData->azSat,
+                                                 found)
+                + satData->lkB * _antex->rcvCorr(OPT->_antNameRover, frqB,
+                                                 satData->eleSat, satData->azSat,
+                                                 found);
+    if (!found) {
+      LOG << "ANTEX: antenna >" << OPT->_antNameRover << "< not found\n";
+    }
+  }
+
+  double antennaOffset = 0.0;
+  double cosa = cos(satData->azSat);
+  double sina = sin(satData->azSat);
+  double cose = cos(satData->eleSat);
+  double sine = sin(satData->eleSat);
+  antennaOffset = -OPT->_neuEccRover(1) * cosa*cose 
+                  -OPT->_neuEccRover(2) * sina*cose 
+                  -OPT->_neuEccRover(3) * sine;
+
+  return satData->rho + phaseCenter + antennaOffset + clk() 
+                      + offset - satData->clk + tropDelay + wind;
+}
+
+// Tropospheric Model (Saastamoinen)
+////////////////////////////////////////////////////////////////////////////
+double t_pppFilter::delay_saast(double Ele) {
+
+  Tracer tracer("t_pppFilter::delay_saast");
+
+  double xyz[3]; 
+  xyz[0] = x();
+  xyz[1] = y();
+  xyz[2] = z();
+  double ell[3]; 
+  xyz2ell(xyz, ell);
+  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 h_km = height / 1000.0;
   
-  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) {
-      _QFlt[ii][ii] = par1->sigma0() * par1->sigma0(); // new parameter
-    }
-    else {
-      _QFlt[ii][ii] = QFltOld[iOld][iOld] + par1->noise() * par1->noise();
-      _xFlt[ii]     = xFltOld[iOld];
-      for (unsigned jj = 0; jj < ii; jj++) {
-        const t_pppParam* par2 = params[jj];
-        int               jOld = par2->indexOld();
-        if (jOld >= 0) {
-          _QFlt[ii][jj] = QFltOld(iOld+1,jOld+1);
+  if (h_km < 0.0) h_km = 0.0;
+  if (h_km > 5.0) h_km = 5.0;
+  int    ii   = int(h_km + 1);
+  double href = ii - 1;
+  
+  double bCor[6]; 
+  bCor[0] = 1.156;
+  bCor[1] = 1.006;
+  bCor[2] = 0.874;
+  bCor[3] = 0.757;
+  bCor[4] = 0.654;
+  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)));
+}
+
+// Prediction Step of the Filter
+////////////////////////////////////////////////////////////////////////////
+void t_pppFilter::predict(int iPhase, t_epoData* epoData) {
+
+  Tracer tracer("t_pppFilter::predict");
+
+  if (iPhase == 0) {
+
+    const double maxSolGap = 60.0;
+
+    bool firstCrd = false;
+    if (!_lastTimeOK.valid() || (maxSolGap > 0.0 && _time - _lastTimeOK > maxSolGap)) {
+      firstCrd = true;
+      _startTime = epoData->tt;
+      reset();
+    }
+    
+    // Use different white noise for Quick-Start mode
+    // ----------------------------------------------
+    double sigCrdP_used = OPT->_noiseCrd(1);
+    if ( OPT->_seedingTime > 0.0 && OPT->_seedingTime > (epoData->tt - _startTime) ) {
+      sigCrdP_used   = 0.0;
+    }
+
+    // Predict Parameter values, add white noise
+    // -----------------------------------------
+    for (int iPar = 1; iPar <= _params.size(); iPar++) {
+      t_pppParam* pp = _params[iPar-1];
+    
+      // Coordinates
+      // -----------
+      if      (pp->type == t_pppParam::CRD_X) {
+        if (firstCrd) {
+          if (OPT->xyzAprRoverSet()) {
+            pp->xx = OPT->_xyzAprRover[0];
+          }
+          else {
+            pp->xx = _xcBanc(1);
+          }
         }
-      }
-    }
-  }
-
-  predictCovCrdPart(QFltOld);
-
-  // Process Satellite Systems separately
-  // ------------------------------------
-  for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
-    char system = OPT->systems()[iSys];
-    vector<t_pppSatObs*> 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) != success ) {
-      return failure;
-    }
-  }
-   
-  cmpDOP(allObs);
-
-  _parlist->printResult(_epoTime, _QFlt, _xFlt);
-
-  return success;
-}
-
-// Process Selected LCs
-////////////////////////////////////////////////////////////////////////////
-t_irc t_pppFilter::processSystem(const vector<t_lc::type>& LCs, 
-                                 const vector<t_pppSatObs*>& obsVector) {
-
-  LOG.setf(ios::fixed);
-
-  // Detect Cycle Slips
-  // ------------------
-  if (detectCycleSlips(LCs, obsVector) != success) {
-    return failure;
-  }
-
-  ColumnVector            xSav       = _xFlt;
-  SymmetricMatrix         QSav       = _QFlt;
-  string                  epoTimeStr = string(_epoTime);
-  const vector<t_pppParam*>& params     = _parlist->params();
-  unsigned                maxObs     = obsVector.size() * LCs.size();
-    
+        _QQ(iPar,iPar) += sigCrdP_used * sigCrdP_used;
+      }
+      else if (pp->type == t_pppParam::CRD_Y) {
+        if (firstCrd) {
+          if (OPT->xyzAprRoverSet()) {
+            pp->xx = OPT->_xyzAprRover[1];
+          }
+          else {
+            pp->xx = _xcBanc(2);
+          }
+        }
+        _QQ(iPar,iPar) += sigCrdP_used * sigCrdP_used;
+      }
+      else if (pp->type == t_pppParam::CRD_Z) {
+        if (firstCrd) {
+          if (OPT->xyzAprRoverSet()) {
+            pp->xx = OPT->_xyzAprRover[2];
+          }
+          else {
+            pp->xx = _xcBanc(3);
+          }
+        }
+        _QQ(iPar,iPar) += sigCrdP_used * sigCrdP_used;
+      }   
+    
+      // Receiver Clocks
+      // ---------------
+      else if (pp->type == t_pppParam::RECCLK) {
+        pp->xx = _xcBanc(4);
+        for (int jj = 1; jj <= _params.size(); jj++) {
+          _QQ(iPar, jj) = 0.0;
+        }
+        _QQ(iPar,iPar) = OPT->_noiseClk * OPT->_noiseClk;
+      }
+    
+      // Tropospheric Delay
+      // ------------------
+      else if (pp->type == t_pppParam::TROPO) {
+        _QQ(iPar,iPar) += OPT->_noiseTrp * OPT->_noiseTrp;
+      }
+    
+      // Glonass Offset
+      // --------------
+      else if (pp->type == t_pppParam::GLONASS_OFFSET) {
+        pp->xx = 0.0;
+        for (int jj = 1; jj <= _params.size(); jj++) {
+          _QQ(iPar, jj) = 0.0;
+        }
+        _QQ(iPar,iPar) = 1000.0 * 1000.0;
+      }
+
+      // Galileo Offset
+      // --------------
+      else if (pp->type == t_pppParam::GALILEO_OFFSET) {
+        _QQ(iPar,iPar) += 0.1 * 0.1;
+      }
+
+      // BDS Offset
+      // ----------
+      else if (pp->type == t_pppParam::BDS_OFFSET) {
+        _QQ(iPar,iPar) = 1000.0 * 1000.0;    //TODO: TEST
+      }
+    }
+  }
+
+  // Add New Ambiguities if necessary
+  // --------------------------------
+  if (OPT->ambLCs('G').size() || OPT->ambLCs('R').size() ||
+      OPT->ambLCs('E').size() || OPT->ambLCs('C').size()) {
+
+    // Make a copy of QQ and xx, set parameter indices
+    // -----------------------------------------------
+    SymmetricMatrix QQ_old = _QQ;
+    
+    for (int iPar = 1; iPar <= _params.size(); iPar++) {
+      _params[iPar-1]->index_old = _params[iPar-1]->index;
+      _params[iPar-1]->index     = 0;
+    }
+    
+    // Remove Ambiguity Parameters without observations
+    // ------------------------------------------------
+    int iPar = 0;
+    QMutableVectorIterator<t_pppParam*> im(_params);
+    while (im.hasNext()) {
+      t_pppParam* par = im.next();
+      bool removed = false;
+      if (par->type == t_pppParam::AMB_L3) {
+        if (epoData->satData.find(par->prn) == epoData->satData.end()) {
+          removed = true;
+          delete par;
+          im.remove();
+        }
+      }
+      if (! removed) {
+        ++iPar;
+        par->index = iPar;
+      }
+    }
+    
+    // Add new ambiguity parameters
+    // ----------------------------
+    QMapIterator<QString, t_satData*> it(epoData->satData);
+    while (it.hasNext()) {
+      it.next();
+      t_satData* satData = it.value();
+      addAmb(satData);
+    }
+    
+    int nPar = _params.size();
+    _QQ.ReSize(nPar); _QQ = 0.0;
+    for (int i1 = 1; i1 <= nPar; i1++) {
+      t_pppParam* p1 = _params[i1-1];
+      if (p1->index_old != 0) {
+        _QQ(p1->index, p1->index) = QQ_old(p1->index_old, p1->index_old);
+        for (int i2 = 1; i2 <= nPar; i2++) {
+          t_pppParam* p2 = _params[i2-1];
+          if (p2->index_old != 0) {
+            _QQ(p1->index, p2->index) = QQ_old(p1->index_old, p2->index_old);
+          }
+        }
+      }
+    }
+    
+    for (int ii = 1; ii <= nPar; ii++) {
+      t_pppParam* par = _params[ii-1];
+      if (par->index_old == 0) {
+        _QQ(par->index, par->index) = OPT->_aprSigAmb * OPT->_aprSigAmb;
+      }
+      par->index_old = par->index;
+    }
+  }
+}
+
+// Update Step of the Filter (currently just a single-epoch solution)
+////////////////////////////////////////////////////////////////////////////
+t_irc t_pppFilter::update(t_epoData* epoData) {
+
+  Tracer tracer("t_pppFilter::update");
+
+  _time = epoData->tt; // current epoch time
+
+  if (OPT->useOrbClkCorr()) {
+    LOG << "Precise Point Positioning of Epoch " << _time.datestr() <<  "_" << _time.timestr(1)
+        << "\n---------------------------------------------------------------\n";
+  }
+  else {
+    LOG << "Single Point Positioning of Epoch " << _time.datestr() <<  "_" << _time.timestr(1)
+        << "\n--------------------------------------------------------------\n";
+  }
+
   // Outlier Detection Loop
   // ----------------------
-  for (unsigned iOutlier = 0; iOutlier < maxObs; iOutlier++) {
-
-    if (iOutlier > 0) {
-      _xFlt = xSav;
-      _QFlt = QSav;
-    }
-
-    // First-Design Matrix, Terms Observed-Computed, Weight Matrix
-    // -----------------------------------------------------------
-    Matrix                AA(maxObs, _parlist->nPar());
-    ColumnVector          ll(maxObs);
-    DiagonalMatrix        PP(maxObs); PP = 0.0;
-    
-    int iObs = -1;
-    vector<t_pppSatObs*> usedObs;
-    vector<t_lc::type>   usedTypes;
-    for (unsigned ii = 0; ii < obsVector.size(); ii++) {
-      t_pppSatObs* obs = obsVector[ii];
-      if (!obs->outlier()) {
-        for (unsigned jj = 0; jj < LCs.size(); jj++) {
-          const t_lc::type tLC = LCs[jj];
-          ++iObs;
-          usedObs.push_back(obs);
-          usedTypes.push_back(tLC);
-          for (unsigned iPar = 0; iPar < params.size(); iPar++) {
-            const t_pppParam* par = params[iPar];
-            AA[iObs][iPar] = par->partial(_epoTime, obs, tLC);
+  if (update_p(epoData) != success) {
+    return failure;
+  }
+
+  // Set Solution Vector
+  // -------------------
+  LOG.setf(ios::fixed);
+  QVectorIterator<t_pppParam*> itPar(_params);
+  while (itPar.hasNext()) {
+    t_pppParam* par = itPar.next();
+    if      (par->type == t_pppParam::RECCLK) {
+      LOG << "\n    clk     = " << setw(10) << setprecision(3) << par->xx 
+          << " +- " << setw(6) << setprecision(3) 
+          << sqrt(_QQ(par->index,par->index));
+    }
+    else if (par->type == t_pppParam::AMB_L3) {
+      ++par->numEpo;
+      LOG << "\n    amb " << par->prn.mid(0,3).toAscii().data() << " = "
+          << setw(10) << setprecision(3) << par->xx 
+          << " +- " << setw(6) << setprecision(3) 
+          << sqrt(_QQ(par->index,par->index))
+          << "   nEpo = " << par->numEpo;
+    }
+    else if (par->type == t_pppParam::TROPO) {
+      double aprTrp = delay_saast(M_PI/2.0);
+      LOG << "\n    trp     = " << par->prn.mid(0,3).toAscii().data()
+          << setw(7) << setprecision(3) << aprTrp << " "
+          << setw(6) << setprecision(3) << showpos << par->xx << noshowpos
+          << " +- " << setw(6) << setprecision(3) 
+          << sqrt(_QQ(par->index,par->index));
+    }
+    else if (par->type == t_pppParam::GLONASS_OFFSET) {
+      LOG << "\n    offGlo  = " << setw(10) << setprecision(3) << par->xx 
+          << " +- " << setw(6) << setprecision(3) 
+          << sqrt(_QQ(par->index,par->index));
+    }
+    else if (par->type == t_pppParam::GALILEO_OFFSET) {
+      LOG << "\n    offGal  = " << setw(10) << setprecision(3) << par->xx 
+          << " +- " << setw(6) << setprecision(3) 
+          << sqrt(_QQ(par->index,par->index));
+    }
+    else if (par->type == t_pppParam::BDS_OFFSET) {
+      LOG << "\n    offBds  = " << setw(10) << setprecision(3) << par->xx
+          << " +- " << setw(6) << setprecision(3)
+          << sqrt(_QQ(par->index,par->index));
+    }
+  }
+
+  LOG << endl << endl;
+
+  // Compute dilution of precision
+  // -----------------------------
+  cmpDOP(epoData);
+
+  // Final Message (both log file and screen)
+  // ----------------------------------------
+  LOG << OPT->_roverName << "  PPP " 
+      << epoData->tt.timestr(1) << " " << epoData->sizeAll() << " "
+      << setw(14) << setprecision(3) << x()                  << " +- "
+      << setw(6)  << setprecision(3) << sqrt(_QQ(1,1))       << " "
+      << setw(14) << setprecision(3) << y()                  << " +- "
+      << setw(6)  << setprecision(3) << sqrt(_QQ(2,2))       << " "
+      << setw(14) << setprecision(3) << z()                  << " +- "
+      << setw(6)  << setprecision(3) << sqrt(_QQ(3,3));
+
+  // NEU Output
+  // ----------
+  if (OPT->xyzAprRoverSet()) {
+    double xyz[3];
+    xyz[0] = x() - OPT->_xyzAprRover[0];
+    xyz[1] = y() - OPT->_xyzAprRover[1];
+    xyz[2] = z() - OPT->_xyzAprRover[2];
+
+    double ellRef[3];
+    xyz2ell(OPT->_xyzAprRover.data(), ellRef);
+    xyz2neu(ellRef, xyz, _neu.data());
+
+    LOG << "  NEU "
+        << setw(8) << setprecision(3) << _neu[0] << " "
+        << setw(8) << setprecision(3) << _neu[1] << " "
+        << setw(8) << setprecision(3) << _neu[2] << endl << endl;
+  }
+  else {
+    LOG << endl << endl;
+  }
+
+  _lastTimeOK = _time; // remember time of last successful update
+  return success;
+}
+
+// Outlier Detection
+////////////////////////////////////////////////////////////////////////////
+QString t_pppFilter::outlierDetection(int iPhase, const ColumnVector& vv,
+                                   QMap<QString, t_satData*>& satData) {
+
+  Tracer tracer("t_pppFilter::outlierDetection");
+
+  QString prnGPS;
+  QString prnGlo;
+  double  maxResGPS = 0.0; // all the other systems except GLONASS
+  double  maxResGlo = 0.0; // GLONASS
+  findMaxRes(vv, satData, prnGPS, prnGlo, maxResGPS, maxResGlo);
+
+  if      (iPhase == 1) {
+    if      (maxResGlo > 2.98 * OPT->_maxResL1) { 
+      LOG << "Outlier Phase " << prnGlo.mid(0,3).toAscii().data() << ' ' << maxResGlo << endl;
+      return prnGlo;
+    }
+    else if (maxResGPS > 2.98 * OPT->_maxResL1) { 
+      LOG << "Outlier Phase " << prnGPS.mid(0,3).toAscii().data() << ' ' << maxResGPS << endl;
+      return prnGPS;
+    }
+  }
+  else if (iPhase == 0 && maxResGPS > 2.98 * OPT->_maxResC1) {
+    LOG << "Outlier Code  " << prnGPS.mid(0,3).toAscii().data() << ' ' << maxResGPS << endl;
+    return prnGPS;
+  }
+
+  return QString();
+}
+
+// Phase Wind-Up Correction
+///////////////////////////////////////////////////////////////////////////
+double t_pppFilter::windUp(const QString& prn, const ColumnVector& rSat,
+                        const ColumnVector& rRec) {
+
+  Tracer tracer("t_pppFilter::windUp");
+
+  double Mjd = _time.mjd() + _time.daysec() / 86400.0;
+
+  // First time - initialize to zero
+  // -------------------------------
+  if (!_windUpTime.contains(prn)) {
+    _windUpSum[prn]  = 0.0;
+  }
+
+  // Compute the correction for new time
+  // -----------------------------------
+  if (!_windUpTime.contains(prn) || _windUpTime[prn] != Mjd) {
+    _windUpTime[prn] = Mjd; 
+
+    // Unit Vector GPS Satellite --> Receiver
+    // --------------------------------------
+    ColumnVector rho = rRec - rSat;
+    rho /= rho.norm_Frobenius();
+    
+    // GPS Satellite unit Vectors sz, sy, sx
+    // -------------------------------------
+    ColumnVector sz = -rSat / rSat.norm_Frobenius();
+
+    ColumnVector xSun = t_astro::Sun(Mjd);
+    xSun /= xSun.norm_Frobenius();
+
+    ColumnVector sy = crossproduct(sz, xSun);
+    ColumnVector sx = crossproduct(sy, sz);
+
+    // Effective Dipole of the GPS Satellite Antenna
+    // ---------------------------------------------
+    ColumnVector dipSat = sx - rho * DotProduct(rho,sx) 
+                                                - crossproduct(rho, sy);
+    
+    // Receiver unit Vectors rx, ry
+    // ----------------------------
+    ColumnVector rx(3);
+    ColumnVector ry(3);
+
+    double recEll[3]; xyz2ell(rRec.data(), recEll) ;
+    double neu[3];
+    
+    neu[0] = 1.0;
+    neu[1] = 0.0;
+    neu[2] = 0.0;
+    neu2xyz(recEll, neu, rx.data());
+    
+    neu[0] =  0.0;
+    neu[1] = -1.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);
+    
+    // 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 dphi = acos(alpha) / 2.0 / M_PI;  // in cycles
+    
+    if ( DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0 ) {
+      dphi = -dphi;
+    }
+
+    _windUpSum[prn] = floor(_windUpSum[prn] - dphi + 0.5) + dphi;
+  }
+
+  return _windUpSum[prn];  
+}
+
+// 
+///////////////////////////////////////////////////////////////////////////
+void t_pppFilter::cmpEle(t_satData* satData) {
+  Tracer tracer("t_pppFilter::cmpEle");
+  ColumnVector rr = satData->xx - _xcBanc.Rows(1,3);
+  double       rho = rr.norm_Frobenius();
+
+  double neu[3];
+  xyz2neu(_ellBanc.data(), rr.data(), neu);
+
+  satData->eleSat = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / rho );
+  if (neu[2] < 0) {
+    satData->eleSat *= -1.0;
+  }
+  satData->azSat  = atan2(neu[1], neu[0]);
+}
+
+// 
+///////////////////////////////////////////////////////////////////////////
+void t_pppFilter::addAmb(t_satData* satData) {
+  Tracer tracer("t_pppFilter::addAmb");
+  if (!OPT->ambLCs(satData->system()).size()){
+    return;
+  }
+  bool    found = false;
+  for (int iPar = 1; iPar <= _params.size(); iPar++) {
+    if (_params[iPar-1]->type == t_pppParam::AMB_L3 && 
+        _params[iPar-1]->prn == satData->prn) {
+      found = true;
+      break;
+    }
+  }
+  if (!found) {
+    t_pppParam* par = new t_pppParam(t_pppParam::AMB_L3, 
+                                 _params.size()+1, satData->prn);
+    _params.push_back(par);
+    par->xx = satData->L3 - cmpValue(satData, true);
+  }
+}
+
+// 
+///////////////////////////////////////////////////////////////////////////
+void t_pppFilter::addObs(int iPhase, unsigned& iObs, t_satData* satData,
+                      Matrix& AA, ColumnVector& ll, DiagonalMatrix& PP) {
+
+  Tracer tracer("t_pppFilter::addObs");
+
+  const double ELEWGHT = 20.0;
+  double ellWgtCoef = 1.0;
+  double eleD = satData->eleSat * 180.0 / M_PI; 
+  if (eleD < ELEWGHT) {
+    ellWgtCoef = 1.5 - 0.5 / (ELEWGHT - 10.0) * (eleD - 10.0);
+  }
+
+  // Remember Observation Index
+  // --------------------------
+  ++iObs;
+  satData->obsIndex = iObs;
+
+  // Phase Observations
+  // ------------------
+
+  if (iPhase == 1) {
+    ll(iObs)      = satData->L3 - cmpValue(satData, true);
+    double sigL3 = 2.98 * OPT->_sigmaL1;
+    if (satData->system() == 'R') {
+      sigL3 *= GLONASS_WEIGHT_FACTOR;
+    }
+    if  (satData->system() == 'C') {
+      sigL3 *= BDS_WEIGHT_FACTOR;
+    }
+    PP(iObs,iObs) = 1.0 / (sigL3 * sigL3) / (ellWgtCoef * ellWgtCoef);
+    for (int iPar = 1; iPar <= _params.size(); iPar++) {
+      if (_params[iPar-1]->type == t_pppParam::AMB_L3 &&
+          _params[iPar-1]->prn  == satData->prn) {
+        ll(iObs) -= _params[iPar-1]->xx;
+      } 
+      AA(iObs, iPar) = _params[iPar-1]->partial(satData, true);
+    }
+  }
+
+  // Code Observations
+  // -----------------
+  else {
+    double sigP3 = 2.98 * OPT->_sigmaC1;
+    if  (satData->system() == 'C') {
+      sigP3 *= BDS_WEIGHT_FACTOR;
+    }
+    ll(iObs)      = satData->P3 - cmpValue(satData, false);
+    PP(iObs,iObs) = 1.0 / (sigP3 * sigP3) / (ellWgtCoef * ellWgtCoef);
+    for (int iPar = 1; iPar <= _params.size(); iPar++) {
+      AA(iObs, iPar) = _params[iPar-1]->partial(satData, false);
+    }
+  }
+}
+
+// 
+///////////////////////////////////////////////////////////////////////////
+QByteArray t_pppFilter::printRes(int iPhase, const ColumnVector& vv, 
+                              const QMap<QString, t_satData*>& satDataMap) {
+
+  Tracer tracer("t_pppFilter::printRes");
+
+  ostringstream str;
+  str.setf(ios::fixed);
+  bool useObs;
+  QMapIterator<QString, t_satData*> it(satDataMap);
+  while (it.hasNext()) {
+    it.next();
+    t_satData* satData = it.value();
+    (iPhase == 0) ? useObs = OPT->codeLCs(satData->system()).size() :
+                    useObs = OPT->ambLCs(satData->system()).size();
+    if (satData->obsIndex != 0 && useObs) {
+      str << _time.timestr(1)
+          << " RES " << satData->prn.mid(0,3).toAscii().data()
+          << (iPhase ? "   L3 " : "   P3 ")
+          << setw(9) << setprecision(4) << vv(satData->obsIndex) << endl;
+    }
+  }
+
+  return QByteArray(str.str().c_str());
+}
+
+// 
+///////////////////////////////////////////////////////////////////////////
+void t_pppFilter::findMaxRes(const ColumnVector& vv,
+                          const QMap<QString, t_satData*>& satData,
+                          QString& prnGPS, QString& prnGlo, 
+                          double& maxResGPS, double& maxResGlo) { 
+
+  Tracer tracer("t_pppFilter::findMaxRes");
+
+  maxResGPS  = 0.0;
+  maxResGlo  = 0.0;
+
+  QMapIterator<QString, t_satData*> it(satData);
+  while (it.hasNext()) {
+    it.next();
+    t_satData* satData = it.value();
+    if (satData->obsIndex != 0) {
+      QString prn = satData->prn;
+      if (prn[0] == 'R') {
+        if (fabs(vv(satData->obsIndex)) > maxResGlo) {
+          maxResGlo = fabs(vv(satData->obsIndex));
+          prnGlo    = prn;
+        }
+      }
+      else {
+        if (fabs(vv(satData->obsIndex)) > maxResGPS) {
+          maxResGPS = fabs(vv(satData->obsIndex));
+          prnGPS    = prn;
+        }
+      }
+    }
+  }
+}
+ 
+// Update Step (private - loop over outliers)
+////////////////////////////////////////////////////////////////////////////
+t_irc t_pppFilter::update_p(t_epoData* epoData) {
+
+  Tracer tracer("t_pppFilter::update_p");
+
+  // Save Variance-Covariance Matrix, and Status Vector
+  // --------------------------------------------------
+  rememberState(epoData);
+
+  QString lastOutlierPrn;
+
+  // Try with all satellites, then with all minus one, etc.
+  // ------------------------------------------------------
+  while (selectSatellites(lastOutlierPrn, epoData->satData) == success) {
+
+    QByteArray strResCode;
+    QByteArray strResPhase;
+
+    // Bancroft Solution
+    // -----------------
+    if (cmpBancroft(epoData) != success) {
+      break;
+    }
+
+    // First update using code observations, then phase observations
+    // -------------------------------------------------------------      
+    bool usePhase = OPT->ambLCs('G').size() || OPT->ambLCs('R').size() ||
+                    OPT->ambLCs('E').size() || OPT->ambLCs('C').size() ;
+
+    for (int iPhase = 0; iPhase <= (usePhase ? 1 : 0); iPhase++) {
+
+      // Status Prediction
+      // -----------------
+      predict(iPhase, epoData);
+      
+      // Create First-Design Matrix
+      // --------------------------
+      unsigned nPar = _params.size();
+      unsigned nObs = 0;
+      nObs = epoData->sizeAll();
+      bool useObs = false;
+      char sys[] ={'G', 'R', 'E', 'C'};
+      for (unsigned ii = 0; ii < sizeof(sys); ii++) {
+        const char s = sys[ii];
+        (iPhase == 0) ? useObs = OPT->codeLCs(s).size() : useObs = OPT->ambLCs(s).size();
+        if (!useObs) {
+          nObs -= epoData->sizeSys(s);
+        }
+      }
+      
+      // Prepare first-design Matrix, vector observed-computed
+      // -----------------------------------------------------
+      Matrix          AA(nObs, nPar);  // first design matrix
+      ColumnVector    ll(nObs);        // tems observed-computed
+      DiagonalMatrix  PP(nObs); PP = 0.0;
+      
+      unsigned iObs = 0;
+      QMapIterator<QString, t_satData*> it(epoData->satData);
+      while (it.hasNext()) {
+        it.next();
+        t_satData* satData = it.value();
+        QString prn = satData->prn;
+        (iPhase == 0) ? useObs = OPT->codeLCs(satData->system()).size() :
+                        useObs = OPT->ambLCs(satData->system()).size();
+        if (useObs) {
+          addObs(iPhase, iObs, satData, AA, ll, PP);
+        } else {
+          satData->obsIndex = 0;
+        }
+      }
+
+      // Compute Filter Update
+      // ---------------------
+      ColumnVector dx(nPar); dx = 0.0;
+      kalman(AA, ll, PP, _QQ, dx);
+      ColumnVector vv = ll - AA * dx;
+      
+      // Print Residuals
+      // ---------------
+      if      (iPhase == 0) {
+        strResCode  = printRes(iPhase, vv, epoData->satData);
+      }
+      else {
+        strResPhase = printRes(iPhase, vv, epoData->satData);
+      }
+
+      // Check the residuals
+      // -------------------
+      lastOutlierPrn = outlierDetection(iPhase, vv, epoData->satData);
+
+      // No Outlier Detected
+      // -------------------
+      if (lastOutlierPrn.isEmpty()) {
+
+        QVectorIterator<t_pppParam*> itPar(_params);
+        while (itPar.hasNext()) {
+          t_pppParam* par = itPar.next();
+          par->xx += dx(par->index);
+        }
+
+        if (!usePhase || iPhase == 1) {
+          if (_outlierGPS.size() > 0 || _outlierGlo.size() > 0) {
+            LOG << "Neglected PRNs: ";
+            if (!_outlierGPS.isEmpty()) {
+              LOG << _outlierGPS.last().mid(0,3).toAscii().data() << ' ';
+            }
+            QStringListIterator itGlo(_outlierGlo);
+            while (itGlo.hasNext()) {
+              QString prn = itGlo.next();
+              LOG << prn.mid(0,3).toAscii().data() << ' ';
+            }
           }
-          ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1));
-          PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC));
+          LOG << strResCode.data() << strResPhase.data();
+
+          return success;
         }
       }
-    }
-
-    // Check number of observations, truncate matrices
-    // -----------------------------------------------
-    if (iObs == -1) {
-      return failure;
-    }
-    AA = AA.Rows(1, iObs+1);
-    ll = ll.Rows(1, iObs+1);
-    PP = PP.SymSubMatrix(1, iObs+1);
-
-    // Kalman update step
-    // ------------------
-    kalman(AA, ll, PP, _QFlt, _xFlt);
-
-    // Check Residuals
-    // ---------------
-    ColumnVector vv = AA * _xFlt - ll;
-    double     maxOutlier      = 0.0;
-    int        maxOutlierIndex = -1;
-    t_lc::type maxOutlierLC    = t_lc::dummy;
-    for (unsigned ii = 0; ii < usedObs.size(); ii++) {
-      const t_lc::type tLC = usedTypes[ii];
-      double res = fabs(vv[ii]);
-      if (res > usedObs[ii]->maxRes(tLC)) {
-        if (res > fabs(maxOutlier)) {
-          maxOutlier      = vv[ii];
-          maxOutlierIndex = ii;
-          maxOutlierLC    = tLC;
-        }
-      }
-    }
-
-    // Mark outlier or break outlier detection loop
-    // --------------------------------------------
-    if (maxOutlierIndex > -1) {
-      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;
-      for (unsigned iPar = 0; iPar < params.size(); iPar++) {
-        t_pppParam* hlp = params[iPar];
-        if (hlp->type() == t_pppParam::amb && hlp->prn()  == obs->prn() && 
-            hlp->tLC()  == usedTypes[maxOutlierIndex]) {
-          par = hlp;
-        }
-      }
-      if (par) {
-        if (par->ambResetCandidate()) {
-          resetAmb(par->prn(), obsVector, &QSav, &xSav);
-        }
-        else {
-          par->setAmbResetCandidate();
-          obs->setOutlier();
-        }
-      }
+
+      // Outlier Found
+      // -------------
       else {
-        obs->setOutlier();
-      }
-    }
-
-    // Print Residuals
-    // ---------------
+        restoreState(epoData);
+        break;
+      }
+
+    } // for iPhase
+
+  } // while selectSatellites
+
+  restoreState(epoData);
+  return failure;
+}
+
+// Remeber Original State Vector and Variance-Covariance Matrix
+////////////////////////////////////////////////////////////////////////////
+void t_pppFilter::rememberState(t_epoData* epoData) {
+
+  _QQ_sav = _QQ;
+
+  QVectorIterator<t_pppParam*> itSav(_params_sav);
+  while (itSav.hasNext()) {
+    t_pppParam* par = itSav.next();
+    delete par;
+  }
+  _params_sav.clear();
+
+  QVectorIterator<t_pppParam*> it(_params);
+  while (it.hasNext()) {
+    t_pppParam* par = it.next();
+    _params_sav.push_back(new t_pppParam(*par));
+  }
+
+  _epoData_sav->deepCopy(epoData);
+}
+
+// Restore Original State Vector and Variance-Covariance Matrix
+////////////////////////////////////////////////////////////////////////////
+void t_pppFilter::restoreState(t_epoData* epoData) {
+
+  _QQ = _QQ_sav;
+
+  QVectorIterator<t_pppParam*> it(_params);
+  while (it.hasNext()) {
+    t_pppParam* par = it.next();
+    delete par;
+  }
+  _params.clear();
+
+  QVectorIterator<t_pppParam*> itSav(_params_sav);
+  while (itSav.hasNext()) {
+    t_pppParam* par = itSav.next();
+    _params.push_back(new t_pppParam(*par));
+  }
+
+  epoData->deepCopy(_epoData_sav);
+}
+
+// 
+////////////////////////////////////////////////////////////////////////////
+t_irc t_pppFilter::selectSatellites(const QString& lastOutlierPrn, 
+                                 QMap<QString, t_satData*>& satData) {
+
+  // First Call 
+  // ----------
+  if (lastOutlierPrn.isEmpty()) {
+    _outlierGPS.clear();
+    _outlierGlo.clear();
+    return success;
+  }
+
+  // Second and next trials
+  // ----------------------
+  else {
+
+    if (lastOutlierPrn[0] == 'R') {
+      _outlierGlo << lastOutlierPrn;
+    }
+
+    // Remove all Glonass Outliers
+    // ---------------------------
+    QStringListIterator it(_outlierGlo);
+    while (it.hasNext()) {
+      QString prn = it.next();
+      if (satData.contains(prn)) {
+        delete satData.take(prn);
+      }
+    }
+
+    if (lastOutlierPrn[0] == 'R') {
+      _outlierGPS.clear();
+      return success;
+    }
+
+    // GPS Outlier appeared for the first time - try to delete it
+    // ----------------------------------------------------------
+    if (_outlierGPS.indexOf(lastOutlierPrn) == -1) {
+      _outlierGPS << lastOutlierPrn;
+      if (satData.contains(lastOutlierPrn)) {
+        delete satData.take(lastOutlierPrn);
+      }
+      return success;
+    }
+
+  }
+
+  return failure;
+}
+
+// 
+////////////////////////////////////////////////////////////////////////////
+double lorentz(const ColumnVector& aa, const ColumnVector& bb) {
+  return aa(1)*bb(1) +  aa(2)*bb(2) +  aa(3)*bb(3) -  aa(4)*bb(4);
+}
+
+// 
+////////////////////////////////////////////////////////////////////////////
+void t_pppFilter::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 {
-      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;
-          }
-        }
-      }
-      break;
-    }
-  }
-
-  return success;
-}
-
-// Cycle-Slip Detection
-////////////////////////////////////////////////////////////////////////////
-t_irc t_pppFilter::detectCycleSlips(const vector<t_lc::type>& LCs, 
-                                    const vector<t_pppSatObs*>& obsVector) {
-
-  const double            SLIP       = 20.0;  // slip threshold
-  string                  epoTimeStr = string(_epoTime);
-  const vector<t_pppParam*>& params     = _parlist->params();
-
-  for (unsigned ii = 0; ii < LCs.size(); ii++) {
-    const t_lc::type& tLC = LCs[ii];
-    if (t_lc::includesPhase(tLC)) {
-      for (unsigned iObs = 0; iObs < obsVector.size(); iObs++) {
-        const t_pppSatObs* obs = obsVector[iObs];
-
-        // Check set Slips and Jump Counters 
-        // ---------------------------------
-        bool slip = false;
-
-        if (obs->slip()) {
-          LOG << "cycle slip set (obs)" << endl;;
-          slip = true;
-        }
-
-        if (_slips[obs->prn()]._obsSlipCounter != -1 &&
-            _slips[obs->prn()]._obsSlipCounter != obs->slipCounter()) {
-          LOG << "cycle slip set (obsSlipCounter)" << endl;
-          slip = true;
-        }
-        _slips[obs->prn()]._obsSlipCounter = obs->slipCounter();
-
-        if (_slips[obs->prn()]._biasJumpCounter != -1 &&
-            _slips[obs->prn()]._biasJumpCounter != obs->biasJumpCounter()) {
-          LOG << "cycle slip set (biasJumpCounter)" << endl;
-          slip = true;
-        }
-        _slips[obs->prn()]._biasJumpCounter = obs->biasJumpCounter();
-
-        // Slip Set
-        // --------  
-        if (slip) {
-          resetAmb(obs->prn(), obsVector);
-        }
-  
-        // Check Pre-Fit Residuals
-        // -----------------------
-        else {
-          ColumnVector AA(params.size());
-          for (unsigned iPar = 0; iPar < params.size(); iPar++) {
-            const t_pppParam* par = params[iPar];
-            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);
-          }
-        }
-      }
-    }
-  }
-
-  return success;
-}
-
-// Reset Ambiguity Parameter (cycle slip)
-////////////////////////////////////////////////////////////////////////////
-t_irc t_pppFilter::resetAmb(t_prn prn, const vector<t_pppSatObs*>& obsVector,
-                            SymmetricMatrix* QSav, ColumnVector* xSav) {
-  t_irc irc = failure;
-  vector<t_pppParam*>& params = _parlist->params();
-  for (unsigned iPar = 0; iPar < params.size(); iPar++) {
-    t_pppParam* par = params[iPar];
-    if (par->type() == t_pppParam::amb && par->prn() == prn) {
-      int ind = par->indexNew();
-      t_lc::type tLC = par->tLC();
-      LOG << string(_epoTime) << " RESET " << par->toString() << endl;
-      delete par; par = new t_pppParam(t_pppParam::amb, prn, tLC, &obsVector);
-      par->setIndex(ind);
-      params[iPar] = par;
-      for (unsigned ii = 1; ii <= params.size(); ii++) {
-        _QFlt(ii, ind+1) = 0.0;
-        if (QSav) {
-          (*QSav)(ii, ind+1) = 0.0;
-        }
-      }
-      _QFlt(ind+1,ind+1) = par->sigma0() * par->sigma0();
-      if (QSav) {
-        (*QSav)(ind+1,ind+1) = _QFlt(ind+1,ind+1);
-      }
-      _xFlt[ind] = 0.0;
-      if (xSav) {
-        (*xSav)[ind] = _xFlt[ind];
-      }
-      _x0[ind] = par->x0();
-      irc = success;
-    }
-  }
-
-  return irc;
-}
-
-// Compute various DOP Values
-////////////////////////////////////////////////////////////////////////////
-void t_pppFilter::cmpDOP(const vector<t_pppSatObs*>& obsVector) {
-
-  _dop.reset();
-
-  try {
-    const unsigned numPar = 4;
-    Matrix AA(obsVector.size(), numPar);
-    _numSat = 0;
-    for (unsigned ii = 0; ii < obsVector.size(); ii++) {
-      t_pppSatObs* obs = obsVector[ii];
-      if (obs->isValid() && !obs->outlier()) {
-        ++_numSat;
-        for (unsigned iPar = 0; iPar < numPar; iPar++) {
-          const t_pppParam* par = _parlist->params()[iPar];
-          AA[_numSat-1][iPar] = par->partial(_epoTime, obs, t_lc::c1);
-        }
-      }
-    }
-    if (_numSat < 4) {
-      return;
-    }
-    AA = AA.Rows(1, _numSat);
-    SymmetricMatrix NN; NN << AA.t() * AA;  
-    SymmetricMatrix QQ = NN.i();
-    
-    _dop.P = sqrt(QQ(1,1) + QQ(2,2) + QQ(3,3));
-    _dop.T = sqrt(QQ(4,4));
-    _dop.G = sqrt(QQ(1,1) + QQ(2,2) + QQ(3,3) + QQ(4,4));
-  }
-  catch (...) {
-  }
-}
-
-// Compute various DOP Values
-////////////////////////////////////////////////////////////////////////////
-void t_pppFilter::predictCovCrdPart(const SymmetricMatrix& QFltOld) {
-
-  const vector<t_pppParam*>& params = _parlist->params();
-  if (params.size() < 3) {
+      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);
+    }
+  }
+}
+
+// 
+////////////////////////////////////////////////////////////////////////////
+void t_pppFilter::cmpDOP(t_epoData* epoData) {
+
+  Tracer tracer("t_pppFilter::cmpDOP");
+
+  _numSat = 0;
+  _pDop   = 0.0;
+
+  if (_params.size() < 4) {
     return;
   }
 
-  bool first = (params[0]->indexOld() < 0);
-
-  SymmetricMatrix Qneu(3); Qneu = 0.0;
-  for (unsigned ii = 0; ii < 3; ii++) {
-    const t_pppParam* par = params[ii];
-    if (first) {
-      Qneu[ii][ii] = par->sigma0() * par->sigma0();
-    }
-    else {
-      Qneu[ii][ii] = par->noise() * par->noise();
-    }
-  }
-
-  const t_pppStation* sta = PPP_CLIENT->staRover();
-  SymmetricMatrix Qxyz(3);
-  covariNEU_XYZ(Qneu, sta->ellApr().data(), Qxyz);
-
-  if (first) {
-    _QFlt.SymSubMatrix(1,3) = Qxyz;
-  }
-  else {
-    double dt = _epoTime - _firstEpoTime;
-    if (dt < OPT->_seedingTime) {
-      _QFlt.SymSubMatrix(1,3) = QFltOld.SymSubMatrix(1,3);
-    }
-    else {
-      _QFlt.SymSubMatrix(1,3) = QFltOld.SymSubMatrix(1,3) + Qxyz;
-    }
-  }
-}
+  const unsigned numPar = 4;
+  Matrix AA(epoData->sizeAll(), numPar);
+  QMapIterator<QString, t_satData*> it(epoData->satData);
+  while (it.hasNext()) {
+    it.next();
+    t_satData* satData = it.value();
+    _numSat += 1;
+    for (unsigned iPar = 0; iPar < numPar; iPar++) {
+      AA[_numSat-1][iPar] = _params[iPar]->partial(satData, false);
+    }
+  }
+  if (_numSat < 4) {
+    return;
+  }
+  AA = AA.Rows(1, _numSat);
+  SymmetricMatrix NN; NN << AA.t() * AA;  
+  SymmetricMatrix QQ = NN.i();
+    
+  _pDop = sqrt(QQ(1,1) + QQ(2,2) + QQ(3,3));
+}

trunk/BNC/src/PPP/pppFilter.h

@@ -1 @@
-#ifndef FILTER_H
-#define FILTER_H
-
-#include <vector>
+// 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.
+
+#ifndef PPPFILTER_H
+#define PPPFILTER_H
+
+#include <QtCore>
+#include <QtNetwork>
 #include <newmat.h>
-#include "pppInclude.h"
-#include "pppParlist.h"
+
+#include "bncconst.h"
 #include "bnctime.h"
-#include "t_prn.h"
+
+class bncAntex;
 
 namespace BNC_PPP {
 
-class t_pppParlist;
-class t_pppObsPool;
-class t_pppSatObs;
+class t_pppClient;
+class t_pppOptions;
+class t_epoData;
+class t_satData;
+class t_tides;
+
+class t_satData {
+ public:
+  t_satData() {
+    obsIndex = 0;
+    P1      = 0.0;
+    P2      = 0.0;
+    P5      = 0.0;
+    P3      = 0.0;
+    L1      = 0.0;
+    L2      = 0.0;
+    L5      = 0.0;
+    L3      = 0.0;
+    lkA     = 0.0;
+    lkB     = 0.0;
+  }
+  ~t_satData() {}
+  bncTime      tt;
+  QString      prn;
+  double       P1;
+  double       P2;
+  double       P5;
+  double       P7;
+  double       P3;
+  double       L1;
+  double       L2;
+  double       L5;
+  double       L7;
+  double       L3;
+  ColumnVector xx;
+  ColumnVector vv;
+  double       clk;
+  double       eleSat;
+  double       azSat;
+  double       rho;
+  bool         slipFlag;
+  double       lambda3;
+  double       lkA;
+  double       lkB;
+  unsigned     obsIndex;
+  char system() const {return prn.toAscii()[0];}
+};
+
+class t_epoData {
+ public:
+  t_epoData() {}
+
+  ~t_epoData() {
+    clear();
+  }
+
+  void clear() {
+    QMapIterator<QString, t_satData*> it(satData);
+    while (it.hasNext()) {
+      it.next();
+      delete it.value();
+    }
+    satData.clear();
+    tt.reset();
+  }
+
+  void deepCopy(const t_epoData* from) {
+    clear();
+    tt = from->tt;
+    QMapIterator<QString, t_satData*> it(from->satData);
+    while (it.hasNext()) {
+      it.next();
+      satData[it.key()] = new t_satData(*it.value());
+    }
+  }
+
+  unsigned sizeSys(char system) const {
+    unsigned ans = 0;
+    QMapIterator<QString, t_satData*> it(satData);
+    while (it.hasNext()) {
+      it.next();
+      if (it.value()->system() == system) {
+        ++ans;
+      }
+    }
+    return ans;
+  }
+  unsigned sizeAll() const {return satData.size();}
+
+  bncTime                   tt;
+  QMap<QString, t_satData*> satData;
+};
+
+class t_pppParam {
+ public:
+  enum parType {CRD_X, CRD_Y, CRD_Z, RECCLK, TROPO, AMB_L3, 
+                GLONASS_OFFSET, GALILEO_OFFSET, BDS_OFFSET};
+  t_pppParam(parType typeIn, int indexIn, const QString& prn);
+  ~t_pppParam();
+  double partial(t_satData* satData, bool phase);
+  bool isCrd() const {
+    return (type == CRD_X || type == CRD_Y || type == CRD_Z);
+  }
+  parType  type;
+  double   xx;
+  int      index;
+  int      index_old;
+  int      numEpo;
+  QString  prn;
+};
 
 class t_pppFilter {
  public:
-  t_pppFilter();
+  t_pppFilter(t_pppClient* pppClient);
   ~t_pppFilter();
-
-  t_irc processEpoch(t_pppObsPool* obsPool);
-
-  const ColumnVector&    x() const {return _xFlt;}
-  const SymmetricMatrix& Q() const {return _QFlt;}
-
+  t_irc update(t_epoData* epoData);
+  bncTime time()  const {return _time;}
+  const SymmetricMatrix& Q() const {return _QQ;}
+  const ColumnVector& neu() const {return _neu;}
   int    numSat() const {return _numSat;}
-  double PDOP() const {return _dop.P;}
-  double GDOP() const {return _dop.G;}
+  double PDOP()   const {return _pDop;}
+  double x()      const {return _params[0]->xx;}
+  double y()      const {return _params[1]->xx;}
+  double z()      const {return _params[2]->xx;}
+  double clk()    const {return _params[3]->xx;}
+  double trp0()   {return delay_saast(M_PI/2.0);}
   double trp() const {
-    const std::vector<t_pppParam*>& par = _parlist->params();
-    for (unsigned ii = 0; ii < par.size(); ++ii) {
-      if (par[ii]->type() == t_pppParam::trp) {
-        return x()[ii];
-      }
-    }
-    return 0.0;
-  };
+    for (int ii = 0; ii < _params.size(); ++ii) {
+      t_pppParam* pp = _params[ii];
+      if (pp->type == t_pppParam::TROPO) {
+        return pp->xx;
+      }
+    }
+    return 0.0;
+  }
   double trpStdev() const {
-    const std::vector<t_pppParam*>& par = _parlist->params();
-    for (unsigned ii = 0; ii < par.size(); ++ii) {
-      if (par[ii]->type() == t_pppParam::trp) {
+    for (int ii = 0; ii < _params.size(); ++ii) {
+      t_pppParam* pp = _params[ii];
+      if (pp->type == t_pppParam::TROPO) {
         return sqrt(Q()[ii][ii]);
       }
     }
     return 0.0;
-  };
-
+  }
+  double Glonass_offset() const {
+    for (int ii = 0; ii < _params.size(); ++ii) {
+      t_pppParam* pp = _params[ii];
+      if (pp->type == t_pppParam::GLONASS_OFFSET) {
+        return pp->xx;
+      }
+    }
+    return 0.0;
+  }
+  double Galileo_offset() const {
+    for (int ii = 0; ii < _params.size(); ++ii) {
+      t_pppParam* pp = _params[ii];
+      if (pp->type == t_pppParam::GALILEO_OFFSET) {
+        return pp->xx;
+      }
+    }
+    return 0.0;
+  }
+  double Bds_offset() const {
+    for (int ii = 0; ii < _params.size(); ++ii) {
+      t_pppParam* pp = _params[ii];
+      if (pp->type == t_pppParam::BDS_OFFSET) {
+        return pp->xx;
+      }
+    }
+    return 0.0;
+  }
  private:
-  class t_slip {
-   public:
-    t_slip() {
-      _slip            = false;
-      _obsSlipCounter  = -1;
-      _biasJumpCounter = -1;
-    }
-    bool _slip;
-    int  _obsSlipCounter;
-    int  _biasJumpCounter;
-  };
-
-  class t_dop {
-   public:
-    t_dop() {reset();}
-    void reset() {P = T = G = 0.0;}
-    double P;
-    double T;
-    double G;
-  };
-
-  t_irc processSystem(const std::vector<t_lc::type>& LCs, 
-                      const std::vector<t_pppSatObs*>& obsVector);
-
-  t_irc detectCycleSlips(const std::vector<t_lc::type>& LCs, 
-                         const std::vector<t_pppSatObs*>& obsVector);
-
-  t_irc resetAmb(t_prn prn, const std::vector<t_pppSatObs*>& obsVector,
-                 SymmetricMatrix* QSav = 0, ColumnVector* xSav = 0);
-
-  void cmpDOP(const std::vector<t_pppSatObs*>& obsVector);
-
-  void predictCovCrdPart(const SymmetricMatrix& QFltOld);
-
-  bncTime         _epoTime;
-  t_pppParlist*   _parlist;
-  SymmetricMatrix _QFlt;
-  ColumnVector    _xFlt;
-  ColumnVector    _x0;
-  t_slip          _slips[t_prn::MAXPRN+1];
-  int             _numSat;
-  t_dop           _dop;
-  bncTime         _firstEpoTime;
+  void   reset();
+  t_irc  cmpBancroft(t_epoData* epoData);
+  void   cmpEle(t_satData* satData);
+  void   addAmb(t_satData* satData);
+  void   addObs(int iPhase, unsigned& iObs, t_satData* satData,
+                Matrix& AA, ColumnVector& ll, DiagonalMatrix& PP);
+  QByteArray printRes(int iPhase, const ColumnVector& vv, 
+                      const QMap<QString, t_satData*>& satDataMap);
+  void   findMaxRes(const ColumnVector& vv,
+                    const QMap<QString, t_satData*>& satData,
+                    QString& prnGPS, QString& prnGlo,  
+                    double& maxResGPS, double& maxResGlo); 
+  double cmpValue(t_satData* satData, bool phase);
+  double delay_saast(double Ele);
+  void   predict(int iPhase, t_epoData* epoData);
+  t_irc  update_p(t_epoData* epoData);
+  QString outlierDetection(int iPhase, const ColumnVector& vv,
+                           QMap<QString, t_satData*>& satData);
+
+  double windUp(const QString& prn, const ColumnVector& rSat,
+                const ColumnVector& rRec);
+
+  bncTime  _startTime;
+
+  void rememberState(t_epoData* epoData);
+  void restoreState(t_epoData* epoData);
+  
+  t_irc selectSatellites(const QString& lastOutlierPrn, 
+                         QMap<QString, t_satData*>& satData);
+
+  void bancroft(const Matrix& BBpass, ColumnVector& pos);
+
+  void cmpDOP(t_epoData* epoData);
+
+  t_pppClient*          _pppClient;
+  bncTime               _time;
+  bncTime               _lastTimeOK;
+  QVector<t_pppParam*>  _params;
+  SymmetricMatrix       _QQ;
+  QVector<t_pppParam*>  _params_sav;
+  SymmetricMatrix       _QQ_sav;
+  t_epoData*            _epoData_sav;
+  ColumnVector          _xcBanc;
+  ColumnVector          _ellBanc;
+  QMap<QString, double> _windUpTime;
+  QMap<QString, double> _windUpSum;
+  QStringList           _outlierGPS;
+  QStringList           _outlierGlo;
+  bncAntex*             _antex;
+  t_tides*              _tides;
+  ColumnVector          _neu;
+  int                   _numSat;
+  double                _pDop;
 };
 

trunk/BNC/src/src.pri

@@ -117 +117 @@
 }
 
-exists(PPP) {
-  INCLUDEPATH += PPP
-  HEADERS += PPP/pppClient.h    PPP/pppObsPool.h   PPP/pppEphPool.h   \
-             PPP/pppStation.h   PPP/pppFilter.h    PPP/pppParlist.h   \
-             PPP/pppSatObs.h           
-  SOURCES += PPP/pppClient.cpp  PPP/pppObsPool.cpp PPP/pppEphPool.cpp \
-             PPP/pppStation.cpp PPP/pppFilter.cpp  PPP/pppParlist.cpp \
-             PPP/pppSatObs.cpp      
+exists(PPP_AR) {
+  INCLUDEPATH += PPP_AR
+  DEFINES += USE_PPP_AR
+  HEADERS += PPP_AR/pppClient.h    PPP_AR/pppObsPool.h   PPP_AR/pppEphPool.h   \
+             PPP_AR/pppStation.h   PPP_AR/pppFilter.h    PPP_AR/pppParlist.h   \
+             PPP_AR/pppSatObs.h           
+  SOURCES += PPP_AR/pppClient.cpp  PPP_AR/pppObsPool.cpp PPP_AR/pppEphPool.cpp \
+             PPP_AR/pppStation.cpp PPP_AR/pppFilter.cpp  PPP_AR/pppParlist.cpp \ 
+             PPP_AR/pppSatObs.cpp      
 }
 else {
-  INCLUDEPATH += PPP_SSR_I
-  DEFINES += USE_PPP_SSR_I
-  HEADERS += PPP_SSR_I/pppClient.h   PPP_SSR_I/pppFilter.h
-  SOURCES += PPP_SSR_I/pppClient.cpp PPP_SSR_I/pppFilter.cpp
+  INCLUDEPATH += PPP
+  DEFINES += USE_PPP
+  HEADERS += PPP/pppClient.h   PPP/pppFilter.h
+  SOURCES += PPP/pppClient.cpp PPP/pppFilter.cpp
 }