Changeset 10946 in ntrip

Timestamp:

Jun 24, 2026, 4:39:53 PM (12 hours ago)

Author:

stuerze

Message:

added GLONASS-M yaw-fixed attitude model for the APC/CoM offset applied when saving SP3 files, reducing along-track/cross-track errors near the orbit noon/midnight points

Location:

trunk/BNC

Files:

• CHANGELOG.md (modified) (1 diff)
• src/bncantex.cpp (modified) (3 diffs)
• src/bncantex.h (modified) (1 diff)
• src/bnchelp.html (modified) (1 diff)
• src/combination/bnccomb.cpp (modified) (1 diff)

trunk/BNC/CHANGELOG.md

@@ -4 +4 @@
 - ADDED: L1C/B (RINEX code: '1E') in QZSS Signal ID Mapping
 - ADDED: Minimum Elevation parameter to BNC's RINEX Editing & QC feature [(#195)](https://software.rtcm-ntrip.org/ticket/195)
+- ADDED: GLONASS-M yaw-fixed attitude model for the APC/CoM offset applied when saving SP3 files, reducing along-track/cross-track errors near the orbit noon/midnight points
 - FIXED: BNC is now able to work with com ports with higher numbers, e.g. COM11 or COM17  [(#205)](https://software.rtcm-ntrip.org/ticket/205)
 - FIXED: Handling of ionospheric constraints [(#220)](https://software.rtcm-ntrip.org/ticket/220)
+- FIXED: SP3/Clock RINEX clock values produced by the Upload Corrections feature no longer depend on whether 'Center of Mass' is ticked, since that setting should only affect the uploaded SSR Broadcast Correction stream
 
 

trunk/BNC/src/bncantex.cpp

@@ -40 +40 @@
 
 #include <iostream>
+#include <cmath>
 #include <newmatio.h>
 
@@ -330 +331 @@
 }
 
+// GLONASS Yaw Angle (Sun-pointing law overridden near noon/midnight when
+// the satellite cannot mechanically keep up, following the GLONASS-M
+// "yaw-fixed" behaviour described in Dilssner, Springer, Flohrer, Dow
+// (2011), "The GLONASS-M satellite yaw-attitude model", Advances in Space
+// Research 47(1), 160-171.
+//
+// Unlike GPS/Galileo/BeiDou, which are commonly approximated by the
+// nominal Sun-pointing yaw-steering law of Bar-Sever (1996) at all times,
+// GLONASS-M has been found to stop tracking that law and hold a constant
+// (frozen) yaw angle whenever the required yaw rate would exceed the
+// satellite's slew capability - which happens close to the orbit
+// noon/midnight points whenever the Sun's elevation above the orbital
+// plane (the "beta" angle) is small. The maximum yaw rate used below
+// (0.25 deg/s) and the general approach follow that paper; satellite
+// telemetry was not available to validate the exact rate against this
+// installation's GLONASS satellites, so it should be checked against
+// independently-known attitude/orbit residuals if high accuracy matters.
+////////////////////////////////////////////////////////////////////////////
+double bncAntex::glonassYawAngle(const QString& prn, const ColumnVector& xSat,
+                                  const ColumnVector& vSat,
+                                  const ColumnVector& xSun) {
+
+  const double MAX_YAW_RATE = 0.25 * M_PI / 180.0; // [rad/s], approximate
+
+  // Inertial-consistent velocity (xSat, vSat are Earth-fixed; remove the
+  // Earth-rotation contribution so that the orbit normal below is not
+  // contaminated by it)
+  // -------------------------------------------------------------------
+  ColumnVector Omega(3); Omega(1) = 0.0; Omega(2) = 0.0; Omega(3) = t_CST::omega;
+  ColumnVector vInert = vSat + crossproduct(Omega, xSat);
+
+  // Orbit normal and instantaneous orbital rate from r x v (exact, valid
+  // for any Keplerian orbit, not just circular ones)
+  // ---------------------------------------------------------------------
+  ColumnVector h = crossproduct(xSat, vInert);
+  double       hNorm = sqrt(DotProduct(h, h));
+  ColumnVector orbNormal = h / hNorm;
+  double       r = sqrt(DotProduct(xSat, xSat));
+  double       nRate = hNorm / (r * r); // [rad/s]
+
+  // Beta angle (Sun elevation above the orbital plane)
+  // ----------------------------------------------------
+  double beta = asin(DotProduct(orbNormal, xSun));
+
+  // Orbit angle mu, measured from the orbit midnight point, increasing in
+  // the direction of satellite motion
+  // -----------------------------------------------------------------------
+  ColumnVector sunProj = xSun - DotProduct(xSun, orbNormal) * orbNormal;
+  sunProj /= sqrt(DotProduct(sunProj, sunProj));
+  ColumnVector eX = -1.0 * sunProj;                 // midnight direction, mu = 0
+  ColumnVector eY = crossproduct(orbNormal, eX);
+  ColumnVector rHat = xSat / r;
+  double mu = atan2(DotProduct(rHat, eY), DotProduct(rHat, eX));
+
+  // Nominal Sun-pointing yaw angle and its rate (Bar-Sever 1996)
+  // ----------------------------------------------------------------
+  double tanBeta = tan(beta);
+  double sinMu   = sin(mu);
+  double psiNom  = atan2(-tanBeta, sinMu);
+  double denom   = tanBeta * tanBeta + sinMu * sinMu;
+  double psiRate = (denom > 1e-12)
+                  ? nRate * fabs(tanBeta * cos(mu)) / denom
+                  : 1e9;
+
+  t_glonassYaw& st = _glonassYaw[prn];
+  double psiEff;
+  if (psiRate > MAX_YAW_RATE) {
+    if (!st.valid) {
+      st.yaw = psiNom; // no prior history - best available estimate
+    }
+    psiEff = st.yaw;    // hold the frozen yaw angle
+  }
+  else {
+    st.yaw   = psiNom;
+    st.valid = true;
+    psiEff   = psiNom;
+  }
+
+  return psiEff;
+}
+
 // Satellite Antenna Offset
 ////////////////////////////////////////////////////////////////////////////
 t_irc bncAntex::satCoMcorrection(const QString& prn, double Mjd,
-                                 const ColumnVector& xSat, ColumnVector& dx) {
+                                 const ColumnVector& xSat,
+                                 const ColumnVector& vSat, ColumnVector& dx) {
 
   t_frequency::type frqType = t_frequency::dummy;
@@ -374 +457 @@
       xSun /= sqrt(DotProduct(xSun,xSun));
 
-      ColumnVector sy = crossproduct(sz, xSun);
-      sy /= sqrt(DotProduct(sy,sy));
-
-      ColumnVector sx = crossproduct(sy, sz);
+      ColumnVector sy, sx;
+
+      // GLONASS: override the nominal Sun-pointing attitude near the
+      // orbit noon/midnight points when the beta angle is small (see
+      // glonassYawAngle() above). Elsewhere GLONASS-M follows the same
+      // nominal law as the other constellations, so the result is
+      // identical to the direct Sun-pointing computation used below.
+      // -----------------------------------------------------------------
+      if (prn[0] == 'R' && vSat.size() == 3) {
+        double psi = glonassYawAngle(prn, xSat, vSat, xSun);
+
+        ColumnVector vInert = vSat;
+        ColumnVector Omega(3); Omega(1) = 0.0; Omega(2) = 0.0; Omega(3) = t_CST::omega;
+        vInert += crossproduct(Omega, xSat);
+
+        ColumnVector sy0 = crossproduct(sz, vInert);
+        sy0 /= sqrt(DotProduct(sy0,sy0));
+        ColumnVector sx0 = crossproduct(sy0, sz);
+
+        // Rodrigues rotation of (sx0, sy0) around sz by angle psi
+        double cosY = cos(psi);
+        double sinY = sin(psi);
+        sx = sx0 * cosY + crossproduct(sz, sx0) * sinY;
+        sy = sy0 * cosY + crossproduct(sz, sy0) * sinY;
+      }
+      else {
+        sy = crossproduct(sz, xSun);
+        sy /= sqrt(DotProduct(sy,sy));
+        sx = crossproduct(sy, sz);
+      }
 
       dx[0] = sx[0] * neu[0] + sy[0] * neu[1] + sz[0] * neu[2];

trunk/BNC/src/bncantex.h

@@ -47 +47 @@
                   double eleSat, double azSat, bool& found) const;
   t_irc   satCoMcorrection(const QString& prn, double Mjd,
-                           const ColumnVector& xSat, ColumnVector& dx);
+                           const ColumnVector& xSat, const ColumnVector& vSat,
+                           ColumnVector& dx);
 
  private:
+  // Per-satellite GLONASS yaw state, used to freeze the yaw angle while
+  // the satellite cannot follow the nominal Sun-pointing attitude law
+  // (see glonassYawAngle() below).
+  class t_glonassYaw {
+   public:
+    t_glonassYaw() {
+      yaw   = 0.0;
+      valid = false;
+    }
+    double yaw;
+    bool   valid;
+  };
+
+  double glonassYawAngle(const QString& prn, const ColumnVector& xSat,
+                          const ColumnVector& vSat, const ColumnVector& xSun);
+
+  QMap<QString, t_glonassYaw> _glonassYaw;
+
   class t_frqMap {
    public:

trunk/BNC/src/bnchelp.html

@@ -5640 +5640 @@
 </p>
 <p>
+For GLONASS satellites, that APC/CoM offset is rotated into the satellite-fixed frame using the satellite's
+yaw attitude. BNC assumes the nominal Sun-pointing yaw-steering law used for the other GNSS systems, except
+close to the orbit noon and midnight points when the Sun's elevation above the orbital plane (the 'beta'
+angle) is small. There, GLONASS-M satellites are known to stop tracking that law and instead hold the yaw
+angle fixed, following the model described in Dilssner, Springer, Flohrer, Dow (2011), 'The GLONASS-M
+satellite yaw-attitude model', Advances in Space Research 47(1), 160-171. Without that correction, the
+converted GLONASS CoM position can be off by several decimeters in the along-track and cross-track
+components in that situation. The maximum yaw rate used to decide when GLONASS-M can no longer follow the
+nominal law (0.25 deg/s) is taken from the literature and not calibrated against any specific satellite, so
+results should be checked against independently known attitude or orbit information where high accuracy on
+GLONASS is required.
+</p>
+<p>
  Note that clocks in the SP3 orbit files are not corrected for the conventional periodic relativistic effect.
 </p>

trunk/BNC/src/combination/bnccomb.cpp

@@ -1205 +1205 @@
       char sys = corr->_eph->prn().system();
       masterIsAPC = _masterIsAPC[sys];
-      if (_antex->satCoMcorrection(corr->_prn, Mjd, xc.Rows(1,3), dx) != success) {
+      if (_antex->satCoMcorrection(corr->_prn, Mjd, xc.Rows(1,3), vv, dx) != success) {
         dx = 0;
         _log += "antenna not found " + corr->_prn.mid(0,3).toLatin1() + '\n';