Changeset 11058 in ntrip

Timestamp:

Sep 30, 2026, 10:52:22 PM (12 hours ago)

Author:

stuerze

Message:

some more fixes

Location:

trunk/BNC

Files:

• Example_Configs/10_PPPQuickStart.bnc (modified) (2 diffs)
• Example_Configs/17_Combi.bnc (modified) (1 diff)
• Example_Configs/18_CombiPPP.bnc (modified) (1 diff)
• src/RTCM3/RTCM3coDecoder.cpp (modified) (4 diffs)
• src/RTCM3/RTCM3coDecoder.h (modified) (1 diff)
• src/bnchelp.html (modified) (2 diffs)
• src/bncoutf.h (modified) (1 diff)
• src/combination/bnccomb.cpp (modified) (5 diffs)
• src/combination/bnccomb.h (modified) (1 diff)
• src/upload/bncrtnetuploadcaster.cpp (modified) (3 diffs)
• src/upload/bncrtnetuploadcaster.h (modified) (1 diff)
• tests/README.md (modified) (1 diff)
• tests/cases/ppp_rt_quickstart/case.json (modified) (2 diffs)
• tests/cases/replay_save_sp3/expected/BKG0MGXRTS_20253210000_01D_05S_CLK.CLK.gz (modified) ( previous)
• tests/cases/replay_save_sp3/expected/BKG0MGXRTS_20253210000_01D_30S_ORB.SP3.gz (modified) ( previous)

trunk/BNC/Example_Configs/10_PPPQuickStart.bnc

@@ -33 +33 @@
 miscPort=
 miscScanRTCM=0
-mountPoints=//Example:Configs@euref-ip.net:2101/WTZR00DEU0 RTCM_3.3 DEU 49.14 12.88 no 2, //Example:Configs@products.igs-ip.net:2101/BCEP00BKG0 RTCM_3 DEU 50.09 8.66 no 2, //Example:Configs@products.igs-ip.net:2101/SSRA01CAS1 RTCM_3.1 CHN 40.07 116.28 no 2
+mountPoints=//Example:Configs@euref-ip.net:2101/WTZR00DEU0 RTCM_3.3 DEU 49.14 12.88 no 2, //Example:Configs@products.igs-ip.net:2101/BCEP00BKG0 RTCM_3 DEU 50.09 8.66 no 2, //Example:Configs@products.igs-ip.net:2101/SSRA00CNE1 RTCM_3.1 FRA 43.56 1.48 no 2
 ntripVersion=2
 onTheFlyInterval=no
@@ -151 +151 @@
 constraints=Ionosphere: pseudo-obs
 corrFile=
-corrMount=SSRA00CAS1
+corrMount=SSRA00CNE1
 corrWaitTime=5
 crdFile=Input/BNC_IGS20.CRD

trunk/BNC/Example_Configs/17_Combi.bnc

@@ -17 +17 @@
 cmbNavic=0
 cmbQzss=0
-cmbSampl=1 sec
+cmbSampl=5 sec
 cmbSbas=0
 cmbStreams="SSRA00BKG1 BKG 1.0  Computed ", "SSRA00CNE1 CNE 1.0  Computed "

trunk/BNC/Example_Configs/18_CombiPPP.bnc

@@ -17 +17 @@
 cmbNavic=0
 cmbQzss=0
-cmbSampl=1 sec
+cmbSampl=5 sec
 cmbSbas=0
 cmbStreams="SSRA00BKG1 BKG 1.0  Computed ", "SSRA00CNE1 CNE 1.0  Computed "

trunk/BNC/src/RTCM3/RTCM3coDecoder.cpp

@@ -1132 +1132 @@
   // Set _lastTime close to currentTime
   // ----------------------------------
