source: ntrip/trunk/BNC/src/ssrQc.cpp

// Part of BNC, a utility for retrieving decoding and
// converting GNSS data streams from NTRIP broadcasters.
//
// Copyright (C) 2007
// German Federal Agency for Cartography and Geodesy (BKG)
// http://www.bkg.bund.de
// Czech Technical University Prague, Department of Geodesy
// http://www.fsv.cvut.cz
//
// Email: euref-ip@bkg.bund.de
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation, version 2.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

/* -------------------------------------------------------------------------
 * BKG NTRIP Client
 * -------------------------------------------------------------------------
 *
 * Class:      t_ssrQc
 *
 * Purpose:    Quality control of decoded SSR streams
 *
 * Created:    05-Oct-2026
 *
 * Changes:
 *
 * -----------------------------------------------------------------------*/

#include <algorithm>
#include <cmath>
#include <iomanip>

#include "ssrQc.h"
#include "ssrRefSig.h"
#include "combination/bncbiassnx.h"
#include "bnccore.h"
#include "bncutils.h"
#include "bncconst.h"
#include "bncsettings.h"
#include "combination/bnccomb.h"

using namespace std;

namespace {
  const double MAX_CLK_JUMP   = 0.5;  // [m] clock change beyond its rate between consecutive epochs
  const double EPH_WARMUP     = 120.0; // [s] no ephemeris checks at the start
  const double CMP_SAMPL      = 30.0;  // [s] sampling of the stream comparison
  const double CMP_DELAY      = 60.0;  // [s] epochs compared when data this much newer arrived
  const double HIST_LENGTH    = 300.0; // [s] corrections kept for the comparison
  const double DEV_MIN        = 0.05;  // [m] deviating satellites: RMS above max(DEV_MIN,
  const double DEV_FACTOR     = 3.0;   //     DEV_FACTOR x median RMS of the system)
  const double VTEC_CMP_SAMPL = 300.0;  // [s] sampling of the VTEC comparison
  const double VTEC_CMP_DELAY = 120.0;  // [s] VTEC epochs compared when data this much newer arrived
  const double VTEC_HIST      = 3600.0; // [s] VTEC models kept for the comparison
  const double VTEC_GRID_LAT  = 10.0;   // [deg] grid of the VTEC evaluation, latitudes -80 ... 80
  const double VTEC_GRID_LON  = 15.0;   // [deg]

  double tKey(const bncTime& tt) {
    return tt.gpsw() * 604800.0 + tt.gpssec();
  }
  const char*  TYPE_NAME[]    = {"ORBIT", "CLOCK", "CODE_BIAS", "PHASE_BIAS", "VTEC"};

  double median(QVector<double> vv) {
    if (vv.isEmpty()) {
      return 0.0;
    }
    std::sort(vv.begin(), vv.end());
    int nn = vv.size();
    return (nn % 2) ? vv[nn/2] : 0.5 * (vv[nn/2-1] + vv[nn/2]);
  }

  bncTime currentTime() {
    int    week = 0;
    double sec  = 0.0;
    currentGPSWeeks(week, sec);
    return bncTime(week, sec);
  }

  // satellite list "G01 G05(12)", counts shown if not all epochs concerned
  QString satList(const QMap<QString, unsigned>& sats, unsigned maxNum = 20) {
    QString str;
    unsigned num = 0;
    QMapIterator<QString, unsigned> it(sats);
    while (it.hasNext()) {
      it.next();
      if (++num > maxNum) {
        str += QString(" ... (%1 satellites)").arg(sats.size());
        break;
      }
      str += " " + it.key();
      if (it.value() > 0) {
        str += QString("(%1)").arg(it.value());
      }
    }
    return str.isEmpty() ? " -" : str;
  }
}

// Constructor
////////////////////////////////////////////////////////////////////////////
t_ssrQc::t_ssrQc(const QString& dirName, double interval, const QStringList& streams,
                 const QString& reference) : _ephUser(true) {
  _dirName   = dirName;
  _reference = reference.isEmpty() ? QString("CONSENSUS") : reference;
  _interval = interval > 0.0 ? interval : 3600.0;
  // "ALL" or a list of streams; streams with a leading '!' are excluded
  // (from the reports and from the consensus)
  bool all = false;
  for (const QString& sta : streams) {
    if (sta.startsWith('!')) {
      _excluded.insert(sta.mid(1));
    }
    else if (sta.toUpper() == "ALL") {
      all = true;
    }
    else {
      _streams.insert(sta);
    }
  }
  if (all) {
    _streams.clear();
  }

  connect(this, SIGNAL(newMessage(QByteArray,bool)), BNC_CORE, SLOT(slotMessage(const QByteArray,bool)));
  for (const QString& err : t_ssrRefSig::init()) {
    emit newMessage(("SSR QC: " + err).toLatin1(), true);
  }
  // reference code biases from a Bias SINEX file (optional)
  _bsxRef = 0;
  bncSettings settings;
  QString bsxFile = settings.value("ssrQcBiasFile").toString();
  expandEnvVar(bsxFile);
  if (!bsxFile.isEmpty()) {
    _bsxRef = new bncBiasSnx(true);
    if (_bsxRef->readFile(bsxFile) == success) {
      _bsxName = QFileInfo(bsxFile).fileName();
    }
    else {
      emit newMessage(("SSR QC: wrong Bias SINEX file " + bsxFile).toLatin1(), true);
      delete _bsxRef;
      _bsxRef = 0;
    }
  }
  connect(BNC_CORE, SIGNAL(newOrbCorrections(QList<t_orbCorr>)),
          this, SLOT(slotNewOrbCorrections(QList<t_orbCorr>)), Qt::DirectConnection);
  connect(BNC_CORE, SIGNAL(newClkCorrections(QList<t_clkCorr>)),
          this, SLOT(slotNewClkCorrections(QList<t_clkCorr>)), Qt::DirectConnection);
  connect(BNC_CORE, SIGNAL(newCodeBiases(QList<t_satCodeBias>)),
          this, SLOT(slotNewCodeBiases(QList<t_satCodeBias>)), Qt::DirectConnection);
  connect(BNC_CORE, SIGNAL(newPhaseBiases(QList<t_satPhaseBias>)),
          this, SLOT(slotNewPhaseBiases(QList<t_satPhaseBias>)), Qt::DirectConnection);
  connect(BNC_CORE, SIGNAL(newTec(t_vTec)),
          this, SLOT(slotNewTec(t_vTec)), Qt::DirectConnection);
  // combined corrections (stream INTERNAL)
  if (bncComb::exists()) {
    connect(BNC_CMB, SIGNAL(newOrbCorrections(QList<t_orbCorr>)),
            this, SLOT(slotNewOrbCorrections(QList<t_orbCorr>)), Qt::DirectConnection);
    connect(BNC_CMB, SIGNAL(newClkCorrections(QList<t_clkCorr>)),
            this, SLOT(slotNewClkCorrections(QList<t_clkCorr>)), Qt::DirectConnection);
  }
}

// Stream selected for the QC
////////////////////////////////////////////////////////////////////////////
bool t_ssrQc::selected(const string& staID) const {
  QString sta = QString::fromStdString(staID);
  return !_excluded.contains(sta) && (_streams.isEmpty() || _streams.contains(sta));
}

// Analysis center of a stream: the 3 characters of the AC in mountpoints
// like SSRA00CNE1 or OSBC00WHU1, otherwise the stream name itself
////////////////////////////////////////////////////////////////////////////
QString t_ssrQc::acOf(const QString& staID) {
  static const QRegularExpression re("^(?:SSR[AC]|OSB[AC]|IONO)[0-9]{2}([A-Z0-9]{3})[0-9]$");
  QRegularExpressionMatch mm = re.match(staID);
  return mm.hasMatch() ? mm.captured(1) : staID;
}

// Solution of an AC: AC and the two digits of the mountpoint, e.g. GFZ00 for
// SSRA00GFZ0/SSRC00GFZ1 and GFZ01 for SSRA01GFZ0; the streams of a solution
// (formats, APC/CoM) have one vote in the consensus
////////////////////////////////////////////////////////////////////////////
QString t_ssrQc::solOf(const QString& staID) {
  static const QRegularExpression re("^(?:SSR[AC]|OSB[AC]|IONO)([0-9]{2})([A-Z0-9]{3})[0-9]$");
  QRegularExpressionMatch mm = re.match(staID);
  return mm.hasMatch() ? mm.captured(2) + mm.captured(1) : staID;
}

