Changeset 3482 in ntrip


Ignore:
Timestamp:
Oct 18, 2011, 4:49:32 PM (10 years ago)
Author:
mervart
Message:
 
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/BNC/combination/bnccomb.cpp

    r3481 r3482  
    783783                                        ColumnVector& dx) {
    784784
    785   // Remove Satellites that are not in Master
    786   // ----------------------------------------
    787   QMutableVectorIterator<cmbCorr*> it(corrs());
    788   while (it.hasNext()) {
    789     cmbCorr* corr = it.next();
    790     QString  prn  = corr->prn;
    791     bool foundMaster = false;
    792     QVectorIterator<cmbCorr*> itHlp(corrs());
    793     while (itHlp.hasNext()) {
    794       cmbCorr* corrHlp = itHlp.next();
    795       QString  prnHlp  = corrHlp->prn;
    796       QString  ACHlp   = corrHlp->acName;
    797       if (ACHlp == _masterOrbitAC && prn == prnHlp) {
    798         foundMaster = true;
    799         break;
    800       }
    801     }
    802     if (!foundMaster) {
    803       it.remove();
    804     }
    805   }
    806 
    807   // Count Number of Observations per Satellite and per AC
    808   // -----------------------------------------------------
    809   QMap<QString, int> numObsPrn;
    810   QMap<QString, int> numObsAC;
    811   QVectorIterator<cmbCorr*> itCorr(corrs());
    812   while (itCorr.hasNext()) {
    813     cmbCorr* corr = itCorr.next();
    814     QString  prn  = corr->prn;
    815     QString  AC   = corr->acName;
    816     if (numObsPrn.find(prn) == numObsPrn.end()) {
    817       numObsPrn[prn]  = 1;
     785  // Outlier Detection Loop
     786  // ----------------------
     787  while (true) {
     788   
     789    // Remove Satellites that are not in Master
     790    // ----------------------------------------
     791    QMutableVectorIterator<cmbCorr*> it(corrs());
     792    while (it.hasNext()) {
     793      cmbCorr* corr = it.next();
     794      QString  prn  = corr->prn;
     795      bool foundMaster = false;
     796      QVectorIterator<cmbCorr*> itHlp(corrs());
     797      while (itHlp.hasNext()) {
     798        cmbCorr* corrHlp = itHlp.next();
     799        QString  prnHlp  = corrHlp->prn;
     800        QString  ACHlp   = corrHlp->acName;
     801        if (ACHlp == _masterOrbitAC && prn == prnHlp) {
     802          foundMaster = true;
     803          break;
     804        }
     805      }
     806      if (!foundMaster) {
     807        it.remove();
     808      }
     809    }
     810   
     811    // Count Number of Observations per Satellite and per AC
     812    // -----------------------------------------------------
     813    QMap<QString, int> numObsPrn;
     814    QMap<QString, int> numObsAC;
     815    QVectorIterator<cmbCorr*> itCorr(corrs());
     816    while (itCorr.hasNext()) {
     817      cmbCorr* corr = itCorr.next();
     818      QString  prn  = corr->prn;
     819      QString  AC   = corr->acName;
     820      if (numObsPrn.find(prn) == numObsPrn.end()) {
     821        numObsPrn[prn]  = 1;
     822      }
     823      else {
     824        numObsPrn[prn] += 1;
     825      }
     826      if (numObsAC.find(AC) == numObsAC.end()) {
     827        numObsAC[AC]  = 1;
     828      }
     829      else {
     830        numObsAC[AC] += 1;
     831      }
     832    }
     833   
     834    // Clean-Up the Paramters
     835    // ----------------------
     836    for (int iPar = 1; iPar <= _params.size(); iPar++) {
     837      delete _params[iPar-1];
     838    }
     839    _params.clear();
     840   
     841    // Set new Parameters
     842    // ------------------
     843    int nextPar = 0;
     844   
     845    QMapIterator<QString, int> itAC(numObsAC);
     846    while (itAC.hasNext()) {
     847      itAC.next();
     848      const QString& AC     = itAC.key();
     849      int            numObs = itAC.value();
     850      if (AC != _masterOrbitAC && numObs > 0) {
     851        _params.push_back(new cmbParam(cmbParam::offAC, ++nextPar, AC, ""));
     852      }
     853    }
     854   
     855    QMapIterator<QString, int> itPrn(numObsPrn);
     856    while (itPrn.hasNext()) {
     857      itPrn.next();
     858      const QString& prn    = itPrn.key();
     859      int            numObs = itPrn.value();
     860      if (numObs > 0) {
     861        _params.push_back(new cmbParam(cmbParam::clkSat, ++nextPar, "", prn));
     862      }
     863    } 
     864   
     865    int nPar = _params.size();
     866    ColumnVector x0(nPar);
     867    x0 = 0.0;
     868   
     869    // Create First-Design Matrix
     870    // --------------------------
     871    Matrix         AA;
     872    ColumnVector   ll;
     873    DiagonalMatrix PP;
     874    if (createAmat(AA, ll, PP, x0, resCorr) != success) {
     875      return failure;
     876    }
     877   
     878    ColumnVector vv;
     879    try {
     880      Matrix          ATP = AA.t() * PP;
     881      SymmetricMatrix NN; NN << ATP * AA;
     882      ColumnVector    bb = ATP * ll;
     883      _QQ = NN.i();
     884      dx  = _QQ * bb;
     885      vv  = ll - AA * dx;
     886    }
     887    catch (Exception& exc) {
     888      out << exc.what() << endl;
     889      return failure;
     890    }
     891
     892    int     maxResIndex;
     893    double  maxRes = vv.maximum_absolute_value1(maxResIndex);   
     894    out.setRealNumberNotation(QTextStream::FixedNotation);
     895    out.setRealNumberPrecision(3); 
     896    out << _resTime.datestr().c_str() << " " << _resTime.timestr().c_str()
     897        << " Maximum Residuum " << maxRes << ' '
     898        << corrs()[maxResIndex-1]->acName << ' ' << corrs()[maxResIndex-1]->prn;
     899
     900    if (maxRes > _MAXRES) {
     901      out << "  Outlier" << endl;
     902      corrs().remove(maxResIndex-1);
    818903    }
    819904    else {
    820       numObsPrn[prn] += 1;
    821     }
    822     if (numObsAC.find(AC) == numObsAC.end()) {
    823       numObsAC[AC]  = 1;
    824     }
    825     else {
    826       numObsAC[AC] += 1;
    827     }
    828   }
    829 
    830   // Clean-Up the Paramters
    831   // ----------------------
    832   for (int iPar = 1; iPar <= _params.size(); iPar++) {
    833     delete _params[iPar-1];
    834   }
    835   _params.clear();
    836 
    837   // Set new Parameters
    838   // ------------------
    839   int nextPar = 0;
    840 
    841   QMapIterator<QString, int> itAC(numObsAC);
    842   while (itAC.hasNext()) {
    843     itAC.next();
    844     const QString& AC     = itAC.key();
    845     int            numObs = itAC.value();
    846     if (AC != _masterOrbitAC && numObs > 0) {
    847       _params.push_back(new cmbParam(cmbParam::offAC, ++nextPar, AC, ""));
    848     }
    849   }
    850 
    851   QMapIterator<QString, int> itPrn(numObsPrn);
    852   while (itPrn.hasNext()) {
    853     itPrn.next();
    854     const QString& prn    = itPrn.key();
    855     int            numObs = itPrn.value();
    856     if (numObs > 0) {
    857       _params.push_back(new cmbParam(cmbParam::clkSat, ++nextPar, "", prn));
    858     }
    859   } 
    860 
    861   int nPar = _params.size();
    862   ColumnVector x0(nPar);
    863   x0 = 0.0;
    864 
    865   // Create First-Design Matrix
    866   // --------------------------
    867   Matrix         AA;
    868   ColumnVector   ll;
    869   DiagonalMatrix PP;
    870   if (createAmat(AA, ll, PP, x0, resCorr) != success) {
    871     return failure;
    872   }
    873 
    874   ColumnVector vv;
    875   try {
    876     Matrix          ATP = AA.t() * PP;
    877     SymmetricMatrix NN; NN << ATP * AA;
    878     ColumnVector    bb = ATP * ll;
    879     _QQ = NN.i();
    880     dx  = _QQ * bb;
    881     vv  = ll - AA * dx;
    882   }
    883   catch (Exception& exc) {
    884     out << exc.what() << endl;
    885     return failure;
    886   }
    887 
    888   out.setRealNumberNotation(QTextStream::FixedNotation);
    889   out.setRealNumberPrecision(3); 
    890   for (int ii = 0; ii < vv.Nrows(); ii++) {
    891     const cmbCorr* corr = corrs()[ii];
    892     out << _resTime.datestr().c_str() << ' '
    893         << _resTime.timestr().c_str() << " "
    894         << corr->acName << ' ' << corr->prn;
    895     out.setFieldWidth(6);
    896     out << " dClk = " << corr->dClk * t_CST::c << " res = " << vv[ii] << endl;
    897     out.setFieldWidth(0);
    898   }
    899 
    900   return success;
    901 }
     905      out << "  OK" << endl;
     906      out.setRealNumberNotation(QTextStream::FixedNotation);
     907      out.setRealNumberPrecision(3); 
     908      for (int ii = 0; ii < vv.Nrows(); ii++) {
     909        const cmbCorr* corr = corrs()[ii];
     910        out << _resTime.datestr().c_str() << ' '
     911            << _resTime.timestr().c_str() << " "
     912            << corr->acName << ' ' << corr->prn;
     913        out.setFieldWidth(6);
     914        out << " dClk = " << corr->dClk * t_CST::c << " res = " << vv[ii] << endl;
     915        out.setFieldWidth(0);
     916      }
     917      return success;
     918    }
     919
     920  }
     921
     922  return failure;
     923}
Note: See TracChangeset for help on using the changeset viewer.