Changeset 10938 in ntrip for trunk/BNC

Timestamp:

Jun 18, 2026, 4:02:38 PM (3 days ago)

Author:

stuerze

Message:

bugfix regarding ionospheric constraints

Location:

trunk/BNC/src

Files:

• PPP/pppClient.cpp (modified) (2 diffs)
• PPP/pppClient.h (modified) (2 diffs)
• PPP/pppFilter.cpp (modified) (2 diffs)
• PPP/pppParlist.cpp (modified) (2 diffs)
• PPP/pppParlist.h (modified) (2 diffs)
• PPP/pppSatObs.cpp (modified) (4 diffs)
• PPP/pppSatObs.h (modified) (2 diffs)
• pppInclude.h (modified) (1 diff)

trunk/BNC/src/PPP/pppClient.cpp

@@ -230 +230 @@
 
   if (_opt->_pseudoObsIono) {
-    vector<t_pppSatObs*>::iterator it = obsVector.begin();
-    while (it != obsVector.end()) {
-      t_pppSatObs* satObs = *it;
-      pseudoObsIono = satObs->setPseudoObsIono(t_frequency::G1);
-      it++;
+
+    // Select reference satellite per system: keep the current one if still
+    // visible; only switch to highest elevation when the current one disappears
+    // -------------------------------------------------------------------------
+    map<char, t_pppSatObs*> refSatMap;
+    map<char, t_pppSatObs*> bestEleMap;
+    for (t_pppSatObs* satObs : obsVector) {
+      satObs->resetReference();
+      double stec = satObs->getIonoCodeDelay(t_frequency::G1);
+      if (stec == 0.0) continue;
+      char sys = satObs->prn().system();
+      if (!bestEleMap.count(sys) || satObs->eleSat() > bestEleMap[sys]->eleSat()) {
+        bestEleMap[sys] = satObs;
+      }
+      if (_refPrnMap.count(sys) && satObs->prn() == _refPrnMap[sys]) {
+        refSatMap[sys] = satObs;
+      }
+    }
+    for (auto& kv : bestEleMap) {
+      if (!refSatMap.count(kv.first)) {
+        refSatMap[kv.first] = kv.second;
+      }
+    }
+    for (auto& kv : refSatMap) {
+      _refPrnMap[kv.first] = kv.second->prn();
+      kv.second->setAsReference();
+    }
+
+    if (bestEleMap.empty()) {
+      LOG << "GIM pseudo-obs unavailable: no valid STEC (vTec data missing?)" << endl;
+    }
+
+    // Set STEC pseudo-observations (single-differenced vs reference satellite)
+    // -------------------------------------------------------------------------
+    int nGIM = 0;
+    for (t_pppSatObs* satObs : obsVector) {
+      char sys = satObs->prn().system();
+      double stecRefSat = refSatMap.count(sys)
+                          ? refSatMap[sys]->getIonoCodeDelay(t_frequency::G1)
+                          : 0.0;
+      if (satObs->setPseudoObsIono(t_frequency::G1, stecRefSat)) {
+        pseudoObsIono = true;
+        nGIM++;
+      }
+    }
+    if (nGIM > 0) {
+      LOG << "GIM pseudo-obs: " << nGIM << " satellites, ref:";
+      for (auto& kv : refSatMap) {
+        LOG << ' ' << kv.second->prn().toString();
+      }
+      LOG << endl;
     }
   }
@@ -598 +644 @@
 
   LOG << "pppClient: reset" << endl;
+  _refPrnMap.clear();
   // to delete old orbit and clock corrections
   delete _ephPool;

trunk/BNC/src/PPP/pppClient.h

@@ -2 +2 @@
 #define PPPCLIENT_H
 
+#include <map>
 #include <sstream>
 #include <vector>
@@ -76 +77 @@
   bool                      _pseudoObsIono;
   bool                      _running;
+  std::map<char, t_prn>     _refPrnMap;
 };
 

trunk/BNC/src/PPP/pppFilter.cpp

@@ -229 +229 @@
             const t_lc LC = LCs[jj];
             if (LC._type == t_lc::GIM) {
+              if (!obs->isValid(LC)) continue;
               ++iObs;
             } else {
@@ -264 +265 @@
     for (unsigned ii = 0; ii < usedObs.size(); ii++) {
       const t_lc LC = usedTypes[ii];
+      if (LC._type == t_lc::GIM) continue;   // GIM is a soft constraint; never reject as outlier
       double res = fabs(vv[ii]);
       if (res > usedObs[ii]->maxRes(LC)) {

trunk/BNC/src/PPP/pppParlist.cpp

@@ -112 +112 @@
   }
 
-  // Special Case - GIM
-  // ------------------
+  // Special Case - GIM (single-differenced: A[ION_ref]=+1, A[ION_i]=-1)
+  // -------------------------------------------------------------------
   if (obsLC._type == t_lc::GIM) {
     if (_type == ion) {
-      return 1.0;
-    }
-    else {
-      return 0.0;
-    }
+      if (obs->prn() == _prn) {
+        return -1.0;
+      }
+      if (_prn == _refPrn && _prn.system() == obs->prn().system()) {
+        return  1.0;
+      }
+    }
+    return 0.0;
   }
 
@@ -333 +336 @@
   }
   