// AC of a solution (GFZ01 -> GFZ)
////////////////////////////////////////////////////////////////////////////
QString t_ssrQc::acOfSol(const QString& sol) {
  static const QRegularExpression re("^([A-Z0-9]{3})[0-9]{2}$");
  QRegularExpressionMatch mm = re.match(sol);
  return mm.hasMatch() ? mm.captured(1) : sol;
}

// Stream part of the consensus: not excluded, and not a combination product
// (INTERNAL, IGS combinations), as these are derived from the other streams
////////////////////////////////////////////////////////////////////////////
bool t_ssrQc::votes(const QString& staID) const {
  return !_excluded.contains(staID) && staID != "INTERNAL" && acOf(staID) != "IGS";
}

// Daily QC file of a stream: <staID>_S_<YYYYDDD>0000_01D_QC.txt
////////////////////////////////////////////////////////////////////////////
ofstream& t_ssrQc::streamOut(const QString& staID, const bncTime& time) {
  unsigned year, month, day;
  time.civil_date(year, month, day);
  QDate date(year, month, day);
  QString name = QDir(_dirName).filePath(QString("%1_S_%2%3%4_01D_QC.txt").arg(staID)
                 .arg(year).arg(date.dayOfYear(), 3, 10, QChar('0')).arg("0000"));
  if (_outNames.value(staID) != name) {
    delete _outs.value(staID, 0);
    bool exists = QFile::exists(name);
    ofstream* out = new ofstream(name.toLocal8Bit().data(), ios_base::out | ios_base::app);
    if (!exists) {
      *out << "SSR QC of stream " << staID.toStdString() << ", report interval " << _interval << " s" << endl
           << "Checks: epochs and update intervals, gaps (interval > 2 x max(declared, median)), latency,"
           << " completeness w.r.t. healthy broadcast ephemerides," << endl
           << "IOD matching, code and phase bias signals, yaw angle, code biases of the reference signals, clock jumps > "
           << MAX_CLK_JUMP << " m, largest orbit corrections, VTEC models (system V)" << endl
           << "Reference signals of the clocks: " << t_ssrRefSig::currentString().toStdString() << endl;
    }
    _outs[staID]     = out;
    _outNames[staID] = name;
  }
  return *_outs[staID];
}

// Destructor
////////////////////////////////////////////////////////////////////////////
t_ssrQc::~t_ssrQc() {
  QMutexLocker locker(&_mutex);
  writeReport(true);
  for (ofstream* out : _outs) {
    delete out;
  }
  delete _bsxRef;
}

// Key of the epoch statistics
////////////////////////////////////////////////////////////////////////////
QString t_ssrQc::key(const string& staID, char sys, e_type type) {
  return QString("%1 %2 %3").arg(staID.c_str()).arg(sys).arg(int(type));
}

// Reference signals of the clocks (as in the combination, option ssrRefSignals)
////////////////////////////////////////////////////////////////////////////
bool t_ssrQc::refSignals(char sys, QString& sig1, QString& sig2,
                         t_frequency::type& frq1, t_frequency::type& frq2) {
  t_ssrRefSig::t_pair pair;
  if (!t_ssrRefSig::get(sys, pair)) {
    return false;
  }
  sig1 = QString::fromStdString(pair.sig1());
  sig2 = QString::fromStdString(pair.sig2());
  frq1 = pair.frq1;
  frq2 = pair.frq2;
  return true;
}

// Start a new report interval if the epoch is beyond the current one
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::checkInterval(const bncTime& time) {
  if (!_intervalStart.valid()) {
    double sec = floor(time.gpssec() / _interval) * _interval;
    _intervalStart = bncTime(time.gpsw(), sec);
  }
  // stream comparison of the epochs already complete for all streams
  if (!_nextCmpEpoch.valid()) {
    double sec = ceil(time.gpssec() / CMP_SAMPL) * CMP_SAMPL;
    _nextCmpEpoch = bncTime(time.gpsw(), sec);
  }
  while (time - _nextCmpEpoch >= CMP_DELAY) {
    compareEpoch(_nextCmpEpoch);
    _nextCmpEpoch = _nextCmpEpoch + CMP_SAMPL;
  }
  // VTEC comparison
  if (!_vtecHist.isEmpty()) {
    if (!_nextVtecCmp.valid()) {
      double sec = ceil(time.gpssec() / VTEC_CMP_SAMPL) * VTEC_CMP_SAMPL;
      _nextVtecCmp = bncTime(time.gpsw(), sec);
    }
    while (time - _nextVtecCmp >= VTEC_CMP_DELAY) {
      compareVtec(_nextVtecCmp);
      _nextVtecCmp = _nextVtecCmp + VTEC_CMP_SAMPL;
    }
  }
  if (time - _intervalStart >= _interval) {
    writeReport(false);
    double sec = floor(time.gpssec() / _interval) * _interval;
    _intervalStart = bncTime(time.gpsw(), sec);
  }
  if (!_lastTime.valid() || _lastTime < time) {
    _lastTime = time;
  }
}

// Epoch statistics
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::newEpoch(const string& staID, char sys, e_type type, const bncTime& time,
                       unsigned updateInt, const QString& prn) {
  QString sta = QString::fromStdString(staID);
  if (!_staStart.contains(sta)) {
    _staStart[sta] = time;
  }
  if (!_staLast.contains(sta) || _staLast[sta] < time) {
    _staLast[sta] = time;
  }
  t_epoStat& stat = _epoStat[key(staID, sys, type)];
  if (updateInt < unsigned(ssrUpdateInt.size())) {
    stat.nominal = ssrUpdateInt[updateInt];
  }
  stat.sats.insert(prn);
  if (!stat.last.valid() || stat.last < time) {
    if (stat.last.valid()) {
      stat.intervals.push_back(time - stat.last);
    }
    stat.last = time;
    stat.numEpo += 1;
    stat.latencies.push_back(currentTime() - time);
  }
  else if (time < stat.last) {
    stat.numOutOfOrder += 1;
  }
}

// Broadcast ephemeris the correction refers to
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::checkEph(t_satStat& stat, const QString& prnInt, unsigned iod, const bncTime& time,
                       const string& staID) {
  // not before all broadcast ephemerides may have been received
  if (time - _staStart[QString::fromStdString(staID)] < EPH_WARMUP) {
    return;
  }
  t_eph* ephLast = _ephUser.ephLast(prnInt);
  t_eph* ephPrev = _ephUser.ephPrev(prnInt);
  t_eph* eph = 0;
  if      (ephLast && ephLast->IOD() == iod) eph = ephLast;
  else if (ephPrev && ephPrev->IOD() == iod) eph = ephPrev;
  if (!ephLast) {
    stat.noEph += 1;
  }
  else if (!eph) {
    stat.iodUnknown += 1;
  }
  else if (eph->checkState() == t_eph::unhealthy || eph->checkState() == t_eph::bad) {
    stat.ephUnhealthy += 1;
  }
}

// Orbit corrections
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::slotNewOrbCorrections(QList<t_orbCorr> orbCorrections) {
  QMutexLocker locker(&_mutex);
  for (const t_orbCorr& corr : orbCorrections) {
    if (!_excluded.contains(QString::fromStdString(corr._staID)) || corr._staID == _reference.toStdString()) {
      storeHist(&corr, 0);
    }
    if (!selected(corr._staID)) continue;
    QString prn = QString::fromStdString(corr._prn.toString());
    checkInterval(corr._time);
    newEpoch(corr._staID, corr._prn.system(), ORB, corr._time, corr._updateInt, prn);
    t_satStat& stat = _satStat[QString::fromStdString(corr._staID) + ' ' + prn];
    stat.numOrb += 1;
    checkEph(stat, QString::fromStdString(corr._prn.toInternalString()), corr._iod, corr._time, corr._staID);
    double radial = fabs(corr._xr[0]);
    double alongCross = max(fabs(corr._xr[1]), fabs(corr._xr[2]));
    stat.maxRadial     = max(stat.maxRadial, radial);
    stat.maxAlongCross = max(stat.maxAlongCross, alongCross);
  }
}

