Changeset 8905 in ntrip for trunk/BNC/src/PPP

Timestamp:

Mar 18, 2020, 11:13:50 AM (4 years ago)

Author:

stuerze

Message:

some developments regarding PPP, not completed!

Location:

trunk/BNC/src/PPP

Files:

• pppClient.cpp (modified) (13 diffs)
• pppClient.h (modified) (4 diffs)
• pppFilter.cpp (modified) (18 diffs)
• pppFilter.h (modified) (3 diffs)
• pppObsPool.cpp (modified) (4 diffs)
• pppObsPool.h (modified) (4 diffs)
• pppParlist.cpp (modified) (8 diffs)
• pppParlist.h (modified) (5 diffs)
• pppRefSat.cpp (added)
• pppRefSat.h (added)
• pppSatObs.cpp (modified) (24 diffs)
• pppSatObs.h (modified) (7 diffs)
• pppStation.h (modified) (3 diffs)

trunk/BNC/src/PPP/pppClient.cpp

@@ -55 +55 @@
   _obsPool  = new t_pppObsPool();
   _staRover = new t_pppStation();
-  _filter   = new t_pppFilter();
+  _filter   = new t_pppFilter(_obsPool);
   _tides    = new t_tides();
-
+  _antex    = 0;
   if (!_opt->_antexFileName.empty()) {
     _antex = new bncAntex(_opt->_antexFileName.c_str());
   }
-  else {
-    _antex = 0;
-  }
-
   if (!_opt->_blqFileName.empty()) {
-    _loading = new t_loading(_opt->_blqFileName.c_str());
-  }
-  else {
-    _loading = 0;
-  }
-
+    if (_tides->readBlqFile(_opt->_blqFileName.c_str()) == success) {
+      //_tides->printAllBlqSets();
+    }
+  }
+
+  if (OPT->_refSatRequired) {
+    for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
+      char system = OPT->systems()[iSys];
+      _obsPool->initRefSatMapElement(system);
+    }
+  }
+  _offGG = 0.0;
   CLIENTS.setLocalData(this);  // CLIENTS takes ownership over "this"
 }
@@ -80 +82 @@
   delete _log;
   delete _opt;
+  delete _filter;
   delete _ephPool;
   delete _obsPool;
@@ -86 +89 @@
     delete _antex;
   }
-  delete _filter;
   delete _tides;
-  if (_loading) {
-    delete _loading;
-  }
   clearObs();
 }
@@ -157 +156 @@
 t_irc t_pppClient::prepareObs(const vector<t_satObs*>& satObs,
                               vector<t_pppSatObs*>& obsVector, bncTime& epoTime) {
+
   // Default
   // -------
@@ -176 +176 @@
   }
 
+  // (re)set reference satellites per system if required
+  // ---------------------------------------------------
+  if (_opt->_refSatRequired) {
+    // reference satellite definition per system
+    for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
+      char system = OPT->systems()[iSys];
+      bool refSatDefined = false;
+      t_pppRefSat* refSat = _obsPool->getRefSatMapElement(system);
+      for (unsigned ii = 0; ii < obsVector.size(); ii++) {
+        const t_pppSatObs* satObs = obsVector.at(ii);
+        // reference satellite is unchanged
+        if      (!_obsPool->epoReProcessing() && refSat->prn() == satObs->prn()) {
+          refSatDefined = true;
+          obsVector[ii]->setAsReference();
+          refSat->setStatus(t_pppRefSat::unchanged);
+        }
+        // reference satellite has changed
+        else if (_obsPool->epoReProcessing() && refSat->prn() != satObs->prn()) {
+          if (satObs->prn().system() == system) {
+            refSatDefined = true;
+            obsVector[ii]->setAsReference();
+            refSat->setStatus(t_pppRefSat::changed);
+            refSat->setPrn(satObs->prn());
+          }
+        }
+        if (refSatDefined) {
+          break;
+        }
+      }
+      // reference satellite has to be initialized
+      if (!refSatDefined) {
+        for (unsigned ii = 0; ii < obsVector.size(); ii++) {
+          const t_pppSatObs* satObs = obsVector.at(ii);
+          if (satObs->prn().system() == system) {
+            obsVector[ii]->setAsReference();
+            refSat->setStatus(t_pppRefSat::initialized);
+            refSat->setPrn(satObs->prn());
+          }
+        }
+      }
+    }
+    _obsPool->setEpoReProcessing(false); //TODO: später erst nach Trafo false setzen
+  }
+
   // Check whether data are synchronized, compute epoTime
   // ----------------------------------------------------
@@ -202 +246 @@
 }
 
+//
+//////////////////////////////////////////////////////////////////////////////
+bool t_pppClient::preparePseudoObs(std::vector<t_pppSatObs*>& obsVector) {
+
+  bool pseudoObsIono = false;
+
+  if (OPT->_pseudoObsIono) {
+    vector<t_pppSatObs*>::iterator it = obsVector.begin();
+    while (it != obsVector.end()) {
+      t_pppSatObs* satObs = *it;
+      char system = satObs->prn().system();
+      t_pppRefSat* refSat = _obsPool->getRefSatMapElement(system);
+      double stecRef = refSat->stecValue();
+      if (stecRef && !satObs->isReference()) {
+        pseudoObsIono = true;
+        satObs->setPseudoObsIono(t_frequency::G1, stecRef);
+      }
+      satObs->printObsMinusComputed();
+      it++;
+    }
+  }
+
+  return pseudoObsIono;
+}
+
 // Compute the Bancroft position, check for blunders
 //////////////////////////////////////////////////////////////////////////////