+  if (_type == IGSssr) {
+    if      (epoSecGPS  != -1) _igsRefEpoSec = epoSecGPS;
+    else if (epoSecGal  != -1) _igsRefEpoSec = epoSecGal;
+    else if (epoSecQzss != -1) _igsRefEpoSec = epoSecQzss;
+    else if (epoSecSbas != -1) _igsRefEpoSec = epoSecSbas;
+  }
+
   if      (epoSecGPS != -1) {
     _lastTime.set(currentWeek, epoSecGPS);
@@ -1137 +1144 @@
   else if (epoSecGlo != -1) {
     QDate date = dateAndTimeFromGPSweek(currentTime.gpsw(), currentTime.gpssec()).date();
+    double leap = gnumleap(date.year(), date.month(), date.day());
     // IGS-SSR: the SSR Epoch Time (IDF003, 20 bit) is given in one continuous
     // time scale for all GNSS, i.e. GPS seconds of week also for GLONASS
-    // (IGS SSR v1.00) - no conversion
+    // (IGS SSR v1.00) - no conversion. Old encoders (e.g. BNC before
+    // 2026-09-30) wrote UTC seconds of day instead: recognised by comparison
+    // with the GPS/Galileo epoch of the same stream and converted, with a
+    // warning, during a transition period.
+    if (_type == IGSssr && _igsRefEpoSec >= 0.0 &&
+        igsLegacyEpoch(epoSecGlo, fmod(_igsRefEpoSec - leap + 7 * 86400.0, 86400.0), 86400.0)) {
+      double sow = floor(_igsRefEpoSec / 86400.0) * 86400.0 + epoSecGlo + leap;
+      if      (sow - _igsRefEpoSec >  43200.0) sow -= 86400.0;
+      else if (_igsRefEpoSec - sow >  43200.0) sow += 86400.0;
+      epoSecGlo = sow;
+      bncTime& warned = _legacyEpochWarned['R'];
+      if (!warned.valid() || currentTime - warned >= 3600.0) {
+        warned = currentTime;
+        emit newMessage(_staID.toLatin1() + ": IGS-SSR GLONASS epoch time in outdated convention "
+                        "(UTC seconds of day instead of GPS seconds of week, IGS SSR v1.00 IDF003) "
+                        "detected and converted - the SSR encoder should be updated", true);
+      }
+    }
     if (_type == RTCMssr || _type == RTCMnewssr) {
       // Same GLONASS epoch wire convention (UTC+3h) as RTCMssr
-      epoSecGlo = epoSecGlo - 3 * 3600 + gnumleap(date.year(), date.month(), date.day());
+      epoSecGlo = epoSecGlo - 3 * 3600 + leap;
     }
     _lastTime.set(currentWeek, epoSecGlo);
@@ -1156 +1181 @@
   }
   else if (epoSecBds != -1) {
-    // RTCM-SSR: BDT; IGS-SSR: GPS time scale for all GNSS (IDF003)
-    if (_type != IGSssr) {
+    // RTCM-SSR: BDT; IGS-SSR: GPS time scale for all GNSS (IDF003). Old
+    // IGS-SSR encoders wrote BDT: recognised and converted as for GLONASS.
+    bool bdt = (_type != IGSssr);
+    if (_type == IGSssr && _igsRefEpoSec >= 0.0 &&
+        igsLegacyEpoch(epoSecBds, fmod(_igsRefEpoSec - 14.0 + 604800.0, 604800.0), 604800.0)) {
+      bdt = true;
+      bncTime& warned = _legacyEpochWarned['C'];
+      if (!warned.valid() || currentTime - warned >= 3600.0) {
+        warned = currentTime;
+        emit newMessage(_staID.toLatin1() + ": IGS-SSR BDS epoch time in outdated convention "
+                        "(BDT instead of GPS seconds of week, IGS SSR v1.00 IDF003) "
+                        "detected and converted - the SSR encoder should be updated", true);
+      }
+    }
+    if (bdt) {
       epoSecBds += 14.0;
       if (epoSecBds > 604800.0) {
@@ -1182 +1220 @@
     }
   }
+}
+
+// IGS-SSR: true if value (seconds, periodic with period) matches the epoch
+// expected in the outdated convention (legacyRef) better than the GPS time
+// reference of the stream, within the spread of epochs of one update cycle
+////////////////////////////////////////////////////////////////////////////
+bool RTCM3coDecoder::igsLegacyEpoch(double value, double legacyRef, double period) const {
+  auto dist = [period](double a, double b) {
+    double d = fmod(fabs(a - b), period);
+    return std::min(d, period - d);
+  };
+  double dSpec   = dist(value, fmod(_igsRefEpoSec, period));
+  double dLegacy = dist(value, legacyRef);
+  return dLegacy < dSpec && dLegacy <= 30.0;
 }
 

trunk/BNC/src/RTCM3/RTCM3coDecoder.h

@@ -82 +82 @@
   void checkPrnRange(const t_prn& prn);
   QSet<QString>                         _prnOutOfRangeLogged;
+  // IGS-SSR: latest GPS/Galileo/QZSS/SBAS epoch (GPS seconds of week) of
+  // this stream, reference to detect GLONASS/BDS epochs in the outdated
+  // convention of old encoders (GLONASS: UTC seconds of day, BDS: BDT)
+  double                                _igsRefEpoSec = -1.0;
+  QMap<QChar, bncTime>                  _legacyEpochWarned; // warned at most once per hour
+  bool   igsLegacyEpoch(double value, double legacyRef, double period) const;
 
 

trunk/BNC/src/bnchelp.html

@@ -8105 +8105 @@
     Specify a combination sampling interval. Orbit and clock corrections will be produced following that interval.
     A value of 10 sec may be an appropriate choice.</p>
+  <p>
+    As the incoming corrections of the individual providers arrive with different latencies, BNC combines an epoch as
+    soon as each provider has delivered corrections for a later epoch, i.e. has completed the epoch. The combined
+    product is thus only delayed by the latency of the slowest provider. If a provider is late or missing, the epoch is
+    combined without it at the latest when the newest incoming corrections are 10 seconds (or the sampling interval,
+    if larger) younger than the epoch. Corrections arriving after an epoch was combined are ignored.</p>
 
   <p>
@@ -9108 +9114 @@
     A value of 10 sec may be appropriate.
     A value of zero '0' tells BNC to store all available samples into Clock RINEX files.
+  </p>
+  <p>
+    SP3 and Clock RINEX epochs are multiples of the sampling interval. Corrections whose epochs are not on this grid,
+    e.g. Galileo HAS corrections at ...:02, ...:12 seconds, are written at the preceding grid epoch if they refer to an
+    epoch at most 5 seconds (and less than half the sampling interval) after it. The satellite clock is then propagated
+    to the grid epoch using its rate, the orbit using the broadcast orbit at the grid epoch plus the difference between
+    precise and broadcast orbit at the correction epoch. Corrections on the grid are always used unchanged.
   </p>
 

trunk/BNC/src/bncoutf.h

@@ -16 +16 @@
   QString agencyFromFileName();
   bool isProductFile();
+  int  sampl() const {return _sampl;}
  protected:
   virtual t_irc reopen(int GPSweek, double GPSweeks);

trunk/BNC/src/combination/bnccomb.cpp

@@ -535 +535 @@
 
 
+// True if every AC has delivered clock corrections for an epoch after epoTime
+////////////////////////////////////////////////////////////////////////////
+bool bncComb::allACsBeyond(const bncTime& epoTime) const {
+  QListIterator<cmbAC*> itAC(_ACs);
+  while (itAC.hasNext()) {
+    const cmbAC* AC = itAC.next();
+    bncTime acLast = _acLastClkTime.value(AC->name);
+    if (!acLast.valid() || !(acLast > epoTime)) {
+      return false;
+    }
+  }
+  return true;
+}
+
 // Process clock corrections
 ////////////////////////////////////////////////////////////////////////////
@@ -552 +566 @@
     if (!_cmbSysPrn.contains(sys)){
       continue;
+    }
+
+    // Remember the latest correction epoch of each AC (whether or not the
+    // correction is used), to know when an AC has finished an epoch
+    // --------------------------------------------------------------------
+    QListIterator<cmbAC*> itAC(_ACs);
+    while (itAC.hasNext()) {
+      cmbAC* AC = itAC.next();
+      if (AC->mountPoint == QString(newClk._staID.c_str())) {
+        bncTime& acLast = _acLastClkTime[AC->name];
+        if (!acLast.valid() || newClk._time > acLast) {
+          acLast = newClk._time;
+        }
+        break;
+      }
     }
 
@@ -665 +694 @@
       }
     }