// Clock corrections
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::slotNewClkCorrections(QList<t_clkCorr> clkCorrections) {
  QMutexLocker locker(&_mutex);
  for (const t_clkCorr& corr : clkCorrections) {
    if (!_excluded.contains(QString::fromStdString(corr._staID)) || corr._staID == _reference.toStdString()) {
      storeHist(0, &corr);
    }
    if (!selected(corr._staID)) continue;
    QString prn = QString::fromStdString(corr._prn.toString());
    checkInterval(corr._time);
    newEpoch(corr._staID, corr._prn.system(), CLK, corr._time, corr._updateInt, prn);
    t_satStat& stat = _satStat[QString::fromStdString(corr._staID) + ' ' + prn];
    stat.numClk += 1;
    checkEph(stat, QString::fromStdString(corr._prn.toInternalString()), corr._iod, corr._time, corr._staID);

    // jump: change between consecutive epochs (same IOD) beyond the clock rate
    double clk = corr._dClk * t_CST::c;
    if (stat.lastClkTime.valid() && stat.lastIod == corr._iod) {
      double dt = corr._time - stat.lastClkTime;
      int nominal = (corr._updateInt < unsigned(ssrUpdateInt.size())) ? ssrUpdateInt[corr._updateInt] : 0;
      if (dt > 0.0 && dt <= max(2.0 * nominal, 60.0)) {
        double jump = fabs(clk - stat.lastClk - corr._dotDClk * t_CST::c * dt);
        stat.maxJump = max(stat.maxJump, jump);
        if (jump > MAX_CLK_JUMP) {
          stat.clkJumps += 1;
        }
      }
    }
    stat.lastClkTime = corr._time;
    stat.lastIod     = corr._iod;
    stat.lastClk     = clk;
  }
}

// Code biases
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::slotNewCodeBiases(QList<t_satCodeBias> codeBiases) {
  QMutexLocker locker(&_mutex);
  for (const t_satCodeBias& bias : codeBiases) {
    bool sel = selected(bias._staID);
    if (_excluded.contains(QString::fromStdString(bias._staID)) && bias._staID != _reference.toStdString()) continue;
    char sys = bias._prn.system();
    QString prn = QString::fromStdString(bias._prn.toString());
    if (sel) {
      checkInterval(bias._time);
      newEpoch(bias._staID, sys, CBIAS, bias._time, bias._updateInt, prn);
    }
    QString staSys = QString("%1 %2").arg(bias._staID.c_str()).arg(sys);
    QMap<QString, double> values;
    for (const t_frqCodeBias& frqBias : bias._bias) {
      QString sig = QString::fromStdString(frqBias._rnxType2ch);
      _codeSignals[staSys].insert(sig);
      values[sig] = frqBias._value;
    }
    if (!values.isEmpty()) {
      _cbLast[QString::fromStdString(bias._staID)][prn] = values;
    }
    t_satStat& stat = _satStat[QString::fromStdString(bias._staID) + ' ' + prn];
    QString sig1, sig2;
    t_frequency::type frq1, frq2;
    if (refSignals(sys, sig1, sig2, frq1, frq2)) {
      stat.missingRefSig.clear();
      if (!values.contains(sig1)) stat.missingRefSig.insert("C" + sig1);
      if (!values.contains(sig2)) stat.missingRefSig.insert("C" + sig2);
      stat.refBiasOK = stat.missingRefSig.isEmpty();
      if (stat.refBiasOK) {
        double channel = 0.0;
        if (sys == 'R') {
          t_eph* eph = _ephUser.ephLast(QString::fromStdString(bias._prn.toInternalString()));
          if (!eph) {
            stat.refBiasOK = false;
            continue;
          }
          channel = eph->slotNum();
        }
        double f1 = t_CST::freq(frq1, channel);
        double f2 = t_CST::freq(frq2, channel);
        stat.biasIF = (f1 * f1 * values[sig1] - f2 * f2 * values[sig2]) / (f1 * f1 - f2 * f2);
      }
    }
  }
}

// Phase biases
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::slotNewPhaseBiases(QList<t_satPhaseBias> phaseBiases) {
  QMutexLocker locker(&_mutex);
  for (const t_satPhaseBias& bias : phaseBiases) {
    if (!selected(bias._staID)) continue;
    checkInterval(bias._time);
    newEpoch(bias._staID, bias._prn.system(), PBIAS, bias._time, bias._updateInt,
             QString::fromStdString(bias._prn.toString()));
    QString staSys = QString("%1 %2").arg(bias._staID.c_str()).arg(bias._prn.system());
    for (const t_frqPhaseBias& frqBias : bias._bias) {
      _phaseSignals[staSys].insert(QString::fromStdString(frqBias._rnxType2ch));
    }
    // Yaw angle: new RTCM-SSR has an indicator, old RTCM-SSR and IGS-SSR always
    // carry the field, providers without attitude information send zeros
    t_satStat& stat = _satStat[QString::fromStdString(bias._staID) + ' ' +
                               QString::fromStdString(bias._prn.toString())];
    stat.numPB++;
    bool yawProvided = (bias._ssrFormat == ssrRtcmNew) ? (bias._satYawInfoInd != 0)
                                                       : (bias._yaw != 0.0 || bias._yawRate != 0.0);
    if (yawProvided) {
      stat.numYaw++;
    }
  }
}

// Reference DCBs (bias(base) - bias(signal)) of a satellite: median of the
// other ACs (one vote per AC, own AC left out, combinations not voting) or
// the selected reference stream, or from the Bias SINEX file
////////////////////////////////////////////////////////////////////////////
bool t_ssrQc::refDcbs(const QString& staID, const QString& prn, const QString& base, bool fromFile,
                      QMap<QString, double>& dcbs, const QString& onlyAC, QMap<QString, int>* numVotes) {
  dcbs.clear();
  if (numVotes) {
    numVotes->clear();
  }
  if (fromFile) {
    if (!_bsxRef) {
      return false;
    }
    char sys = prn[0].toLatin1();
    t_ssrRefSig::t_pair pair;
    if (!t_ssrRefSig::get(sys, pair)) {
      return false;
    }
    // the DCBs do not depend on the coefficients (as long as they add up to 1)
    double f1 = t_CST::freq(pair.frq1, 0.0);
    double f2 = t_CST::freq(pair.frq2, 0.0);
    double aIF1 =  f1 * f1 / (f1 * f1 - f2 * f2);
    double aIF2 = -f2 * f2 / (f1 * f1 - f2 * f2);
    t_satCodeBias satCodeBias;
    _bsxRef->determineSsrSatCodeBiases(prn, aIF1, aIF2, satCodeBias);
    QMap<QString, double> osb;
    for (const t_frqCodeBias& bb : satCodeBias._bias) {
      osb[QString::fromStdString(bb._rnxType2ch)] = bb._value;
    }
    if (!osb.contains(base)) {
      return false;
    }
    for (auto it = osb.constBegin(); it != osb.constEnd(); ++it) {
      if (it.key() != base) {
        dcbs[it.key()] = osb[base] - it.value();
      }
    }
    return !dcbs.isEmpty();
  }

  // per AC: DCBs of its streams
  QMap<QString, QMap<QString, QVector<double> > > acDcbs;  // AC, signal
  for (auto itS = _cbLast.constBegin(); itS != _cbLast.constEnd(); ++itS) {
    const QString& sta = itS.key();
    if (sta == staID) {
      continue;
    }
    if (_reference == "CONSENSUS") {
      if (!votes(sta)) {
        continue;
      }
      if (onlyAC.isEmpty() ? (acOf(sta) == acOf(staID))
                           : (solOf(sta) != onlyAC || solOf(sta) == solOf(staID))) {
        continue;
      }
    }
    else if (sta != _reference) {
      continue;
    }
    if (!itS.value().contains(prn)) {
      continue;
    }
    const QMap<QString, double>& bb = itS.value()[prn];
    if (!bb.contains(base)) {
      continue;
    }
    for (auto it = bb.constBegin(); it != bb.constEnd(); ++it) {
      if (it.key() != base) {
        acDcbs[solOf(sta)][it.key()].append(bb[base] - it.value());
      }
    }
  }
  QMap<QString, QVector<double> > votesPerSig;
  for (auto itA = acDcbs.constBegin(); itA != acDcbs.constEnd(); ++itA) {
    for (auto it = itA.value().constBegin(); it != itA.value().constEnd(); ++it) {
      votesPerSig[it.key()].append(median(it.value()));
    }
  }
  for (auto it = votesPerSig.constBegin(); it != votesPerSig.constEnd(); ++it) {
    dcbs[it.key()] = median(it.value());
    if (numVotes) {
      (*numVotes)[it.key()] = it.value().size();
    }
  }
  return !dcbs.isEmpty();
}