@@ -250 +319 @@
     for (unsigned ii = 0; ii < obsVector.size(); ii++) {
       const t_pppSatObs* satObs = obsVector.at(ii);
-      if ( satObs->isValid() && satObs->prn().system() == 'G' &&
-           (!satObs->modelSet() || satObs->eleSat() >= OPT->_minEle) ) {
+      if (satObs->isValid() &&
+          satObs->prn().system() == 'G' &&
+          (!satObs->modelSet() || satObs->eleSat() >= OPT->_minEle) ) {
         ColumnVector rr = satObs->xc().Rows(1,3) - xyzc.Rows(1,3);
-        double res = rr.norm_Frobenius() - satObs->obsValue(tLC)
+        double res = rr.NormFrobenius() - satObs->obsValue(tLC)
                    - (satObs->xc()[3] - xyzc[3]) * t_CST::c;
         if (std::isnan(res) || fabs(res) > maxRes) {
@@ -393 +463 @@
 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);
+  station->setEpochTime(time);
   if (OPT->xyzAprRoverSet()) {
     station->setXyzApr(OPT->_xyzAprRover);
@@ -412 +482 @@
   // Tides
   // -----
-  station->setTideDspl( _tides->displacement(time, station->xyzApr()) );
-
-  // Ionosphere
-  // ----------
-  station->setIonoEpochTime(time);
+  station->setTideDsplEarth(_tides->earth(time, station->xyzApr()));
+  station->setTideDsplOcean(_tides->ocean(time, station->xyzApr(), station->name()));
 
   // Observation model
@@ -430 +497 @@
         satObs->eleSat() >= OPT->_minEle &&
         modelSetup == success) {
+      if (satObs->isReference() && OPT->_pseudoObsIono) {
+        char system = satObs->prn().system();
+        t_pppRefSat* refSat = _obsPool->getRefSatMapElement(system);
+        refSat->setStecValue(satObs->getIonoCodeDelay(t_frequency::G1));
+      }
       ++it;
     }
@@ -447 +519 @@
   try {
     initOutput(output);
-
-    // Prepare Observations of the Rover
-    // ---------------------------------
-    if (prepareObs(satObs, _obsRover, _epoTimeRover) != success) {
-      return finish(failure);
-    }
-
-    LOG << "\nPPP of Epoch ";
-    if (!_epoTimeRover.undef()) LOG << string(_epoTimeRover);
-    LOG << "\n---------------------------------------------------------------\n";
-
-    for (int iter = 1; iter <= 2; iter++) {
-      ColumnVector xyzc(4); xyzc = 0.0;
-      bool print = (iter == 2);
-      if (cmpBancroft(_epoTimeRover, _obsRover, xyzc, print) != success) {
+    _num = 0;
+    _obsPool->setEpoReProcessing(false); // initialize for epoch
+
+    do {
+      _num++;
+
+      // Prepare Observations of the Rover
+      // ---------------------------------
+      if (prepareObs(satObs, _obsRover, _epoTimeRover) != success) {
         return finish(failure);
       }
-      if (cmpModel(_staRover, xyzc, _obsRover) != success) {
+
+      LOG << "\nPPP of Epoch ";
+      if (!_epoTimeRover.undef()) LOG << string(_epoTimeRover);
+      LOG << "\n---------------------------------------------------------------\n";
+
+      for (int iter = 1; iter <= 2; iter++) {
+        ColumnVector xyzc(4); xyzc = 0.0;
+        bool print = (iter == 2);
+        if (cmpBancroft(_epoTimeRover, _obsRover, xyzc, print) != success) {
+          return finish(failure);
+        }
+        if (cmpModel(_staRover, xyzc, _obsRover) != success) {
+          return finish(failure);
+        }
+      }
+
+      _offGG = cmpOffGG(_obsRover);
+
+      // Prepare Pseudo Observations of the Rover
+      // ----------------------------------------
+      _pseudoObsIono = preparePseudoObs(_obsRover);//qDebug() << "_pseudoObsIonoAvailable: " << _pseudoObsIono;
+
+      if (int(_obsRover.size()) < OPT->_minObs) {
+        LOG << "t_pppClient::processEpoch not enough observations" << endl;
         return finish(failure);
       }
-    }
-
-    _offGG = cmpOffGG(_obsRover);
-
-    if (int(_obsRover.size()) < OPT->_minObs) {
-      LOG << "t_pppClient::processEpoch not enough observations" << endl;
-      return finish(failure);
-    }
-
-    // Store last epoch of data
-    // ------------------------
-    _obsPool->putEpoch(_epoTimeRover, _obsRover);
-
-    // Process Epoch in Filter
-    // -----------------------
-    if (_filter->processEpoch(_obsPool) != success) {
-      return finish(failure);
-    }
+
+      if (OPT->_refSatRequired) {
+        LOG.setf(ios::fixed);
+        QMapIterator<char, t_pppRefSat*> it(_obsPool->getRefSatMap());
+        while (it.hasNext()) {
+          it.next();
+          char   sys = it.key();
+          string prn = it.value()->prn().toString();
+          LOG << string(_epoTimeRover) << " REFSAT  " << sys << ": " << prn << endl;
+        }
+      }
+
+      // Store last epoch of data
+      // ------------------------
+      //_obsRover.resize(2);
+      _obsPool->putEpoch(_epoTimeRover, _obsRover, _pseudoObsIono);
+
+      // Process Epoch in Filter
+      // -----------------------
+      if (_filter->processEpoch(_num) != success) {
+        return finish(failure);
+      }
+      // if num > 1 und !obsPool->epoReProcessing() => filter ->datumTransformation
+    } while (_obsPool->epoReProcessing());
   }
   catch (Exception& exc) {
@@ -588 +683 @@
 void t_pppClient::reset() {
 
-  // to delete all parameters
-  delete _filter;
-  _filter   = new t_pppFilter();
-
   // to delete old orbit and clock corrections
   delete _ephPool;
@@ -599 +690 @@
   delete _obsPool;
   _obsPool  = new t_pppObsPool();
-}
+
+  // to delete all parameters
+  delete _filter;
+  _filter   = new t_pppFilter(_obsPool);
+
+}

trunk/BNC/src/PPP/pppClient.h

@@ -10 +10 @@
 
 class bncAntex;
+class t_pppRefSat;
 
 namespace BNC_PPP {
@@ -34 +35 @@
   const t_pppEphPool* ephPool() const {return _ephPool;}
   const t_pppObsPool* obsPool() const {return _obsPool;}
-  const bncAntex*  antex() const {return _antex;}
+  const bncAntex*     antex() const {return _antex;}
   const t_pppStation* staRover() const {return _staRover;}
-  double           offGG() const {return _offGG;}
+  double              offGG() const {return _offGG;}
 
   std::ostringstream& log() {return *_log;}
@@ -53 +54 @@
   t_irc prepareObs(const std::vector<t_satObs*>& satObs,
                    std::vector<t_pppSatObs*>& obsVector, bncTime& epoTime);
+  bool  preparePseudoObs(std::vector<t_pppSatObs*>& obsVector);
   t_irc cmpModel(t_pppStation* station, const ColumnVector& xyzc,
                  std::vector<t_pppSatObs*>& obsVector);
@@ -71 +73 @@
   t_pppOptions*             _opt;
   t_tides*                  _tides;
-  t_loading*                _loading;
+  bool                      _pseudoObsIono;
+  int                       _num;
 };
 

trunk/BNC/src/PPP/pppFilter.cpp

@@ -34 +34 @@
 // Constructor
 ////////////////////////////////////////////////////////////////////////////
-t_pppFilter::t_pppFilter() {
+t_pppFilter::t_pppFilter(t_pppObsPool* obsPool) {
   _parlist = 0;
+  _numSat = 0;
+  _obsPool = obsPool;
 }
 
@@ -46 +48 @@
 // Process Single Epoch
 ////////////////////////////////////////////////////////////////////////////
-t_irc t_pppFilter::processEpoch(t_pppObsPool* obsPool) {
-
+t_irc t_pppFilter::processEpoch(int num) {
   _numSat     = 0;
   const double maxSolGap = 60.0;
@@ -57 +58 @@
   // Vector of all Observations
   // --------------------------
-  t_pppObsPool::t_epoch* epoch = obsPool->lastEpoch();
+  t_pppObsPool::t_epoch* epoch = _obsPool->lastEpoch();
   if (!epoch) {
     return failure;
@@ -75 +76 @@
   string epoTimeStr = string(_epoTime);
 
+  //--
   // Set Parameters
   // --------------
@@ -91 +93 @@
   for (unsigned ii = 0; ii < params.size(); ii++) {
     const t_pppParam* par1 = params[ii];
-
     _x0[ii] = par1->x0();
-
     int iOld = par1->indexOld();
     if (iOld < 0) {
@@ -113 +113 @@
   predictCovCrdPart(QFltOld);
 
+
+  // Pre-Process Satellite Systems separately
+  // ----------------------------------------
+  bool preProcessing = false;
+  if (OPT->_obsModelType == OPT->DCMcodeBias ||
+      OPT->_obsModelType == OPT->DCMphaseBias) {
+    preProcessing = true;
+    for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
+      char system = OPT->systems()[iSys];
+      t_prn refPrn = (_obsPool->getRefSatMapElement(system))->prn();
+      vector<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, refPrn,
+                        epoch->pseudoObsIono(), preProcessing) != success) {
+        return failure;
+      }
+    }
+  }
+
   // Process Satellite Systems separately
   // ------------------------------------
   for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
     char system = OPT->systems()[iSys];
+    t_prn refPrn = (_obsPool->getRefSatMapElement(system))->prn();
     unsigned int num = 0;
     vector<t_pppSatObs*> obsVector;
@@ -126 +150 @@
     }
     LOG << epoTimeStr << " SATNUM " << system << ' ' << right << setw(2) << num << endl;
-    if ( processSystem(OPT->LCs(system), obsVector) != success ) {
+    if (processSystem(OPT->LCs(system), obsVector, refPrn,
+                      epoch->pseudoObsIono(), preProcessing) != success) {
       return failure;
     }
   }
-
-  cmpDOP(allObs);
-
-  _parlist->printResult(_epoTime, _QFlt, _xFlt);
-  _lastEpoTimeOK = _epoTime;  // remember time of last successful epoch processing
+  if (_obsPool->epoReProcessing()) {
+    // set A1 und A2 abhängig von num
+    // if num == 1 => A1
+    // if num >  1 => A2
+  }
+  else {
+    cmpDOP(allObs);
+    _parlist->printResult(_epoTime, _QFlt, _xFlt);
+    _lastEpoTimeOK = _epoTime;  // remember time of last successful epoch processing
+  }
+
   return success;
 }
@@ -141 +172 @@
 ////////////////////////////////////////////////////////////////////////////
 t_irc t_pppFilter::processSystem(const vector<t_lc::type>& LCs,
-                                 const vector<t_pppSatObs*>& obsVector) {
-
+                                 const vector<t_pppSatObs*>& obsVector,
+                                 const t_prn& refPrn,
+                                 bool pseudoObsIonoAvailable,
+                                 bool preProcessing) {
+  qDebug() << "======t_pppFilter::processSystem=======";
   LOG.setf(ios::fixed);
 
   // Detect Cycle Slips
   // ------------------
-  if (detectCycleSlips(LCs, obsVector) != success) {
+  if (detectCycleSlips(LCs, obsVector, refPrn, preProcessing) != success) {
     return failure;
   }
 
-  ColumnVector            xSav       = _xFlt;
-  SymmetricMatrix         QSav       = _QFlt;
-  string                  epoTimeStr = string(_epoTime);
+  unsigned usedLCs = LCs.size(); //qDebug() << "usedLCs: " << usedLCs;
+  if (OPT->_pseudoObsIono && !pseudoObsIonoAvailable) {
+    usedLCs -= 1;  // GIM not used
+  }
+  ColumnVector               xSav       = _xFlt;
+  SymmetricMatrix            QSav       = _QFlt;
+  string                     epoTimeStr = string(_epoTime);
   const vector<t_pppParam*>& params     = _parlist->params();
-  unsigned                maxObs     = obsVector.size() * LCs.size();
+  unsigned                   maxObs     = obsVector.size() * usedLCs;
+  //unsigned                   maxObs     = 2 * usedLCs;
+
+  if (OPT->_pseudoObsIono && pseudoObsIonoAvailable) { // stecDiff w.r.t refSat
+    maxObs -= 1;
+  }
+  qDebug() << "par.size()      : " << _parlist->nPar() << " LCs.size(): " << usedLCs;
+  qDebug() << "obsVector.size(): " << obsVector.size() << " maxObs    : " << maxObs;
 
   // Outlier Detection Loop
   // ----------------------
   for (unsigned iOutlier = 0; iOutlier < maxObs; iOutlier++) {
+    qDebug() << "iOutlier: " << iOutlier;
 
     if (iOutlier > 0) {
@@ -168 +214 @@
     // First-Design Matrix, Terms Observed-Computed, Weight Matrix
     // -----------------------------------------------------------
+    qDebug() << "A(" << maxObs << "," << _parlist->nPar() << ")";
     Matrix                AA(maxObs, _parlist->nPar());
     ColumnVector          ll(maxObs);
     DiagonalMatrix        PP(maxObs); PP = 0.0;
+    //TETSPLOT
+    for (unsigned iPar = 0; iPar < params.size(); iPar++) {
+      const t_pppParam* par = params[iPar];
+      cout << "  " << par->toString() <<endl;
+    }
+    cout << endl;
 
     int iObs = -1;
@@ -176 +229 @@
     vector<t_lc::type>   usedTypes;
     for (unsigned ii = 0; ii < obsVector.size(); ii++) {
-      t_pppSatObs* obs = obsVector[ii];
+      t_pppSatObs* obs = obsVector[ii]; qDebug() << "SATELLITE: " << obs->prn().toString().c_str() << "isRef: " << obs->isReference();
       if (!obs->outlier()) {
-        for (unsigned jj = 0; jj < LCs.size(); jj++) {
+        for (unsigned jj = 0; jj < usedLCs; jj++) {
           const t_lc::type tLC = LCs[jj];
+          if (tLC == t_lc::GIM &&  obs->isReference()) {continue;}
           ++iObs;
           usedObs.push_back(obs);
@@ -186 +240 @@
             const t_pppParam* par = params[iPar];
             AA[iObs][iPar] = par->partial(_epoTime, obs, tLC);
+            cout << setw(10) << par->partial(_epoTime, obs, tLC);
           }
+          cout << endl;
           ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1));
           PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC));
+          //qDebug() << "ll[iObs]: " << ll[iObs];
+          //qDebug() << "PP[iObs]: " << PP[iObs];
         }
       }
@@ -229 +287 @@
       t_pppSatObs* obs = usedObs[maxOutlierIndex];
       t_pppParam* par = 0;
-      LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' '
-          << obs->prn().toString()                        << ' '
-          << setw(8) << setprecision(4) << maxOutlier << endl;
+      if (!preProcessing) {
+        LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' '
+            << obs->prn().toString()                        << ' '
+            << setw(8) << setprecision(4) << maxOutlier << endl;
+      }
       for (unsigned iPar = 0; iPar < params.size(); iPar++) {
         t_pppParam* hlp = params[iPar];
-        if (hlp->type() == t_pppParam::amb && hlp->prn()  == obs->prn() &&
+        if (hlp->type() == t_pppParam::amb &&
+            hlp->prn()  == obs->prn() &&
             hlp->tLC()  == usedTypes[maxOutlierIndex]) {
           par = hlp;
         }
       }
-      if (par) {
+      if      (par && preProcessing) {
+        if (par->prn() == refPrn) {
+          _obsPool->setEpoReProcessing(true);
+        }
+      }
+      else if (par && !preProcessing) {
         if (par->ambResetCandidate()) {
           resetAmb(par->prn(), obsVector, &QSav, &xSav);
@@ -248 +314 @@
         }
       }
-      else {
+      else if (!par && !preProcessing){
         obs->setOutlier();
       }
@@ -256 +322 @@
     // ---------------
     else {
-      for (unsigned jj = 0; jj < LCs.size(); jj++) {
-        for (unsigned ii = 0; ii < usedObs.size(); ii++) {
-          const t_lc::type tLC = usedTypes[ii];
-          t_pppSatObs*        obs = usedObs[ii];
-          if (tLC == LCs[jj]) {
-            obs->setRes(tLC, vv[ii]);
-            LOG << epoTimeStr << " RES "
-                << left << setw(3) << t_lc::toString(tLC) << right << ' '
-                << obs->prn().toString().substr(0,3) << ' '
-                << setw(8) << setprecision(4) << vv[ii] << endl;
+      if (!preProcessing) {
+        for (unsigned jj = 0; jj < LCs.size(); jj++) {
+          for (unsigned ii = 0; ii < usedObs.size(); ii++) {
+            const t_lc::type tLC = usedTypes[ii];
+            t_pppSatObs*        obs = usedObs[ii];
+            if (tLC == LCs[jj]) {
+              obs->setRes(tLC, vv[ii]);
+              LOG << epoTimeStr << " RES "
+                  << left << setw(3) << t_lc::toString(tLC) << right << ' '
+                  << obs->prn().toString().substr(0,3) << ' '
+                  << setw(8) << setprecision(4) << vv[ii] << endl;
+            }
           }
         }
@@ -279 +347 @@
 ////////////////////////////////////////////////////////////////////////////
 t_irc t_pppFilter::detectCycleSlips(const vector<t_lc::type>& LCs,
-                                    const vector<t_pppSatObs*>& obsVector) {
+                                    const vector<t_pppSatObs*>& obsVector,
+                                    const t_prn& refPrn,
+                                    bool preProcessing) {
 
   const double            SLIP       = 20.0;  // slip threshold
   string                  epoTimeStr = string(_epoTime);
-  const vector<t_pppParam*>& params     = _parlist->params();
+  const vector<t_pppParam*>& params  = _parlist->params();
 
   for (unsigned ii = 0; ii < LCs.size(); ii++) {
@@ -294 +364 @@
         // ---------------------------------
         bool slip = false;
-
         if (obs->slip()) {
           LOG << epoTimeStr  << "cycle slip set (obs) " << obs->prn().toString()  << endl;
@@ -312 +381 @@
           slip = true;
         }
-        _slips[obs->prn()]._biasJumpCounter = obs->biasJumpCounter();
 
         // Slip Set
         // --------
         if (slip) {
-          resetAmb(obs->prn(), obsVector);
-        }
-
+          if (preProcessing) {
+            if (obs->prn() == refPrn) {
+              _obsPool->setEpoReProcessing(true);
+            }
+          }
+          else {
+            resetAmb(obs->prn(), obsVector);
+          }
+        }
         // Check Pre-Fit Residuals
         // -----------------------
@@ -328 +402 @@
             AA[iPar] = par->partial(_epoTime, obs, tLC);
           }
-
           double ll = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA);
           double vv = DotProduct(AA, _xFlt) - ll;
-
           if (fabs(vv) > SLIP) {
-            LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' '
-                << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl;
-            resetAmb(obs->prn(), obsVector);
+            if (preProcessing) {
+              if (obs->prn() == refPrn) {
+                _obsPool->setEpoReProcessing(true);
+              }
+            }
+            else {
+              LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' '
+                  << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl;
+              resetAmb(obs->prn(), obsVector);
+            }
           }
         }

trunk/BNC/src/PPP/pppFilter.h

@@ -17 +17 @@
 class t_pppFilter {
  public:
-  t_pppFilter();
+  t_pppFilter(t_pppObsPool* obsPool);
   ~t_pppFilter();
 
-  t_irc processEpoch(t_pppObsPool* obsPool);
+  t_irc processEpoch(int num);
 
   const ColumnVector&    x() const {return _xFlt;}
@@ -74 +74 @@
 
   t_irc processSystem(const std::vector<t_lc::type>& LCs,
-                      const std::vector<t_pppSatObs*>& obsVector);
+                      const std::vector<t_pppSatObs*>& obsVector,
+                      const t_prn& refPrn,
+                      bool pseudoObsIonoAvailable,
+                      bool preProcessing);
 
   t_irc detectCycleSlips(const std::vector<t_lc::type>& LCs,
-                         const std::vector<t_pppSatObs*>& obsVector);
+                         const std::vector<t_pppSatObs*>& obsVector,
+                         const t_prn& refPrn,
+                         bool preProcessing);
 
   t_irc resetAmb(t_prn prn, const std::vector<t_pppSatObs*>& obsVector,
@@ -88 +93 @@
   bncTime         _epoTime;
   t_pppParlist*   _parlist;
+  t_pppObsPool*   _obsPool;
   SymmetricMatrix _QFlt;
   ColumnVector    _xFlt;

trunk/BNC/src/PPP/pppObsPool.cpp

@@ -22 +22 @@
 // Constructor
 /////////////////////////////////////////////////////////////////////////////
-t_pppObsPool::t_epoch::t_epoch(const bncTime& epoTime, vector<t_pppSatObs*>& obsVector) {
-  _epoTime   = epoTime;
+t_pppObsPool::t_epoch::t_epoch(const bncTime& epoTime, vector<t_pppSatObs*>& obsVector,
+                               bool pseudoObsIono) {
+  _epoTime        = epoTime;
+  _pseudoObsIono  = pseudoObsIono;
   for (unsigned ii = 0; ii < obsVector.size(); ii++) {
     _obsVector.push_back(obsVector[ii]);
@@ -48 +50 @@
   }
   _vTec = 0;
+  _epoReProcessing = false;
 }
 
@@ -64 +67 @@
     _epochs.pop_front();
   }
+  clearRefSatMap();
 }
 
@@ -92 +96 @@
 //
 /////////////////////////////////////////////////////////////////////////////
-void t_pppObsPool::putEpoch(const bncTime& epoTime, vector<t_pppSatObs*>& obsVector) {
+void t_pppObsPool::putEpoch(const bncTime& epoTime, vector<t_pppSatObs*>& obsVector,
+    bool pseudoObsIono) {
   const unsigned MAXSIZE = 2;
-  _epochs.push_back(new t_epoch(epoTime, obsVector));
+
+  if (epoReProcessing()) {
+    delete _epochs.back();
+    _epochs.pop_back();
+  }
+  _epochs.push_back(new t_epoch(epoTime, obsVector,pseudoObsIono));
   if (_epochs.size() > MAXSIZE) {
     delete _epochs.front();

trunk/BNC/src/PPP/pppObsPool.h

@@ -4 +4 @@
 #include <vector>
 #include <deque>
+#include <QMap>
 #include "pppSatObs.h"
 #include "bnctime.h"
+#include "pppRefSat.h"
 
 namespace BNC_PPP {
@@ -14 +16 @@
   class t_epoch {
    public:
-    t_epoch(const bncTime& epoTime, std::vector<t_pppSatObs*>& obsVector);
+    t_epoch(const bncTime& epoTime, std::vector<t_pppSatObs*>& obsVector,
+            bool pseudoObsIono);
     ~t_epoch();
-    std::vector<t_pppSatObs*>& obsVector() {return _obsVector;}
+          std::vector<t_pppSatObs*>& obsVector() {return _obsVector;}
     const std::vector<t_pppSatObs*>& obsVector() const {return _obsVector;}
     const bncTime& epoTime() const {return _epoTime;}
+    bool pseudoObsIono() const {return _pseudoObsIono;}
    private:
-    bncTime                _epoTime;
+    bncTime                   _epoTime;
     std::vector<t_pppSatObs*> _obsVector;
+    bool                      _pseudoObsIono;
   };
 
@@ -30 +35 @@
   void putTec(t_vTec* _vTec);
 
-  void putEpoch(const bncTime& epoTime, std::vector<t_pppSatObs*>& obsVector);
+  void putEpoch(const bncTime& epoTime, std::vector<t_pppSatObs*>& obsVector,
+                bool pseudoObs);
 
   const t_satCodeBias* satCodeBias(const t_prn& prn) const {
@@ -49 +55 @@
   }
 
+  void initRefSatMapElement(char system) {
+    _refSatMap[system] = new t_pppRefSat();
+  }
+  void clearRefSatMap() {
+    QMapIterator<char, t_pppRefSat*> it(_refSatMap);
+    while (it.hasNext()) {
+      it.next();
+      delete it.value();
+    }
+    _refSatMap.clear();
+  }
+  t_pppRefSat* getRefSatMapElement(char system) {
+    return _refSatMap[system];
+  }
+  QMap<char, t_pppRefSat*> getRefSatMap() {return _refSatMap;}
+
+  void setEpoReProcessing(bool epoReProcessing) {
+    _epoReProcessing = epoReProcessing;
+  }
+  bool epoReProcessing() {return _epoReProcessing;}
+
  private:
-  t_satCodeBias*       _satCodeBiases[t_prn::MAXPRN+1];
-  t_satPhaseBias*      _satPhaseBiases[t_prn::MAXPRN+1];
-  t_vTec*              _vTec;
-  std::deque<t_epoch*> _epochs;
+  t_satCodeBias*           _satCodeBiases[t_prn::MAXPRN+1];
+  t_satPhaseBias*          _satPhaseBiases[t_prn::MAXPRN+1];
+  t_vTec*                  _vTec;
+  std::deque<t_epoch*>     _epochs;
+  QMap<char, t_pppRefSat*> _refSatMap;
+  bool                     _epoReProcessing;
 };
 

trunk/BNC/src/PPP/pppParlist.cpp

@@ -94 +94 @@
      _noise   = OPT->_noiseTrp;
      break;
+    case ion:
+     _epoSpec = false;
+     _sigma0  = OPT->_aprSigIon;
+     _noise   = OPT->_noiseIon;
+     break;
+   case cBias1:
+   case cBias2:
+     _epoSpec = false;
+     _sigma0  = OPT->_aprSigCodeBias;
+     _noise   = OPT->_noiseCodeBias;
+     break;
+   case pBias1:
+   case pBias2:
+     _epoSpec = false;
+     _sigma0  = OPT->_aprSigPhaseBias;
+     _noise   = OPT->_noisePhaseBias;
+     break;
   }
 }
@@ -106 +123 @@
 ////////////////////////////////////////////////////////////////////////////
 double t_pppParam::partial(const bncTime& /* epoTime */, const t_pppSatObs* obs,
-                        const t_lc::type& tLC) const {
+                           const t_lc::type& tLC) const {//qDebug() << "t_pppParam::partial: " << tLC;
 
   // Special Case - Melbourne-Wuebbena
@@ -115 +132 @@
 
   const t_pppStation* sta  = PPP_CLIENT->staRover();
-  ColumnVector     rhoV = sta->xyzApr() - obs->xc().Rows(1,3);
+  ColumnVector        rhoV = sta->xyzApr() - obs->xc().Rows(1,3);
+
+  map<t_frequency::type, double> codeCoeff;
+  map<t_frequency::type, double> phaseCoeff;
+  map<t_frequency::type, double> ionoCoeff;
+  obs->lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
 
   switch (_type) {
   case crdX:
-    return (sta->xyzApr()[0] - obs->xc()[0]) / rhoV.norm_Frobenius();
+    if (tLC == t_lc::GIM) {return 0.0;}
+    return (sta->xyzApr()[0] - obs->xc()[0]) / rhoV.NormFrobenius();
   case crdY:
-    return (sta->xyzApr()[1] - obs->xc()[1]) / rhoV.norm_Frobenius();
+    if (tLC == t_lc::GIM) {return 0.0;}
+    return (sta->xyzApr()[1] - obs->xc()[1]) / rhoV.NormFrobenius();
   case crdZ:
-    return (sta->xyzApr()[2] - obs->xc()[2]) / rhoV.norm_Frobenius();
+    if (tLC == t_lc::GIM) {return 0.0;}
+    return (sta->xyzApr()[2] - obs->xc()[2]) / rhoV.NormFrobenius();
   case clkR:
+    if (tLC == t_lc::GIM) {return 0.0;}
     return 1.0;
   case offGG:
+    if (tLC == t_lc::GIM) {return 0.0;}
     return (obs->prn().system() == 'R') ? 1.0 : 0.0;
   case amb:
+    if      (tLC == t_lc::GIM) {return 0.0;}
+    else if ((OPT->_obsModelType == OPT->IF)     ||
+             (OPT->_obsModelType == OPT->PPPRTK) ||
+             (OPT->_obsModelType == OPT->UncombPPP) ||
+             (OPT->_obsModelType == OPT->DCMcodeBias  && !obs->isReference()) ||
+             (OPT->_obsModelType == OPT->DCMphaseBias && !obs->isReference())   ) {
+      if (obs->prn() == _prn) {
+        if      (tLC == _tLC) {
+          return (obs->lambda(tLC));
+        }
+        else if (tLC == t_lc::lIF && _tLC == t_lc::MW) {
+          return obs->lambda(t_lc::lIF) * obs->lambda(t_lc::MW) / obs->lambda(t_lc::l2);
+        }
+        else {
+          if      (_tLC == t_lc::l1) {
+            return obs->lambda(t_lc::l1) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
+          }
+          else if (_tLC == t_lc::l2) {
+            return obs->lambda(t_lc::l2) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
+          }
+        }
+      }
+    }
+    break;
+  case trp:
+    if (tLC == t_lc::GIM) {return 0.0;}
+    return  1.0 / sin(obs->eleSat());
+  case ion:
+//    qDebug() << "refPrn: " << _refPrn.toString().c_str();
     if (obs->prn() == _prn) {
-      if      (tLC == _tLC) {
-        return (obs->lambda(tLC));
-      }
-      else if (tLC == t_lc::lIF && _tLC == t_lc::MW) {
-        return obs->lambda(t_lc::lIF) * obs->lambda(t_lc::MW) / obs->lambda(t_lc::l2);
-      }
-      else {
-        map<t_frequency::type, double> codeCoeff;
-        map<t_frequency::type, double> phaseCoeff;
-        obs->lcCoeff(tLC, codeCoeff, phaseCoeff);
-        if      (_tLC == t_lc::l1) {
-          return obs->lambda(t_lc::l1) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
-        }
-        else if (_tLC == t_lc::l2) {
-          return obs->lambda(t_lc::l2) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
-        }
-      }
-    }
-    return 0.0;
-  case trp:
-    return 1.0 / sin(obs->eleSat());
-  }
-
+      if      (tLC == t_lc::c1) {
+        return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::c1)];
+      }
+      else if (tLC == t_lc::c2) {
+        return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::c2)];
+      }
+      else if (tLC == t_lc::l1) {
+        return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
+      }
+      else if (tLC == t_lc::l2) {
+        return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
+      }
+      else if (tLC == t_lc::GIM) {
+        return -1.0;
+      }
+    }
+    if (tLC == t_lc::GIM && _prn == _refPrn) {
+      return  1.0;
+    }
+    break;
+  case cBias1:
+    if  (tLC == t_lc::c1) {
+      return  1.0;
+    }
+    else {
+      return 0.0;
+    }
+    break;
+  case cBias2:
+     if (tLC == t_lc::c2) {
+      return  1.0;
+    }
+    else {
+      return 0.0;
+    }
+    break;
+  case pBias1:
+    if  (tLC == t_lc::l1) {
+      return  1.0;
+    }
+    else {
+      return 0.0;
+    }
+    break;
+  case pBias2:
+    if  (tLC == t_lc::l2) {
+      return  1.0;
+    }
+    else {
+      return 0.0;
+    }
+  }
   return 0.0;
 }
@@ -171 +256 @@
     break;
   case clkR:
-    ss << "CLK        ";
-    break;
-  case amb:
-    ss << "AMB " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
+    ss << "REC_CLK    ";
     break;
   case offGG:
@@ -181 +263 @@
   case trp:
     ss << "TRP        ";
+    break;
+  case amb:
+    ss << "AMB  " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
+    break;
+  case ion:
+    ss << "ION  " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
+    break;
+  case cBias1:
+  case cBias2:
+  case pBias1:
+  case pBias2:
+    ss << "BIAS " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << "REC";
     break;
   }
@@ -201 +295 @@
 //
 ////////////////////////////////////////////////////////////////////////////
-t_irc t_pppParlist::set(const bncTime& epoTime, const std::vector<t_pppSatObs*>& obsVector) {
+t_irc t_pppParlist::set(const bncTime& epoTime,
+    const std::vector<t_pppSatObs*>& obsVector) {
 
   // Remove some Parameters
@@ -271 +366 @@
   required.push_back(new t_pppParam(t_pppParam::clkR, t_prn(), t_lc::dummy));
 
-  // GPS-Glonass Clock Offset
+  // GPS-GLONASS Clock Offset
   // ------------------------
   if (OPT->useSystem('R')) {
@@ -283 +378 @@
   }
 
+  // Ionosphere
+  // ----------
+  if (OPT->_obsModelType == OPT->UncombPPP    ||
+      OPT->_obsModelType == OPT->DCMcodeBias  ||
+      OPT->_obsModelType == OPT->DCMphaseBias   ) {
+    for (unsigned jj = 0; jj < obsVector.size(); jj++) {
+      const t_pppSatObs*        satObs = obsVector[jj];
+      required.push_back(new t_pppParam(t_pppParam::ion, satObs->prn(), t_lc::dummy));
+    }
+  }
+
   // Ambiguities
   // -----------
   for (unsigned jj = 0; jj < obsVector.size(); jj++) {
     const t_pppSatObs*        satObs = obsVector[jj];
-    const vector<t_lc::type>& ambLCs = OPT->ambLCs(satObs->prn().system());
-    for (unsigned ii = 0; ii < ambLCs.size(); ii++) {
-      required.push_back(new t_pppParam(t_pppParam::amb, satObs->prn(), ambLCs[ii], &obsVector));
-    }
+    if ((OPT->_obsModelType == OPT->IF)        ||
+        (OPT->_obsModelType == OPT->PPPRTK)    ||
+        (OPT->_obsModelType == OPT->UncombPPP) ||
+        (OPT->_obsModelType == OPT->DCMcodeBias  && !satObs->isReference()) ||
+        (OPT->_obsModelType == OPT->DCMphaseBias && !satObs->isReference())   ) {
+      const vector<t_lc::type>& ambLCs = OPT->ambLCs(satObs->prn().system());
+      for (unsigned ii = 0; ii < ambLCs.size(); ii++) {
+        required.push_back(new t_pppParam(t_pppParam::amb, satObs->prn(), ambLCs[ii], &obsVector));
+      }
+    }
+  }
+
+  // Receiver Code Biases
+  // --------------------
+  if (OPT->_obsModelType == OPT->DCMcodeBias) {
+    required.push_back(new t_pppParam(t_pppParam::cBias1, t_prn(), t_lc::c1));
+    required.push_back(new t_pppParam(t_pppParam::cBias2, t_prn(), t_lc::c2));
+  }
+
+  // Receiver Phase Biases
+  // ---------------------
+  if ((OPT->_obsModelType == OPT->DCMphaseBias) ||
+      (OPT->_obsModelType == OPT->PPPRTK)     ) {
+    required.push_back(new t_pppParam(t_pppParam::pBias1, t_prn(), t_lc::l1));
+    required.push_back(new t_pppParam(t_pppParam::pBias2, t_prn(), t_lc::l2));
   }
 

trunk/BNC/src/PPP/pppParlist.h

@@ -14 +14 @@
 class t_pppParam {
  public:
-  enum e_type {crdX, crdY, crdZ, clkR, amb, offGG, trp};
+  enum e_type {crdX, crdY, crdZ, clkR, amb, offGG, trp, ion, cBias1, cBias2, pBias1, pBias2};
 
   t_pppParam(e_type type, const t_prn& prn, t_lc::type tLC,
@@ -20 +20 @@
 
   ~t_pppParam();
+
   e_type type() const {return _type;}
   double x0()  const {return _x0;}
-  double partial(const bncTime& epoTime, const t_pppSatObs* obs, 
-                 const t_lc::type& tLC) const;
+  double partial(const bncTime& epoTime, const t_pppSatObs* obs, const t_lc::type& tLC) const;
   bool   epoSpec() const {return _epoSpec;}
   bool   isEqual(const t_pppParam* par2) const {
@@ -51 +51 @@
   unsigned ambNumEpo() const           {return _ambInfo ? _ambInfo->_numEpo : 0;}
   void     stepAmbNumEpo()             {if (_ambInfo) _ambInfo->_numEpo += 1;}
+  void     setRefPrn(t_prn prn)        {_refPrn = prn;}
 
   static bool sortFunction(const t_pppParam* p1, const t_pppParam* p2) {
@@ -78 +79 @@
     unsigned _numEpo;
   };
-  e_type     _type;
-  t_prn      _prn;
-  t_lc::type _tLC;
-  double     _x0;
-  bool       _epoSpec;
-  int        _indexOld;
-  int        _indexNew;
-  double     _sigma0;
-  double     _noise;
-  t_ambInfo* _ambInfo;
-  bncTime    _lastObsTime;
-  bncTime    _firstObsTime;
+  e_type       _type;
+  t_prn        _prn;
+  t_prn        _refPrn;
+  t_lc::type   _tLC;
+  double       _x0;
+  bool         _epoSpec;
+  int          _indexOld;
+  int          _indexNew;
+  double       _sigma0;
+  double       _noise;
+  t_ambInfo*   _ambInfo;
+  bncTime      _lastObsTime;
+  bncTime      _firstObsTime;
 };
 
@@ -100 +102 @@
   unsigned nPar() const {return _params.size();}
   const std::vector<t_pppParam*>& params() const {return _params;}
-  std::vector<t_pppParam*>& params() {return _params;}
-  void printResult(const bncTime& epoTime, const SymmetricMatrix& QQ, 
+        std::vector<t_pppParam*>& params()       {return _params;}
+  void printResult(const bncTime& epoTime, const SymmetricMatrix& QQ,
                    const ColumnVector& xx) const;
+
  private:
   std::vector<t_pppParam*> _params;

trunk/BNC/src/PPP/pppSatObs.cpp

@@ -35 +35 @@
 ////////////////////////////////////////////////////////////////////////////
 t_pppSatObs::t_pppSatObs(const t_satObs& pppSatObs) {
-  _prn     = pppSatObs._prn;
-  _time    = pppSatObs._time;
-  _outlier = false;
-  _valid   = true;
+  _prn        = pppSatObs._prn;
+  _time       = pppSatObs._time;
+  _outlier    = false;
+  _valid      = true;
+  _reference  = false;
+  _stecRefSat = 0.0;
+  _stecSat    = 0.0;
   for (unsigned ii = 0; ii < t_frequency::max; ii++) {
     _obs[ii] = 0;
@@ -61 +64 @@
   // Select pseudoranges and phase observations
   // ------------------------------------------
-  const string preferredAttrib = "CWPXI_";
-  //const string preferredAttrib = "G:12&PWCSLXYN G:5&IQX R:12&PC R:3&IQX E:16&BCX E:578&IQX J:1&SLXCZ J:26&SLX J:5&IQX C:IQX I:ABCX S:1&C S:5&IQX";
+  const string preferredAttrib = "G:12&PWCSLXYN G:5&IQX R:12&PC R:3&IQX E:16&BCX E:578&IQX J:1&SLXCZ J:26&SLX J:5&IQX C:IQX I:ABCX S:1&C S:5&IQX";
 
   for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
@@ -90 +92 @@
   for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) {
     t_lc::type tLC = OPT->LCs(_prn.system())[ii];
+    if (tLC == t_lc::GIM) {continue;}
     if (!isValid(tLC)) {
       _valid = false;
@@ -121 +124 @@
     ColumnVector dx = _xcSat - satPosOld;
     dx[3] *= t_CST::c;
-    if (dx.norm_Frobenius() < 1.e-4) {
+    if (dx.NormFrobenius() < 1.e-4) {
       totOK = true;
       break;
@@ -140 +143 @@
 void t_pppSatObs::lcCoeff(t_lc::type tLC,
                           map<t_frequency::type, double>& codeCoeff,
-                          map<t_frequency::type, double>& phaseCoeff) const {
+                          map<t_frequency::type, double>& phaseCoeff,
+                          map<t_frequency::type, double>& ionoCoeff) const {
 
   codeCoeff.clear();
   phaseCoeff.clear();
+  ionoCoeff.clear();
 
   double f1 = t_CST::freq(_fType1, _channel);
   double f2 = t_CST::freq(_fType2, _channel);
+  double f1GPS = t_CST::freq(t_frequency::G1, 0);
 
   switch (tLC) {
   case t_lc::l1:
-    phaseCoeff[_fType1] = 1.0;
+    phaseCoeff[_fType1] =  1.0;
+    ionoCoeff [_fType1] = -1.0 * pow(f1GPS, 2) / pow(f1, 2);
     return;
   case t_lc::l2:
-    phaseCoeff[_fType2] = 1.0;
+    phaseCoeff[_fType2] =  1.0;
+    ionoCoeff [_fType2] = -1.0 * pow(f1GPS, 2) / pow(f2, 2);
     return;
   case t_lc::lIF:
@@ -167 +175 @@
   case t_lc::CL:
     phaseCoeff[_fType1] =  0.5;
-    codeCoeff[_fType1]  =  0.5;
+    codeCoeff [_fType1] =  0.5;
     return;
   case t_lc::c1:
     codeCoeff[_fType1] = 1.0;
+    ionoCoeff[_fType1] = pow(f1GPS, 2) / pow(f1, 2);
     return;
   case t_lc::c2:
     codeCoeff[_fType2] = 1.0;
+    ionoCoeff[_fType2] = pow(f1GPS, 2) / pow(f2, 2);
     return;
   case t_lc::cIF:
@@ -179 +189 @@
     codeCoeff[_fType2] = -f2 * f2 / (f1 * f1 - f2 * f2);
     return;
+  case t_lc::GIM:
   case t_lc::dummy:
   case t_lc::maxLc:
@@ -190 +201 @@
   bool valid = true;
   obsValue(tLC, &valid);
+  //qDebug() << "tLC: " << tLC << "  valid: " << valid;
   return valid;
 }
@@ -195 +207 @@
 ////////////////////////////////////////////////////////////////////////////
 double t_pppSatObs::obsValue(t_lc::type tLC, bool* valid) const {
+
+  double retVal = 0.0;
+  if (valid) *valid = true;
+
+  // Pseudo observations
+  if (tLC == t_lc::GIM) {
+    if (_stecRefSat == 0.0 || _stecSat == 0.0) {
+      if (valid) *valid = false;
+      return 0.0;
+    }
+    else {
+      return _stecRefSat - _stecSat;
+    }
+  }
 
   map<t_frequency::type, double> codeCoeff;
   map<t_frequency::type, double> phaseCoeff;
-  lcCoeff(tLC, codeCoeff, phaseCoeff);
-
-  double retVal = 0.0;
-  if (valid) *valid = true;
+  map<t_frequency::type, double> ionoCoeff;
+  lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
 
   map<t_frequency::type, double>::const_iterator it;
+
+  // Code observations
   for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
     t_frequency::type tFreq = it->first;
@@ -214 +240 @@
     }
   }
+  // Phase observations
   for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
     t_frequency::type tFreq = it->first;
@@ -224 +251 @@
     }
   }
-
   return retVal;
 }
@@ -258 +284 @@
 double t_pppSatObs::sigma(t_lc::type tLC) const {
 
+  double retVal = 0.0;
   map<t_frequency::type, double> codeCoeff;
   map<t_frequency::type, double> phaseCoeff;
-  lcCoeff(tLC, codeCoeff, phaseCoeff);
-
-  double retVal = 0.0;
+  map<t_frequency::type, double> ionoCoeff;
+  lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
+
+  if (tLC == t_lc::GIM) {
+    retVal = OPT->_sigmaGIMdiff * OPT->_sigmaGIMdiff;
+  }
 
   map<t_frequency::type, double>::const_iterator it;
-  for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
+  for (it = codeCoeff.begin(); it != codeCoeff.end(); it++)   {//qDebug() << "codeCoeff : " << t_frequency::toString(it->first).c_str() << ": " << it->second;
     retVal += it->second * it->second * OPT->_sigmaC1 * OPT->_sigmaC1;
   }
-  for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
+
+  for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {//qDebug() << "phaseCoeff: " << t_frequency::toString(it->first).c_str()  << ": " << it->second;
     retVal += it->second * it->second * OPT->_sigmaL1 * OPT->_sigmaL1;
   }
@@ -295 +326 @@
 //
 ////////////////////////////////////////////////////////////////////////////
-double t_pppSatObs::maxRes(t_lc::type tLC) const {
+double t_pppSatObs::maxRes(t_lc::type tLC) const {//qDebug() << "t_pppSatObs::maxRes(t_lc::type tLC)";
+  double retVal = 0.0;
 
   map<t_frequency::type, double> codeCoeff;
   map<t_frequency::type, double> phaseCoeff;
-  lcCoeff(tLC, codeCoeff, phaseCoeff);
-
-  double retVal = 0.0;
+  map<t_frequency::type, double> ionoCoeff;
+  lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
 
   map<t_frequency::type, double>::const_iterator it;
-  for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
+  for (it = codeCoeff.begin(); it != codeCoeff.end(); it++)   {//qDebug() << "codeCoeff: " << it->first << ": " << it->second;
     retVal += it->second * it->second * OPT->_maxResC1 * OPT->_maxResC1;
   }
-  for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
+  for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {//qDebug() << "phaseCoeff: " << it->first << ": " << it->second;
     retVal += it->second * it->second * OPT->_maxResL1 * OPT->_maxResL1;
   }
-
+  if (tLC == t_lc::GIM) {
+    retVal = 3 * (OPT->_sigmaGIMdiff * OPT->_sigmaGIMdiff);
+  }
   return sqrt(retVal);
 }
@@ -326 +359 @@
   // ------------------------------
   ColumnVector rSat = _xcSat.Rows(1,3);
-  ColumnVector rhoV = rSat - station->xyzApr();
-  _model._rho = rhoV.norm_Frobenius();
+  ColumnVector rRec = station->xyzApr();
+  ColumnVector rhoV = rSat - rRec;
+  _model._rho = rhoV.NormFrobenius();
 
   ColumnVector vSat = _vvSat;
@@ -334 +368 @@
   xyz2neu(station->ellApr().data(), rhoV.data(), neu.data());
 
-  _model._eleSat = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / _model._rho );
+  _model._eleSat = acos(sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / _model._rho);
   if (neu[2] < 0) {
     _model._eleSat *= -1.0;
@@ -354 +388 @@
   Omega[1] = 0.0;
   Omega[2] = t_CST::omega / t_CST::c;
-  _model._sagnac = DotProduct(Omega, crossproduct(rSat, station->xyzApr()));
+  _model._sagnac = DotProduct(Omega, crossproduct(rSat, rRec));
 
   // Antenna Eccentricity
@@ -373 +407 @@
   // Tropospheric Delay
   // ------------------
-  _model._tropo = t_tropo::delay_saast(station->xyzApr(), _model._eleSat);
+  _model._tropo = t_tropo::delay_saast(rRec, _model._eleSat);
 
   // Code Biases
@@ -392 +426 @@
   // Phase Biases
   // -----------
-  // TODO: consideration of fix indicators and jump counter
   const t_satPhaseBias* satPhaseBias = PPP_CLIENT->obsPool()->satPhaseBias(_prn);
   double yaw = 0.0;
   bool ssr = false;
   if (satPhaseBias) {
-    yaw = satPhaseBias->_yaw;
+    double dt = station->epochTime() - satPhaseBias->_time;
+    if (satPhaseBias->_updateInt) {
+      dt -= (0.5 * ssrUpdateInt[satPhaseBias->_updateInt]);
+    }
+    yaw = satPhaseBias->_yaw + satPhaseBias->_yawRate * dt;
     ssr = true;
     for (unsigned ii = 0; ii < satPhaseBias->_bias.size(); ii++) {
@@ -414 +451 @@
   _model._windUp = station->windUp(_time, _prn, rSat, ssr, yaw, vSat) ;
 
+  // Relativistic effect due to earth gravity
+  // ----------------------------------------
+  double a = rSat.NormFrobenius() + rRec.NormFrobenius();
+  double b = (rSat - rRec).NormFrobenius();
+  double gm = 3.986004418e14; // m3/s2
+  _model._rel = 2 * gm / t_CST::c / t_CST::c * log((a + b) / (a - b));
 
   // Tidal Correction
   // ----------------
-  _model._tide = -DotProduct(station->tideDspl(), rhoV) / _model._rho;
+  _model._tideEarth = -DotProduct(station->tideDsplEarth(), rhoV) / _model._rho;
+  _model._tideOcean = -DotProduct(station->tideDsplOcean(), rhoV) / _model._rho;
 
   // Ionospheric Delay
@@ -429 +473 @@
     }
   }
+
   if (vTecUsage && vTec) {
-    double stec = station->stec(vTec, _signalPropagationTime, rSat);
+    double stec  = station->stec(vTec, _signalPropagationTime, rSat);
+    double f1GPS = t_CST::freq(t_frequency::G1, 0);
     for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
-      t_frequency::type frqType = static_cast<t_frequency::type>(iFreq);
-      _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(t_CST::freq(frqType, _channel), 2) * stec;
-    }
-  }
-
-  // Relativistic effect due to earth gravity
-  // ----------------------------------------
-  // TODO
-
-  // Ocean Loading
-  // -------------
-  // TODO
+      if (OPT->_pseudoObsIono) { // DCMcodeBias, DCMphaseBias
+        // For scaling the slant ionospheric delays the trick is to be consistent with units!
+        // The conversion of TECU into meters requires the frequency of the signal.
+        // Hence, GPS L1 frequency is used for all systems. The same is true for mu_i in lcCoeff().
+        _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(f1GPS, 2) * stec;
+      }
+      else { // PPP-RTK
+        t_frequency::type frqType = static_cast<t_frequency::type>(iFreq);
+        _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(t_CST::freq(frqType, _channel), 2) * stec;
+      }
+    }
+  }
 
   // Set Model Set Flag
@@ -449 +495 @@
   _model._set = true;
 
-  //printModel();
+  printModel();
 
   return success;
@@ -457 +503 @@
 ////////////////////////////////////////////////////////////////////////////
 void t_pppSatObs::printModel() const {
-
-  LOG.setf(ios::fixed);
-  LOG << "MODEL for Satellite " << _prn.toString() << endl
-      << "RHO:        " << setw(12) << setprecision(3) << _model._rho     << endl
-      << "ELE:        " << setw(12) << setprecision(3) << _model._eleSat * 180.0 / M_PI << endl
-      << "AZI:        " << setw(12) << setprecision(3) << _model._azSat  * 180.0 / M_PI << endl
-      << "SATCLK:     " << setw(12) << setprecision(3) << _model._satClkM << endl
-      << "RECCLK:     " << setw(12) << setprecision(3) << _model._recClkM << endl
-      << "SAGNAC:     " << setw(12) << setprecision(3) << _model._sagnac  << endl
-      << "ANTECC:     " << setw(12) << setprecision(3) << _model._antEcc  << endl
-      << "TROPO:      " << setw(12) << setprecision(3) << _model._tropo   << endl
-      << "WINDUP:     " << setw(12) << setprecision(3) << _model._windUp  << endl
-      << "TIDES:      " << setw(12) << setprecision(3) << _model._tide    << endl;
+// TODO: cout should be LOG
+  cout.setf(ios::fixed);
+  cout << "\nMODEL for Satellite " << _prn.toString() << (isReference() ? " (Reference Satellite)" : "") << endl
+       << "======================= " << endl
+       << "PPP STRATEGY  : " <<  OPT->_obsmodelTypeStr.at((int)OPT->_obsModelType).toLocal8Bit().constData()
+       <<  ((OPT->_pseudoObsIono) ? " with pseudo-observations for STEC" : "")  <<  endl
+      << "RHO           : " << setw(12) << setprecision(3) << _model._rho              << endl
+      << "ELE           : " << setw(12) << setprecision(3) << _model._eleSat * RHO_DEG << endl
+      << "AZI           : " << setw(12) << setprecision(3) << _model._azSat  * RHO_DEG << endl
+      << "SATCLK        : " << setw(12) << setprecision(3) << _model._satClkM          << endl
+      << "RECCLK        : " << setw(12) << setprecision(3) << _model._recClkM          << endl
+      << "SAGNAC        : " << setw(12) << setprecision(3) << _model._sagnac           << endl
+      << "ANTECC        : " << setw(12) << setprecision(3) << _model._antEcc           << endl
+      << "TROPO         : " << setw(12) << setprecision(3) << _model._tropo            << endl
+      << "WINDUP        : " << setw(12) << setprecision(3) << _model._windUp           << endl
+      << "REL           : " << setw(12) << setprecision(3) << _model._rel              << endl
+      << "EARTH TIDES   : " << setw(12) << setprecision(3) << _model._tideEarth        << endl
+      << "OCEAN TIDES   : " << setw(12) << setprecision(3) << _model._tideOcean        << endl
+      << endl
+      << "FREQUENCY DEPENDENT CORRECTIONS:" << endl
+      << "-------------------------------" << endl;
   for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
     if (_obs[iFreq]) {
       string frqStr = t_frequency::toString(t_frequency::type(iFreq));
       if (_prn.system() == frqStr[0]) {
-      LOG << "PCO           : " << frqStr << setw(12) << setprecision(3) << _model._antPCO[iFreq]    << endl
-          << "BIAS CODE     : " << frqStr << setw(12) << setprecision(3) << _model._codeBias[iFreq]  << endl
-          << "BIAS PHASE    : " << frqStr << setw(12) << setprecision(3) << _model._phaseBias[iFreq]  << endl
-          << "IONO CODEDELAY: " << frqStr << setw(12) << setprecision(3) << _model._ionoCodeDelay[iFreq] << endl;
-      }
-    }
-  }
+      cout << "PCO           : " << frqStr << setw(12) << setprecision(3) << _model._antPCO[iFreq]       << endl
+           << "BIAS CODE     : " << frqStr << setw(12) << setprecision(3) << _model._codeBias[iFreq]     << endl
+           << "BIAS PHASE    : " << frqStr << setw(12) << setprecision(3) << _model._phaseBias[iFreq]    << endl
+           << "IONO CODEDELAY: " << frqStr << setw(12) << setprecision(3) << _model._ionoCodeDelay[iFreq]<< endl;
+      }
+    }
+  }
+}
+
+//
+////////////////////////////////////////////////////////////////////////////
+void t_pppSatObs::printObsMinusComputed() const {
+// TODO: cout should be LOG
+  cout.setf(ios::fixed);
+  cout << "\nOBS-COMP for Satellite " << _prn.toString() << (isReference() ? " (Reference Satellite)" : "") << endl
+       << "========================== " << endl;
   for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) {
     t_lc::type tLC = OPT->LCs(_prn.system())[ii];
-    LOG << "OBS-CMP " << t_lc::toString(tLC) << ": " << _prn.toString() << " "
+    cout << "OBS-CMP " << setw(4) << t_lc::toString(tLC) << ": " << _prn.toString() << " "
         << setw(12) << setprecision(3) << obsValue(tLC) << " "
         << setw(12) << setprecision(3) << cmpValue(tLC) << " "
         << setw(12) << setprecision(3) << obsValue(tLC) - cmpValue(tLC) << endl;
-
-  }
-  LOG << "OBS-CMP MW: " << _prn.toString() << " "
-      << setw(12) << setprecision(3) << obsValue(t_lc::MW) << " "
-      << setw(12) << setprecision(3) << cmpValue(t_lc::MW) << " "
-      << setw(12) << setprecision(3) << obsValue(t_lc::MW) - cmpValue(t_lc::MW) << endl;
-}
+  }
+}
+
 
 //
@@ -504 +563 @@
 ////////////////////////////////////////////////////////////////////////////
 double t_pppSatObs::cmpValue(t_lc::type tLC) const {
-
-  if (!isValid(tLC)) {
-    return 0.0;
-  }
-
-  // Non-Dispersive Part
-  // -------------------
-  double nonDisp = _model._rho    + _model._recClkM - _model._satClkM
-                 + _model._sagnac + _model._antEcc  + _model._tropo
-                 + _model._tide;
-
-  // Add Dispersive Part
-  // -------------------
-  map<t_frequency::type, double> codeCoeff;
-  map<t_frequency::type, double> phaseCoeff;
-  lcCoeff(tLC, codeCoeff, phaseCoeff);
-
-  double dispPart = 0.0;
-
-  map<t_frequency::type, double>::const_iterator it;
-  for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
-    t_frequency::type tFreq = it->first;
-    dispPart += it->second * (_model._antPCO[tFreq] - _model._codeBias[tFreq] +
-                              _model._ionoCodeDelay[tFreq]);
-  }
-  for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
-    t_frequency::type tFreq = it->first;
-    dispPart += it->second * (_model._antPCO[tFreq] - _model._phaseBias[tFreq] +
-                              _model._windUp * t_CST::lambda(tFreq, _channel)  -
-                              _model._ionoCodeDelay[tFreq]);
-  }
-
-    return nonDisp + dispPart;
+  double cmpValue;
+
+  if      (!isValid(tLC)) {
+    cmpValue =  0.0;
+  }
+  else if (tLC == t_lc::GIM) {
+    cmpValue = 0.0;
+  }
+  else {
+    // Non-Dispersive Part
+    // -------------------
+    double nonDisp = _model._rho
+                   + _model._recClkM   - _model._satClkM
+                   + _model._sagnac    + _model._antEcc    + _model._tropo
+                   + _model._tideEarth + _model._tideOcean + _model._rel;
+
+    // Add Dispersive Part
+    // -------------------
+    double dispPart = 0.0;
+    map<t_frequency::type, double> codeCoeff;
+    map<t_frequency::type, double> phaseCoeff;
+    map<t_frequency::type, double> ionoCoeff;
+    lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
+    map<t_frequency::type, double>::const_iterator it;
+    for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
+      t_frequency::type tFreq = it->first;
+      dispPart += it->second * (_model._antPCO[tFreq] - _model._codeBias[tFreq]);
+      if (OPT->PPPRTK) {
+        dispPart += it->second * (_model._ionoCodeDelay[tFreq]);
+      }
+    }
+    for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
+      t_frequency::type tFreq = it->first;
+      dispPart += it->second * (_model._antPCO[tFreq] - _model._phaseBias[tFreq] +
+                                _model._windUp * t_CST::lambda(tFreq, _channel));
+      if (OPT->PPPRTK) {
+        dispPart += it->second * (- _model._ionoCodeDelay[tFreq]);
+      }
+    }
+    cmpValue = nonDisp + dispPart;
+  }
+
+  return cmpValue;
 }
 
@@ -556 +625 @@
   }
 }
