Changeset 10955 in ntrip


Ignore:
Timestamp:
Jul 3, 2026, 11:10:42 AM (12 hours ago)
Author:
stuerze
Message:

Possibility to select how satellite attitude is modelled when converting Antenna Phase Center (APC) corrections to
Center-of-Mass (CoM) positions required for SP3 output

Location:
trunk/BNC
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/BNC/CHANGELOG.md

    r10947 r10955  
    44- ADDED: L1C/B (RINEX code: '1E') in QZSS Signal ID Mapping
    55- ADDED: Minimum Elevation parameter to BNC's RINEX Editing & QC feature [(#195)](https://software.rtcm-ntrip.org/ticket/195)
    6 - 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 [(#210)](https://software.rtcm-ntrip.org/ticket/210)
     6- ADDED: Possibility to select how satellite attitude is modelled when converting Antenna Phase Center (APC) corrections to
     7Center-of-Mass (CoM) positions required for SP3 output. Three options are available: (1) Computed (default): BNC applies its own kinematic attitude model: GPS noon/midnight turn manoeuvres (Kouba 2009/2015, Bar-Sever 1996), GLONASS yaw-fixed mode (Dilssner et al. 2011), and
     8Galileo / BDS orbit-normal mode switching (Kouba 2017, Dai et al. 2015, Steigenberger et al. 2018). (2) Nominal: a simplified, continuous Sun-pointing model is used without any manoeuvre modelling. (3) SSR: the yaw angle transmitted in the SSR phase bias message is used directly, if present for the satellite and epoch.  [(#210)](https://software.rtcm-ntrip.org/ticket/210)
    79- 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)
    810- FIXED: Handling of ionospheric constraints [(#220)](https://software.rtcm-ntrip.org/ticket/220)

  • trunk/BNC/src/bncantex.cpp

    r10951 r10955  
    192192            line.indexOf("NavIC") == 0 ){
    193193          newAntMap->antName = line.mid(20,3);
     194          if (line.indexOf("BLOCK I") == 0) {
     195            // Extract GPS block type: "BLOCK IIF   " → "IIF"
     196            QString bt = line.mid(0, 20).trimmed();
     197            if (bt.startsWith("BLOCK "))
     198              newAntMap->blockType = bt.mid(6);
     199          }
     200          else if (line.indexOf("GALILEO") == 0) {
     201            // "GALILEO-1" → "1" (IOV), "GALILEO-2" → "2" (FOC)
     202            QString bt = line.mid(0, 20).trimmed();
     203            if (bt.startsWith("GALILEO-"))
     204              newAntMap->blockType = bt.mid(8);
     205          }
     206          else if (line.indexOf("BEIDOU") == 0) {
     207            // "BEIDOU-2M" → "2M", "BEIDOU-3M-CAST" → "3M-CAST", etc.
     208            QString bt = line.mid(0, 20).trimmed();
     209            if (bt.startsWith("BEIDOU-"))
     210              newAntMap->blockType = bt.mid(7);
     211          }
    194212        }
    195213        else {
    … …  
    448466}
    449467
     468// GPS satellite yaw angle during nominal tracking and noon/midnight turns.
     469//
     470// GPS satellites perform a "noon turn" (and, for some blocks, a "midnight
     471// turn") when the Sun's elevation angle above the orbital plane (beta) is
     472// small: the required nominal yaw rate then exceeds the satellite's
     473// mechanical maximum, so the satellite yaws at that maximum rate until it
     474// catches back up to the nominal Sun-pointing orientation (Kouba 2009/2015,
     475// Bar-Sever 1996).
     476//
     477// The max yaw rates below are the best published estimates per block type:
     478//   IIA:   0.12 °/s (Kouba 2009)
     479//   IIR:   0.20 °/s (Bar-Sever 1996)
     480//   IIR-M: 0.20 °/s
     481//   IIF:   0.11 °/s (Kouba 2015)
     482//   IIIA:  0.15 °/s (tentative)
     483//
     484// Returns the effective yaw angle [rad] in the velocity-referenced frame,
     485// for use with the same Rodrigues rotation as the GLONASS model. During
     486// nominal tracking this equals psiNom and the result is identical to the
     487// simple sz×xSun formula.
     488////////////////////////////////////////////////////////////////////////////
     489double bncAntex::gpsYawAngle(const QString& prn, const QString& blockType,
     490                              double Mjd,
     491                              const ColumnVector& xSat,
     492                              const ColumnVector& vSat,
     493                              const ColumnVector& xSun) {
     494
     495  // Max yaw rate [rad/s] by GPS block type
     496  double psiDotMax;
     497  if      (blockType == "IIA")   psiDotMax = 0.12 * M_PI / 180.0;
     498  else if (blockType == "IIR")   psiDotMax = 0.20 * M_PI / 180.0;
     499  else if (blockType == "IIR-M") psiDotMax = 0.20 * M_PI / 180.0;
     500  else if (blockType == "IIF")   psiDotMax = 0.11 * M_PI / 180.0;
     501  else if (blockType == "IIIA")  psiDotMax = 0.15 * M_PI / 180.0;
     502  else return 0.0; // unknown block: caller uses simple Sun-pointing
     503
     504  const double MAX_CALL_GAP = 1800.0 / 86400.0; // 30 min in days
     505
     506  // Inertial velocity
     507  ColumnVector Omega(3); Omega(1) = 0.0; Omega(2) = 0.0; Omega(3) = t_CST::omega;
     508  ColumnVector vInert = vSat + crossproduct(Omega, xSat);
     509
     510  // Orbital angular momentum vector → orbital rate
     511  ColumnVector h     = crossproduct(xSat, vInert);
     512  double       hNorm = sqrt(DotProduct(h, h));
     513  ColumnVector orbNormal = h / hNorm;
     514  double       r    = sqrt(DotProduct(xSat, xSat));
     515  double       nRate = hNorm / (r * r); // [rad/s]
     516
     517  // Beta angle
     518  double beta = asin(DotProduct(orbNormal, xSun));
     519
     520  // Mu: orbit angle from midnight (same geometry as GLONASS)
     521  ColumnVector sunProj = xSun - DotProduct(xSun, orbNormal) * orbNormal;
     522  sunProj /= sqrt(DotProduct(sunProj, sunProj));
     523  ColumnVector eX = -1.0 * sunProj; // midnight direction
     524  ColumnVector eY = crossproduct(orbNormal, eX);
     525  ColumnVector rHat = xSat / r;
     526  double mu = atan2(DotProduct(rHat, eY), DotProduct(rHat, eX));
     527
     528  // Nominal yaw and its rate
     529  double tanBeta = tan(beta);
     530  double sinMu   = sin(mu);
     531  double psiNom  = atan2(-tanBeta, sinMu);
     532  double denom   = tanBeta * tanBeta + sinMu * sinMu;
     533  // |dPsi/dt| (always non-negative, sign comes from sign of tanBeta*cos(mu))
     534  double psiDotNomAbs = (denom > 1e-12)
     535                        ? nRate * fabs(tanBeta * cos(mu)) / denom
     536                        : 1e9;
     537  // Sign of the nominal yaw rate: dPsi/dt = nRate * tanBeta * cos(mu) / denom
     538  double psiDotNomSign = (tanBeta * cos(mu) >= 0.0) ? 1.0 : -1.0;
     539
     540  t_gpsYaw& st      = _gpsYaw[prn];
     541  double    callGap = st.valid ? (Mjd - st.lastCallMjd) : 0.0;
     542
     543  // Stale state: reset if there has been a gap in calls
     544  if (st.valid && callGap > MAX_CALL_GAP) {
     545    _gpsYawLog += QString().asprintf(
     546      "%s gpsYaw STALE-RESET  Mjd=%.6f beta=%6.3f mu=%7.2f"
     547      " old=%7.2f new=%7.2f gap=%.1fmin\n",
     548      prn.toLatin1().data(), Mjd, beta*180.0/M_PI, mu*180.0/M_PI,
     549      st.yaw*180.0/M_PI, psiNom*180.0/M_PI, callGap*1440.0);
     550    st.valid  = false;
     551    st.inTurn = false;
     552  }
     553
     554  double psiEff;
     555  if (psiDotNomAbs > psiDotMax) {
     556    // Constrained: satellite yaws at psiDotMax in the nominal direction
     557    if (!st.valid || !st.inTurn) {
     558      // Entering a noon/midnight turn
     559      psiEff = st.valid ? st.yaw : psiNom;
     560      _gpsYawLog += QString().asprintf(
     561        "%s gpsYaw TURN-ENTER   Mjd=%.6f beta=%6.3f mu=%7.2f"
     562        " psiNom=%7.2f psiEff=%7.2f psiDotNom=%6.3f max=%5.3f\n",
     563        prn.toLatin1().data(), Mjd, beta*180.0/M_PI, mu*180.0/M_PI,
     564        psiNom*180.0/M_PI, psiEff*180.0/M_PI,
     565        psiDotNomAbs*psiDotNomSign*180.0/M_PI, psiDotMax*180.0/M_PI);
     566      st.inTurn = true;
     567    } else {
     568      // Continuing the turn: integrate at constrained rate
     569      double dt  = callGap * 86400.0; // [s]
     570      psiEff = st.yaw + psiDotNomSign * psiDotMax * dt;
     571    }
     572    // Wrap to [-pi, pi]
     573    while (psiEff >  M_PI) psiEff -= 2.0 * M_PI;
     574    while (psiEff < -M_PI) psiEff += 2.0 * M_PI;
     575  }
     576  else {
     577    // Nominal tracking
     578    psiEff = psiNom;
     579    if (st.inTurn) {
     580      _gpsYawLog += QString().asprintf(
     581        "%s gpsYaw TURN-EXIT    Mjd=%.6f beta=%6.3f mu=%7.2f"
     582        " frozen=%7.2f nominal=%7.2f\n",
     583        prn.toLatin1().data(), Mjd, beta*180.0/M_PI, mu*180.0/M_PI,
     584        st.yaw*180.0/M_PI, psiNom*180.0/M_PI);
     585      st.inTurn = false;
     586    }
     587    st.yaw = psiNom;
     588  }
     589
     590  st.yaw         = psiEff;
     591  st.lastCallMjd = Mjd;
     592  st.valid       = true;
     593
     594  return psiEff;
     595}
     596
     597//
     598// Orbit-Normal Mode Yaw Angle — shared model for Galileo and BDS
     599//
     600// When |beta| < betaThr the satellite rotates toward yaw = 0 (orbit-normal)
     601// at the block-specific maximum yaw rate. When |beta| >= betaThr it returns
     602// to nominal yaw-steering (psiNom). Transitions are rate-limited in both
     603// directions to match physical satellite behaviour.
     604//
     605// References:
     606//   Galileo IOV/FOC: Kouba (2017), Steigenberger et al. (2018)
     607//   BDS MEO/IGSO:    Dai et al. (2015), Wang et al. (2018)
     608////////////////////////////////////////////////////////////////////////////
     609double bncAntex::onModeYawAngle(const QString& prn,
     610                                 double betaThr, double psiDotMax, double Mjd,
     611                                 const ColumnVector& xSat,
     612                                 const ColumnVector& vSat,
     613                                 const ColumnVector& xSun) {
     614
     615  const double MAX_CALL_GAP = 1800.0 / 86400.0;  // 30 min [days]
     616
     617  // Inertial velocity
     618  ColumnVector Omega(3); Omega(1) = 0.0; Omega(2) = 0.0; Omega(3) = t_CST::omega;
     619  ColumnVector vInert = vSat + crossproduct(Omega, xSat);
     620
     621  // Orbit geometry
     622  ColumnVector h        = crossproduct(xSat, vInert);
     623  double       hNorm    = sqrt(DotProduct(h, h));
     624  ColumnVector orbNormal = h / hNorm;
     625  double       r        = sqrt(DotProduct(xSat, xSat));
     626
     627  double beta = asin(DotProduct(orbNormal, xSun));
     628
     629  ColumnVector sunProj = xSun - DotProduct(xSun, orbNormal) * orbNormal;
     630  sunProj /= sqrt(DotProduct(sunProj, sunProj));
     631  ColumnVector eX   = -1.0 * sunProj;
     632  ColumnVector eY   = crossproduct(orbNormal, eX);
     633  ColumnVector rHat = xSat / r;
     634  double mu     = atan2(DotProduct(rHat, eY), DotProduct(rHat, eX));
     635  double psiNom = atan2(-tan(beta), sin(mu));
     636
     637  // Target yaw: orbit-normal (0) when below threshold, nominal otherwise
     638  bool   wantsON   = (fabs(beta) < betaThr);
     639  double psiTarget = wantsON ? 0.0 : psiNom;
     640
     641  t_onYaw& st      = _onYaw[prn];
     642  double   callGap = st.valid ? (Mjd - st.lastCallMjd) : 0.0;
     643
     644  if (st.valid && callGap > MAX_CALL_GAP) {
     645    _onYawLog += QString().asprintf(
     646      "%s onYaw STALE-RESET  Mjd=%.6f beta=%5.2f gap=%.1fmin\n",
     647      prn.toLatin1().data(), Mjd, beta*180.0/M_PI, callGap*1440.0);
     648    st.valid = false;
     649  }
     650
     651  double psiEff;
     652  if (!st.valid) {
     653    // Cold start: initialise at psiNom regardless of mode
     654    psiEff    = psiNom;
     655    st.inON   = false;
     656  }
     657  else {
     658    double dt      = callGap * 86400.0;       // [s]
     659    double dpsi    = psiTarget - st.yaw;
     660    while (dpsi >  M_PI) dpsi -= 2.0 * M_PI;
     661    while (dpsi < -M_PI) dpsi += 2.0 * M_PI;
     662    double maxDpsi = psiDotMax * dt;
     663
     664    if (fabs(dpsi) <= maxDpsi) {
     665      // Target reached this step
     666      psiEff = psiTarget;
     667      bool wasON = st.inON;
     668      st.inON    = wantsON;
     669      if (!wasON && wantsON && fabs(st.yaw) < 0.5*M_PI/180.0) {
     670        _onYawLog += QString().asprintf(
     671          "%s onYaw ON-ENTER   Mjd=%.6f beta=%5.2f psiNom=%7.2f\n",
     672          prn.toLatin1().data(), Mjd, beta*180.0/M_PI, psiNom*180.0/M_PI);
     673      }
     674      else if (wasON && !wantsON) {
     675        _onYawLog += QString().asprintf(
     676          "%s onYaw ON-EXIT    Mjd=%.6f beta=%5.2f psiNom=%7.2f\n",
     677          prn.toLatin1().data(), Mjd, beta*180.0/M_PI, psiNom*180.0/M_PI);
     678        st.inON = false;
     679      }
     680    }
     681    else {
     682      // Still rotating toward target
     683      psiEff = st.yaw + (dpsi > 0.0 ? 1.0 : -1.0) * maxDpsi;
     684    }
     685  }
     686
     687  while (psiEff >  M_PI) psiEff -= 2.0 * M_PI;
     688  while (psiEff < -M_PI) psiEff += 2.0 * M_PI;
     689
     690  st.yaw         = psiEff;
     691  st.lastCallMjd = Mjd;
     692  st.valid       = true;
     693
     694  return psiEff;
     695}
     696
     697//
     698////////////////////////////////////////////////////////////////////////////
     699double bncAntex::galileoYawAngle(const QString& prn, const QString& blockType,
     700                                  double Mjd,
     701                                  const ColumnVector& xSat,
     702                                  const ColumnVector& vSat,
     703                                  const ColumnVector& xSun) {
     704  // IOV (type "1") and FOC (type "2"): orbit-normal for |beta| < 2 deg.
     705  // Max yaw rate 0.20 deg/s applies to both generations.
     706  // Kouba (2017), Steigenberger et al. (2018)
     707  Q_UNUSED(blockType);
     708  const double betaThr   = 2.0 * M_PI / 180.0;
     709  const double psiDotMax = 0.20 * M_PI / 180.0;
     710  return onModeYawAngle(prn, betaThr, psiDotMax, Mjd, xSat, vSat, xSun);
     711}
     712
     713//
     714////////////////////////////////////////////////////////////////////////////
     715double bncAntex::bdsYawAngle(const QString& prn, const QString& blockType,
     716                              double Mjd,
     717                              const ColumnVector& xSat,
     718                              const ColumnVector& vSat,
     719                              const ColumnVector& xSun) {
     720  // GEO satellites are always in orbit-normal mode (yaw = 0).
     721  if (blockType == "2G" || blockType == "3G-CAST")
     722    return 0.0;
     723
     724  double betaThr, psiDotMax;
     725  if (blockType == "2M" || blockType == "2I") {
     726    // BDS-2 MEO / IGSO: orbit-normal for |beta| < 4 deg (Dai et al. 2015)
     727    betaThr   = 4.0 * M_PI / 180.0;
     728    psiDotMax = 0.10 * M_PI / 180.0;
     729  }
     730  else {
     731    // BDS-3 (CAST and SECM MEO/IGSO): orbit-normal for |beta| < 3 deg
     732    betaThr   = 3.0 * M_PI / 180.0;
     733    psiDotMax = 0.15 * M_PI / 180.0;
     734  }
     735  return onModeYawAngle(prn, betaThr, psiDotMax, Mjd, xSat, vSat, xSun);
     736}
     737
     738//
    450739// Satellite Antenna Offset
    451740////////////////////////////////////////////////////////////////////////////
    452741t_irc bncAntex::satCoMcorrection(const QString& prn, double Mjd,
    453742                                 const ColumnVector& xSat,
    454                                  const ColumnVector& vSat, ColumnVector& dx) {
     743                                 const ColumnVector& vSat, ColumnVector& dx,
     744                                 e_attMode mode, double externalYaw) {
    455745
    456746  t_frequency::type frqType = t_frequency::dummy;
    … …  
    495785      ColumnVector sy, sx;
    496786
    497       // GLONASS: override the nominal Sun-pointing attitude near the
    498       // orbit noon/midnight points when the beta angle is small (see
    499       // glonassYawAngle() above). Elsewhere GLONASS-M follows the same
    500       // nominal law as the other constellations, so the result is
    501       // identical to the direct Sun-pointing computation used below.
    502       // -----------------------------------------------------------------
    503       if (prn[0] == 'R' && vSat.size() == 3) {
    504         double psi = glonassYawAngle(prn, Mjd, xSat, vSat, xSun);
    505 
    506         ColumnVector vInert = vSat;
     787      // Determine the satellite body frame orientation.
     788      //
     789      // ATT_NOMINAL: simple nominal Sun-pointing for all systems.
     790      // ATT_COMPUTED or ATT_EXTERNAL: use per-system attitude models.
     791      //   GLONASS  → yaw-fixed model (Dilssner et al. 2011)
     792      //   GPS      → noon/midnight turn model (Kouba 2009/2015)
     793      //   others   → simple nominal Sun-pointing
     794      // ATT_EXTERNAL: caller supplies the yaw angle [rad] directly
     795      //   (velocity-referenced frame, same convention as GLONASS/GPS models).
     796      // -----------------------------------------------------------------------
     797      bool useVelocityFrame = false;
     798      double psiEff = 0.0;
     799
     800      if (mode == ATT_COMPUTED && prn[0] == 'R' && vSat.size() == 3) {
     801        psiEff = glonassYawAngle(prn, Mjd, xSat, vSat, xSun);
     802        useVelocityFrame = true;
     803      }
     804      else if (mode == ATT_COMPUTED && prn[0] == 'G' && vSat.size() == 3
     805               && !map->blockType.isEmpty()) {
     806        psiEff = gpsYawAngle(prn, map->blockType, Mjd, xSat, vSat, xSun);
     807        useVelocityFrame = true;
     808      }
     809      else if (mode == ATT_COMPUTED && prn[0] == 'E' && vSat.size() == 3
     810               && !map->blockType.isEmpty()) {
     811        psiEff = galileoYawAngle(prn, map->blockType, Mjd, xSat, vSat, xSun);
     812        useVelocityFrame = true;
     813      }
     814      else if (mode == ATT_COMPUTED && prn[0] == 'C' && vSat.size() == 3
     815               && !map->blockType.isEmpty()) {
     816        psiEff = bdsYawAngle(prn, map->blockType, Mjd, xSat, vSat, xSun);
     817        useVelocityFrame = true;
     818      }
     819      else if (mode == ATT_EXTERNAL && vSat.size() == 3) {
     820        psiEff = externalYaw;
     821        useVelocityFrame = true;
     822      }
     823
     824      if (useVelocityFrame) {
    507825        ColumnVector Omega(3); Omega(1) = 0.0; Omega(2) = 0.0; Omega(3) = t_CST::omega;
    508         vInert += crossproduct(Omega, xSat);
     826        ColumnVector vInert = vSat + crossproduct(Omega, xSat);
    509827
    510828        ColumnVector sy0 = crossproduct(sz, vInert);
    … …  
    512830        ColumnVector sx0 = crossproduct(sy0, sz);
    513831
    514         // Rodrigues rotation of (sx0, sy0) around sz by angle psi
    515         double cosY = cos(psi);
    516         double sinY = sin(psi);
     832        // Rodrigues rotation of (sx0, sy0) around sz by psiEff
     833        double cosY = cos(psiEff);
     834        double sinY = sin(psiEff);
    517835        sx = sx0 * cosY + crossproduct(sz, sx0) * sinY;
    518836        sy = sy0 * cosY + crossproduct(sz, sy0) * sinY;
    519837      }
    520838      else {
     839        // Nominal Sun-pointing: direct formula (ATT_NOMINAL, or computed
     840        // with no block-type-specific model available)
    521841        sy = crossproduct(sz, xSun);
    522842        sy /= sqrt(DotProduct(sy,sy));

  • trunk/BNC/src/bncantex.h

    r10951 r10955  
    4646  double  rcvCorr(const std::string& antName, t_frequency::type frqType,
    4747                  double eleSat, double azSat, bool& found) const;
     48
     49  // Attitude model selection for satCoMcorrection().
     50  enum e_attMode {
     51    ATT_COMPUTED = 0,  // full model: GLONASS yaw-fixed + GPS noon/midnight turn
     52    ATT_NOMINAL  = 1,  // simple nominal Sun-pointing only (no maneuver model)
     53    ATT_EXTERNAL = 2,  // use caller-supplied externalYaw angle [rad]
     54  };
     55
    4856  t_irc   satCoMcorrection(const QString& prn, double Mjd,
    4957                           const ColumnVector& xSat, const ColumnVector& vSat,
    50                            ColumnVector& dx);
     58                           ColumnVector& dx,
     59                           e_attMode mode = ATT_COMPUTED,
     60                           double externalYaw = 0.0);
    5161
    52   // Drains and returns the diagnostic log accumulated by the GLONASS
    53   // yaw-fixed model (mode transitions, with beta/mu/rate at the time).
     62  // Drain and return diagnostic logs accumulated by the yaw models.
    5463  // Empty most of the time; only non-empty right after a transition.
    5564  QString takeGlonassYawLog() {
    … …  
    5867    return s;
    5968  }
     69  QString takeGpsYawLog() {
     70    QString s = _gpsYawLog;
     71    _gpsYawLog.clear();
     72    return s;
     73  }
     74  QString takeOnYawLog() {
     75    QString s = _onYawLog;
     76    _onYawLog.clear();
     77    return s;
     78  }
    6079
    6180 private:
    62   // Per-satellite GLONASS yaw state, used to freeze the yaw angle while
    63   // the satellite cannot follow the nominal Sun-pointing attitude law
    64   // (see glonassYawAngle() below).
     81  // Per-satellite GLONASS yaw state: freezes the yaw angle while the
     82  // satellite cannot follow the nominal Sun-pointing law (Dilssner 2011).
    6583  class t_glonassYaw {
    6684   public:
    … …  
    7795  };
    7896
     97  // Per-satellite GPS yaw state: rate-limits the yaw during noon/midnight
     98  // turns when the required yaw rate exceeds the block's mechanical maximum
     99  // (Kouba 2009/2015, Bar-Sever 1996).
     100  class t_gpsYaw {
     101   public:
     102    t_gpsYaw() {
     103      yaw         = 0.0;
     104      lastCallMjd = 0.0;
     105      valid       = false;
     106      inTurn      = false;
     107    }
     108    double yaw;
     109    double lastCallMjd;
     110    bool   valid;
     111    bool   inTurn;  // true while in a constrained noon/midnight turn
     112  };
     113
     114  // Per-satellite Galileo/BDS yaw state: tracks rate-limited rotation between
     115  // yaw-steering (psiNom) and orbit-normal mode (psi=0) when |beta| crosses
     116  // the constellation-specific threshold (Kouba 2017, Dai et al. 2015).
     117  class t_onYaw {
     118   public:
     119    t_onYaw() {
     120      yaw         = 0.0;
     121      lastCallMjd = 0.0;
     122      valid       = false;
     123      inON        = false;
     124    }
     125    double yaw;
     126    double lastCallMjd;
     127    bool   valid;
     128    bool   inON;  // true while satellite is in (or transitioning to) orbit-normal mode
     129  };
     130
    79131  QString _glonassYawLog;
     132  QString _gpsYawLog;
     133  QString _onYawLog;
    80134
    81135  double glonassYawAngle(const QString& prn, double Mjd, const ColumnVector& xSat,
    82136                          const ColumnVector& vSat, const ColumnVector& xSun);
    83137
     138  double gpsYawAngle(const QString& prn, const QString& blockType, double Mjd,
     139                     const ColumnVector& xSat, const ColumnVector& vSat,
     140                     const ColumnVector& xSun);
     141
     142  double onModeYawAngle(const QString& prn, double betaThr, double psiDotMax,
     143                        double Mjd, const ColumnVector& xSat,
     144                        const ColumnVector& vSat, const ColumnVector& xSun);
     145
     146  double galileoYawAngle(const QString& prn, const QString& blockType, double Mjd,
     147                         const ColumnVector& xSat, const ColumnVector& vSat,
     148                         const ColumnVector& xSun);
     149
     150  double bdsYawAngle(const QString& prn, const QString& blockType, double Mjd,
     151                     const ColumnVector& xSat, const ColumnVector& vSat,
     152                     const ColumnVector& xSun);
     153
    84154  QMap<QString, t_glonassYaw> _glonassYaw;
     155  QMap<QString, t_gpsYaw>     _gpsYaw;
     156  QMap<QString, t_onYaw>      _onYaw;
    85157
    86158  class t_frqMap {
    … …  
    111183    }
    112184    QString                            antName;
     185    QString                            blockType;  // e.g. "IIF", "IIR-M", "IIA", "IIIA" (GPS only)
    113186    double                             zen1;
    114187    double                             zen2;

  • trunk/BNC/src/bnchelp.html

    r10946 r10955  
    50965096Default is an empty option field, meaning that no satellite is excluded from this individual AC.</p>
    50975097<p>
     5098Use the 'Attitude' field to select how satellite attitude is modelled when converting Antenna Phase Center (APC) corrections to
     5099Center-of-Mass (CoM) positions required for SP3 output. Three options are available:
     5100<ul>
     5101<li><b>Computed</b> (default): BNC applies its own kinematic attitude model:
     5102GPS noon/midnight turn manoeuvres (Kouba 2009/2015, Bar-Sever 1996),
     5103GLONASS yaw-fixed mode (Dilssner et al. 2011), and
     5104Galileo / BDS orbit-normal mode switching (Kouba 2017, Dai et al. 2015, Steigenberger et al. 2018).</li>
     5105<li><b>Nominal</b>: a simplified, continuous Sun-pointing model is used without any manoeuvre modelling.</li>
     5106<li><b>SSR</b>: the yaw angle transmitted in the SSR phase bias message is used directly, if present for the satellite and epoch.
     5107If no yaw angle is available for a particular satellite in a given epoch, BNC falls back to 'Computed'.
     5108Select this option only if you trust the yaw values provided by the Analysis Center.</li>
     5109</ul>
     5110Note that the attitude model affects APC-referenced correction streams only.
     5111For CoM-referenced streams (SSRC) the Analysis Center has already applied its own attitude model before encoding.</p>
     5112<p>
    50985113Note that the orbit information in the resulting combination stream is just copied from one of the incoming streams.
    50995114The stream used for providing the orbits may vary over time: if the orbit providing stream has an outage

  • trunk/BNC/src/bncwindow.cpp

    r10945 r10955  
    472472  // Combine Corrections
    473473  // -------------------
    474   _cmbTable = new QTableWidget(0, 4);
    475   _cmbTable->setHorizontalHeaderLabels(QString("Mountpoint, AC Name, Weight Factor, Exclude Satellites").split(","));
     474  _cmbTable = new QTableWidget(0, 5);
     475  _cmbTable->setHorizontalHeaderLabels(QString("Mountpoint, AC Name, Weight Factor, Exclude Satellites, Attitude").split(","));
    476476  _cmbTable->setSelectionMode(QAbstractItemView::ExtendedSelection);
    477477  _cmbTable->setSelectionBehavior(QAbstractItemView::SelectRows);
    478   _cmbTable->setMaximumWidth(40 * ww);
     478  _cmbTable->setMaximumWidth(50 * ww);
    479479  _cmbTable->horizontalHeader()->resizeSection(0, 10 * ww);
    480480  _cmbTable->horizontalHeader()->resizeSection(1, 6 * ww);
    481481  _cmbTable->horizontalHeader()->resizeSection(2, 9 * ww);
    482482  _cmbTable->horizontalHeader()->resizeSection(3, 9 * ww);
     483  _cmbTable->horizontalHeader()->resizeSection(4, 8 * ww);
    483484#if QT_VERSION < 0x050000
    484485  _cmbTable->horizontalHeader()->setResizeMode(QHeaderView::Interactive);
    … …  
    17011702  // WhatsThis, Combine Corrections
    17021703  // ------------------------------
    1703   _cmbTable->setWhatsThis(tr("<p>BNC allows to process several orbit and clock correction streams in real-time to produce, encode, upload and save a combination of correctors coming from different providers. </p><p>To add a line to the 'Combine Corrections' table hit the 'Add Row' button, double click on the 'Mountpoint' field to specify a Broadcast Ephemeris Correction mountpoint from the 'Streams' section below and hit Enter. Then double click on the 'AC Name' field to enter your choice of an abbreviation for the Analysis Center (AC) providing the stream. Double click on the 'Weight Factor' field to enter a weight factor to be applied for this stream in the combination. A Factor greater than 1 will enlarge the sigma of the clock pseudo-observations and with it down-weight its contribution. Finally, double click on the 'Exclude Satellites' field and specify satellites, to exclude them for an individual AC. An entry 'G04,G31,R' means to excludes GPS satellites PRN 4 and 31 as well as all GLONASS satellites from one individual AC. Default is an empty option field, meaning that no satellite is excluded from this individual AC.</p><p>Note that the orbit information in the resulting combination stream is just copied from one of the incoming streams. The stream used for providing the orbits may vary over time: if the orbit providing stream has an outage then BNC switches to the next remaining stream for getting hold of the orbit information.</p><p>The combination process requires Broadcast Ephemeris. Besides orbit and clock correction streams BNC should therefore pull a stream carrying Broadcast Ephemeris in the form of RTCM Version 3 messages.</p><p>It is possible to specify only one Broadcast Ephemeris Correction stream in the 'Combine Corrections' table. Instead of combining corrections BNC will then add the corrections to the Broadcast Ephemeris with the possibility to save final orbit and clock results in SP3 and/or Clock RINEX format. <i>[key: cmbStreams]</i></p>"));
     1704  _cmbTable->setWhatsThis(tr("<p>BNC allows to process several orbit and clock correction streams in real-time to produce, encode, upload and save a combination of correctors coming from different providers. </p><p>To add a line to the 'Combine Corrections' table hit the 'Add Row' button, double click on the 'Mountpoint' field to specify a Broadcast Ephemeris Correction mountpoint from the 'Streams' section below and hit Enter. Then double click on the 'AC Name' field to enter your choice of an abbreviation for the Analysis Center (AC) providing the stream. Double click on the 'Weight Factor' field to enter a weight factor to be applied for this stream in the combination. A Factor greater than 1 will enlarge the sigma of the clock pseudo-observations and with it down-weight its contribution. Finally, double click on the 'Exclude Satellites' field and specify satellites, to exclude them for an individual AC. An entry 'G04,G31,R' means to excludes GPS satellites PRN 4 and 31 as well as all GLONASS satellites from one individual AC. Default is an empty option field, meaning that no satellite is excluded from this individual AC.</p><p>Use the 'Attitude' field to select how satellite attitude is modelled when converting Antenna Phase Center (APC) corrections to Center-of-Mass (CoM) positions for SP3 output. 'Computed' (default) applies BNC's kinematic attitude model, including GPS noon/midnight turn manoeuvres, GLONASS yaw-fixed mode, and Galileo/BDS orbit-normal mode switching. 'Nominal' uses a simplified continuous Sun-pointing model without manoeuvre modelling. 'SSR' uses the yaw angle transmitted in the SSR phase bias message; select this option only if you trust the yaw values provided by the AC. If no yaw angle is available for a particular satellite in a given epoch, BNC falls back to 'Computed'. Note that the attitude model affects APC-referenced streams only; for CoM-referenced streams the AC has already applied attitude before encoding.</p><p>Note that the orbit information in the resulting combination stream is just copied from one of the incoming streams. The stream used for providing the orbits may vary over time: if the orbit providing stream has an outage then BNC switches to the next remaining stream for getting hold of the orbit information.</p><p>The combination process requires Broadcast Ephemeris. Besides orbit and clock correction streams BNC should therefore pull a stream carrying Broadcast Ephemeris in the form of RTCM Version 3 messages.</p><p>It is possible to specify only one Broadcast Ephemeris Correction stream in the 'Combine Corrections' table. Instead of combining corrections BNC will then add the corrections to the Broadcast Ephemeris with the possibility to save final orbit and clock results in SP3 and/or Clock RINEX format. <i>[key: cmbStreams]</i></p>"));
    17041705  addCmbRowButton->setWhatsThis(tr("<p>Hit 'Add Row' button to add another line to the 'Combine Corrections' table.</p>"));
    17051706  delCmbRowButton->setWhatsThis(tr("<p>Hit 'Delete' button to delete the highlighted line(s) from the 'Combine Corrections' table.</p>"));
    … …  
    22182219    QString hlp;
    22192220    for (int iCol = 0; iCol < _cmbTable->columnCount(); iCol++) {
    2220       if (_cmbTable->item(iRow, iCol)) {
     2221      if (iCol == 4) {
     2222        QComboBox* attCombo = qobject_cast<QComboBox*>(_cmbTable->cellWidget(iRow, iCol));
     2223        hlp += (attCombo ? attCombo->currentText() : "Computed") + " ";
     2224      }
     2225      else if (_cmbTable->item(iRow, iCol)) {
    22212226        hlp += _cmbTable->item(iRow, iCol)->text() + " ";
    22222227      }
    … …  
    30563061  _cmbTable->insertRow(iRow);
    30573062  for (int iCol = 0; iCol < _cmbTable->columnCount(); iCol++) {
    3058     _cmbTable->setItem(iRow, iCol, new QTableWidgetItem(""));
     3063    if (iCol == 4) {
     3064      QComboBox* attCombo = new QComboBox();
     3065      attCombo->setEditable(false);
     3066      attCombo->addItems(QString("Computed,SSR,Nominal").split(","));
     3067      attCombo->setFrame(false);
     3068      _cmbTable->setCellWidget(iRow, iCol, attCombo);
     3069    }
     3070    else {
     3071      _cmbTable->setItem(iRow, iCol, new QTableWidgetItem(""));
     3072    }
    30593073  }
    30603074}
    … …  
    31103124      _cmbTable->insertRow(iRow);
    31113125    }
    3112     for (int iCol = 0; iCol < hlp.size(); iCol++) {
    3113       _cmbTable->setItem(iRow, iCol, new QTableWidgetItem(hlp[iCol]));
     3126    for (int iCol = 0; iCol < hlp.size() && iCol < _cmbTable->columnCount(); iCol++) {
     3127      if (iCol == 4) {
     3128        QComboBox* attCombo = new QComboBox();
     3129        attCombo->setEditable(false);
     3130        attCombo->addItems(QString("Computed,SSR,Nominal").split(","));
     3131        attCombo->setFrame(false);
     3132        int idx = attCombo->findText(hlp[iCol]);
     3133        if (idx != -1) attCombo->setCurrentIndex(idx);
     3134        _cmbTable->setCellWidget(iRow, iCol, attCombo);
     3135      }
     3136      else {
     3137        _cmbTable->setItem(iRow, iCol, new QTableWidgetItem(hlp[iCol]));
     3138      }
     3139    }
     3140    // Ensure the Attitude combo exists even for old configs without column 4
     3141    if (hlp.size() <= 4 && _cmbTable->columnCount() > 4 && iRow >= 0
     3142        && !_cmbTable->cellWidget(iRow, 4)) {
     3143      QComboBox* attCombo = new QComboBox();
     3144      attCombo->setEditable(false);
     3145      attCombo->addItems(QString("Computed,SSR,Nominal").split(","));
     3146      attCombo->setFrame(false);
     3147      _cmbTable->setCellWidget(iRow, 4, attCombo);
    31143148    }
    31153149  }

  • trunk/BNC/src/combination/bnccomb.cpp

    r10954 r10955  
    181181      QStringList hlp = it.next().split(" ");
    182182      cmbAC* newAC = new cmbAC();
    183       newAC->mountPoint   = hlp[0];
    184       newAC->name         = hlp[1];
    185       newAC->weightFactor = hlp[2].toDouble();
    186       newAC->excludeSats  = hlp[3].split(QRegExp("[ ,]"), Qt::SkipEmptyParts);
    187       newAC->isAPC        = bool(newAC->mountPoint.mid(0,4) == "SSRA");
     183      newAC->mountPoint     = hlp[0];
     184      newAC->name           = hlp[1];
     185      newAC->weightFactor   = hlp[2].toDouble();
     186      newAC->excludeSats    = hlp[3].split(QRegExp("[ ,]"), Qt::SkipEmptyParts);
     187      newAC->isAPC          = bool(newAC->mountPoint.mid(0,4) == "SSRA");
     188      newAC->attitudeSource = (hlp.size() > 4) ? hlp[4] : "Computed";
    188189      QMapIterator<char, unsigned> itSys(_cmbSysPrn);
    189190      // init
    … …  
    859860          dt -= (0.5 * ssrUpdateInt[pb._updateInt]);
    860861        }
    861         _newCorr->_satYawAngle = pb._yaw + pb._yawRate * dt;
     862        _newCorr->_satYawAngle      = pb._yaw + pb._yawRate * dt;
     863        _newCorr->_satYawAngleValid = true;
    862864
    863865        // _lambdaIF
    … …  
    12131215      if (corr->_prn.startsWith("R08")) {
    12141216        emit newMessage(("bncComb: DIAG " + corr->_prn.mid(0,3) + " getCrd FAILED").toLatin1(), false);
    1215       }
    1216       */
     1217      } */
    12171218      delete corr;
    12181219      it.remove();
    … …  
    12221223    else if (corr->_prn.startsWith("R08")) {
    12231224      emit newMessage(("bncComb: DIAG " + corr->_prn.mid(0,3) + " getCrd OK").toLatin1(), false);
    1224     }
    1225    
    1226     */
     1225    }*/
    12271226    // TEMPORARY DIAGNOSTIC: how far is this epoch from the GLONASS broadcast
    12281227    // ephemeris reference time? t_ephGlo::position() numerically integrates
    … …  
    12511250      char sys = corr->_eph->prn().system();
    12521251      masterIsAPC = _masterIsAPC[sys];
    1253       if (_antex->satCoMcorrection(corr->_prn, Mjd, xc.Rows(1,3), vv, dx) != success) {
     1252
     1253      // Determine attitude mode from the master orbit AC for this system
     1254      bncAntex::e_attMode attMode   = bncAntex::ATT_COMPUTED;
     1255      double              extYaw    = 0.0;
     1256      const QString&      masterName = _masterOrbitAC[sys];
     1257      for (const cmbAC* ac : _ACs) {
     1258        if (ac->name == masterName) {
     1259          if (ac->attitudeSource == "SSR" && corr->_satYawAngleValid) {
     1260            attMode = bncAntex::ATT_EXTERNAL;
     1261            extYaw  = corr->_satYawAngle;
     1262          }
     1263          else if (ac->attitudeSource == "Nominal") {
     1264            attMode = bncAntex::ATT_NOMINAL;
     1265          }
     1266          break;
     1267        }
     1268      }
     1269
     1270      if (_antex->satCoMcorrection(corr->_prn, Mjd, xc.Rows(1,3), vv, dx,
     1271                                   attMode, extYaw) != success) {
    12541272        dx = 0;
    12551273        emit newMessage("bncComb: antenna not found " + corr->_prn.mid(0,3).toLatin1(), false);
    … …  
    12581276      if (!glonassYawLog.isEmpty()) {
    12591277        emit newMessage(("bncComb: " + glonassYawLog.trimmed()).toLatin1(), false);
     1278      }
     1279      QString gpsYawLog = _antex->takeGpsYawLog();
     1280      if (!gpsYawLog.isEmpty()) {
     1281        emit newMessage(("bncComb: " + gpsYawLog.trimmed()).toLatin1(), false);
     1282      }
     1283      QString onYawLog = _antex->takeOnYawLog();
     1284      if (!onYawLog.isEmpty()) {
     1285        emit newMessage(("bncComb: " + onYawLog.trimmed()).toLatin1(), false);
    12601286      }
    12611287    }

  • trunk/BNC/src/combination/bnccomb.h

    r10914 r10955  
    9595      numObs['S']  = 0;
    9696      numObs['I']  = 0;
    97       isAPC = false;
     97      isAPC          = false;
     98      attitudeSource = "Computed";
    9899    }
    99100    ~cmbAC() {
    … …  
    105106    QStringList          excludeSats;
    106107    bool                 isAPC;
     108    QString              attitudeSource; // "Computed", "SSR", or "Nominal"
    107109    QMap<char, unsigned> numObs;
    108110  };
    … …  
    117119      _lambdaIF                    = 0.0;
    118120      _satYawAngle                 = 0.0;
     121      _satYawAngleValid            = false;
    119122      _weightFactor                = 1.0;
    120123      _satPos.ReSize(3); _satPos   = 0.0;
    … …  
    135138    double         _lambdaIF;
    136139    double         _satYawAngle;
     140    bool           _satYawAngleValid;
    137141    double         _dClkResult;
    138142    ColumnVector   _satPos;
Note: See TracChangeset for help on using the changeset viewer.