Changeset 10946 in ntrip for trunk/BNC/src/bncantex.cpp

Timestamp:

Jun 24, 2026, 4:39:53 PM (13 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

File:

• trunk/BNC/src/bncantex.cpp (modified) (3 diffs)

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];