+
+//
+////////////////////////////////////////////////////////////////////////////
+void  t_pppSatObs::setPseudoObsIono(t_frequency::type freq, double stecRefSat) {
+  _stecSat    = _model._ionoCodeDelay[freq];
+  _stecRefSat = stecRefSat;
+}

trunk/BNC/src/PPP/pppSatObs.h

@@ -19 +19 @@
   bool                isValid() const {return _valid;};
   bool                isValid(t_lc::type tLC) const;
+  bool                isReference() const {return _reference;};
+  void                setAsReference() {_reference = true;};
   const t_prn&        prn() const {return _prn;}
   const ColumnVector& xc() const {return _xcSat;}
@@ -31 +33 @@
   bool                modelSet() const {return _model._set;}
   void                printModel() const;
+  void                printObsMinusComputed() const;
   void                lcCoeff(t_lc::type tLC,
                               std::map<t_frequency::type, double>& codeCoeff,
-                              std::map<t_frequency::type, double>& phaseCoeff) const;
+                              std::map<t_frequency::type, double>& phaseCoeff,
+                              std::map<t_frequency::type, double>& ionoCoeff) const;
   double              lambda(t_lc::type tLC) const;
   double              sigma(t_lc::type tLC) const;
@@ -41 +45 @@
   void                setRes(t_lc::type tLC, double res);
   double              getRes(t_lc::type tLC) const;