-    // Process the front epoch
-    // -----------------------
-    if (epoTime < (_lastClkCorrTime - outWait)) {
+    // Process the front epoch as soon as all ACs have delivered corrections
+    // for a later epoch (i.e. have finished this one), at the latest when
+    // the newest correction is maxWait newer (an AC may be late or missing).
+    // Waiting only outWait after the newest correction lost the corrections
+    // of ACs with a few seconds more latency ("old correction").
+    // -----------------------------------------------------------------------
+    const double maxWait = max(10.0, outWait);
+    if (allACsBeyond(epoTime) || epoTime < (_lastClkCorrTime - maxWait)) {
       _resTime = epoTime;
       processEpoch(_resTime, clkCorrVec);
@@ -834 +868 @@
               codeBiasesRefSig[frqType] = bias._value;
             }
+          }
+        }
+        // Code biases present, but not for both reference signals (e.g.
+        // Galileo HAS): the AC clocks cannot be referred to the reference
+        // signals and are used as they are - warn (at most once per hour
+        // per AC and system)
+        if (codeBiasesRefSig.size() < 2) {
+          QString key = acName + sys;
+          bncTime& warned = _refBiasWarned[key];
+          if (!warned.valid() || _newCorr->_time - warned >= 3600.0) {
+            warned = _newCorr->_time;
+            QString refSig;
+            for (unsigned ii = 1; ii < cmbRefSig::cIF; ii++) {
+              t_frequency::type frqType = cmbRefSig::toFreq(sys, static_cast<cmbRefSig::type>(ii));
+              refSig += QString(" C%1%2").arg(t_frequency::toString(frqType)[1])
+                                         .arg(cmbRefSig::toAttrib(sys, static_cast<cmbRefSig::type>(ii)));
+            }
+            emit newMessage("bncComb: " + acName.toLatin1() + " provides no code biases for the " +
+                            QByteArray(1, sys) + " reference signals" + refSig.toLatin1() +
+                            " - its clocks are combined without referring them to these signals", false);
           }
         }