// Code bias check: the DCBs of the stream (relative to the first reference
// signal, otherwise to another common signal) compared with a reference;
// DCBs do not depend on the clock datum. A sign opposite to the reference
// for most DCBs indicates the SINEX instead of the SSR sign convention.
// References: the Bias SINEX file, the selected reference stream, or each
// other AC separately (a median of two ACs cannot decide which one is wrong)
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::compareCodeBiases(const QString& staID, char sys, QStringList& lines) {
  const double MIN_DCB  = 0.3;   // [m] DCBs used for the sign check
  const double MIN_DEV  = 0.5;   // [m] deviating DCBs: > max(MIN_DEV, DEV_FACTOR x median)
  const int    MIN_VOTE = 3;     // ACs needed for a consensus DCB (deviating satellites)
  if (!_cbLast.contains(staID)) {
    return;
  }
  t_ssrRefSig::t_pair pair;
  if (!t_ssrRefSig::get(sys, pair)) {
    return;
  }
  const QString ref1 = QString::fromStdString(pair.sig1());

  class t_cmp {
   public:
    unsigned numSat = 0, numSign = 0, numOpp = 0;
    double   sumSqr = 0.0;
    QVector<double> absDiffs;
    QMap<QString, QPair<QString, double> > maxDiff;   // PRN: signal, difference
    QSet<QString> bases;
    bool opposite() const {return numSign >= 5 && numOpp > numSign / 2;}
    double rms() const {return absDiffs.isEmpty() ? 0.0 : sqrt(sumSqr / absDiffs.size());}
  };
  // mode 0: file, 1: reference stream or median of the other ACs, 2: one other AC
  auto compare = [&](int mode, const QString& onlyAC) {
    t_cmp cc;
    for (auto itP = _cbLast[staID].constBegin(); itP != _cbLast[staID].constEnd(); ++itP) {
      const QString& prn = itP.key();
      if (prn[0].toLatin1() != sys) {
        continue;
      }
      const QMap<QString, double>& bb = itP.value();
      QStringList candidates = QStringList() << ref1 << bb.keys();
      QMap<QString, double> dcbs;
      QMap<QString, int>    numVotes;
      QString base;
      for (const QString& cand : candidates) {
        if (bb.contains(cand) && refDcbs(staID, prn, cand, mode == 0, dcbs, onlyAC, &numVotes)) {
          base = cand;
          break;
        }
      }
      if (base.isEmpty()) {
        continue;
      }
      bool used = false;
      for (auto it = dcbs.constBegin(); it != dcbs.constEnd(); ++it) {
        if (!bb.contains(it.key())) {
          continue;
        }
        if (mode == 1 && _reference == "CONSENSUS" && numVotes.value(it.key()) < MIN_VOTE) {
          continue;
        }
        double dcb  = bb[base] - bb[it.key()];
        double diff = dcb - it.value();
        used = true;
        cc.sumSqr += diff * diff;
        cc.absDiffs.append(fabs(diff));
        if (fabs(diff) > fabs(cc.maxDiff.value(prn).second)) {
          cc.maxDiff[prn] = qMakePair("C" + it.key(), diff);
        }
        if (fabs(it.value()) >= MIN_DCB) {
          cc.numSign++;
          if (dcb * it.value() < 0.0) {
            cc.numOpp++;
          }
        }
      }
      if (used) {
        cc.numSat++;
        cc.bases.insert("C" + base);
      }
    }
    return cc;
  };
  auto baseStr = [](const t_cmp& cc) {
    QStringList ll = cc.bases.values();
    ll.sort();
    return ll.join("/");
  };
  auto deviating = [&](const t_cmp& cc, const QString& refName) {
    if (cc.absDiffs.isEmpty() || cc.opposite()) {
      return;
    }
    double limit = std::max(MIN_DEV, DEV_FACTOR * median(cc.absDiffs));
    QStringList devs;
    for (auto it = cc.maxDiff.constBegin(); it != cc.maxDiff.constEnd(); ++it) {
      if (fabs(it.value().second) > limit) {
        devs << QString("%1(%2 %3)").arg(it.key()).arg(it.value().first).arg(it.value().second, 0, 'f', 2);
      }
    }
    if (!devs.isEmpty()) {
      lines << QString("satellites with DCBs deviating from %1 [m]: %2").arg(refName).arg(devs.join(" "));
    }
  };

  // Bias SINEX file or the selected reference stream
  for (int mode = 0; mode < 2; mode++) {
    if ((mode == 0 && !_bsxRef) || (mode == 1 && _reference == "CONSENSUS")) {
      continue;
    }
    QString refName = (mode == 0) ? QString("Bias SINEX file %1").arg(_bsxName) : _reference;
    t_cmp cc = compare(mode, QString());
    if (cc.absDiffs.isEmpty()) {
      continue;
    }
    lines << QString("code biases compared with %1 (DCBs to %2, %3 satellites, %4 DCBs): RMS of the differences "
                     "%5 m, opposite sign for %6 of %7 DCBs (|DCB| >= %8 m)")
             .arg(refName).arg(baseStr(cc)).arg(cc.numSat).arg(cc.absDiffs.size())
             .arg(cc.rms(), 0, 'f', 3).arg(cc.numOpp).arg(cc.numSign).arg(MIN_DCB, 0, 'f', 1);
    if (cc.opposite()) {
      lines << QString("code biases with the opposite sign of %1 for %2 of %3 DCBs - "
                       "SINEX sign convention instead of the SSR one?").arg(refName).arg(cc.numOpp).arg(cc.numSign);
    }
    deviating(cc, refName);
  }

  // Each other AC separately, and the median of the other ACs for deviating satellites
  if (_reference == "CONSENSUS") {
    QSet<QString> acs;
    for (auto itS = _cbLast.constBegin(); itS != _cbLast.constEnd(); ++itS) {
      if (itS.key() != staID && votes(itS.key()) && solOf(itS.key()) != solOf(staID)) {
        acs.insert(solOf(itS.key()));
      }
    }
    QStringList acList = acs.values();
    acList.sort();
    QStringList items;
    QStringList oppACs;
    int numACs = 0;
    QSet<QString> bases;
    for (const QString& ac : acList) {
      t_cmp cc = compare(2, ac);
      if (cc.absDiffs.isEmpty()) {
        continue;
      }
      bases.unite(cc.bases);
      items << QString("%1 RMS %2 m, opposite sign %3/%4").arg(ac).arg(cc.rms(), 0, 'f', 2)
                                                         .arg(cc.numOpp).arg(cc.numSign);
      // the sign rule counts the solutions of other ACs only
      if (acOfSol(ac) != acOf(staID)) {
        numACs++;
        if (cc.opposite()) {
          oppACs << ac;
        }
      }
    }
    if (!items.isEmpty()) {
      QStringList ll = bases.values();
      ll.sort();
      lines << QString("code biases compared with the other solutions (DCBs to %1): %2").arg(ll.join("/")).arg(items.join("; "));
      if (numACs >= 2 && oppACs.size() > numACs / 2) {
        lines << QString("code biases with the opposite sign of %1 of %2 solutions of other ACs (%3) - "
                         "SINEX sign convention instead of the SSR one?")
                 .arg(oppACs.size()).arg(numACs).arg(oppACs.join(" "));
      }
      deviating(compare(1, QString()), QString("the median of the other ACs (at least %1 solutions)").arg(MIN_VOTE));
    }
  }
}

// Reference point of the orbits of a stream: APC (SSRA...), CoM (SSRC...);
// the combined stream INTERNAL refers to the APC
////////////////////////////////////////////////////////////////////////////
char t_ssrQc::refPoint(const QString& staID) const {
  if (staID == "INTERNAL") {
    return 'A';  // the combined corrections always refer to the APC
  }
  const QString& sta = staID;
  if (sta.startsWith("SSRA")) return 'A';
  if (sta.startsWith("SSRC")) return 'C';
  return '?';
}

// Keep the recent orbit and clock corrections of a stream
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::storeHist(const t_orbCorr* orb, const t_clkCorr* clk) {
  if (orb) {
    QMap<double, t_orbCorr>& hist = _orbHist[QString::fromStdString(orb->_staID)]
                                            [QString::fromStdString(orb->_prn.toInternalString())];
    hist[tKey(orb->_time)] = *orb;
    while (!hist.isEmpty() && hist.firstKey() < tKey(orb->_time) - HIST_LENGTH) {
      hist.erase(hist.begin());
    }
  }
  if (clk) {
    QMap<double, t_clkCorr>& hist = _clkHist[QString::fromStdString(clk->_staID)]
                                            [QString::fromStdString(clk->_prn.toInternalString())];
    hist[tKey(clk->_time)] = *clk;
    while (!hist.isEmpty() && hist.firstKey() < tKey(clk->_time) - HIST_LENGTH) {
      hist.erase(hist.begin());
    }
  }
}

