Changeset 9419 in ntrip for trunk/BNC/src/PPP/pppFilter.cpp


Ignore:
Timestamp:
May 3, 2021, 2:18:39 PM (3 years ago)
Author:
stuerze
Message:

update regarding PPP

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r9398 r9419  
    8181  // Set Parameters
    8282  // --------------
    83   _parlist->set(_epoTime, allObs, _obsPool->getRefSatMap());
     83  if (_parlist->set(_epoTime, allObs, _obsPool->getRefSatMap()) != success) {
     84    return failure;
     85  }
    8486  const vector<t_pppParam*>& params = _parlist->params();
    8587#ifdef BNC_DEBUG_PPP
    … …  
    8890  }
    8991#endif
     92
    9093  // Status Vector, Variance-Covariance Matrix
    9194  // -----------------------------------------
    … …  
    127130      OPT->_obsModelType == OPT->DCMphaseBias) {
    128131    preProcessing = true;
     132    unsigned usableSys = 0;
    129133    for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
    130134      char sys = OPT->systems()[iSys];
    … …  
    138142        }
    139143      }
    140       if (!obsVector.size()) {continue;}
     144      if (!obsVector.size()) {
     145        continue;
     146      }
     147      else {
     148        ++usableSys;
     149        if (usableSys == 1) {
     150          _datumTrafo->setFirstSystem(sys);
     151        }
     152      }
    141153      if (processSystem(OPT->LCs(sys), obsVector, _refPrn,
    142154                        epoch->pseudoObsIono(), preProcessing) != success) {
    … …  
    154166    }
    155167    else if (!_obsPool->refSatChangeRequired()) {
    156       initDatumTransformation(allObs);
     168      initDatumTransformation(allObs, epoch->pseudoObsIono());
    157169    }
    158170  }
    … …  
    162174  preProcessing = false;
    163175  for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
    164     (iSys) ? _datumTrafo->setFirstSystem(false) : _datumTrafo->setFirstSystem(true);
    165     char system = OPT->systems()[iSys];
     176    char sys = OPT->systems()[iSys];
    166177    if (OPT->_refSatRequired) {
    167       _refPrn = (_obsPool->getRefSatMapElement(system))->prn();
     178      _refPrn = (_obsPool->getRefSatMapElement(sys))->prn();
    168179    }
    169180    unsigned int num = 0;
    170181    vector<t_pppSatObs*> obsVector;
    171182    for (unsigned jj = 0; jj < allObs.size(); jj++) {
    172       if (allObs[jj]->prn().system() == system) {
     183      if (allObs[jj]->prn().system() == sys) {
    173184        obsVector.push_back(allObs[jj]);
    174         num++;
    175       }
    176     }
    177     if (!num) {continue;}
    178     LOG << epoTimeStr << " SATNUM " << system << ' ' << right << setw(2) << num << endl;
    179     if (processSystem(OPT->LCs(system), obsVector, _refPrn,
     185        ++num;
     186      }
     187    }
     188    if (!num) {
     189      continue;
     190    }
     191    LOG << epoTimeStr << " SATNUM " << sys << ' ' << right << setw(2) << num << endl;
     192    if (processSystem(OPT->LCs(sys), obsVector, _refPrn,
    180193                      epoch->pseudoObsIono(), preProcessing) != success) {
    181194      return failure;
    … …  
    190203  if (OPT->_refSatRequired) {
    191204    _obsPool->saveLastEpoRefSats();
     205    _datumTrafo->setLastEpoParlist(_parlist);
    192206  }
    193207  return success;
    … …  
    219233
    220234  unsigned usedLCs     = LCs.size();
    221   unsigned realUsedLCs = usedLCs;
    222235  if (OPT->_pseudoObsIono && !pseudoObsIonoAvailable) {
    223236      usedLCs -= 1;  // GIM not used
    … …  
    226239  if (OPT->_pseudoObsTropo) {
    227240    hlpLCs = -1;
    228     realUsedLCs -= 1;
    229241  }
    230242  // max Obs
    231243  unsigned maxObs = obsVector.size() * (usedLCs + hlpLCs);
    232   if (OPT->_pseudoObsTropo) {
     244  if (OPT->_pseudoObsTropo && _datumTrafo->firstSystem(sys)) {
    233245    maxObs += 1;
    234246  }
    … …  
    260272    for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    261273      t_pppSatObs* obs = obsVector[ii];
     274      if (iOutlier == 0 && !preProcessing) {
     275        obs->resetOutlier();
     276      }
    262277      if (!obs->outlier()) {
    263278        for (unsigned jj = 0; jj < usedLCs; jj++) {
    … …  
    280295    // pseudo Obs Tropo
    281296    // ================
    282     if (OPT->_pseudoObsTropo) {
     297    if (OPT->_pseudoObsTropo && _datumTrafo->firstSystem(sys)) {
    283298      for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    284299        t_pppSatObs* obs = obsVector[ii];
    … …  
    302317    if ((!iOutlier) &&
    303318        (OPT->_obsModelType == OPT->DCMcodeBias ||
    304          OPT->_obsModelType == OPT->DCMphaseBias  ) &&  (!preProcessing)) {
    305       _datumTrafo->updateIndices(iObs+1);
     319         OPT->_obsModelType == OPT->DCMphaseBias) &&
     320        (!preProcessing)) {
     321      _datumTrafo->updateIndices(sys, iObs+1);
    306322      _datumTrafo->prepareAA(AA.SubMatrix(1, iObs+1 , 1, _parlist->nPar()), 1);
    307323    }
    … …  
    499515        // -----------------------
    500516        else {
     517          if (refPrn != t_prn()) {return success;}
    501518          ColumnVector AA(params.size());
    502519          for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    … …  
    666683////////////////////////////////////////////////////////////////////////////
    667684t_irc t_pppFilter::datumTransformation() {
     685  // get last epoch
    668686  t_pppObsPool::t_epoch* epoch = _obsPool->lastEpoch();
    669   if (!epoch) {LOG << "!epoch" << endl;
     687  if (!epoch) {
     688    LOG << "!lastEpoch" << endl;
    670689    return failure;
    671690  }
    … …  
    674693    LOG << string(epoch->epoTime()) << " DATUM TRANSFORMATION " << endl;
    675694  }
     695
    676696  vector<t_pppSatObs*>& allObs = epoch->obsVector();
    677697
    678   // reset old and set new refSats in last epoch (ambiguities)
    679   // ========================================================
     698  // reset old and set new refSats in last epoch (ambiguities/GIM)
     699  // =============================================================
    680700  if (resetRefSatellitesLastEpoch(allObs) != true) {
     701    LOG  << "resetRefSatellitesLastEpoch = failure" << endl;
    681702    return failure;
     703  }
     704
     705  if (OPT->_obsModelType == OPT->UncombPPP) {
     706    return success;
    682707  }
    683708
    684709  // set AA2
    685710  // =======
    686   _parlist->set(epoch->epoTime(), allObs, _obsPool->getRefSatMap());
    687   const vector<t_pppParam*>& params = _parlist->params();
     711  t_pppParlist* parlist = _datumTrafo->lastEpoParlist();
     712  if (parlist->set(epoch->epoTime(), allObs, _obsPool->getRefSatMap()) != success) {
     713    return failure;
     714  }
     715  vector<t_pppParam*>& params = parlist->params();
    688716#ifdef BNC_DEBUG_PPP
     717  LOG << " parameters of last epoch" << endl;
    689718  for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    690719    LOG << params[iPar]->toString() << "\t\t" << endl;
    691720  }
    692721#endif
     722  unsigned nPar = parlist->nPar();
     723  unsigned usableSys = 0;
    693724  for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
    694     (iSys) ? _datumTrafo->setFirstSystem(false) : _datumTrafo->setFirstSystem(true);
    695725    char sys = OPT->systems()[iSys];
    696726    t_prn refPrn = (_obsPool->getRefSatMapElement(sys))->prn();
    697 #ifdef BNC_DEBUG_PPP
    698     LOG << "refPrn: " << refPrn.toString() << endl;
    699 #endif
    700727    vector<t_pppSatObs*> obsVector;
    701728    for (unsigned jj = 0; jj < allObs.size(); jj++) {
    … …  
    704731      }
    705732    }
    706 
     733    if (!obsVector.size()) {
     734      continue;
     735    }
     736    else {
     737      ++usableSys;
     738      if (usableSys == 1) {
     739        _datumTrafo->setFirstSystem(sys);
     740      }
     741    }
    707742    vector<t_lc::type> LCs = OPT->LCs(sys);
    708743    unsigned usedLCs = LCs.size();
    709     unsigned realUsedLCs = usedLCs;
    710744    if (OPT->_pseudoObsIono && !epoch->pseudoObsIono()) {
    711745        usedLCs -= 1;  // GIM not used
    … …  
    714748    if (OPT->_pseudoObsTropo) {
    715749      hlpLCs = -1;
    716       realUsedLCs -= 1;
    717750    }
    718751    // max Obs
    719752    unsigned maxObs = obsVector.size() * (usedLCs + hlpLCs);
    720     if (OPT->_pseudoObsTropo) {
     753
     754    if (OPT->_pseudoObsTropo && _datumTrafo->firstSystem(sys)) {
    721755      maxObs += 1;
    722756    }
    … …  
    724758      maxObs -= 1; // pseudo obs iono with respect to refSat
    725759    }
    726 
    727     Matrix  AA(maxObs, _parlist->nPar());
     760    Matrix  AA(maxObs, nPar);
    728761
    729762    // Real Observations
    … …  
    732765    for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    733766      t_pppSatObs* obs = obsVector[ii];
    734       if (!obs->outlier()) {
    735         for (unsigned jj = 0; jj < usedLCs; jj++) {
    736           const t_lc::type tLC = LCs[jj];
    737           if (tLC == t_lc::GIM) {continue;}
    738           if (tLC == t_lc::Tz0) {continue;}
    739           ++iObs;
    740           for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    741             const t_pppParam* par = params[iPar];
    742             AA[iObs][iPar] = par->partial(_epoTime, obs, tLC, refPrn);
    743           }
     767      for (unsigned jj = 0; jj < usedLCs; jj++) {
     768        const t_lc::type tLC = LCs[jj];
     769        if (tLC == t_lc::GIM) {continue;}
     770        if (tLC == t_lc::Tz0) {continue;}
     771        ++iObs;
     772        for (unsigned iPar = 0; iPar < params.size(); iPar++) {
     773          const t_pppParam* par = params[iPar];
     774          AA[iObs][iPar] = par->partial(_epoTime, obs, tLC, refPrn);
    744775        }
    745776      }
    … …  
    747778    // pseudo Obs Tropo
    748779    // ================
    749     if (OPT->_pseudoObsTropo) {
     780    if (OPT->_pseudoObsTropo && _datumTrafo->firstSystem(sys)) {
    750781      for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    751782        t_pppSatObs* obs = obsVector[ii];
    … …  
    762793      }
    763794    }
    764     _datumTrafo->updateIndices(iObs+1);
    765     _datumTrafo->prepareAA(AA.SubMatrix(1, iObs+1 , 1, _parlist->nPar()), 2);
    766   }
     795    if (!iObs) {
     796      continue;
     797    }
     798    _datumTrafo->updateIndices(sys, iObs+1);
     799    _datumTrafo->prepareAA(AA.SubMatrix(1, iObs+1 , 1, nPar), 2);
     800  }
     801  _datumTrafo->updateNumObs();
    767802
    768803  // Datum Transformation
    769804  // ====================
    770805#ifdef BNC_DEBUG_PPP
    771       LOG << "AA1\n"; _datumTrafo->printMatrix(_datumTrafo->AA1(), _datumTrafo->obsNum(), _datumTrafo->parNum());
    772       LOG << "AA2\n"; _datumTrafo->printMatrix(_datumTrafo->AA2(), _datumTrafo->obsNum(), _datumTrafo->parNum());
     806      //LOG << "AA1\n"; _datumTrafo->printMatrix(_datumTrafo->AA1(), _datumTrafo->numObs(), _datumTrafo->numPar());
     807      //LOG << "AA2\n"; _datumTrafo->printMatrix(_datumTrafo->AA2(), _datumTrafo->numObs(), _datumTrafo->numPar());
    773808#endif
    774   Matrix D21 = _datumTrafo->computeTrafoMatrix();
     809  if(_datumTrafo->computeTrafoMatrix() != success) {
     810    return failure;
     811  }
    775812#ifdef BNC_DEBUG_PPP
    776       LOG << "D21" << endl; _datumTrafo->printMatrix(D21, _datumTrafo->parNum(), _datumTrafo->parNum());
     813      //LOG << "D21" << endl; _datumTrafo->printMatrix(_datumTrafo->D21(), _datumTrafo->numObs(), _datumTrafo->numPar());
    777814#endif
    778815  ColumnVector    xFltOld = _xFlt;
    779816  SymmetricMatrix QFltOld = _QFlt;
    780817
    781   _QFlt << D21 * QFltOld * D21.t();
    782   _xFlt =  D21 * xFltOld;
     818  _QFlt << _datumTrafo->D21() * QFltOld * _datumTrafo->D21().t();
     819  _xFlt =  _datumTrafo->D21() * xFltOld;
    783820
    784821#ifdef BNC_DEBUG_PPP
    … …  
    791828  for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
    792829    char sys = OPT->systems()[iSys];
    793     t_irc irc = resetAmb(_obsPool->getRefSatMapElementLastEpoch(sys), allObs);
    794     if (OPT->_obsModelType == OPT->DCMcodeBias) {
    795       if (irc == success) {
    796         addNoiseToIono(sys);}
     830    t_prn refPrnOld = _obsPool->getRefSatMapElementLastEpoch(sys);
     831    t_prn refPrnNew = (_obsPool->getRefSatMapElement(sys))->prn();
     832    if (refPrnNew != refPrnOld) {
     833      t_irc irc = resetAmb(_obsPool->getRefSatMapElementLastEpoch(sys), allObs);
     834      if (OPT->_obsModelType == OPT->DCMcodeBias) {
     835        if (irc == success) {
     836          addNoiseToIono(sys);
     837        }
     838      }
    797839    }
    798840  }
    … …  
    802844  _datumTrafo->switchAA();
    803845
     846  // save parameter list
     847  // ====================
     848  _datumTrafo->setLastEpoParlist(_parlist);
    804849  return success;
    805850}
    … …  
    807852// Init datum transformation
    808853////////////////////////////////////////////////////////////////////////////
    809 void t_pppFilter::initDatumTransformation(const std::vector<t_pppSatObs*>& allObs) {
     854void t_pppFilter::initDatumTransformation(const std::vector<t_pppSatObs*>& allObs,
     855                                          bool pseudoObsIono) {
    810856  unsigned trafoObs = 0;
    811857  for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
    812     char system = OPT->systems()[iSys];
     858    char sys = OPT->systems()[iSys];
    813859    int satNum = 0;
    814860    for (unsigned jj = 0; jj < allObs.size(); jj++) {
    815       if (allObs[jj]->prn().system() == system) {
     861      if (allObs[jj]->prn().system() == sys) {
    816862        satNum++;
    817863      }
    818864    }
     865    if (!satNum) {
     866      continue;
     867    }
    819868    // all LCs
    820     unsigned realUsedLCs = OPT->LCs(system).size();
     869    unsigned realUsedLCs = OPT->LCs(sys).size();
    821870    // exclude pseudo obs GIM
    822     if (OPT->_pseudoObsIono) {
     871    if (OPT->_pseudoObsIono && !pseudoObsIono) {
    823872      realUsedLCs -= 1;
    824873    }
    … …  
    828877    trafoObs += satNum * realUsedLCs;
    829878
    830     if (OPT->_pseudoObsTropo) {
     879    if (OPT->_pseudoObsTropo && _datumTrafo->firstSystem(sys)) {
    831880      trafoObs += 1;
    832881    }
    833 
    834   }
    835   _datumTrafo->setObsNum(trafoObs);
    836   _datumTrafo->setParNum(_parlist->nPar());
     882    if (OPT->_pseudoObsIono && pseudoObsIono) {
     883      trafoObs -= 1; // pseudo obs iono with respect to refSat
     884    }
     885  }
     886  _datumTrafo->setNumObs(trafoObs);
     887  _datumTrafo->setNumPar(_parlist->nPar());
    837888  _datumTrafo->initAA();
    838889}
Note: See TracChangeset for help on using the changeset viewer.