@@ -1398 +1452 @@
     double dC_apriori_corr = dC_orb;// + dC_att; can be added if useful values are provided by the ACs
 
-    ll(iObs) = ((corr->_clkCorr._dClk * t_CST::c) - dC_apriori_corr - corr->_satCodeBiasIF) - DotProduct(AA.Row(iObs), x0);
+    // Clock referred to the ionosphere-free combination of the reference
+    // signals: dt - (alpha b1 + beta b2) with biases in Bias-SINEX sign
+    // (subtracted from observations, Banville et al. 2020, Eq. 1). The code
+    // biases taken from the SSR streams have the opposite sign (added to
+    // observations), hence + _satCodeBiasIF.
+    ll(iObs) = ((corr->_clkCorr._dClk * t_CST::c) - dC_apriori_corr + corr->_satCodeBiasIF) - DotProduct(AA.Row(iObs), x0);
 
     PP(iObs, iObs) *= 1.0 / (corr->_weightFactor * corr->_weightFactor);

trunk/BNC/src/combination/bnccomb.h

@@ -267 +267 @@
   std::deque<epoClkData*>                    _epoClkData;
   bncTime                                    _lastClkCorrTime;
+  QMap<QString, bncTime>                     _acLastClkTime; // latest correction epoch per AC
+  QMap<QString, bncTime>                     _refBiasWarned; // last warning per AC and system
+  bool allACsBeyond(const bncTime& epoTime) const;
   bncTime                                    _resTime;
   cmbCorr*                                   _newCorr;

trunk/BNC/src/upload/bncrtnetuploadcaster.cpp

@@ -548 +548 @@
             // Make sure the clock messages refer to same IOD as orbit messages
             // ----------------------------------------------------------------