// Compare the orbits and clocks of all streams at epoch tt
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::compareEpoch(const bncTime& tt) {
  double tk = tKey(tt);

  // Satellite positions and clocks of all streams (latest corrections <= tt)
  QMap<QString, QMap<QString, t_satCrd> > crd;   // stream, PRN (external)
  QMapIterator<QString, QMap<QString, QMap<double, t_clkCorr> > > itSta(_clkHist);
  while (itSta.hasNext()) {
    itSta.next();
    const QString& staID = itSta.key();
    QMapIterator<QString, QMap<double, t_clkCorr> > itPrn(itSta.value());
    while (itPrn.hasNext()) {
      itPrn.next();
      const QString& prnInt = itPrn.key();
      QMap<double, t_clkCorr>::const_iterator itC = itPrn.value().upperBound(tk);
      if (itC == itPrn.value().constBegin()) {
        continue;
      }
      const t_clkCorr& clk = (--itC).value();
      const QMap<double, t_orbCorr>& orbHist = _orbHist[staID][prnInt];
      QMap<double, t_orbCorr>::const_iterator itO = orbHist.upperBound(tk);
      if (itO == orbHist.constBegin()) {
        continue;
      }
      const t_orbCorr& orb = (--itO).value();
      t_eph* eph = 0;
      t_eph* ephLast = _ephUser.ephLast(prnInt);
      t_eph* ephPrev = _ephUser.ephPrev(prnInt);
      if      (ephLast && ephLast->IOD() == orb._iod) eph = ephLast;
      else if (ephPrev && ephPrev->IOD() == orb._iod) eph = ephPrev;
      if (!eph || clk._iod != orb._iod) {
        continue;
      }
      eph->setOrbCorr(&orb);
      eph->setClkCorr(&clk);
      ColumnVector xc(6), vv(3);
      if (eph->getCrd(tt, xc, vv, true) != success) {
        continue;
      }
      t_satCrd& sc = crd[staID][QString::fromStdString(clk._prn.toString())];
      sc.xx  = xc.Rows(1,3);
      sc.vv  = vv;
      sc.clk = xc(4) * t_CST::c;
      sc.refPoint = refPoint(staID);
      // clock referred to the reference signals (as in the combination):
      // + B_IF of the stream's code biases; the combination (INTERNAL)
      // refers to them already
      sc.clkOK = (staID == "INTERNAL");
      QString satKey = staID + ' ' + QString::fromStdString(clk._prn.toString());
      if (!sc.clkOK && _satStat.contains(satKey) && _satStat[satKey].refBiasOK) {
        sc.clk  += _satStat[satKey].biasIF;
        sc.clkOK = true;
      }
    }
  }

  // Satellites per system
  QMap<char, QSet<QString> > satsSys;
  QMapIterator<QString, QMap<QString, t_satCrd> > itCrd(crd);
  while (itCrd.hasNext()) {
    itCrd.next();
    for (const QString& prn : itCrd.value().keys()) {
      satsSys[prn[0].toLatin1()].insert(prn);
    }
  }

  QMapIterator<char, QSet<QString> > itSys(satsSys);
  while (itSys.hasNext()) {
    itSys.next();
    const QSet<QString>& sats = itSys.value();
    QStringList staList;
    for (const QString& staID : crd.keys()) {
      for (const QString& prn : sats) {
        if (crd[staID].contains(prn)) {
          staList << staID;
          break;
        }
      }
    }

    // Consensus: one value per AC (median of its streams), so that an AC
    // with several streams (formats, APC/CoM) has one vote; the streams of
    // the AC compared are left out (leave-one-AC-out)
    auto acClocks = [&](const QString& prn, const QMap<QString, double>& off, bool useOff) {
      QMap<QString, QVector<double> > perAC;
      for (const QString& oo : staList) {
        if (!votes(oo) || !crd[oo].contains(prn) || !crd[oo][prn].clkOK) continue;
        if (useOff && !off.contains(oo)) continue;
        perAC[solOf(oo)] << crd[oo][prn].clk - (useOff ? off[oo] : 0.0);
      }
      QMap<QString, double> res;
      for (const QString& ac : perAC.keys()) res[ac] = median(perAC[ac]);
      return res;
    };

    // Clock datum: offset of each stream w.r.t. the per-satellite consensus
    // of all ACs (or w.r.t. the reference stream)
    QMap<QString, double> offset;
    QMap<QString, double> noOff;
    QMap<QString, QMap<QString, double> > rawAC;  // PRN, AC
    if (_reference == "CONSENSUS") {
      for (const QString& prn : sats) rawAC[prn] = acClocks(prn, noOff, false);
    }
    for (const QString& staID : staList) {
      QVector<double> diffs;
      for (const QString& prn : sats) {
        if (!crd[staID].contains(prn) || !crd[staID][prn].clkOK) continue;
        if (_reference == "CONSENSUS") {
          if (rawAC[prn].size() >= 3) diffs << crd[staID][prn].clk - median(rawAC[prn].values().toVector());
        }
        else if (crd.value(_reference).contains(prn) && crd[_reference][prn].clkOK) {
          diffs << crd[staID][prn].clk - crd[_reference][prn].clk;
        }
      }
      if (!diffs.isEmpty()) {
        offset[staID] = median(diffs);
      }
    }
    QMap<QString, QMap<QString, double> > alignedAC;  // PRN, AC
    if (_reference == "CONSENSUS") {
      for (const QString& prn : sats) alignedAC[prn] = acClocks(prn, offset, true);
    }

    // Differences of each selected stream to its reference
    for (const QString& staID : staList) {
      if (!selected(staID.toStdString()) || staID == _reference) {
        continue;
      }
      const QString ownAC = acOf(staID);
      QMap<QString, double> dClk;
      QMap<QString, ColumnVector> dRac;
      for (const QString& prn : sats) {
        if (!crd[staID].contains(prn)) continue;
        const t_satCrd& sc = crd[staID][prn];
        QVector<double> clks, xs[3];
        if (_reference == "CONSENSUS") {
          QMapIterator<QString, double> itA(alignedAC[prn]);
          while (itA.hasNext()) {
            itA.next();
            if (acOfSol(itA.key()) != ownAC) clks << itA.value();
          }
          // orbits: per AC the median of its streams with the same reference point
          QMap<QString, QVector<double> > perAC[3];
          for (const QString& oo : staList) {
            if (!votes(oo) || acOf(oo) == ownAC || !crd[oo].contains(prn)) continue;
            const t_satCrd& rc = crd[oo][prn];
            if (rc.refPoint != sc.refPoint || sc.refPoint == '?') continue;
            for (int ii = 0; ii < 3; ii++) perAC[ii][solOf(oo)] << rc.xx[ii];
          }
          for (const QString& ac : perAC[0].keys()) {
            for (int ii = 0; ii < 3; ii++) xs[ii] << median(perAC[ii][ac]);
          }
          if (clks.size() < 2) clks.clear();
          if (xs[0].size() < 2) for (int ii = 0; ii < 3; ii++) xs[ii].clear();
        }
        else {
          if (!crd.value(_reference).contains(prn)) continue;
          const t_satCrd& rc = crd[_reference][prn];
          if (rc.clkOK) clks << rc.clk;
          if (rc.refPoint == sc.refPoint && sc.refPoint != '?') {
            for (int ii = 0; ii < 3; ii++) xs[ii] << rc.xx[ii];
          }
        }
        if (sc.clkOK && offset.contains(staID) && !clks.isEmpty()) {
          dClk[prn] = sc.clk - offset[staID] - median(clks);
        }
        else if (!sc.clkOK) {
          _cmpStat[staID][prn].clkNotReferred += 1;
        }
        if (!xs[0].isEmpty()) {
          ColumnVector dx(3);
          for (int ii = 0; ii < 3; ii++) dx[ii] = sc.xx[ii] - median(xs[ii]);
          ColumnVector er = sc.xx / sc.xx.NormFrobenius();
          ColumnVector ec = crossproduct(sc.xx, sc.vv);
          ec /= ec.NormFrobenius();
          ColumnVector ea = crossproduct(ec, er);
          ColumnVector rac(3);
          rac[0] = DotProduct(dx, er);
          rac[1] = DotProduct(dx, ea);
          rac[2] = DotProduct(dx, ec);
          dRac[prn] = rac;
        }
      }
      if (dClk.isEmpty() && dRac.isEmpty()) {
        continue;
      }
      // remaining common clock offset of this stream and epoch
      double off = dClk.isEmpty() ? 0.0 : median(dClk.values().toVector());
      _cmpEpochs[staID] += 1;
      QSet<QString> prns;
      for (const QString& prn : dClk.keys()) prns.insert(prn);
      for (const QString& prn : dRac.keys()) prns.insert(prn);
      for (const QString& prn : prns) {
        t_cmpStat& stat = _cmpStat[staID][prn];
        if (dRac.contains(prn)) {
          const ColumnVector& rac = dRac[prn];
          stat.numOrb += 1;
          stat.sR   += rac[0] * rac[0];
          stat.sA   += rac[1] * rac[1];
          stat.sC   += rac[2] * rac[2];
        }
        if (dClk.contains(prn)) {
          double dc = dClk[prn] - off;
          stat.numClk += 1;
          stat.sClk   += dc * dc;
          stat.mClk   += dc;
          if (dRac.contains(prn)) {
            double sis = dRac[prn][0] - dc;
            stat.numSis += 1;
            stat.sSis   += sis * sis;
            stat.mSis   += sis;
          }
        }
      }
    }
  }
}

