Index: trunk/BNC/src/rinex/reqcanalyze.cpp =================================================================== --- trunk/BNC/src/rinex/reqcanalyze.cpp (revision 11037) +++ trunk/BNC/src/rinex/reqcanalyze.cpp (revision 11057) @@ -131,4 +131,5 @@ // ---------------- t_reqcEdit::readEphemerides(_navFileNames, _ephs, _logStream, _checkEph); + t_reqcEdit::indexEphemerides(_ephs, _ephsByPrn); // Loop over all RINEX Files @@ -207,4 +208,9 @@ continue; } + // observations of unhealthy satellites are excluded, as those + // below the elevation mask (see also setExpectedObs) + if (_navFileNames.size() && unhealthyAt(rnxSat.prn, _currEpo->tt)) { + continue; + } t_satObs satObs; t_rnxObsFile::setObsFromRnx(obsFile, _currEpo, rnxSat, satObs); @@ -266,11 +272,8 @@ } - t_eph* eph = 0; - for (int ie = 0; ie < _ephs.size(); ie++) { - if (_ephs[ie]->prn() == prn) { - eph = _ephs[ie]; - break; - } - } + if (unhealthyAt(prn, _currEpo->tt)) { + continue; + } + t_eph* eph = ephForTime(prn, _currEpo->tt); if (eph) { ColumnVector xSat(6); @@ -328,11 +331,5 @@ t_qcSat& qcSat) { - t_eph* eph = 0; - for (int ie = 0; ie < _ephs.size(); ie++) { - if (_ephs[ie]->prn() == satObs._prn) { - eph = _ephs[ie]; - break; - } - } + t_eph* eph = ephForTime(satObs._prn, epoTime); if (eph) { ColumnVector xc(6); @@ -976,15 +973,26 @@ double interval, const ColumnVector& xyzSta) { for(QMap::iterator it = _numExpObs.begin(); it != _numExpObs.end(); it++) { - t_eph* eph = 0; + bool navAvail = false; for (int ie = 0; ie < _ephs.size(); ie++) { if (_ephs[ie]->prn() == it.key()) { - eph = _ephs[ie]; + navAvail = true; break; } } - if (eph) { + if (navAvail) { int numExpObs = 0; bncTime epoTime; for (epoTime = startTime - interval; epoTime < endTime; epoTime = epoTime + interval) { + // an observation is only expected if a valid ephemeris exists, + // the satellite is healthy (unhealthy satellite-epochs are also + // excluded from the observations, see analyzeFile) and above + // the elevation mask + if (unhealthyAt(it.key(), epoTime)) { + continue; + } + t_eph* eph = ephForTime(it.key(), epoTime, true); + if (!eph) { + continue; + } ColumnVector xc(6); ColumnVector vv(3); @@ -1007,2 +1015,18 @@ } } + +// A satellite is considered unhealthy at epoch tt if the valid navigation +// data set closest to tt is flagged unhealthy. Its orbit is not used then +// (the geometry may be erroneous) and its observations are excluded. +//////////////////////////////////////////////////////////////////////////// +bool t_reqcAnalyze::unhealthyAt(const t_prn& prn, const bncTime& tt) const { + t_eph* eph = ephForTime(prn, tt, true, true); + return eph && eph->checkState() == t_eph::unhealthy; +} + +// Ephemeris of the given satellite for epoch tt, see t_reqcEdit::ephForTime +//////////////////////////////////////////////////////////////////////////// +t_eph* t_reqcAnalyze::ephForTime(const t_prn& prn, const bncTime& tt, + bool validOnly, bool allowUnhealthy) const { + return t_reqcEdit::ephForTime(_ephsByPrn.value(prn), prn, tt, validOnly, allowUnhealthy); +} Index: trunk/BNC/src/rinex/reqcanalyze.h =================================================================== --- trunk/BNC/src/rinex/reqcanalyze.h (revision 11037) +++ trunk/BNC/src/rinex/reqcanalyze.h (revision 11057) @@ -194,4 +194,9 @@ void setExpectedObs(const bncTime& startTime, const bncTime& endTime, double interval, const ColumnVector& xyzSta); + + t_eph* ephForTime(const t_prn& prn, const bncTime& tt, bool validOnly = false, + bool allowUnhealthy = false) const; + + bool unhealthyAt(const t_prn& prn, const bncTime& tt) const; void analyzeMultipath(); @@ -219,4 +224,5 @@ QStringList _defaultSignalTypes; QVector _ephs; + QMap > _ephsByPrn; t_rnxObsFile::t_rnxEpo* _currEpo; t_qcFile _qcFile; Index: trunk/BNC/src/rinex/reqcedit.cpp =================================================================== --- trunk/BNC/src/rinex/reqcedit.cpp (revision 11037) +++ trunk/BNC/src/rinex/reqcedit.cpp (revision 11057) @@ -230,4 +230,5 @@ if (_minEle > 0.0 && _ephs.isEmpty()) { t_reqcEdit::readEphemerides(_navFileNames, _ephs, _logStream, _checkEph); + t_reqcEdit::indexEphemerides(_ephs, _ephsByPrn); } @@ -546,11 +547,5 @@ const t_rnxObsFile::t_rnxSat& rnxSat = epo->rnxSat[iSat]; - t_eph* eph = 0; - for (int ie = 0; ie < _ephs.size(); ie++) { - if (_ephs[ie]->prn() == rnxSat.prn) { - eph = _ephs[ie]; - break; - } - } + t_eph* eph = ephForTime(_ephsByPrn.value(rnxSat.prn), rnxSat.prn, epo->tt); if (eph) { ColumnVector xc(6); @@ -567,4 +562,57 @@ } epo->rnxSat = keptSats; +} + +// Ephemeris of the given satellite for epoch tt: the one with TOC closest +// to tt among those valid at tt (see outDatedBcep). If none is valid, +// return 0 when validOnly is set, otherwise the one with TOC closest to tt. +// Only ephemerides accepted by t_eph::getCrd are considered (not bad or +// outdated; unhealthy ones only if allowUnhealthy is set, to be used with +// getCrd(..., ignoreHealth = true)). +// Using the first ephemeris of the day for all epochs gives wrong +// geometry late in the day and, for GLONASS, integrations over many +// hours (t_ephGlo::position integrates from TOC on every call). +//////////////////////////////////////////////////////////////////////////// +t_eph* t_reqcEdit::ephForTime(const QVector& ephs, const t_prn& prn, + const bncTime& tt, bool validOnly, + bool allowUnhealthy) { + t_eph* bestValid = 0; + t_eph* bestAny = 0; + double bestValidDt = 0.0; + double bestAnyDt = 0.0; + for (int ie = 0; ie < ephs.size(); ie++) { + t_eph* eph = ephs[ie]; + if (eph->prn() != prn) { + continue; + } + if (eph->checkState() == t_eph::bad || + eph->checkState() == t_eph::outdated || + (eph->checkState() == t_eph::unhealthy && !allowUnhealthy)) { + continue; + } + double dt = fabs(tt - eph->TOC()); + if (bestAny == 0 || dt < bestAnyDt) { + bestAny = eph; + bestAnyDt = dt; + } + if (!outDatedBcep(eph, tt) && (bestValid == 0 || dt < bestValidDt)) { + bestValid = eph; + bestValidDt = dt; + } + } + if (bestValid || validOnly) { + return bestValid; + } + return bestAny; +} + +// Group ephemerides by satellite (speeds up ephForTime) +//////////////////////////////////////////////////////////////////////////// +void t_reqcEdit::indexEphemerides(const QVector& ephs, + QMap >& ephsByPrn) { + ephsByPrn.clear(); + for (int ie = 0; ie < ephs.size(); ie++) { + ephsByPrn[ephs[ie]->prn()].append(ephs[ie]); + } } Index: trunk/BNC/src/rinex/reqcedit.h =================================================================== --- trunk/BNC/src/rinex/reqcedit.h (revision 11037) +++ trunk/BNC/src/rinex/reqcedit.h (revision 11057) @@ -56,4 +56,9 @@ QVector& ephs, QTextStream* log, bool checkEph); + static void indexEphemerides(const QVector& ephs, + QMap >& ephsByPrn); + static t_eph* ephForTime(const QVector& ephs, const t_prn& prn, + const bncTime& tt, bool validOnly = false, + bool allowUnhealthy = false); private: @@ -82,4 +87,5 @@ QMap > _lli; QVector _ephs; + QMap > _ephsByPrn; }; Index: trunk/BNC/src/rinex/rnxobsfile.cpp =================================================================== --- trunk/BNC/src/rinex/rnxobsfile.cpp (revision 11037) +++ trunk/BNC/src/rinex/rnxobsfile.cpp (revision 11057) @@ -349,55 +349,37 @@ } else { - _obsTypes['G'] << "C1C" << "L1C" << "S1C" - << "C1W" << "L1W" << "S1W" - << "C2X" << "L2X" << "S2X" - << "C2W" << "L2W" << "S2W" - << "C5X" << "L5X" << "S5X"; - - _obsTypes['J'] << "C1C" << "L1C" << "S1C" - << "C1Z" << "L1Z" << "S1Z" - << "C1X" << "L1X" << "S1X" - << "C2L" << "L2L" << "S2L" - << "C2X" << "L2X" << "S2X" - << "C5Q" << "L5Q" << "S5Q" - << "C5X" << "L5X" << "S5X" - << "C6L" << "L6L" << "S6L"; - - _obsTypes['R'] << "C1C" << "L1C" << "S1C" - << "C1P" << "L1P" << "S1P" - << "C2C" << "L2C" << "S2C" - << "C2P" << "L2P" << "S2P" - << "C3I" << "L3I" << "S3I" - << "C4X" << "L4X" << "S4X" - << "C6X" << "L6X" << "S6X"; - - _obsTypes['E'] << "C1C" << "L1C" << "S1C" - << "C1X" << "L1X" << "S1X" - << "C5Q" << "L5Q" << "S5Q" - << "C5X" << "L5X" << "S5X" - << "C6C" << "L6C" << "S6C" - << "C6X" << "L6X" << "S6X" - << "C7Q" << "L7Q" << "S7Q" - << "C7X" << "L7X" << "S7X" - << "C8Q" << "L8Q" << "S8Q" - << "C8X" << "L8X" << "S8X"; - - _obsTypes['S'] << "C1C" << "L1C" << "S1C" - << "C5I" << "L5I" << "S5I" - << "C5Q" << "L5Q" << "S5Q" - << "C5X" << "L5X" << "S5X"; - - _obsTypes['C'] << "C2I" << "L2I" << "S2I" - << "C2Q" << "L2Q" << "S2Q" - << "C2X" << "L2X" << "S2X" - << "C6I" << "L6I" << "S6I" - << "C6Q" << "L6Q" << "S6Q" - << "C6X" << "L6X" << "S6X" - << "C7I" << "L7I" << "S7I" - << "C7Q" << "L7Q" << "S7Q" - << "C7X" << "L7X" << "S7X"; - - _obsTypes['I'] << "C5A" << "L5A" << "S5A" - << "C9A" << "L9A" << "S9A"; + // All signals the RTCM3 MSM decoder can deliver (see RTCM3Decoder.cpp), + // each with code, phase, Doppler and SNR. Empty observations are blank + // padded up to the last non-empty one of a record, so the order keeps + // records short: per system first the signals most receivers track + // (Q/C/L/P tracking modes), then the combined 'X' modes, then rarely + // used ones. + QMap sigCodes; + sigCodes['G'] << "1C" << "1W" << "2W" << "2L" << "2S" << "5Q" << "1L" + << "2X" << "5X" << "1X" + << "5I" << "1S" << "1P" << "2C" << "2P"; + sigCodes['R'] << "1C" << "1P" << "2C" << "2P" << "3Q" + << "3X" << "4X" << "6X" + << "3I" << "4A" << "4B" << "6A" << "6B"; + sigCodes['E'] << "1C" << "5Q" << "7Q" << "8Q" << "6C" + << "1X" << "5X" << "7X" << "8X" << "6X" + << "1B" << "5I" << "7I" << "8I" << "6B" << "1A" << "1Z" << "6A" << "6Z"; + sigCodes['J'] << "1C" << "2L" << "5Q" << "1L" + << "2X" << "5X" << "1X" + << "2S" << "5I" << "1S" << "1E" << "1Z" << "6L" << "6S" << "6X"; + sigCodes['C'] << "2I" << "7I" << "6I" << "5P" << "1P" << "7D" + << "2X" << "7X" << "6X" << "5X" << "1X" + << "2Q" << "7Q" << "6Q" << "5D" << "1D"; + sigCodes['S'] << "1C" << "5I" << "5Q" << "5X"; + sigCodes['I'] << "5A" << "9A"; + + QMapIterator it(sigCodes); + while (it.hasNext()) { + it.next(); + for (int ii = 0; ii < it.value().size(); ii++) { + const QString& sig = it.value()[ii]; + _obsTypes[it.key()] << "C" + sig << "L" + sig << "D" + sig << "S" + sig; + } + } } } @@ -1528,8 +1510,27 @@ char sys = rnxSat.prn.system(); + // If all observation types are given as v3 types, a direct lookup is + // equivalent to the (slow) matching loops below + // -------------------------------------------------------------------- + bool allTypesV3 = true; + QMapIterator itChk(rnxSat.obs); + while (itChk.hasNext() && allTypesV3) { + itChk.next(); + if (itChk.key().length() != 3 || type2to3(sys, itChk.key()) != itChk.key()) { + allTypesV3 = false; + } + } + std::vector hlp(header.nTypes(sys)); for (int iTypeV3 = 0; iTypeV3 < header.nTypes(sys); iTypeV3++) { hlp[iTypeV3] = 0; QString typeV3 = header.obsType(sys, iTypeV3); + if (allTypesV3) { + QMap::const_iterator itV3 = rnxSat.obs.constFind(typeV3); + if (itV3 != rnxSat.obs.constEnd() && itV3.value().value != 0.0) { + hlp[iTypeV3] = &itV3.value(); + } + continue; + } QMapIterator itObs(rnxSat.obs); @@ -1559,27 +1560,33 @@ if (header.nTypes(sys)) { - *stream << rnxSat.prn.toString().c_str(); + QString record = rnxSat.prn.toString().c_str(); for (int iTypeV3 = 0; iTypeV3 < header.nTypes(sys); iTypeV3++) { const t_rnxObs* rnxObs = hlp[iTypeV3]; if (rnxObs == 0) { - *stream << QString().leftJustified(16); + record += QString().leftJustified(16); } else { - *stream << QString("%1").arg(rnxObs->value, 14, 'f', 3); + record += QString("%1").arg(rnxObs->value, 14, 'f', 3); if (rnxObs->lli != 0.0) { - *stream << QString("%1").arg(rnxObs->lli, 1); + record += QString("%1").arg(rnxObs->lli, 1); } else { - *stream << ' '; + record += ' '; } if (rnxObs->snr != 0.0) { - *stream << QString("%1").arg(rnxObs->snr, 1); + record += QString("%1").arg(rnxObs->snr, 1); } else { - *stream << ' '; + record += ' '; } } } - *stream << Qt::endl; + // Trailing blanks can be removed from the records (RINEX 4.02, 4.) + int len = record.length(); + while (len > 3 && record[len-1] == ' ') { + len--; + } + record.truncate(len); + *stream << record << Qt::endl; } }