+  // Reference satellite PRN per system (from satellite marked by preparePseudoObs)
+  // -------------------------------------------------------------------------------
+  if (OPT->_pseudoObsIono) {
+    map<char, t_prn> refPrnMap;
+    for (unsigned jj = 0; jj < obsVector.size(); jj++) {
+      if (obsVector[jj]->isReference()) {
+        refPrnMap[obsVector[jj]->prn().system()] = obsVector[jj]->prn();
+      }
+    }
+    for (unsigned ii = 0; ii < required.size(); ii++) {
+      if (required[ii]->type() == t_pppParam::ion) {
+        char sys = required[ii]->prn().system();
+        if (refPrnMap.count(sys)) {
+          required[ii]->setRefPrn(refPrnMap[sys]);
+        }
+      }
+    }
+  }
+
   // Ambiguities
   // -----------

trunk/BNC/src/PPP/pppParlist.h

@@ -64 +64 @@
   }
   void     setSystem(char sys) {_sys = sys;}
+  t_prn    refPrn() const {return _refPrn;}
+  void     setRefPrn(const t_prn& prn) {_refPrn = prn;}
 
   static bool sortFunction(const t_pppParam* p1, const t_pppParam* p2) {
@@ -93 +95 @@
   e_type       _type;
   t_prn        _prn;
+  t_prn        _refPrn;
   char         _sys;
   t_lc         _LC;

trunk/BNC/src/PPP/pppSatObs.cpp

@@ -45 +45 @@
   _reference  = false;
   _stecSat    = 0.0;
+  _stecRefSat = 0.0;
   for (unsigned ii = 0; ii < t_frequency::max; ii++) {
     _obs[ii] = 0;
@@ -265 +266 @@
   if (valid) *valid = true;
 
-  // Pseudo observations
+  // Pseudo observations (single-differenced: reference minus satellite)
   if (LC._type == t_lc::GIM) {
-    if (_stecSat == 0.0) {
+    if (_stecSat == 0.0 || _stecRefSat == 0.0) {
       if (valid) *valid = false;
       return 0.0;
     }
     else {
-      return _stecSat;
+      return _stecRefSat - _stecSat;
     }
   }
@@ -686 +687 @@
   }
   else if (LC._type == t_lc::GIM) {
-    cmpValue =  _stecSat;
+    cmpValue = 0.0;
   }
   else {
@@ -739 +740 @@
 //
 ////////////////////////////////////////////////////////////////////////////
-bool  t_pppSatObs::setPseudoObsIono(t_frequency::type freq) {
-  bool pseudoObsIono = false;
+bool  t_pppSatObs::setPseudoObsIono(t_frequency::type freq, double stecRefSat) {
   _stecSat    = _model._ionoCodeDelay[freq];
-  if (_stecSat) {
-    pseudoObsIono = true;
-  }
-  return pseudoObsIono;
+  _stecRefSat = stecRefSat;
+  return (_stecSat != 0.0 && _stecRefSat != 0.0);
 }
 

trunk/BNC/src/PPP/pppSatObs.h

@@ -48 +48 @@
     void                setRes(t_lc LC, double res);
     double              getRes(t_lc LC) const;
-    bool                setPseudoObsIono(t_frequency::type freq);
+    bool                setPseudoObsIono(t_frequency::type freq, double stecRefSat);
+    void                resetGIM() {_stecSat = 0.0;}
     double              getIonoCodeDelay(t_frequency::type freq) {return _model._ionoCodeDelay[freq];}
     double              getCodeBias(t_frequency::type freq) {return _model._codeBias[freq];}
@@ -153 +154 @@
     double                 _signalPropagationTime;
     double                 _stecSat;
+    double                 _stecRefSat;
     double                 _tropo0;
   };

trunk/BNC/src/pppInclude.h

@@ -151 +151 @@
     }
     else if (_type == GIM) {
-      out << t_frequency::toSystem(_frq1) << "GIM" ;
+      out << "GIM" ;
     }
     return out.str();