+  void                setPseudoObsIono(t_frequency::type freq, double stecRefSat);
+  double              getIonoCodeDelay(t_frequency::type freq) {return _model._ionoCodeDelay[freq];}
 
   // RINEX
@@ -79 +85 @@
     ~t_model() {}
     void reset() {
-      _set     = false;
-      _rho     = 0.0;
-      _eleSat  = 0.0;
-      _azSat   = 0.0;
-      _recClkM = 0.0;
-      _satClkM = 0.0;
-      _sagnac  = 0.0;
-      _antEcc  = 0.0;
-      _tropo   = 0.0;
-      _tide    = 0.0;
-      _windUp  = 0.0;
-      _rel    = 0.0;
+      _set       = false;
+      _rho       = 0.0;
+      _eleSat    = 0.0;
+      _azSat     = 0.0;
+      _recClkM   = 0.0;
+      _satClkM   = 0.0;
+      _sagnac    = 0.0;
+      _antEcc    = 0.0;
+      _tropo     = 0.0;
+      _tideEarth = 0.0;
+      _tideOcean = 0.0;
+      _windUp    = 0.0;
+      _rel       = 0.0;
       for (unsigned ii = 0; ii < t_frequency::max; ii++) {
         _antPCO[ii]        = 0.0;
@@ -107 +114 @@
     double _antEcc;
     double _tropo;
-    double _tide;
+    double _tideEarth;
+    double _tideOcean;
     double _windUp;
     double _rel;
@@ -119 +127 @@
 
   bool                         _valid;
+  bool                         _reference;
   t_frequency::type            _fType1;
   t_frequency::type            _fType2;
@@ -131 +140 @@
   std::map<t_lc::type, double> _res;
   double                       _signalPropagationTime;
+  double                       _stecRefSat;
+  double                       _stecSat;
 };
 

trunk/BNC/src/PPP/pppStation.h

@@ -21 +21 @@
   void setNeuEcc(const ColumnVector& neuEcc);
   void setDClk(double dClk) {_dClk = dClk;}
-  void setTideDspl(const ColumnVector& tideDspl) {_tideDspl = tideDspl;}
-  void setIonoEpochTime(const bncTime& epochTime) {_epochTime = epochTime;}
+  void setTideDsplEarth(const ColumnVector& tideDsplEarth) {_tideDsplEarth = tideDsplEarth;}
+  void setTideDsplOcean(const ColumnVector& tideDsplOcean) {_tideDsplOcean = tideDsplOcean;}
+  void setEpochTime(const bncTime& epochTime) {_epochTime = epochTime;}
   const std::string&  name()      const {return _name;}
   const std::string&  antName()   const {return _antName;}
@@ -29 +30 @@
   const ColumnVector& neuEcc()    const {return _neuEcc;}
   const ColumnVector& xyzEcc()    const {return _xyzEcc;}
-  const ColumnVector& tideDspl()  const {return _tideDspl;}
+  const ColumnVector& tideDsplEarth()  const {return _tideDsplEarth;}
+  const ColumnVector& tideDsplOcean()  const {return _tideDsplOcean;}
+  const bncTime       epochTime() const {return _epochTime;}
   double dClk() const {return _dClk;}
   double windUp(const bncTime& time, t_prn prn, const ColumnVector& rSat, bool ssr,
@@ -42 +45 @@
   ColumnVector      _neuEcc;
   ColumnVector      _xyzEcc;
-  ColumnVector      _tideDspl;
+  ColumnVector      _tideDsplEarth;
+  ColumnVector      _tideDsplOcean;
   double            _dClk;
   mutable t_windUp* _windUp;