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

Timestamp:

Apr 9, 2020, 10:31:37 PM (5 years ago)

Author:

stuerze

Message:

minor changes regarding PPP

Location:

trunk/BNC/src/PPP

Files:

• pppFilter.cpp (modified) (8 diffs)
• pppFilter.h (modified) (2 diffs)

trunk/BNC/src/PPP/pppFilter.cpp

@@ -39 +39 @@
   _obsPool = obsPool;
   _refPrn  = t_prn();
+  _datumTrafo = new t_datumTrafo();
 }
 
@@ -45 +46 @@
 t_pppFilter::~t_pppFilter() {
   delete _parlist;
+  delete _datumTrafo;
 }
 
@@ -51 +53 @@
 t_irc t_pppFilter::processEpoch(int num) {
   _numSat     = 0;
-  _numOfEpochProcessing = num;
+  _numEpoProcessing = num;
   const double maxSolGap = 60.0;
 
@@ -119 +121 @@
       OPT->_obsModelType == OPT->DCMphaseBias) {
     preProcessing = true;
+    _numAllUsedLCs = 0;
     for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
       char system = OPT->systems()[iSys];
+      _numAllUsedLCs += OPT->LCs(system).size();
+      if (OPT->_pseudoObsIono && !epoch->pseudoObsIono()) {
+        _numAllUsedLCs -= 1;  // GIM not used
+      }
       if (OPT->_refSatRequired) {
         _refPrn = (_obsPool->getRefSatMapElement(system))->prn();
@@ -137 +144 @@
   }
 
+  if (_numEpoProcessing == 1) {
+    int maxObs = allObs.size() * _numAllUsedLCs;
+    _datumTrafo->initAA(maxObs, _parlist->nPar());
+  }
+
   // Process Satellite Systems separately
   // ------------------------------------
-  int numOfAllUsedLCs = 0;
   preProcessing = false;
   for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
     char system = OPT->systems()[iSys];
-    numOfAllUsedLCs += OPT->LCs(system).size();
+    (iSys) ? _datumTrafo->setFirstSystem(false) :
+                 _datumTrafo->setFirstSystem(true);
     if (OPT->_refSatRequired) {
       _refPrn = (_obsPool->getRefSatMapElement(system))->prn();
@@ -167 +179 @@
     _xFlt = xFltOld;
     _QFlt = QFltOld;
-    initDatumTransformation(numOfAllUsedLCs);
   }
   // close epoch processing
@@ -264 +275 @@
     }
 
+    if ((!preProcessing) &&
+       (OPT->_obsModelType == OPT->DCMcodeBias ||
+        OPT->_obsModelType == OPT->DCMphaseBias)) {
+          _datumTrafo->updateIndices(maxObs);
+      _datumTrafo->prepareAA(AA, _numEpoProcessing, _parlist->nPar());
+    }
+
     // Check number of observations, truncate matrices
     // -----------------------------------------------
@@ -551 +569 @@
 }
 
-//
-////////////////////////////////////////////////////////////////////////////
-void t_pppFilter::initDatumTransformation(int numOfUsedLCs) {
-
-  if (_numOfEpochProcessing == 1) {
-           _AA1.ReSize(numOfUsedLCs, _parlist->nPar());
-           _AA1 = 0.0;
-     _AA2.ReSize(numOfUsedLCs, _parlist->nPar());
-     _AA2 = 0.0;
-  }
-}
-
-//
-////////////////////////////////////////////////////////////////////////////
-void t_pppFilter::datumTransformation() {
-
-  Matrix D21 = (_AA2.t() * _AA2).i() * _AA2.t() * _AA1;
-
-  _QFlt = D21 * _QFlt * D21.t();
-}
-
-
-
-
+// Compute datum transformation
+////////////////////////////////////////////////////////////////////////////
+void t_pppFilter::datumTransformation(void) {
+  _QFlt = _datumTrafo->varCov(Q());
+}
+

trunk/BNC/src/PPP/pppFilter.h

@@ -75 +75 @@
   };
 
+  class t_datumTrafo{
+  public:
+    t_datumTrafo () {initIndices();}
+    void initIndices() {_firstRow = 1; _lastRow = 0;}
+    void setFirstSystem(bool firstSys) {_firstSys = firstSys;}
+    bool firstSystem() {return _firstSys;}
+    void updateIndices(int maxObs) {
+      if (_firstSys) {
+        initIndices();
+      }
+      else {
+        _firstRow += maxObs;
+      }
+      _lastRow += maxObs;
+    };
+    void initAA(int maxObs, int numPar) {
+      _AA1.ReSize(maxObs, numPar); _AA1 = 0.0;
+      _AA2.ReSize(maxObs, numPar); _AA2 = 0.0;
+    }
+    void prepareAA(Matrix& AA, int _numEpoProcessing, int nPar) {
+      Matrix& Prep = _AA2;
+      if (_numEpoProcessing == 1) {
+        Prep = _AA1;
+      }
+      Prep.SubMatrix(_firstRow, _lastRow, 1, nPar) = AA;
+    }
+    Matrix varCov(const SymmetricMatrix& QFlt) {
+      Matrix D21 = (_AA2.t() * _AA2).i() * _AA2.t() * _AA1;
+      return D21 * QFlt * D21.t();
+    }
+    int     _firstRow;
+    int     _lastRow;
+    Matrix  _AA1;
+    Matrix  _AA2;
+    bool    _firstSys;
+  };
+
   t_irc processSystem(const std::vector<t_lc::type>& LCs,
                       const std::vector<t_pppSatObs*>& obsVector,
@@ -93 +130 @@
   void predictCovCrdPart(const SymmetricMatrix& QFltOld);
 
-  void initDatumTransformation(int numOfAllUsedLCs);
-
   bncTime         _epoTime;
   t_pppParlist*   _parlist;
   t_pppObsPool*   _obsPool;
+  t_datumTrafo*   _datumTrafo;
   SymmetricMatrix _QFlt;
-  Matrix          _AA1;
-  Matrix          _AA2;
   ColumnVector    _xFlt;
   ColumnVector    _x0;
   t_slip          _slips[t_prn::MAXPRN+1];
   int             _numSat;
-  int             _numOfEpochProcessing;
+  int             _numEpoProcessing;
+  int             _numAllUsedLCs;
   t_dop           _dop;
   bncTime         _firstEpoTime;