// Ionospheric model (VTEC)
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::slotNewTec(t_vTec vTec) {
  QMutexLocker locker(&_mutex);
  if (vTec._layers.empty()) {
    return;
  }
  QString sta = QString::fromStdString(vTec._staID);
  if (!_excluded.contains(sta) || sta == _reference) {
    QMap<double, t_vTec>& hist = _vtecHist[sta];
    hist[tKey(vTec._time)] = vTec;
    while (!hist.isEmpty() && tKey(vTec._time) - hist.firstKey() > VTEC_HIST) {
      hist.erase(hist.begin());
    }
  }
  if (selected(vTec._staID)) {
    newEpoch(vTec._staID, 'V', VTEC, vTec._time, vTec._updateInt, "VTEC");
    QStringList layers;
    for (const t_vTecLayer& layer : vTec._layers) {
      layers << QString("%1 km (degree/order %2/%3)").arg(layer._height / 1000.0, 0, 'f', 0)
                .arg(layer._C.Nrows() - 1).arg(layer._C.Ncols() - 1);
    }
    _vtecStat[sta].model = QString("%1 layer(s): %2").arg(vTec._layers.size()).arg(layers.join(", "));
  }
  checkInterval(vTec._time);
}

// VTEC [TECU] of a model at a geocentric latitude and longitude [rad] and
// epoch: spherical harmonics in the sun-fixed frame as used in PPP, sum of
// the layers (negative layer values set to zero as in PPP)
////////////////////////////////////////////////////////////////////////////
double t_ssrQc::vtecAt(const t_vTec& vTec, double lat, double lon, const bncTime& tt, bool& negative) {
  double epoch = fmod(tt.gpssec(), 86400.0);
  double lonS  = fmod(lon + (epoch - 50400.0) * M_PI / 43200.0, 2.0 * M_PI);
  double vtec  = 0.0;
  negative = false;
  for (const t_vTecLayer& layer : vTec._layers) {
    int    N  = layer._C.Nrows() - 1;
    int    M  = layer._C.Ncols() - 1;
    double vv = 0.0;
    for (int n = 0; n <= N; n++) {
      for (int m = 0; m <= min(n, M); m++) {
        double pnm = associatedLegendreFunction(n, m, sin(lat));
        double fac = (m == 0) ? sqrt(2.0 * n + 1.0)
                              : sqrt(2.0 * (2.0 * n + 1.0) * factorial(n - m) / factorial(n + m));
        vv += (layer._C[n][m] * cos(m * lonS) + layer._S[n][m] * sin(m * lonS)) * pnm * fac;
      }
    }
    if (vv < 0.0) {
      negative = true;
      vv = 0.0;
    }
    vtec += vv;
  }
  return vtec;
}

// VTEC comparison: each model (latest one not older than twice its update
// interval) evaluated on a global grid and compared with the median of the
// other AC solutions (one vote per solution, own AC left out, combination
// products not voting) or with the reference stream
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::compareVtec(const bncTime& tt) {
  double tk = tKey(tt);

  QMap<QString, const t_vTec*> models;
  QMapIterator<QString, QMap<double, t_vTec> > itH(_vtecHist);
  while (itH.hasNext()) {
    itH.next();
    QMap<double, t_vTec>::const_iterator it = itH.value().upperBound(tk);
    if (it == itH.value().constBegin()) {
      continue;
    }
    --it;
    const t_vTec& vTec = it.value();
    double interval = (vTec._updateInt < unsigned(ssrUpdateInt.size()) && ssrUpdateInt[vTec._updateInt] > 0)
                    ? ssrUpdateInt[vTec._updateInt] : VTEC_CMP_SAMPL;
    double maxAge = min(3600.0, max(120.0, 2.0 * interval));
    if (tk - it.key() <= maxAge) {
      models[itH.key()] = &vTec;
    }
  }
  if (models.isEmpty()) {
    return;
  }

  // grid values of all models
  QVector<QPair<double, double> > grid;
  for (double lat = -80.0; lat <= 80.0 + 1e-6; lat += VTEC_GRID_LAT) {
    for (double lon = -180.0; lon < 180.0 - 1e-6; lon += VTEC_GRID_LON) {
      grid << QPair<double, double>(lat * M_PI / 180.0, lon * M_PI / 180.0);
    }
  }
  QMap<QString, QVector<double> > val;
  QMapIterator<QString, const t_vTec*> itM(models);
  while (itM.hasNext()) {
    itM.next();
    QVector<double>& vv = val[itM.key()];
    vv.resize(grid.size());
    bool selectedSta = selected(itM.key().toStdString());
    t_vtecStat* st = selectedSta ? &_vtecStat[itM.key()] : 0;
    for (int ii = 0; ii < grid.size(); ii++) {
      bool negative = false;
      vv[ii] = vtecAt(*itM.value(), grid[ii].first, grid[ii].second, tt, negative);
      if (st) {
        if (st->numGrid == 0 || vv[ii] < st->minVal) st->minVal = vv[ii];
        if (st->numGrid == 0 || vv[ii] > st->maxVal) st->maxVal = vv[ii];
        st->numGrid += 1;
        if (negative) st->numNeg += 1;
      }
    }
  }

  // comparison
  for (const QString& staID : models.keys()) {
    if (!selected(staID.toStdString())) {
      continue;
    }
    QMap<QString, QStringList> perSol;   // solution: streams
    for (const QString& oo : models.keys()) {
      if (oo == staID) continue;
      if (_reference == "CONSENSUS") {
        if (!votes(oo) || acOf(oo) == acOf(staID)) continue;
      }
      else if (oo != _reference) {
        continue;
      }
      perSol[solOf(oo)] << oo;
    }
    if (perSol.isEmpty()) {
      continue;
    }
    t_vtecStat& st = _vtecStat[staID];
    st.numCmp += 1;
    st.numSol += perSol.size();
    for (int ii = 0; ii < grid.size(); ii++) {
      QVector<double> solVal;
      for (const QStringList& streams : perSol) {
        QVector<double> vv;
        for (const QString& oo : streams) vv << val[oo][ii];
        solVal << median(vv);
      }
      double diff = val[staID][ii] - median(solVal);
      st.numDiff += 1;
      st.sDiff   += diff;
      st.sDiff2  += diff * diff;
      st.maxAbs   = max(st.maxAbs, fabs(diff));
    }
  }
}


