Changeset 9642 in ntrip for trunk/BNC/src/PPP


Ignore:
Timestamp:
Mar 4, 2022, 10:38:38 AM (3 years ago)
Author:
stuerze
Message:

minor changes regarding PPP

Location:
trunk/BNC/src/PPP
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/BNC/src/PPP/pppFilter.cpp

    r9641 r9642  
    3333// Constructor
    3434////////////////////////////////////////////////////////////////////////////
    35 t_pppFilter::t_pppFilter(t_pppObsPool* obsPool) {
     35t_pppFilter::t_pppFilter(t_pppObsPool *obsPool) {
    3636  _numSat = 0;
    3737  _obsPool = obsPool;
    38   _refPrn  = t_prn();
     38  _refPrn = t_prn();
    3939  _datumTrafo = new t_datumTrafo();
    4040}
    … …  
    4949////////////////////////////////////////////////////////////////////////////
    5050t_irc t_pppFilter::processEpoch() {
    51   _numSat     = 0;
     51  _numSat = 0;
    5252  const double maxSolGap = 60.0;
    5353
    5454  // Vector of all Observations
    5555  // --------------------------
    56   t_pppObsPool::t_epoch* epoch = _obsPool->lastEpoch();
     56  t_pppObsPool::t_epoch *epoch = _obsPool->lastEpoch();
    5757  if (!epoch) {
    5858    return failure;
    5959  }
    60   vector<t_pppSatObs*>& allObs = epoch->obsVector();
     60  vector<t_pppSatObs*> &allObs = epoch->obsVector();
    6161
    6262  // Time of the Epoch
    … …  
    6464  _epoTime = epoch->epoTime();
    6565
    66   if (!_firstEpoTime.valid() ||
    67       !_lastEpoTimeOK.valid() ||
    68       (maxSolGap > 0.0 && _epoTime - _lastEpoTimeOK > maxSolGap)) {
     66  if (!_firstEpoTime.valid() || !_lastEpoTimeOK.valid()
     67      || (maxSolGap > 0.0 && _epoTime - _lastEpoTimeOK > maxSolGap)) {
    6968    _firstEpoTime = _epoTime;
    7069  }
    … …  
    7271  string epoTimeStr = string(_epoTime);
    7372
    74   const QMap<char, t_pppRefSat*>& refSatMap = epoch->refSatMap();
     73  const QMap<char, t_pppRefSat*> &refSatMap = epoch->refSatMap();
    7574
    7675  //LOG << "processEpoch: printParams before set" << endl;
    … …  
    8483
    8584  if (OPT->_obsModelType == OPT->DCMcodeBias ||
    86       OPT->_obsModelType == OPT->DCMphaseBias) {
     85  OPT->_obsModelType == OPT->DCMphaseBias) {
    8786#ifdef BNC_DEBUG_PPP
    8887  _parlist.printParams(_epoTime);
    … …  
    9190  // Status Vector, Variance-Covariance Matrix
    9291  // -----------------------------------------
    93   ColumnVector    xFltOld = _xFlt;
     92  ColumnVector xFltOld = _xFlt;
    9493  SymmetricMatrix QFltOld = _QFlt;
    9594  setStateVectorAndVarCovMatrix(xFltOld, QFltOld);
    … …  
    9998  bool preProcessing = false;
    10099  if (OPT->_obsModelType == OPT->DCMcodeBias ||
    101       OPT->_obsModelType == OPT->DCMphaseBias) {
     100  OPT->_obsModelType == OPT->DCMphaseBias) {
    102101    preProcessing = true;
    103     const QList<char>& usedSystems = _parlist.usedSystems();
     102    const QList<char> &usedSystems = _parlist.usedSystems();
    104103    for (int iSys = 0; iSys < usedSystems.size(); iSys++) {
    105104      char sys = usedSystems[iSys];
    … …  
    108107        _refPrn = refSatMap[sys]->prn();
    109108      }
    110      vector<t_pppSatObs*> obsVector;
     109      vector<t_pppSatObs*> obsVector;
    111110      for (unsigned jj = 0; jj < allObs.size(); jj++) {
    112111        if (allObs[jj]->prn().system() == sys) {
    … …  
    118117      }
    119118      if (processSystem(OPT->LCs(sys), obsVector, _refPrn,
    120                         epoch->pseudoObsIono(), preProcessing) != success) {
     119          epoch->pseudoObsIono(), preProcessing) != success) {
    121120        LOG << sys << ": processSystem !=  success (pre-processing)" << endl;
    122121        _xFlt = xFltOld;
    … …  
    133132      _QFlt = QFltOld;
    134133      return success;
    135     }
    136     else if (!_obsPool->refSatChangeRequired()) {
     134    } else if (!_obsPool->refSatChangeRequired()) {
    137135      initDatumTransformation(allObs, epoch->pseudoObsIono());
    138136    }
    … …  
    142140  // ------------------------------------
    143141  preProcessing = false;
    144   const QList<char>& usedSystems = _parlist. usedSystems();
     142  const QList<char> &usedSystems = _parlist.usedSystems();
    145143  for (int iSys = 0; iSys < usedSystems.size(); iSys++) {
    146144    char sys = usedSystems[iSys];
    … …  
    160158      _datumTrafo->setFirstSystem(sys);
    161159    }
    162     LOG << epoTimeStr << " SATNUM " << sys << ' ' << right << setw(2) << num << endl;
    163     if (processSystem(OPT->LCs(sys), obsVector, _refPrn,
    164                       epoch->pseudoObsIono(), preProcessing) != success) {
     160    LOG << epoTimeStr << " SATNUM " << sys << ' ' << right << setw(2) << num
     161        << endl;
     162    if (processSystem(OPT->LCs(sys), obsVector, _refPrn, epoch->pseudoObsIono(),
     163        preProcessing) != success) {
    165164      LOG << "processSystem !=  success (fin-processing)" << endl;
    166165      if (OPT->_obsModelType == OPT->DCMcodeBias ||
    167           OPT->_obsModelType == OPT->DCMphaseBias) {
     166      OPT->_obsModelType == OPT->DCMphaseBias) {
    168167        _xFlt = xFltOld;
    169168        _QFlt = QFltOld;
    … …  
    179178  cmpDOP(allObs, refSatMap);
    180179  _parlist.printResult(_epoTime, _QFlt, _xFlt);
    181   _lastEpoTimeOK = _epoTime;  // remember time of last successful epoch processing
     180  _lastEpoTimeOK = _epoTime; // remember time of last successful epoch processing
    182181  return success;
    183182}
    … …  
    185184// Process Selected LCs
    186185////////////////////////////////////////////////////////////////////////////
    187 t_irc t_pppFilter::processSystem(const vector<t_lc::type>& LCs,
    188                                  const vector<t_pppSatObs*>& obsVector,
    189                                  const t_prn& refPrn,
    190                                  bool pseudoObsIonoAvailable,
    191                                  bool preProcessing) {
     186t_irc t_pppFilter::processSystem(const vector<t_lc::type> &LCs,
     187    const vector<t_pppSatObs*> &obsVector, const t_prn &refPrn,
     188    bool pseudoObsIonoAvailable, bool preProcessing) {
    192189  LOG.setf(ios::fixed);
    193190  char sys = refPrn.system();
    … …  
    202199  }
    203200
    204   ColumnVector               xSav      = _xFlt;
    205   SymmetricMatrix            QSav      = _QFlt;
    206   string                     epoTimeStr = string(_epoTime);
    207   const vector<t_pppParam*>& params    = _parlist.params();
     201  ColumnVector xSav = _xFlt;
     202  SymmetricMatrix QSav = _QFlt;
     203  string epoTimeStr = string(_epoTime);
     204  const vector<t_pppParam*> &params = _parlist.params();
    208205  unsigned nPar = _parlist.nPar();
    209206
    210   unsigned usedLCs     = LCs.size();
     207  unsigned usedLCs = LCs.size();
    211208  if (OPT->_pseudoObsIono && !pseudoObsIonoAvailable) {
    212       usedLCs -= 1;  // GIM not used
     209    usedLCs -= 1;  // GIM not used
    213210  }
    214211  // max Obs
    … …  
    230227    // First-Design Matrix, Terms Observed-Computed, Weight Matrix
    231228    // -----------------------------------------------------------
    232     Matrix                AA(maxObs, nPar);
    233     ColumnVector          ll(maxObs);
    234     DiagonalMatrix        PP(maxObs); PP = 0.0;
     229    Matrix AA(maxObs, nPar);
     230    ColumnVector ll(maxObs);
     231    DiagonalMatrix PP(maxObs);
     232    PP = 0.0;
    235233
    236234    int iObs = -1;
    237235    vector<t_pppSatObs*> usedObs;
    238     vector<t_lc::type>   usedTypes;
     236    vector<t_lc::type> usedTypes;
    239237
    240238    // Real Observations
    … …  
    242240    double nSat = 0;
    243241    for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    244       t_pppSatObs* obs = obsVector[ii];
     242      t_pppSatObs *obs = obsVector[ii];
    245243      if (iOutlier == 0 && !preProcessing) {
    246244        obs->resetOutlier();
    … …  
    250248        for (unsigned jj = 0; jj < usedLCs; jj++) {
    251249          const t_lc::type tLC = LCs[jj];
    252           if (tLC == t_lc::GIM) {continue;}
     250          if (tLC == t_lc::GIM) {
     251            continue;
     252          }
    253253          ++iObs;
    254254          usedObs.push_back(obs);
    255255          usedTypes.push_back(tLC);
    256256          for (unsigned iPar = 0; iPar < nPar; iPar++) {
    257             const t_pppParam* par = params[iPar];
     257            const t_pppParam *par = params[iPar];
    258258            AA[iObs][iPar] = par->partial(_epoTime, obs, tLC, refPrn);
    259259          }
    260           ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1));
     260          ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC)
     261              - DotProduct(_x0, AA.Row(iObs + 1));
    261262          PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC));
    262263        }
    … …  
    271272    }
    272273
    273     if ((!iOutlier) &&
    274         (OPT->_obsModelType == OPT->DCMcodeBias ||
    275          OPT->_obsModelType == OPT->DCMphaseBias) &&
    276         (!preProcessing)) {
    277       _datumTrafo->updateIndices(sys, iObs+1);
    278       _datumTrafo->prepareAA(AA.SubMatrix(1, iObs+1 , 1, _parlist.nPar()), 1);
     274    if ((!iOutlier) && (OPT->_obsModelType == OPT->DCMcodeBias ||
     275    OPT->_obsModelType == OPT->DCMphaseBias) && (!preProcessing)) {
     276      _datumTrafo->updateIndices(sys, iObs + 1);
     277      _datumTrafo->prepareAA(AA.SubMatrix(1, iObs + 1, 1, _parlist.nPar()), 1);
    279278    }
    280279
    … …  
    283282    if (OPT->_pseudoObsIono && pseudoObsIonoAvailable) {
    284283      for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    285         t_pppSatObs* obs = obsVector[ii];
     284        t_pppSatObs *obs = obsVector[ii];
    286285        if (!obs->outlier()) {
    287286          for (unsigned jj = 0; jj < usedLCs; jj++) {
    … …  
    289288            if (tLC == t_lc::GIM && !obs->isReference()) {
    290289              ++iObs;
    291             } else {continue;}
     290            } else {
     291              continue;
     292            }
    292293            usedObs.push_back(obs);
    293294            usedTypes.push_back(tLC);
    294295            for (unsigned iPar = 0; iPar < nPar; iPar++) {
    295               const t_pppParam* par = params[iPar];
     296              const t_pppParam *par = params[iPar];
    296297              AA[iObs][iPar] = par->partial(_epoTime, obs, tLC, refPrn);
    297298            }
    298             ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1));
     299            ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC)
     300                - DotProduct(_x0, AA.Row(iObs + 1));
    299301            PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC));
    300302          }
    … …  
    305307    // Truncate matrices
    306308    // -----------------
    307     AA = AA.Rows(1, iObs+1);
    308     ll = ll.Rows(1, iObs+1);
    309     PP = PP.SymSubMatrix(1, iObs+1);
     309    AA = AA.Rows(1, iObs + 1);
     310    ll = ll.Rows(1, iObs + 1);
     311    PP = PP.SymSubMatrix(1, iObs + 1);
    310312
    311313    // Kalman update step
    … …  
    316318    // ---------------
    317319    ColumnVector vv = AA * _xFlt - ll;
    318     double     maxOutlier      = 0.0;
    319     int        maxOutlierIndex = -1;
    320     t_lc::type maxOutlierLC    = t_lc::dummy;
     320    double maxOutlier = 0.0;
     321    int maxOutlierIndex = -1;
     322    t_lc::type maxOutlierLC = t_lc::dummy;
    321323    for (unsigned ii = 0; ii < usedObs.size(); ii++) {
    322324      const t_lc::type tLC = usedTypes[ii];
    … …  
    324326      if (res > usedObs[ii]->maxRes(tLC)) {
    325327        if (res > fabs(maxOutlier)) {
    326           maxOutlier      = vv[ii];
     328          maxOutlier = vv[ii];
    327329          maxOutlierIndex = ii;
    328           maxOutlierLC    = tLC;
     330          maxOutlierLC = tLC;
    329331        }
    330332      }
    … …  
    334336    // --------------------------------------------
    335337    if (maxOutlierIndex > -1) {
    336       t_pppSatObs* obs = usedObs[maxOutlierIndex];
    337       t_pppParam* par = 0;
     338      t_pppSatObs *obs = usedObs[maxOutlierIndex];
     339      t_pppParam *par = 0;
    338340      for (unsigned iPar = 0; iPar < nPar; iPar++) {
    339         t_pppParam* hlp = params[iPar];
    340         if (hlp->type() == t_pppParam::amb &&
    341             hlp->prn()  == obs->prn() &&
    342             hlp->tLC()  == usedTypes[maxOutlierIndex]) {
     341        t_pppParam *hlp = params[iPar];
     342        if (hlp->type() == t_pppParam::amb && hlp->prn() == obs->prn()
     343            && hlp->tLC() == usedTypes[maxOutlierIndex]) {
    343344          par = hlp;
    344345        }
    … …  
    346347      if (preProcessing) {
    347348        // for refSats no ambiguity parameter exists
    348          if ((obs->prn() == refPrn) &&
    349              (t_lc::toString(maxOutlierLC) == "l1" ||
    350               t_lc::toString(maxOutlierLC) == "l2" )) {
    351            _obsPool->setRefSatChangeRequired(sys, true);
    352            LOG << epoTimeStr << " Outlier ("
    353                << ((preProcessing) ? "pre-processing) " : "fin-processing) ") << t_lc::toString(maxOutlierLC) << ' '
    354                << obs->prn().toString()  << ' ' << setw(8) << setprecision(4) << maxOutlier << endl;
    355            break;
    356          }
    357          else {
    358            obs->setOutlier();
    359         }
    360       }
    361       else {// fin-processing
     349        if ((obs->prn() == refPrn)
     350            && (t_lc::toString(maxOutlierLC) == "l1"
     351                || t_lc::toString(maxOutlierLC) == "l2")) {
     352          _obsPool->setRefSatChangeRequired(sys, true);
     353          LOG << epoTimeStr << " Outlier ("
     354              << ((preProcessing) ? "pre-processing) " : "fin-processing) ")
     355              << t_lc::toString(maxOutlierLC) << ' ' << obs->prn().toString()
     356              << ' ' << setw(8) << setprecision(4) << maxOutlier << endl;
     357          break;
     358        } else {
     359          obs->setOutlier();
     360        }
     361      } else {    // fin-processing
    362362        LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' '
    363             << obs->prn().toString()  << ' ' << setw(8) << setprecision(4) << maxOutlier << endl;
     363            << obs->prn().toString() << ' ' << setw(8) << setprecision(4)
     364            << maxOutlier << endl;
    364365        if (par) {
    365366          resetAmb(par->prn(), obsVector, &QSav, &xSav);
    366         }
    367         else {
     367        } else {
    368368          obs->setOutlier();
    369369        }
    … …  
    377377          for (unsigned ii = 0; ii < usedObs.size(); ii++) {
    378378            const t_lc::type tLC = usedTypes[ii];
    379             t_pppSatObs* obs = usedObs[ii];
     379            t_pppSatObs *obs = usedObs[ii];
    380380            if (tLC == LCs[jj]) {
    381381              obs->setRes(tLC, vv[ii]);
    382               LOG << epoTimeStr << " RES "
    383                   << left << setw(3) << t_lc::toString(tLC) << right << ' ';
    384               if (t_lc::toString(tLC) == "Tz0") {LOG << sys << "  ";}
    385               else {LOG << obs->prn().toString();}
     382              LOG << epoTimeStr << " RES " << left << setw(3)
     383                  << t_lc::toString(tLC) << right << ' ';
     384              if (t_lc::toString(tLC) == "Tz0") {
     385                LOG << sys << "  ";
     386              } else {
     387                LOG << obs->prn().toString();
     388              }
    386389              LOG << setw(9) << setprecision(4) << vv[ii] << endl;
    387390            }
    … …  
    397400// Cycle-Slip Detection
    398401////////////////////////////////////////////////////////////////////////////
    399 t_irc t_pppFilter::detectCycleSlips(const vector<t_lc::type>& LCs,
    400                                     const vector<t_pppSatObs*>& obsVector,
    401                                     const t_prn& refPrn,
    402                                     bool preProcessing) {
     402t_irc t_pppFilter::detectCycleSlips(const vector<t_lc::type> &LCs,
     403    const vector<t_pppSatObs*> &obsVector, const t_prn &refPrn,
     404    bool preProcessing) {
    403405  const double SLIP = 20.0;
    404406  char sys = refPrn.system();
    405407  string epoTimeStr = string(_epoTime);
    406   const vector<t_pppParam*>& params = _parlist.params();
     408  const vector<t_pppParam*> &params = _parlist.params();
    407409  unsigned nPar = _parlist.nPar();
    408410
    409411  for (unsigned ii = 0; ii < LCs.size(); ii++) {
    410     const t_lc::type& tLC = LCs[ii];
     412    const t_lc::type &tLC = LCs[ii];
    411413    if (t_lc::includesPhase(tLC)) {
    412414      for (unsigned iObs = 0; iObs < obsVector.size(); iObs++) {
    413         const t_pppSatObs* obs = obsVector[iObs];
     415        const t_pppSatObs *obs = obsVector[iObs];
    414416
    415417        // Check set Slips and Jump Counters
    … …  
    423425
    424426        if (obs->slip()) {
    425           LOG << epoTimeStr  << "cycle slip set (obs) " << obs->prn().toString()  << endl;
     427          LOG << epoTimeStr << "cycle slip set (obs) " << obs->prn().toString()
     428              << endl;
    426429          slip = true;
    427430        }
    428431
    429         if (_slips[obs->prn()]._obsSlipCounter != -1 &&
    430             _slips[obs->prn()]._obsSlipCounter != obs->slipCounter()) {
    431           LOG << epoTimeStr  << "cycle slip set (obsSlipCounter) " << obs->prn().toString()  << endl;
     432        if (_slips[obs->prn()]._obsSlipCounter != -1
     433            && _slips[obs->prn()]._obsSlipCounter != obs->slipCounter()) {
     434          LOG << epoTimeStr << "cycle slip set (obsSlipCounter) "
     435              << obs->prn().toString() << endl;
    432436          slip = true;
    433437        }
    … …  
    435439          _slips[obs->prn()]._obsSlipCounter = obs->slipCounter();
    436440        }
    437         if (_slips[obs->prn()]._biasJumpCounter != -1 &&
    438             _slips[obs->prn()]._biasJumpCounter != obs->biasJumpCounter()) {
    439           LOG << epoTimeStr  << "cycle slip set (biasJumpCounter) " << obs->prn().toString() << endl;
     441        if (_slips[obs->prn()]._biasJumpCounter != -1
     442            && _slips[obs->prn()]._biasJumpCounter != obs->biasJumpCounter()) {
     443          LOG << epoTimeStr << "cycle slip set (biasJumpCounter) "
     444              << obs->prn().toString() << endl;
    440445          slip = true;
    441446        }
    … …  
    448453          if (preProcessing) {
    449454            _obsPool->setRefSatChangeRequired(sys, true);
    450           }
    451           else {
     455          } else {
    452456            resetAmb(obs->prn(), obsVector);
    453457          }
    … …  
    456460        // -----------------------
    457461        else {
    458           if (refPrn != t_prn()) {return success;}
     462          if (refPrn != t_prn()) {
     463            return success;
     464          }
    459465          ColumnVector AA(nPar);
    460466          for (unsigned iPar = 0; iPar < nPar; iPar++) {
    461             const t_pppParam* par = params[iPar];
     467            const t_pppParam *par = params[iPar];
    462468            AA[iPar] = par->partial(_epoTime, obs, tLC, refPrn);
    463469          }
    464           double ll = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA);
     470          double ll = obs->obsValue(tLC) - obs->cmpValue(tLC)
     471              - DotProduct(_x0, AA);
    465472          double vv = DotProduct(AA, _xFlt) - ll;
    466473          if (fabs(vv) > SLIP) {
    467             LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' '
    468               << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl;
     474            LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC)
     475                << ' ' << obs->prn().toString() << ' ' << setw(8)
     476                << setprecision(4) << vv << endl;
    469477            if (preProcessing) {
    470478              _obsPool->setRefSatChangeRequired(sys, true);
    471             }
    472             else {
     479            } else {
    473480              resetAmb(obs->prn(), obsVector);
    474481            }
    … …  
    483490// Reset Ambiguity Parameter (cycle slip)
    484491////////////////////////////////////////////////////////////////////////////
    485 t_irc t_pppFilter::resetAmb(t_prn prn, const vector<t_pppSatObs*>& obsVector,
    486                             SymmetricMatrix* QSav, ColumnVector* xSav) {
     492t_irc t_pppFilter::resetAmb(t_prn prn, const vector<t_pppSatObs*> &obsVector,
     493    SymmetricMatrix *QSav, ColumnVector *xSav) {
    487494
    488495  t_irc irc = failure;
    489   vector<t_pppParam*>& params = _parlist.params();
     496  vector<t_pppParam*> &params = _parlist.params();
    490497  for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    491     t_pppParam* par = params[iPar];
     498    t_pppParam *par = params[iPar];
    492499    if (par->type() == t_pppParam::amb && par->prn() == prn) {
    493500      int ind = par->indexNew();
    494501      t_lc::type tLC = par->tLC();
    495502      LOG << string(_epoTime) << " RESET " << par->toString() << endl;
    496       delete par; par = new t_pppParam(t_pppParam::amb, prn, tLC, &obsVector);
     503      delete par;
     504      par = new t_pppParam(t_pppParam::amb, prn, tLC, &obsVector);
    497505      par->setIndex(ind);
    498506      params[iPar] = par;
    499507      for (unsigned ii = 1; ii <= params.size(); ii++) {
    500         _QFlt(ii, ind+1) = 0.0;
     508        _QFlt(ii, ind + 1) = 0.0;
    501509        if (QSav) {
    502           (*QSav)(ii, ind+1) = 0.0;
    503         }
    504       }
    505       _QFlt(ind+1,ind+1) = par->sigma0() * par->sigma0();
     510          (*QSav)(ii, ind + 1) = 0.0;
     511        }
     512      }
     513      _QFlt(ind + 1, ind + 1) = par->sigma0() * par->sigma0();
    506514      if (QSav) {
    507         (*QSav)(ind+1,ind+1) = _QFlt(ind+1,ind+1);
     515        (*QSav)(ind + 1, ind + 1) = _QFlt(ind + 1, ind + 1);
    508516      }
    509517      _xFlt[ind] = 0.0;
    … …  
    521529// Add infinite noise to iono
    522530////////////////////////////////////////////////////////////////////////////
    523 t_irc t_pppFilter::addNoiseToIono(char sys) {
     531t_irc t_pppFilter::addNoiseToPar(t_pppParam::e_type parType, char sys) {
    524532
    525533  t_irc irc = failure;
    526   vector<t_pppParam*>& params = _parlist.params();
     534  vector<t_pppParam*> &params = _parlist.params();
    527535  for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    528     t_pppParam* par = params[iPar];
    529     if (par->type() == t_pppParam::ion && par->prn().system() == sys) {
     536    t_pppParam *par = params[iPar];
     537    if (par->type() == parType && par->prn().system() == sys) {
    530538      int ind = par->indexNew();
    531       LOG << string(_epoTime) << " ADD IONO_NOISE TO "  <<  par->prn().toString() << endl;
     539      LOG << string(_epoTime) << " ADD IONO_NOISE TO " << par->prn().toString()
     540          << endl;
    532541      par->setIndex(ind);
    533       _QFlt(ind+1,ind+1) = par->sigma0() * par->sigma0();
     542      _QFlt(ind + 1, ind + 1) = par->sigma0() * par->sigma0();
    534543      irc = success;
    535544    }
    … …  
    541550// Compute various DOP Values
    542551////////////////////////////////////////////////////////////////////////////
    543 void t_pppFilter::cmpDOP(const vector<t_pppSatObs*>& obsVector,
    544                          const QMap<char, t_pppRefSat*>& refSatMap) {
     552void t_pppFilter::cmpDOP(const vector<t_pppSatObs*> &obsVector,
     553    const QMap<char, t_pppRefSat*> &refSatMap) {
    545554
    546555  _dop.reset();
    … …  
    551560    _numSat = 0;
    552561    for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    553       t_pppSatObs* obs = obsVector[ii];
     562      t_pppSatObs *obs = obsVector[ii];
    554563      char sys = obs->prn().system();
    555564      t_prn refPrn = t_prn();
    … …  
    560569        ++_numSat;
    561570        for (unsigned iPar = 0; iPar < numPar; iPar++) {
    562           const t_pppParam* par = _parlist.params()[iPar];
    563           AA[_numSat-1][iPar] = par->partial(_epoTime, obs, t_lc::c1, refPrn);
     571          const t_pppParam *par = _parlist.params()[iPar];
     572          AA[_numSat - 1][iPar] = par->partial(_epoTime, obs, t_lc::c1, refPrn);
    564573        }
    565574      }
    … …  
    569578    }
    570579    AA = AA.Rows(1, _numSat);
    571     SymmetricMatrix NN; NN << AA.t() * AA;
     580    SymmetricMatrix NN;
     581    NN << AA.t() * AA;
    572582    SymmetricMatrix QQ = NN.i();
    573583
    574     _dop.H = sqrt(QQ(1,1) + QQ(2,2));
    575     _dop.V = sqrt(QQ(3,3));
    576     _dop.P = sqrt(QQ(1,1) + QQ(2,2) + QQ(3,3));
    577     _dop.T = sqrt(QQ(4,4));
    578     _dop.G = sqrt(QQ(1,1) + QQ(2,2) + QQ(3,3) + QQ(4,4));
    579   }
    580   catch (...) {
     584    _dop.H = sqrt(QQ(1, 1) + QQ(2, 2));
     585    _dop.V = sqrt(QQ(3, 3));
     586    _dop.P = sqrt(QQ(1, 1) + QQ(2, 2) + QQ(3, 3));
     587    _dop.T = sqrt(QQ(4, 4));
     588    _dop.G = sqrt(QQ(1, 1) + QQ(2, 2) + QQ(3, 3) + QQ(4, 4));
     589  } catch (...) {
    581590  }
    582591}
    … …  
    584593//
    585594////////////////////////////////////////////////////////////////////////////
    586 void t_pppFilter::predictCovCrdPart(const SymmetricMatrix& QFltOld) {
    587 
    588   const vector<t_pppParam*>& params = _parlist.params();
     595void t_pppFilter::predictCovCrdPart(const SymmetricMatrix &QFltOld) {
     596
     597  const vector<t_pppParam*> &params = _parlist.params();
    589598  if (params.size() < 3) {
    590599    return;
    … …  
    593602  bool first = (params[0]->indexOld() < 0);
    594603
    595   SymmetricMatrix Qneu(3); Qneu = 0.0;
     604  SymmetricMatrix Qneu(3);
     605  Qneu = 0.0;
    596606  for (unsigned ii = 0; ii < 3; ii++) {
    597     const t_pppParam* par = params[ii];
     607    const t_pppParam *par = params[ii];
    598608    if (first) {
    599609      Qneu[ii][ii] = par->sigma0() * par->sigma0();
    600     }
    601     else {
     610    } else {
    602611      Qneu[ii][ii] = par->noise() * par->noise();
    603612    }
    604613  }
    605614
    606   const t_pppStation* sta = PPP_CLIENT->staRover();
     615  const t_pppStation *sta = PPP_CLIENT->staRover();
    607616  SymmetricMatrix Qxyz(3);
    608617  covariNEU_XYZ(Qneu, sta->ellApr().data(), Qxyz);
    609618
    610619  if (first) {
    611     _QFlt.SymSubMatrix(1,3) = Qxyz;
    612   }
    613   else {
     620    _QFlt.SymSubMatrix(1, 3) = Qxyz;
     621  } else {
    614622    double dt = _epoTime - _firstEpoTime;
    615623    if (dt < OPT->_seedingTime) {
    616       _QFlt.SymSubMatrix(1,3) = QFltOld.SymSubMatrix(1,3);
    617     }
    618     else {
    619       _QFlt.SymSubMatrix(1,3) = QFltOld.SymSubMatrix(1,3) + Qxyz;
     624      _QFlt.SymSubMatrix(1, 3) = QFltOld.SymSubMatrix(1, 3);
     625    } else {
     626      _QFlt.SymSubMatrix(1, 3) = QFltOld.SymSubMatrix(1, 3) + Qxyz;
    620627    }
    621628  }
    … …  
    624631//
    625632////////////////////////////////////////////////////////////////////////////
    626 void t_pppFilter::setStateVectorAndVarCovMatrix(const ColumnVector& xFltOld,
    627                                                 const SymmetricMatrix& QFltOld) {
    628 
    629   const vector<t_pppParam*>& params = _parlist.params();
     633void t_pppFilter::setStateVectorAndVarCovMatrix(const ColumnVector &xFltOld,
     634    const SymmetricMatrix &QFltOld) {
     635
     636  const vector<t_pppParam*> &params = _parlist.params();
    630637  unsigned nPar = params.size();
    631638
    632   _QFlt.ReSize(nPar); _QFlt = 0.0;
    633   _xFlt.ReSize(nPar); _xFlt = 0.0;
    634   _x0.ReSize(nPar);   _x0   = 0.0;
     639  _QFlt.ReSize(nPar);
     640  _QFlt = 0.0;
     641  _xFlt.ReSize(nPar);
     642  _xFlt = 0.0;
     643  _x0.ReSize(nPar);
     644  _x0 = 0.0;
    635645
    636646  for (unsigned ii = 0; ii < nPar; ii++) {
    637     t_pppParam* par1 = params[ii];
     647    t_pppParam *par1 = params[ii];
    638648    if (QFltOld.size() == 0) {
    639649      par1->resetIndex();
    … …  
    643653    if (iOld < 0) {
    644654      _QFlt[ii][ii] = par1->sigma0() * par1->sigma0(); // new parameter
    645     }
    646     else {
     655    } else {
    647656      _QFlt[ii][ii] = QFltOld[iOld][iOld] + par1->noise() * par1->noise();
    648       _xFlt[ii]     = xFltOld[iOld];
     657      _xFlt[ii] = xFltOld[iOld];
    649658      for (unsigned jj = 0; jj < ii; jj++) {
    650         t_pppParam* par2 = params[jj];
    651         int  jOld = par2->indexOld();
     659        t_pppParam *par2 = params[jj];
     660        int jOld = par2->indexOld();
    652661        if (jOld >= 0) {
    653           _QFlt[ii][jj] = QFltOld(iOld+1,jOld+1);
     662          _QFlt[ii][jj] = QFltOld(iOld + 1, jOld + 1);
    654663        }
    655664      }
    … …  
    661670// Compute datum transformation
    662671////////////////////////////////////////////////////////////////////////////
    663 t_irc t_pppFilter::datumTransformation(const QMap<char, t_pppRefSat*>& refSatMap) {
     672t_irc t_pppFilter::datumTransformation(
     673    const QMap<char, t_pppRefSat*> &refSatMap) {
    664674
    665675  // get last epoch
    666   t_pppObsPool::t_epoch* epoch = _obsPool->lastEpoch();
     676  t_pppObsPool::t_epoch *epoch = _obsPool->lastEpoch();
    667677  if (!epoch) {
    668678    LOG << "t_pppFilter::datumTransformation: !lastEpoch" << endl;
    … …  
    673683  LOG << string(_epoTime) << " DATUM TRANSFORMATION " << endl;
    674684
    675   vector<t_pppSatObs*>& allObs = epoch->obsVector();
    676 
    677   const QMap<char, t_pppRefSat*>& refSatMapLastEpoch = epoch->refSatMap();
     685  vector<t_pppSatObs*> &allObs = epoch->obsVector();
     686
     687  const QMap<char, t_pppRefSat*> &refSatMapLastEpoch = epoch->refSatMap();
    678688
    679689  bool pseudoObsIono = epoch->pseudoObsIono();
    … …  
    681691  // reset old and set new refSats in last epoch (ambiguities/GIM)
    682692  // =============================================================
    683   if (resetRefSatellitesLastEpoch(allObs, refSatMap, refSatMapLastEpoch) != true) {
    684     LOG  << "datumTransformation: resetRefSatellitesLastEpoch != success" << endl;
     693  if (resetRefSatellitesLastEpoch(allObs, refSatMap, refSatMapLastEpoch)
     694      != true) {
     695    LOG << "datumTransformation: resetRefSatellitesLastEpoch != success"
     696        << endl;
    685697    return failure;
    686698  }
    … …  
    703715#endif
    704716
    705   const QList<char>& usedSystems = _parlist.usedSystems();
     717  const QList<char> &usedSystems = _parlist.usedSystems();
    706718  for (int iSys = 0; iSys < usedSystems.size(); iSys++) {
    707719    char sys = usedSystems[iSys];
    … …  
    720732    unsigned usedLCs = LCs.size();
    721733    if (OPT->_pseudoObsIono && !pseudoObsIono) {
    722         usedLCs -= 1;  // GIM not used
     734      usedLCs -= 1;  // GIM not used
    723735    }
    724736    // max Obs
    … …  
    729741    }
    730742
    731     const vector<t_pppParam*>& params = _parlist.params();
     743    const vector<t_pppParam*> &params = _parlist.params();
    732744    unsigned nPar = _parlist.nPar();
    733745
    734     Matrix  AA(maxObs, nPar);
     746    Matrix AA(maxObs, nPar);
    735747
    736748    // Real Observations
    … …  
    738750    int iObs = -1;
    739751    for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    740       t_pppSatObs* obs = obsVector[ii];
     752      t_pppSatObs *obs = obsVector[ii];
    741753      for (unsigned jj = 0; jj < usedLCs; jj++) {
    742754        const t_lc::type tLC = LCs[jj];
    743         if (tLC == t_lc::GIM) {continue;}
     755        if (tLC == t_lc::GIM) {
     756          continue;
     757        }
    744758        ++iObs;
    745759        for (unsigned iPar = 0; iPar < nPar; iPar++) {
    746           const t_pppParam* par = params[iPar];
     760          const t_pppParam *par = params[iPar];
    747761          AA[iObs][iPar] = par->partial(_epoTime, obs, tLC, refPrn);
    748762        }
    … …  
    753767      continue;
    754768    }
    755     _datumTrafo->updateIndices(sys, iObs+1);    //LOG << "AA Ncols/Nrows: " << AA.Ncols() << "/" << AA.Nrows() << "  nPar: "  << nPar << endl;    //LOG << "AA.SubMatrix(1 .. " << iObs+1 << " , 1 .. " <<  nPar << ")" << endl;
    756     if(_datumTrafo->prepareAA(AA.SubMatrix(1, iObs+1 , 1, nPar), 2) != success) {
     769    _datumTrafo->updateIndices(sys, iObs + 1); //LOG << "AA Ncols/Nrows: " << AA.Ncols() << "/" << AA.Nrows() << "  nPar: "  << nPar << endl;    //LOG << "AA.SubMatrix(1 .. " << iObs+1 << " , 1 .. " <<  nPar << ")" << endl;
     770    if (_datumTrafo->prepareAA(AA.SubMatrix(1, iObs + 1, 1, nPar), 2)
     771        != success) {
    757772      return failure;
    758773    }
    … …  
    766781  //LOG << "AA2\n"; _datumTrafo->printMatrix(_datumTrafo->AA2(), _datumTrafo->numObs(), _datumTrafo->numPar());
    767782#endif
    768   if(_datumTrafo->computeTrafoMatrix() != success) {
     783  if (_datumTrafo->computeTrafoMatrix() != success) {
    769784    return failure;
    770785  }
    … …  
    772787  //LOG << "D21" << endl; _datumTrafo->printMatrix(_datumTrafo->D21(), _datumTrafo->numObs(), _datumTrafo->numPar());
    773788#endif
    774   ColumnVector    xFltOld = _xFlt;
     789  ColumnVector xFltOld = _xFlt;
    775790  SymmetricMatrix QFltOld = _QFlt;
    776791
    777792  _QFlt << _datumTrafo->D21() * QFltOld * _datumTrafo->D21().t();
    778   _xFlt =  _datumTrafo->D21() * xFltOld;
     793  _xFlt = _datumTrafo->D21() * xFltOld;
    779794
    780795#ifdef BNC_DEBUG_PPP
    … …  
    790805    t_prn refPrnNew = refSatMap[sys]->prn();
    791806    if (refPrnNew != refPrnOld) {
    792       t_irc irc = resetAmb(refPrnOld, allObs);
    793       if (OPT->_obsModelType == OPT->DCMcodeBias) {if (irc == success) {addNoiseToIono(sys);}}
     807      if (resetAmb(refPrnOld, allObs) == success) {
     808        if (OPT->_obsModelType == OPT->DCMcodeBias) {
     809          addNoiseToPar(t_pppParam::ion, sys);
     810        } else if (OPT->_obsModelType == OPT->DCMphaseBias) {
     811          if (sys == 'G') {
     812            addNoiseToPar(t_pppParam::pBiasG1, sys);
     813            addNoiseToPar(t_pppParam::pBiasG2, sys);
     814          } else if (sys == 'R') {
     815            addNoiseToPar(t_pppParam::pBiasR1, sys);
     816            addNoiseToPar(t_pppParam::pBiasR2, sys);
     817          } else if (sys == 'E') {
     818            addNoiseToPar(t_pppParam::pBiasE1, sys);
     819            addNoiseToPar(t_pppParam::pBiasE2, sys);
     820          } else if (sys == 'C') {
     821            addNoiseToPar(t_pppParam::pBiasC1, sys);
     822            addNoiseToPar(t_pppParam::pBiasC2, sys);
     823          }
     824        }
     825      }
    794826    }
    795827  }
    … …  
    806838// Init datum transformation
    807839////////////////////////////////////////////////////////////////////////////
    808 void t_pppFilter::initDatumTransformation(const std::vector<t_pppSatObs*>& allObs,
    809                                           bool pseudoObsIono) {
     840void t_pppFilter::initDatumTransformation(
     841    const std::vector<t_pppSatObs*> &allObs, bool pseudoObsIono) {
    810842  unsigned trafoObs = 0;
    811   const QList<char>& usedSystems = _parlist. usedSystems();
     843  const QList<char> &usedSystems = _parlist.usedSystems();
    812844  for (int iSys = 0; iSys < usedSystems.size(); iSys++) {
    813845    char sys = usedSystems[iSys];
    … …  
    838870//
    839871//////////////////////////////////////////////////////////////////////////////
    840 bool t_pppFilter::resetRefSatellitesLastEpoch(std::vector<t_pppSatObs*>& obsVector,
    841     const QMap<char, t_pppRefSat*>& refSatMap,
    842     const QMap<char, t_pppRefSat*>& refSatMapLastEpoch) {
     872bool t_pppFilter::resetRefSatellitesLastEpoch(
     873    std::vector<t_pppSatObs*> &obsVector,
     874    const QMap<char, t_pppRefSat*> &refSatMap,
     875    const QMap<char, t_pppRefSat*> &refSatMapLastEpoch) {
    843876  bool resetRefSat;
    844877  // reference satellite definition per system
    845   const QList<char>& usedSystems = _parlist.usedSystems();
     878  const QList<char> &usedSystems = _parlist.usedSystems();
    846879  for (int iSys = 0; iSys < usedSystems.size(); iSys++) {
    847880    resetRefSat = false;
    … …  
    856889    vector<t_pppSatObs*>::iterator it = obsVector.begin();
    857890    while (it != obsVector.end()) {
    858       t_pppSatObs* satObs = *it;
    859       if      (satObs->prn() == newPrn) {
     891      t_pppSatObs *satObs = *it;
     892      if (satObs->prn() == newPrn) {
    860893        resetRefSat = true;
    861894        satObs->setAsReference();
    … …  
    863896        satObs->resetReference();
    864897      }
    865      it++;
     898      it++;
    866899    }
    867900    if (!resetRefSat) {

  • trunk/BNC/src/PPP/pppFilter.h

    r9552 r9642  
    195195  void predictCovCrdPart(const SymmetricMatrix& QFltOld);
    196196
    197   t_irc addNoiseToIono(char sys);
     197  t_irc addNoiseToPar(t_pppParam::e_type parType, char sys);
    198198
    199199  bool resetRefSatellitesLastEpoch(std::vector<t_pppSatObs*>& obsVector,

  • trunk/BNC/src/PPP/pppParlist.cpp

    r9631 r9642  
    102102      _noise   = OPT->_noiseIon;
    103103     break;
    104    case cBiasG1:   case cBiasR1:   case cBiasE1:   case cBiasC1:
    105    case cBiasG2:   case cBiasR2:   case cBiasE2:   case cBiasC2:
     104   case cBiasG1:   case cBiasE1:   case cBiasC1:
     105   case cBiasG2:   case cBiasE2:   case cBiasC2:
    106106     _epoSpec = false;
    107107     _sigma0  = OPT->_aprSigCodeBias;
    108108     _noise   = OPT->_noiseCodeBias;
    109109     break;
    110    case pBiasG1:   case pBiasR1:   case pBiasE1:   case pBiasC1:
    111    case pBiasG2:   case pBiasR2:   case pBiasE2:   case pBiasC2:
     110   case cBiasR1:
     111   case cBiasR2:
     112     _epoSpec = false;
     113     _sigma0  = OPT->_aprSigCodeBias;
     114     _noise   = 10.0 * OPT->_noiseCodeBias;
     115     break;
     116   case pBiasG1:   case pBiasE1:   case pBiasC1:
     117   case pBiasG2:   case pBiasE2:   case pBiasC2:
    112118     _epoSpec = false;
    113119     _sigma0  = OPT->_aprSigPhaseBias;
    114120     _noise   = OPT->_noisePhaseBias;
     121     break;
     122   case pBiasR1:
     123   case pBiasR2:
     124     _epoSpec = false;
     125     _sigma0  = OPT->_aprSigPhaseBias;
     126     _noise   = 10.0 * OPT->_noisePhaseBias;
    115127     break;
    116128  }
Note: See TracChangeset for help on using the changeset viewer.