+            // The ephemeris is selected at the first correction epoch of each
+            // sampling interval: the grid epoch itself for corrections on the
+            // sampling grid, otherwise (e.g. Galileo HAS at ...:02, ...:12)
+            // the first epoch after it, which previously never selected one
             if (_usedEph) {
-                if (fmod(epoTime.gpssec(), _samplRtcmEphCorr) == 0.0) {
+                bool select = (fmod(epoTime.gpssec(), _samplRtcmEphCorr) == 0.0);
+                if (!select && _samplRtcmEphCorr > 0.0) {
+                    select = !_usedEphEpoch.valid() || epoTime == _usedEphEpoch ||
+                             epoTime.gpsw() != _usedEphEpoch.gpsw() ||
+                             floor(epoTime.gpssec() / _samplRtcmEphCorr) !=
+                             floor(_usedEphEpoch.gpssec() / _samplRtcmEphCorr);
+                }
+                if (select) {
                     (*_usedEph)[prnInternalStr] = eph;
+                    _usedEphEpoch = epoTime;
                 }
                 else {
@@ -1170 +1182 @@
     double clkRnxAcc = rtnClk[2] / t_CST::c;         // [s/s² = -/s]
 
-    if (_rnx) {
+    // Corrections on epochs off the output sampling grid (e.g. Galileo HAS
+    // at ...:02, ...:12) are written at the nearest grid epoch, see gridEpoch;
+    // the clock is propagated with its rate, the orbit via the broadcast
+    // orbit at the grid epoch plus the (slowly varying) difference
+    // precise - broadcast at the correction epoch
+    // -----------------------------------------------------------------------
+    bncTime tt(GPSweek, GPSweeks);
+    bncTime tg;
+    double  clkRate = (clkRnxRate != 0.0) ? clkRnxRate : xB(5);
+
+    if (_rnx && gridEpoch(_rnx->sampl(), tt, tg, _rnxGridEpo, _rnxGridSrc)) {
         double clkRnxSig, clkRnxRateSig, clkRnxAccSig;
         int s = rtnClkSig.size();
@@ -1190 +1212 @@
             break;
         }
-        _rnx->write(GPSweek, GPSweeks, prn, clkRnx, clkRnxRate, clkRnxAcc,
-            clkRnxSig, clkRnxRateSig, clkRnxAccSig);
-    }
-    if (_sp3) {
-        _sp3->write(GPSweek, GPSweeks, prn, rtnCoM, clkRnx, rtnVel, clkRnxRate);
+        _rnx->write(tg.gpsw(), tg.gpssec(), prn, clkRnx + clkRate * (tg - tt),
+            clkRnxRate, clkRnxAcc, clkRnxSig, clkRnxRateSig, clkRnxAccSig);
+    }
+    if (_sp3 && gridEpoch(_sp3->sampl(), tt, tg, _sp3GridEpo, _sp3GridSrc)) {
+        ColumnVector xyz = rtnCoM;
+        ColumnVector vel = rtnVel;
+        if (!(tg == tt)) {
+            ColumnVector xBg(6);
+            ColumnVector vBg(3);
+            if (eph->getCrd(tg, xBg, vBg, false) != success) {
+                return success;
+            }
+            xyz = xBg.Rows(1, 3) + (rtnCoM - xB.Rows(1, 3));
+            vel = vBg + (rtnVel - vB);
+        }
+        _sp3->write(tg.gpsw(), tg.gpssec(), prn, xyz, clkRnx + clkRate * (tg - tt),
+            vel, clkRnxRate);
     }
     return success;
+}
+
+// Output epoch for SP3/CLK files with sampling interval sampl [s]: the
+// correction epoch itself if it is on the sampling grid, otherwise the
+// preceding grid epoch if at most 5 s (and less than half the sampling
+// interval) earlier. Only preceding grid epochs are used, so a correction
+// epoch on the grid always arrives first and takes precedence. Each grid
+// epoch is derived from one correction epoch only.
+////////////////////////////////////////////////////////////////////////////
+bool bncRtnetUploadCaster::gridEpoch(int sampl, const bncTime& tt, bncTime& tg,
+                                     bncTime& lastGrid, bncTime& lastSrc) const {
+    if (sampl <= 0) {
+        tg = tt;
+        return true;
+    }
+    double sec  = tt.gpssec();
+    double back = sec - floor(sec / sampl) * sampl;   // [0, sampl)
+    if (back > 1.e-6 && (back > 5.0 || back >= 0.5 * sampl)) {
+        return false;
+    }
+    tg = (back <= 1.e-6) ? tt : tt - back;
+    if (lastGrid.valid() && tg == lastGrid && !(tt == lastSrc)) {
+        return false;
+    }
+    lastGrid = tg;
+    lastSrc  = tt;
+    return true;
 }
 

trunk/BNC/src/upload/bncrtnetuploadcaster.h

@@ -103 +103 @@
     bncClockRinex* _rnx;
     bncSP3* _sp3;
+    // grid epoch last written to SP3/CLK and the epoch it was derived from
+    bncTime        _sp3GridEpo, _sp3GridSrc, _rnxGridEpo, _rnxGridSrc;
+    bncTime        _usedEphEpoch; // epoch the used ephemerides were last selected
+    bool gridEpoch(int sampl, const bncTime& tt, bncTime& tg,
+                   bncTime& lastGrid, bncTime& lastSrc) const;
     bncBiasSinex* _bsx;
     QMap<QString, const t_eph*>* _usedEph;

trunk/BNC/tests/README.md

@@ -46 +46 @@
   `replay_rinex_eph` (02), `replay_feed_engine` (07, IP port disabled) and
   `replay_save_sp3` (14, with SSRA00CNE1 instead of SSRA00BKG1).
+* `replay_save_sp3_has`: 14 with Galileo HAS (SSRA00EUH0) only; its epochs are
+  off the SP3/CLK sampling grid and are written at the preceding grid epoch.
 * Real-time PPP replays (`ppp_rt_*`) from 11 min slices of the CNES and WHU
   captures (`data/wtzr_cnes_11min.raw`, `data/wtzr_whu_11min.raw`):
-  `ppp_rt_float` (08), `ppp_rt_quickstart` (10, with SSRA00CNE1 instead of
-  SSRA01CAS1, NMEA port disabled), `ppp_rt_ar_cnes` (29) and `ppp_rt_ar_whu`
+  `ppp_rt_float` (08), `ppp_rt_quickstart` (10, NMEA port disabled), `ppp_rt_ar_cnes` (29) and `ppp_rt_ar_whu`
   (30), compared exactly.
+* `ppp_rt_net` (09, two stations) and `ppp_rt_has` (13, Galileo HAS via
+  SSRA00EUH0) from `data/merged_260930_11min.raw`, a slice of two simultaneous
+  recordings merged by time (BCEP00BKG0, WTZR00DEU0, FFMJ00DEU0, SSRA00CNE1,
+  SSRA00EUH0).
+* Combination (`combi`, `combi_ppp`; 17 and 18) from `data/combi_260930_11min.raw`
+  (SSRA00BKG1, SSRA00CNE1, BCEP00BKG0, FFMJ00DEU0), upload to the caster
+  disabled, SP3/CLK written to files.
 * SSR decoding (`ssr_*`): 3 min slices of real streams replayed through
   `03_BrdcCorr.bnc`. They compare the `.ssr` files and the `BNC_DEBUG_SSR`

trunk/BNC/tests/cases/ppp_rt_quickstart/case.json

@@ -1 +1 @@
 {
-  "description": "10_PPPQuickStart replay: Quick-Start PPP with SSRA00CNE1 instead of SSRA01CAS1, NMEA port disabled, 11 min",
+  "description": "10_PPPQuickStart replay: Quick-Start PPP with SSRA00CNE1, NMEA port disabled, 11 min",
   "config": "Example_Configs/10_PPPQuickStart.bnc",
   "raw": "data/wtzr_cnes_11min.raw",
@@ -9 +9 @@
     "rawOutFile",
     "",
-    "-key",
-    "PPP/corrMount",
-    "SSRA00CNE1",
     "-key",
     "PPP/staTable",