Index: trunk/BNC/src/bncantex.cpp =================================================================== --- trunk/BNC/src/bncantex.cpp (revision 10945) +++ trunk/BNC/src/bncantex.cpp (revision 10946) @@ -40,4 +40,5 @@ #include +#include #include @@ -330,8 +331,90 @@ } +// 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,8 +457,34 @@ 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]; Index: trunk/BNC/src/bncantex.h =================================================================== --- trunk/BNC/src/bncantex.h (revision 10945) +++ trunk/BNC/src/bncantex.h (revision 10946) @@ -47,7 +47,26 @@ 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 _glonassYaw; + class t_frqMap { public: Index: trunk/BNC/src/bnchelp.html =================================================================== --- trunk/BNC/src/bnchelp.html (revision 10945) +++ trunk/BNC/src/bnchelp.html (revision 10946) @@ -5640,4 +5640,17 @@

+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. +

+

Note that clocks in the SP3 orbit files are not corrected for the conventional periodic relativistic effect.

Index: trunk/BNC/src/combination/bnccomb.cpp =================================================================== --- trunk/BNC/src/combination/bnccomb.cpp (revision 10945) +++ trunk/BNC/src/combination/bnccomb.cpp (revision 10946) @@ -1205,5 +1205,5 @@ 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';