Index: trunk/BNC/src/RTCM3/RTCM3coDecoder.cpp
===================================================================
--- trunk/BNC/src/RTCM3/RTCM3coDecoder.cpp (revision 11081)
+++ trunk/BNC/src/RTCM3/RTCM3coDecoder.cpp (revision 11083)
@@ -1275,5 +1275,5 @@
if (bdt) {
epoSecBds += 14.0;
- if (epoSecBds > 604800.0) {
+ if (epoSecBds >= 604800.0) {
epoSecBds -= 7.0*24.0*60.0*60.0;
}
@@ -1289,11 +1289,17 @@
}
+ // The epoch is given as seconds of week (or of day): move it by whole days
+ // to the one closest to the current time. Folding in 12 h steps (as
+ // before) stopped at the edge of the +-12 h window, e.g. at the week
+ // change a correction of Saturday 23:59:55 received after midnight ended
+ // up at Sunday 11:59:55.
if (_lastTime.valid()) {
- double maxDiff = 12 * 3600.0;
- while (_lastTime < currentTime - maxDiff) {
- _lastTime = _lastTime + maxDiff;
- }
- while (_lastTime > currentTime + maxDiff) {
- _lastTime = _lastTime - maxDiff;
+ const double halfDay = 12 * 3600.0;
+ const double day = 24 * 3600.0;
+ while (_lastTime < currentTime - halfDay) {
+ _lastTime = _lastTime + day;
+ }
+ while (_lastTime > currentTime + halfDay) {
+ _lastTime = _lastTime - day;
}
}
Index: trunk/BNC/src/RTCM3/clock_and_orbit/clock_orbit_rtcm.cpp
===================================================================
--- trunk/BNC/src/RTCM3/clock_and_orbit/clock_orbit_rtcm.cpp (revision 11081)
+++ trunk/BNC/src/RTCM3/clock_and_orbit/clock_orbit_rtcm.cpp (revision 11083)
@@ -97,5 +97,5 @@
break;
case CLOCKORBIT_SATGLONASS:
- E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
+ E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
E_GLONASS_IOD(co->Sat[i].IOD)
break;
@@ -161,5 +161,5 @@
break;
case CLOCKORBIT_SATGLONASS:
- E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
+ E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
break;
}
@@ -221,5 +221,5 @@
break;
case CLOCKORBIT_SATGLONASS:
- E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
+ E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
E_GLONASS_IOD(co->Sat[i].IOD)
break;
@@ -294,5 +294,5 @@
break;
case CLOCKORBIT_SATGLONASS:
- E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
+ E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
break;
}
@@ -334,5 +334,5 @@
break;
case CLOCKORBIT_SATGLONASS:
- E_GLONASS_SATELLITE_ID(co->Sat[i].ID)
+ E_GLONASS_SATELLITE_ID_OLD(co->Sat[i].ID)
break;
}
@@ -385,5 +385,5 @@
break;
case CLOCKORBIT_SATGLONASS:
- E_GLONASS_SATELLITE_ID(cb->Sat[i].ID)
+ E_GLONASS_SATELLITE_ID_OLD(cb->Sat[i].ID)
break;
}
@@ -442,5 +442,5 @@
break;
case CLOCKORBIT_SATGLONASS:
- E_GLONASS_SATELLITE_ID(pb->Sat[i].ID)
+ E_GLONASS_SATELLITE_ID_OLD(pb->Sat[i].ID)
break;
}
Index: trunk/BNC/src/bnchelp.html
===================================================================
--- trunk/BNC/src/bnchelp.html (revision 11081)
+++ trunk/BNC/src/bnchelp.html (revision 11083)
@@ -4337,10 +4337,21 @@
- To compare satellite clocks provided by the two files, BNC first converts coordinate differences dX,dY,dZ
- into along track, out-of-plane, and radial components. It then corrects the clock differences for the radial
- components
- of coordinate differences. RMS values of clock differences are finally calculated after introducing at first one
- offset
- 'per epoch for all satellites' and secondly one offset 'per satellite for all epochs'.
+ To compare the orbits, BNC converts the coordinate differences dX,dY,dZ into radial, along-track and out-of-plane
+ components, using the satellite velocity derived from neighbouring epochs.
+
+
+ The clocks of two products usually differ by an offset common to all satellites of a system, which changes from
+ epoch to epoch (clock datum), and by a constant offset per satellite (e.g. different realizations of the satellite
+ code biases). For each system, BNC therefore estimates from all clock differences in a least-squares adjustment
+ one offset 'per epoch for all satellites' and one offset 'per satellite for all epochs' (column 'Offset'), with
+ the sum of the satellite offsets constrained to zero. The clock residuals of this adjustment ('clk') are then
+ corrected for the radial component of the orbit differences ('clkRed'). The signal-in-space difference ('SIS',
+ radial orbit difference minus clock difference) removes only the epoch offset, i.e. it contains the satellite
+ offsets.
+
+
+ The epoch offsets are eliminated from the normal equations before the satellite offsets are solved for; the
+ result is identical to that of the complete adjustment. This keeps the comparison fast also for long or densely
+ sampled files: comparing two daily files with 30-second sampling (2880 epochs) takes a few seconds.
Index: trunk/BNC/src/combination/bnccomb.cpp
===================================================================
--- trunk/BNC/src/combination/bnccomb.cpp (revision 11081)
+++ trunk/BNC/src/combination/bnccomb.cpp (revision 11083)
@@ -1525,9 +1525,11 @@
// 30 min and are valid for ±900 s around the reference time (TOC).
// t_ephGlo::position() uses Runge-Kutta integration from TOC; beyond
- // ±900 s the integration error grows unchecked and propagates into both
- // the combined clock estimate and the SP3 position.
+ // that the integration error grows and propagates into both the
+ // combined clock estimate and the SP3 position. A margin of 60 s is
+ // allowed, as ACs switch to the next ephemeris a few seconds before its
+ // window begins (otherwise satellites drop out at every switch).
if (corr->_eph->prn().system() == 'R') {
double age = epoTime - corr->_eph->TOC();
- if (fabs(age) > 900.0) {
+ if (fabs(age) > 960.0) {
emit newMessage(QString().asprintf(
"bncComb: skip %s, GLONASS eph age %.0fs",
Index: trunk/BNC/src/orbComp/sp3Comp.cpp
===================================================================
--- trunk/BNC/src/orbComp/sp3Comp.cpp (revision 11081)
+++ trunk/BNC/src/orbComp/sp3Comp.cpp (revision 11083)
@@ -178,42 +178,62 @@
}
- int nPar = epochs.size() + clkSats.size();
- SymmetricMatrix NN(nPar); NN = 0.0;
- ColumnVector bb(nPar); bb = 0.0;
-
- // Create Matrix A'A and vector b
- // ------------------------------
- for (unsigned ie = 0; ie < epochs.size(); ie++) {
+ // Least-squares adjustment of an epoch-specific offset (common to all
+ // satellites) and a satellite-specific offset (common to all epochs), the
+ // sum of the satellite offsets constrained to zero (weight 1). The normal
+ // equations of the epoch offsets are diagonal, so these offsets are
+ // eliminated (Schur complement) and only the satellite offsets are solved
+ // for; the epoch offsets follow from them. The solution is identical to
+ // that of the full system, but the effort grows only linearly with the
+ // number of epochs (before: cubic).
+ // -----------------------------------------------------------------------
+ int nEpo = epochs.size();
+ int nSat = clkSats.size();
+ int nPar = nEpo + nSat;
+
+ vector > satIdx(nEpo); // satellites (index) of each epoch
+ ColumnVector nEp(nEpo); nEp = 0.0; // number of satellites per epoch
+ ColumnVector bEp(nEpo); bEp = 0.0; // sum of the differences per epoch
+ SymmetricMatrix NS(nSat); NS = 1.0; // regularization: sum of satellite offsets
+ ColumnVector bS(nSat); bS = 0.0;
+ for (int ie = 0; ie < nEpo; ie++) {
const map& dc = epochs[ie]->_dc;
- Matrix AA(dc.size(), nPar); AA = 0.0;
- ColumnVector ll(dc.size()); ll = 0.0;
- map::const_iterator it;
- int ii = -1;
- for (it = dc.begin(); it != dc.end(); it++) {
- const t_prn& prn = it->first;
- if (satIndex(clkSats, prn) != -1) {
- ++ii;
- int index = epochs.size() + satIndex(clkSats, prn);
- AA[ii][ie] = 1.0; // epoch-specfic offset (common for all satellites)
- AA[ii][index] = 1.0; // satellite-specific offset (common for all epochs)
- ll[ii] = it->second;
- }
- }
- SymmetricMatrix dN; dN << AA.t() * AA;
- NN += dN;
- bb += AA.t() * ll;
- }
-
- // Regularize NN
- // -------------
- RowVector HH(nPar);
- HH.columns(1, epochs.size()) = 0.0;
- HH.columns(epochs.size()+1, nPar) = 1.0;
- SymmetricMatrix dN; dN << HH.t() * HH;
- NN += dN;
-
- // Estimate Parameters
- // -------------------
- ColumnVector xx = NN.i() * bb;
+ for (map::const_iterator it = dc.begin(); it != dc.end(); it++) {
+ int is = satIndex(clkSats, it->first);
+ if (is != -1) {
+ satIdx[ie].push_back(is);
+ nEp[ie] += 1.0;
+ bEp[ie] += it->second;
+ NS(is+1, is+1) += 1.0;
+ bS[is] += it->second;
+ }
+ }
+ }
+
+ // Eliminate the epoch offsets
+ // ---------------------------
+ for (int ie = 0; ie < nEpo; ie++) {
+ const vector& idx = satIdx[ie];
+ for (unsigned ia = 0; ia < idx.size(); ia++) {
+ bS[idx[ia]] -= bEp[ie] / nEp[ie];
+ for (unsigned ib = 0; ib <= ia; ib++) {
+ NS(idx[ia]+1, idx[ib]+1) -= 1.0 / nEp[ie];
+ }
+ }
+ }
+
+ // Estimate Parameters (epoch offsets first, then satellite offsets)
+ // ---------------------------------------------------------------
+ ColumnVector yy = NS.i() * bS;
+ ColumnVector xx(nPar);
+ for (int ie = 0; ie < nEpo; ie++) {
+ double sum = bEp[ie];
+ for (unsigned ia = 0; ia < satIdx[ie].size(); ia++) {
+ sum -= yy[satIdx[ie][ia]];
+ }
+ xx[ie] = sum / nEp[ie];
+ }
+ for (int is = 0; is < nSat; is++) {
+ xx[nEpo + is] = yy[is];
+ }
// Compute clock residuals