// Write the report of the current interval and reset the statistics
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::writeReport(bool final) {
  if (_epoStat.isEmpty()) {
    return;
  }
  bncTime endTime = final ? _lastTime : _intervalStart + _interval;

  // Healthy and current broadcast ephemerides per system
  QMap<char, QSet<QString> > ephSats;
  for (const QString& prnInt : _ephUser.prnList()) {
    t_eph* eph = _ephUser.ephLast(prnInt);
    if (eph && eph->checkState() == t_eph::ok && !outDatedBcep(eph, endTime)) {
      ephSats[eph->prn().system()].insert(QString::fromStdString(eph->prn().toString()));
    }
  }

  // Streams and systems
  QMap<QString, QSet<char> > streams;
  QMapIterator<QString, t_epoStat> itEpo(_epoStat);
  while (itEpo.hasNext()) {
    itEpo.next();
    QStringList hlp = itEpo.key().split(' ');
    streams[hlp[0]].insert(hlp[1][0].toLatin1());
  }

  QMapIterator<QString, QSet<char> > itSta(streams);
  while (itSta.hasNext()) {
    itSta.next();
    const QString& staID = itSta.key();
    QList<char> systems = itSta.value().values();
    std::sort(systems.begin(), systems.end());
    ofstream& out = streamOut(staID, _intervalStart);
    if (final) {
      endTime = _staLast.value(staID, _lastTime);
      if (endTime < _intervalStart) { // stream without data in the last interval
        endTime = _lastTime;
      }
    }
    out << endl << "=== SSR QC " << staID.toStdString() << ' '
        << _intervalStart.datestr() << ' ' << _intervalStart.timestr(0)
        << " - " << endTime.datestr() << ' ' << endTime.timestr(0)
        << (final ? " (end of data)" : "") << " ===" << endl
        << "  Sys Type        Epochs  Interval nom/med/max [s]  Gaps  Out-of-order"
          << "  Latency med/max [s]  Sats" << endl;
    for (char sys : systems) {
      for (int type = ORB; type <= VTEC; type++) {
        QString kk = QString("%1 %2 %3").arg(staID).arg(sys).arg(type);
        if (!_epoStat.contains(kk)) {
          continue;
        }
        const t_epoStat& stat = _epoStat[kk];
        double medInt = median(stat.intervals);
        double maxInt = stat.intervals.isEmpty() ? 0.0 :
                        *std::max_element(stat.intervals.begin(), stat.intervals.end());
        double ref = max(double(stat.nominal), medInt);
        unsigned gaps = 0;
        for (double dt : stat.intervals) {
          if (ref > 0.0 && dt > 2.0 * ref) {
            ++gaps;
          }
        }
        double maxLat = stat.latencies.isEmpty() ? 0.0 :
                        *std::max_element(stat.latencies.begin(), stat.latencies.end());
        out << "  " << sys << "   " << left << setw(10) << TYPE_NAME[type] << right
              << setw(8) << stat.numEpo
              << setw(10) << stat.nominal << " /" << setw(6) << setprecision(1) << fixed << medInt
              << " /" << setw(7) << maxInt
              << setw(6) << gaps << setw(14) << stat.numOutOfOrder
              << setw(10) << median(stat.latencies) << " /" << setw(7) << maxLat
              << setw(7) << stat.sats.size() << endl;
      }
    }

    // Per system: completeness, ephemeris, biases, plausibility
    for (char sys : systems) {
      // Ionosphere (VTEC): model, values on the grid, comparison
      if (sys == 'V') {
        QStringList lines;
        QString kk = QString("%1 %2 %3").arg(staID).arg(sys).arg(int(VTEC));
        if (_epoStat.contains(kk)) {
          const t_epoStat& stat = _epoStat[kk];
          double medInt = median(stat.intervals);
          if (stat.nominal > 0 && stat.intervals.size() > 1 &&
              (medInt > 1.5 * stat.nominal || medInt < stat.nominal / 1.5)) {
            lines << QString("VTEC: declared update interval %1 s, actual %2 s")
                     .arg(stat.nominal).arg(medInt, 0, 'f', 1);
          }
        }
        if (_vtecStat.contains(staID)) {
          const t_vtecStat& vs = _vtecStat[staID];
          if (!vs.model.isEmpty()) {
            lines << "VTEC model: " + vs.model;
          }
          if (vs.numGrid) {
            lines << QString("VTEC on a global %1 x %2 deg grid: %3 ... %4 TECU, negative (set to zero) %5 %")
                     .arg(VTEC_GRID_LAT, 0, 'f', 0).arg(VTEC_GRID_LON, 0, 'f', 0)
                     .arg(vs.minVal, 0, 'f', 1).arg(vs.maxVal, 0, 'f', 1)
                     .arg(100.0 * vs.numNeg / vs.numGrid, 0, 'f', 1);
          }
          if (vs.numDiff) {
            QString refText = (_reference == "CONSENSUS")
                            ? QString("consensus of the other ACs (%1 solutions on average)")
                              .arg(double(vs.numSol) / vs.numCmp, 0, 'f', 1)
                            : QString("stream %1").arg(_reference);
            double mean = vs.sDiff / vs.numDiff;
            double rms  = sqrt(vs.sDiff2 / vs.numDiff);
            lines << QString("VTEC comparison with %1, %2 epochs: mean %3 TECU, RMS %4 TECU, max |diff| %5 TECU")
                     .arg(refText).arg(vs.numCmp).arg(mean, 0, 'f', 2).arg(rms, 0, 'f', 2)
                     .arg(vs.maxAbs, 0, 'f', 1);
          }
          else {
            lines << "VTEC not compared: no other model available";
          }
        }
        for (const QString& line : lines) {
          out << "  " << sys << " " << line.toStdString() << endl;
        }
        continue;
      }
      QMap<QString, unsigned> missOrb, missClk, missCB, missPB, iodUnknown, unhealthy, noEph,
                              jumps, missRef, noYaw, yawPartly;
      unsigned numPB = 0, numYaw = 0, yawSats = 0;
      double maxRadial = 0.0, maxAlongCross = 0.0;
      QString maxRadialPrn, maxAlongCrossPrn;
      QVector<double> biasIF;
      double  maxBiasIF = 0.0;
      QString biasIFPrn;
      const QSet<QString>& ephs = ephSats[sys];
      for (int type = ORB; type <= PBIAS; type++) {
        QString kk = QString("%1 %2 %3").arg(staID).arg(sys).arg(type);
        if (!_epoStat.contains(kk)) {
          continue;
        }
        QMap<QString, unsigned>& miss = (type == ORB) ? missOrb : (type == CLK) ? missClk :
                                        (type == CBIAS) ? missCB : missPB;
        for (const QString& prn : ephs) {
          if (!_epoStat[kk].sats.contains(prn)) {
            miss[prn] = 0;
          }
        }
      }
      QMapIterator<QString, t_satStat> itSat(_satStat);
      while (itSat.hasNext()) {
        itSat.next();
        QStringList hlp = itSat.key().split(' ');
        if (hlp[0] != staID || hlp[1][0].toLatin1() != sys) {
          continue;
        }
        const QString& prn = hlp[1];
        const t_satStat& stat = itSat.value();
        if (stat.iodUnknown)   iodUnknown[prn] = stat.iodUnknown;
        if (stat.ephUnhealthy) unhealthy[prn]  = stat.ephUnhealthy;
        if (stat.noEph)        noEph[prn]      = stat.noEph;
        if (stat.clkJumps)     jumps[prn]      = stat.clkJumps;
        if (stat.maxRadial > maxRadial) {
          maxRadial    = stat.maxRadial;
          maxRadialPrn = prn;
        }
        if (stat.maxAlongCross > maxAlongCross) {
          maxAlongCross    = stat.maxAlongCross;
          maxAlongCrossPrn = prn;
        }
        if (!stat.missingRefSig.isEmpty()) {
          QStringList sigs = stat.missingRefSig.values();
          sigs.sort();
          missRef[prn + "[" + sigs.join("/") + "]"] = 0;
        }
        if (stat.refBiasOK) {
          biasIF.push_back(stat.biasIF);
          if (biasIFPrn.isEmpty() || fabs(stat.biasIF) > maxBiasIF) {
            maxBiasIF  = fabs(stat.biasIF);
            biasIFPrn  = prn;
          }
        }
        if (stat.numPB) {
          numPB  += stat.numPB;
          numYaw += stat.numYaw;
          if (stat.numYaw) {
            yawSats++;
            if (stat.numYaw < stat.numPB) {
              yawPartly[prn] = stat.numPB - stat.numYaw;
            }
          }
          else {
            noYaw[prn] = 0;
          }
        }
      }
      // only findings are listed
      QStringList lines;
      for (int type = ORB; type <= PBIAS; type++) {
        QString kk = QString("%1 %2 %3").arg(staID).arg(sys).arg(type);
        if (_epoStat.contains(kk)) {
          const t_epoStat& stat = _epoStat[kk];
          double medInt = median(stat.intervals);
          if (stat.nominal > 0 && stat.intervals.size() > 1 &&
              (medInt > 1.5 * stat.nominal || medInt < stat.nominal / 1.5)) {
            lines << QString("%1: declared update interval %2 s, actual %3 s")
                     .arg(TYPE_NAME[type]).arg(stat.nominal).arg(medInt, 0, 'f', 1);
          }
        }
      }
      auto addSats = [&lines](const QString& text, const QMap<QString, unsigned>& sats) {
        if (!sats.isEmpty()) {
          lines << text + ":" + satList(sats);
        }
      };
      addSats("orbit corrections missing (healthy broadcast ephemeris available)", missOrb);
      addSats("clock corrections missing (healthy broadcast ephemeris available)", missClk);
      addSats("code biases missing (healthy broadcast ephemeris available)", missCB);
      addSats("phase biases missing (healthy broadcast ephemeris available)", missPB);
      addSats("corrections referring to an IOD of no received broadcast ephemeris", iodUnknown);
      addSats("corrections for unhealthy broadcast ephemeris", unhealthy);
      addSats("corrections for satellites without broadcast ephemeris", noEph);
      QString staSys = QString("%1 %2").arg(staID).arg(sys);
      if (_codeSignals.contains(staSys)) {
        QStringList sigs = _codeSignals[staSys].values();
        sigs.sort();
        lines << "code bias signals: " + sigs.join(" ");
        addSats("code biases of the reference signals missing", missRef);
        if (biasIF.size() > 1) {
          double mean = 0.0;
          for (double bb : biasIF) mean += bb / biasIF.size();
          double var = 0.0;
          for (double bb : biasIF) var += (bb - mean) * (bb - mean) / biasIF.size();
          // the median shows whether the clocks refer to other signals (all
          // satellites) or single satellites have inconsistent biases
          QVector<double> absIF;
          for (double bb : biasIF) absIF << fabs(bb);
          lines << QString("ionosphere-free combination B_IF of the reference signal biases: "
                           "std. dev. %1 m, median |B_IF| %2 m, max |B_IF| %3 m (%4) (%5 satellites)")
                   .arg(sqrt(var), 0, 'f', 3).arg(median(absIF), 0, 'f', 3).arg(maxBiasIF, 0, 'f', 3)
                   .arg(biasIFPrn).arg(biasIF.size());
        }
      }
      compareCodeBiases(staID, sys, lines);
      if (_phaseSignals.contains(staSys)) {
        QStringList sigs = _phaseSignals[staSys].values();
        sigs.sort();
        lines << "phase bias signals: " + sigs.join(" ");
      }
      if (numPB) {
        if (numYaw == 0) {
          lines << "yaw angle: not provided";
        }
        else {
          lines << QString("yaw angle: provided for %1 of %2 satellites (%3 % of the phase bias records)")
                   .arg(yawSats).arg(yawSats + noYaw.size()).arg(100.0 * numYaw / numPB, 0, 'f', 1);
          addSats("phase biases without yaw angle", noYaw);
          addSats("phase biases partly without yaw angle (records)", yawPartly);
        }
      }
      addSats(QString("clock jumps > %1 m").arg(MAX_CLK_JUMP), jumps);
      if (!maxRadialPrn.isEmpty()) {
        lines << QString("largest orbit corrections: radial %1 m (%2), along/cross %3 m (%4)")
                 .arg(maxRadial, 0, 'f', 2).arg(maxRadialPrn).arg(maxAlongCross, 0, 'f', 2).arg(maxAlongCrossPrn);
      }
      // stream comparison
      if (_cmpStat.contains(staID)) {
        unsigned nSat = 0, nClk = 0, nOrb = 0, nSis = 0;
        double sClk = 0.0, sR = 0.0, sA = 0.0, sC = 0.0, sSis = 0.0;
        double vClk = 0.0, vSis = 0.0;  // sums of squares about the satellite means
        QMap<QString, double> rmsSis, rmsClk, rmsOrb;
        QMap<QString, unsigned> notReferred;
        QMapIterator<QString, t_cmpStat> itC(_cmpStat[staID]);
        while (itC.hasNext()) {
          itC.next();
          if (itC.key()[0].toLatin1() != sys) continue;
          const t_cmpStat& st = itC.value();
          if (st.numClk == 0 && st.numOrb == 0) {
            if (st.clkNotReferred) notReferred[itC.key()] = 0;
            continue;
          }
          ++nSat;
          if (st.clkNotReferred) notReferred[itC.key()] = 0;
          nClk += st.numClk; sClk += st.sClk;
          nOrb += st.numOrb; sR += st.sR; sA += st.sA; sC += st.sC;
          nSis += st.numSis; sSis += st.sSis;
          if (st.numClk) vClk += st.sClk - st.mClk * st.mClk / st.numClk;
          if (st.numSis) vSis += st.sSis - st.mSis * st.mSis / st.numSis;
          if (st.numSis) rmsSis[itC.key()] = sqrt(st.sSis / st.numSis);
          else if (st.numClk) rmsClk[itC.key()] = sqrt(st.sClk / st.numClk);
          if (st.numOrb) rmsOrb[itC.key()] = sqrt((st.sR + st.sA + st.sC) / st.numOrb);
        }
        if (nSat) {
          QString refText = (_reference == "CONSENSUS") ? QString("consensus of the other ACs")
                                                        : QString("stream %1").arg(_reference);
          QStringList parts;
          if (nClk) parts << QString("clock %1 (STD %2)").arg(sqrt(sClk / nClk), 0, 'f', 3)
                                                          .arg(sqrt(max(0.0, vClk) / nClk), 0, 'f', 3);
          if (nOrb) parts << QString("radial %1, along %2, cross %3")
                             .arg(sqrt(sR / nOrb), 0, 'f', 3).arg(sqrt(sA / nOrb), 0, 'f', 3)
                             .arg(sqrt(sC / nOrb), 0, 'f', 3);
          if (nSis) parts << QString("SIS (radial - clock) %1 (STD %2)").arg(sqrt(sSis / nSis), 0, 'f', 3)
                                                                         .arg(sqrt(max(0.0, vSis) / nSis), 0, 'f', 3);
          lines << QString("comparison with %1 (%2 satellites), RMS (STD about the satellite means) [m]: %3")
                   .arg(refText).arg(nSat).arg(parts.join(", "));
          if (!nOrb) {
            lines << "orbits not compared: no other stream with the same reference point (APC/CoM)";
          }
          addSats("clocks not compared (code biases of the reference signals missing)", notReferred);
          auto deviating = [&lines](const QString& what, const QMap<QString, double>& rms) {
            if (rms.size() < 2) return;
            double med = median(rms.values().toVector());
            QStringList dev;
            QMapIterator<QString, double> itR(rms);
            while (itR.hasNext()) {
              itR.next();
              if (itR.value() > max(DEV_MIN, DEV_FACTOR * med)) {
                dev << QString("%1(%2)").arg(itR.key()).arg(itR.value(), 0, 'f', 2);
              }
            }
            if (!dev.isEmpty()) {
              lines << QString("satellites deviating from the reference, RMS %1 [m]: %2").arg(what).arg(dev.join(" "));
            }
          };
          deviating("SIS", rmsSis);
          deviating("clock", rmsClk);
          deviating("orbit (3D)", rmsOrb);
        }
      }
      for (const QString& line : lines) {
        out << "  " << sys << " " << line.toStdString() << endl;
      }
    }
    out.flush();
  }

  emit newMessage(QString("SSR QC: report %1 %2 of %3 stream(s) written to %4")
                  .arg(QString::fromStdString(_intervalStart.datestr()))
                  .arg(QString::fromStdString(_intervalStart.timestr(0)))
                  .arg(streams.size()).arg(_dirName).toLatin1(), false);

  _cmpStat.clear();
  _cmpEpochs.clear();
  for (QMap<QString, t_vtecStat>::iterator it = _vtecStat.begin(); it != _vtecStat.end(); ++it) {
    t_vtecStat fresh;
    fresh.model = it.value().model;
    it.value() = fresh;
  }

  // reset, keeping the state needed across intervals
  for (QMap<QString, t_epoStat>::iterator it = _epoStat.begin(); it != _epoStat.end(); ++it) {
    t_epoStat& stat = it.value();
    bncTime last = stat.last;
    int nominal = stat.nominal;
    stat = t_epoStat();
    stat.last = last;
    stat.nominal = nominal;
  }
  for (QMap<QString, t_satStat>::iterator it = _satStat.begin(); it != _satStat.end(); ++it) {
    t_satStat& stat = it.value();
    t_satStat fresh;
    fresh.lastClkTime   = stat.lastClkTime;
    fresh.lastIod       = stat.lastIod;
    fresh.lastClk       = stat.lastClk;
    fresh.missingRefSig = stat.missingRefSig;
    fresh.refBiasOK     = stat.refBiasOK;
    fresh.biasIF        = stat.biasIF;
    stat = fresh;
  }
}