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


Ignore:
Timestamp:
Mar 25, 2021, 3:17:35 PM (3 years ago)
Author:
stuerze
Message:

update regarding PPP

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r9088 r9386  
    2727#include "pppObsPool.h"
    2828#include "pppStation.h"
    29 #include "pppClient.h"
    3029
    3130using namespace BNC_PPP;
    … …  
    5150// Process Single Epoch
    5251////////////////////////////////////////////////////////////////////////////
    53 t_irc t_pppFilter::processEpoch(int num) {
     52t_irc t_pppFilter::processEpoch() {
    5453  _numSat     = 0;
    55   _numEpoProcessing = num;
    5654  const double maxSolGap = 60.0;
    5755
    … …  
    8583  _parlist->set(_epoTime, allObs, _obsPool->getRefSatMap());
    8684  const vector<t_pppParam*>& params = _parlist->params();
    87 
     85#ifdef BNC_DEBUG_PPP
     86  for (unsigned iPar = 0; iPar < params.size(); iPar++) {
     87    LOG << params[iPar]->toString() << endl;
     88  }
     89#endif
    8890  // Status Vector, Variance-Covariance Matrix
    8991  // -----------------------------------------
    … …  
    9496  _xFlt.ReSize(_parlist->nPar()); _xFlt = 0.0;
    9597  _x0.ReSize(_parlist->nPar());   _x0   = 0.0;
     98
    9699  for (unsigned ii = 0; ii < params.size(); ii++) {
    97     const t_pppParam* par1 = params[ii];
     100    t_pppParam* par1 = params[ii];
     101    if (QFltOld.size() == 0) {
     102      par1->resetIndex();
     103    }
    98104    _x0[ii] = par1->x0();
    99105    int iOld = par1->indexOld();
    … …  
    105111      _xFlt[ii]     = xFltOld[iOld];
    106112      for (unsigned jj = 0; jj < ii; jj++) {
    107         const t_pppParam* par2 = params[jj];
    108         int               jOld = par2->indexOld();
     113        t_pppParam* par2 = params[jj];
     114        int  jOld = par2->indexOld();
    109115        if (jOld >= 0) {
    110116          _QFlt[ii][jj] = QFltOld(iOld+1,jOld+1);
    … …  
    114120  }
    115121  predictCovCrdPart(QFltOld);
    116 
    117   // Init Datum Trafo
    118   // ----------------------------------------
    119   if ((OPT->_obsModelType == OPT->DCMcodeBias   ||
    120        OPT->_obsModelType == OPT->DCMphaseBias) &&(_numEpoProcessing == 1))  {
    121     _numAllUsedLCs = 0;
    122     for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
    123       char system = OPT->systems()[iSys];
    124       _numAllUsedLCs += OPT->LCs(system).size();
    125       if (OPT->_pseudoObsIono && epoch->pseudoObsIono() == false) {
    126         _numAllUsedLCs -= 1;  // GIM not used
    127         if (OPT->_pseudoObsTropo) { // if no VTEC values /pseudo obs iono then no pseudo obs tropo
    128           _numAllUsedLCs -= 1;
    129         }
    130       }
    131       int modifyLCs = 0;
    132       if (OPT->_pseudoObsTropo &&
    133           OPT->_pseudoObsIono && epoch->pseudoObsIono() == true) {
    134         modifyLCs  = -1;
    135       }
    136       // max Obs
    137       int maxObs = allObs.size() * (_numAllUsedLCs + modifyLCs);
    138       if (OPT->_pseudoObsIono && epoch->pseudoObsIono() == true) {
    139         maxObs -= 1; // pseudo obs iono with respect to refSat
    140       }
    141       if (OPT->_pseudoObsIono && epoch->pseudoObsIono() == true &&
    142           OPT->_pseudoObsTropo) {
    143         maxObs += 1;   // pseudo obs tropo once per station
    144       }
    145       _datumTrafo->initAA(maxObs, _parlist->nPar());
    146     }
    147   }
    148122
    149123  // Pre-Process Satellite Systems separately
    … …  
    154128    preProcessing = true;
    155129    for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
    156       (iSys) ? _datumTrafo->setFirstSystem(false) : _datumTrafo->setFirstSystem(true);
    157       char system = OPT->systems()[iSys];
     130      char sys = OPT->systems()[iSys];
    158131      if (OPT->_refSatRequired) {
    159         _refPrn = (_obsPool->getRefSatMapElement(system))->prn();
    160         if (_obsPool->hasHistoricalRefSat(_refPrn)) {
    161           LOG << epoTimeStr << " Warning: prevent to process erroneous refSat again!";
    162           return failure;
    163         }
    164       }
    165       vector<t_pppSatObs*> obsVector;
     132        _refPrn = (_obsPool->getRefSatMapElement(sys))->prn();
     133      }
     134     vector<t_pppSatObs*> obsVector;
    166135      for (unsigned jj = 0; jj < allObs.size(); jj++) {
    167         if (allObs[jj]->prn().system() == system) {
     136        if (allObs[jj]->prn().system() == sys) {
    168137          obsVector.push_back(allObs[jj]);
    169138        }
    170139      }
    171       if (processSystem(OPT->LCs(system), obsVector, _refPrn,
     140      if (processSystem(OPT->LCs(sys), obsVector, _refPrn,
    172141                        epoch->pseudoObsIono(), preProcessing) != success) {
     142        _xFlt = xFltOld;
     143        _QFlt = QFltOld;
    173144        return failure;
    174145      }
    … …  
    181152      return success;
    182153    }
    183     else if (!_obsPool->refSatChangeRequired() &&
    184              _numEpoProcessing > 1) {
    185       datumTransformation(xFltOld, QFltOld);
     154    else if (!_obsPool->refSatChangeRequired()) {
     155      initDatumTransformation(allObs);
    186156    }
    187157  }
    … …  
    191161  preProcessing = false;
    192162  for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
     163    (iSys) ? _datumTrafo->setFirstSystem(false) : _datumTrafo->setFirstSystem(true);
    193164    char system = OPT->systems()[iSys];
    194165    if (OPT->_refSatRequired) {
    … …  
    215186  _parlist->printResult(_epoTime, _QFlt, _xFlt);
    216187  _lastEpoTimeOK = _epoTime;  // remember time of last successful epoch processing
    217 
     188  if (OPT->_refSatRequired) {
     189    _obsPool->saveLastEpoRefSats();
     190  }
    218191  return success;
    219192}
    … …  
    227200                                 bool preProcessing) {
    228201  LOG.setf(ios::fixed);
     202  char sys = refPrn.system();
    229203
    230204  // Detect Cycle Slips
    … …  
    233207    return failure;
    234208  }
    235 
    236   unsigned usedLCs = LCs.size();
     209  if (preProcessing && _obsPool->refSatChangeRequired(sys)) { // from detectCycleSlips()
     210#ifdef BNC_DEBUG_PPP
     211    LOG << "_obsPool->refSatChangeRequired(" << sys << ") from detectCycleSlips()" << endl;
     212#endif
     213    return success;
     214  }
    237215
    238216  ColumnVector               xSav       = _xFlt;
    … …  
    241219  const vector<t_pppParam*>& params     = _parlist->params();
    242220
     221  unsigned usedLCs     = LCs.size();
     222  unsigned realUsedLCs = usedLCs;
    243223  if (OPT->_pseudoObsIono && !pseudoObsIonoAvailable) {
    244     usedLCs -= 1;  // GIM not used
    245     if (OPT->_pseudoObsTropo) { // if no VTEC values /pseudo obs iono then no pseudo obs tropo
    246       usedLCs -= 1;
    247     }
    248   }
    249 
    250   int modifyLCs = 0;
    251   if (OPT->_pseudoObsTropo &&
    252       OPT->_pseudoObsIono && pseudoObsIonoAvailable) {
    253     modifyLCs  = -1;
     224      usedLCs -= 1;  // GIM not used
     225  }
     226  int hlpLCs = 0;
     227  if (OPT->_pseudoObsTropo) {
     228    hlpLCs = -1;
     229    realUsedLCs -= 1;
    254230  }
    255231  // max Obs
    256   unsigned maxObs = obsVector.size() * (usedLCs + modifyLCs);
     232  unsigned maxObs = obsVector.size() * (usedLCs + hlpLCs);
     233  if (OPT->_pseudoObsTropo) {
     234    maxObs += 1;
     235  }
    257236  if (OPT->_pseudoObsIono && pseudoObsIonoAvailable) {
    258237    maxObs -= 1; // pseudo obs iono with respect to refSat
    259   }
    260   if (OPT->_pseudoObsIono && pseudoObsIonoAvailable &&
    261       OPT->_pseudoObsTropo) {
    262     maxObs += 1;   // pseudo obs tropo once per station
    263238  }
    264239
    … …  
    281256    vector<t_pppSatObs*> usedObs;
    282257    vector<t_lc::type>   usedTypes;
     258
     259    // Real Observations
     260    // =================
    283261    for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    284262      t_pppSatObs* obs = obsVector[ii];
    285       unsigned hlp = ii+1;
    286263      if (!obs->outlier()) {
    287264        for (unsigned jj = 0; jj < usedLCs; jj++) {
    288265          const t_lc::type tLC = LCs[jj];
    289           if (tLC == t_lc::GIM &&  obs->isReference()) {continue;}
    290           if (tLC == t_lc::Tz0 && hlp != obsVector.size()) {continue;}
     266          if (tLC == t_lc::GIM) {continue;}
     267          if (tLC == t_lc::Tz0) {continue;}
    291268          ++iObs;
    292269          usedObs.push_back(obs);
    … …  
    302279    }
    303280
    304     if ((preProcessing && ! iOutlier) &&
    305         (OPT->_obsModelType == OPT->DCMcodeBias || OPT->_obsModelType == OPT->DCMphaseBias)) {
    306           _datumTrafo->updateIndices(maxObs);
    307       _datumTrafo->prepareAA(AA, _numEpoProcessing, _parlist->nPar());
    308       if (_obsPool->refSatChangeRequired()) { // from detectCycleSlips()
    309         return success;
     281    if ((!iOutlier) &&
     282        (OPT->_obsModelType == OPT->DCMcodeBias ||
     283         OPT->_obsModelType == OPT->DCMphaseBias  ) &&  (!preProcessing)) {
     284      _datumTrafo->updateIndices(iObs+1);
     285      _datumTrafo->prepareAA(AA.SubMatrix(1, iObs+1 , 1, _parlist->nPar()), 1);
     286    }
     287
     288    // Pseudo Obs Iono
     289    // ================
     290    if (OPT->_pseudoObsIono && pseudoObsIonoAvailable) {
     291      for (unsigned ii = 0; ii < obsVector.size(); ii++) {
     292        t_pppSatObs* obs = obsVector[ii];
     293        if (!obs->outlier()) {
     294          for (unsigned jj = 0; jj < usedLCs; jj++) {
     295            const t_lc::type tLC = LCs[jj];
     296            if (tLC == t_lc::GIM && !obs->isReference()) {
     297              ++iObs;
     298            } else {continue;}
     299            usedObs.push_back(obs);
     300            usedTypes.push_back(tLC);
     301            for (unsigned iPar = 0; iPar < params.size(); iPar++) {
     302              const t_pppParam* par = params[iPar];
     303              AA[iObs][iPar] = par->partial(_epoTime, obs, tLC, refPrn);
     304            }
     305            ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1));
     306            PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC));
     307          }
     308        }
     309      }
     310    }
     311    // pseudo Obs Tropo
     312    // ================
     313    if (OPT->_pseudoObsTropo) {
     314      for (unsigned ii = 0; ii < obsVector.size(); ii++) {
     315        t_pppSatObs* obs = obsVector[ii];
     316        if (!obs->isReference()) {continue;}
     317        for (unsigned jj = 0; jj < usedLCs; jj++) {
     318          const t_lc::type tLC = LCs[jj];
     319          if (tLC != t_lc::Tz0) {continue;}
     320          ++iObs;
     321          usedObs.push_back(obs);
     322          usedTypes.push_back(tLC);
     323          for (unsigned iPar = 0; iPar < params.size(); iPar++) {
     324            const t_pppParam* par = params[iPar];
     325            AA[iObs][iPar] = par->partial(_epoTime, obs, tLC, refPrn);
     326          }
     327          ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1));
     328          PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC));
     329        }
    310330      }
    311331    }
    … …  
    347367      t_pppSatObs* obs = usedObs[maxOutlierIndex];
    348368      t_pppParam* par = 0;
    349       if (!preProcessing) {
    350         LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' '
    351             << obs->prn().toString()                        << ' '
    352             << setw(8) << setprecision(4) << maxOutlier << endl;
    353       }
     369#ifdef BNC_DEBUG_PPP
     370      LOG << epoTimeStr << " Outlier ("
     371          << ((preProcessing) ? "pre-processing) " : "fin-processing) ") << t_lc::toString(maxOutlierLC) << ' '
     372          << obs->prn().toString()                        << ' '
     373          << setw(8) << setprecision(4) << maxOutlier << endl;
     374#endif
    354375      for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    355376        t_pppParam* hlp = params[iPar];
    … …  
    361382      }
    362383      if (preProcessing) {
    363         if      (par && obs->prn() == refPrn) {
    364           _obsPool->setRefSatChangeRequired(true);
    365           break;
    366         }
    367         else if (par && obs->prn() != refPrn) {
    368           par->setAmbResetCandidate();
    369         }
    370       }
    371       else {
     384        // for refSats  no ambiguity parameter exists
     385         if ((obs->prn() == refPrn) &&
     386             (t_lc::toString(maxOutlierLC) == "l1" ||
     387              t_lc::toString(maxOutlierLC) == "l2") ) {
     388           _obsPool->setRefSatChangeRequired(sys, true);
     389         }
     390         else {
     391          if (par) {
     392            par->setAmbResetCandidate();
     393            obs->setOutlier();
     394          }
     395          else {
     396            obs->setOutlier();
     397          }
     398        }
     399      }
     400      else {// fin-processing
    372401        if (par) {
    373           if (par->ambResetCandidate() || _obsPool->hasHistoricalRefSat(par->prn())) {
     402          if (par->ambResetCandidate()) {
    374403            resetAmb(par->prn(), obsVector, &QSav, &xSav);
    375404          }
    … …  
    391420          for (unsigned ii = 0; ii < usedObs.size(); ii++) {
    392421            const t_lc::type tLC = usedTypes[ii];
    393             t_pppSatObs*        obs = usedObs[ii];
     422            t_pppSatObs* obs = usedObs[ii];
    394423            if (tLC == LCs[jj]) {
    395424              obs->setRes(tLC, vv[ii]);
    396425              LOG << epoTimeStr << " RES "
    397                   << left << setw(3) << t_lc::toString(tLC) << right << ' '
    398                   << obs->prn().toString().substr(0,3) << ' '
    399                   << setw(8) << setprecision(4) << vv[ii] << endl;
     426                  << left << setw(3) << t_lc::toString(tLC) << right << ' ';
     427              if (t_lc::toString(tLC) == "Tz0") {LOG << sys << "  ";}
     428              else {LOG << obs->prn().toString();}
     429              LOG << setw(8) << setprecision(4) << vv[ii] << endl;
    400430            }
    401431          }
    … …  
    405435    }
    406436  }
    407 
    408437  return success;
    409438}
    … …  
    415444                                    const t_prn& refPrn,
    416445                                    bool preProcessing) {
    417   if (_obsPool->hasHistoricalRefSat(refPrn)) {
    418     return success;
    419   }
    420 
    421   //double                  SLIP       = 20.0;  // slip threshold
    422   double                  SLIP       = 200.0;  // slip threshold deactivated
    423   string                  epoTimeStr = string(_epoTime);
     446
     447  const double SLIP = 20.0;
     448  char sys = refPrn.system();
     449  string epoTimeStr = string(_epoTime);
    424450  const vector<t_pppParam*>& params  = _parlist->params();
    425451
    … …  
    430456        const t_pppSatObs* obs = obsVector[iObs];
    431457
    432         if (preProcessing && obs->prn() != refPrn) {continue;}
    433 
    434458        // Check set Slips and Jump Counters
    435459        // ---------------------------------
    436460        bool slip = false;
     461
     462        // in pre-processing only the reference satellite has to be checked
     463        if (preProcessing && obs->prn() != refPrn) {
     464          continue;
     465        }
     466
    437467        if (obs->slip()) {
    438468          LOG << epoTimeStr  << "cycle slip set (obs) " << obs->prn().toString()  << endl;
    … …  
    460490        if (slip) {
    461491          if (preProcessing) {
    462             _obsPool->setRefSatChangeRequired(true);
     492            _obsPool->setRefSatChangeRequired(sys, true);
    463493          }
    464494          else {
    … …  
    474504            AA[iPar] = par->partial(_epoTime, obs, tLC, refPrn);
    475505          }
    476 
    477506          double ll = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA);
    478507          double vv = DotProduct(AA, _xFlt) - ll;
    479 
    480508          if (fabs(vv) > SLIP) {
     509            LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' '
     510              << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl;
    481511            if (preProcessing) {
    482               _obsPool->setRefSatChangeRequired(true);
     512              _obsPool->setRefSatChangeRequired(sys, true);
    483513            }
    484514            else {
    485               LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' '
    486                   << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl;
    487515              resetAmb(obs->prn(), obsVector);
    488516            }
    … …  
    533561}
    534562
     563// Add innfinite noise to iono
     564////////////////////////////////////////////////////////////////////////////
     565t_irc t_pppFilter::addNoiseToIono(char sys) {
     566
     567  t_irc irc = failure;
     568  vector<t_pppParam*>& params = _parlist->params();
     569  for (unsigned iPar = 0; iPar < params.size(); iPar++) {
     570    t_pppParam* par = params[iPar];
     571    if (par->type() == t_pppParam::ion && par->prn().system() == sys) {
     572      int ind = par->indexNew();
     573      LOG << string(_epoTime) << " ADD IONO_NOISE TO "  <<  par->prn().toString() << endl;
     574      par->setIndex(ind);
     575      _QFlt(ind+1,ind+1) += par->sigma0() * par->sigma0();
     576      irc = success;
     577    }
     578  }
     579
     580  return irc;
     581}
     582
    535583// Compute various DOP Values
    536584////////////////////////////////////////////////////////////////////////////
    … …  
    617665// Compute datum transformation
    618666////////////////////////////////////////////////////////////////////////////
    619 void t_pppFilter::datumTransformation(const ColumnVector& xFltOld, const SymmetricMatrix& QFltOld) {
     667t_irc t_pppFilter::datumTransformation() {
     668  t_pppObsPool::t_epoch* epoch = _obsPool->lastEpoch();
     669  if (!epoch) {LOG << "!epoch" << endl;
     670    return failure;
     671  }
     672  else {
     673    LOG.setf(ios::fixed);
     674    LOG << string(epoch->epoTime()) << " DATUM TRANSFORMATION " << endl;
     675  }
     676  vector<t_pppSatObs*>& allObs = epoch->obsVector();
     677
     678  // reset old and set new refSats in last epoch (ambiguities)
     679  // ========================================================
     680  if (resetRefSatellitesLastEpoch(allObs) != true) {
     681    return failure;
     682  }
     683
     684  // set AA2
     685  // =======
     686  _parlist->set(epoch->epoTime(), allObs, _obsPool->getRefSatMap());
     687  const vector<t_pppParam*>& params = _parlist->params();
     688#ifdef BNC_DEBUG_PPP
     689  for (unsigned iPar = 0; iPar < params.size(); iPar++) {
     690    LOG << params[iPar]->toString() << "\t\t" << endl;
     691  }
     692#endif
     693  for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
     694    (iSys) ? _datumTrafo->setFirstSystem(false) : _datumTrafo->setFirstSystem(true);
     695    char sys = OPT->systems()[iSys];
     696    t_prn refPrn = (_obsPool->getRefSatMapElement(sys))->prn();
     697#ifdef BNC_DEBUG_PPP
     698    LOG << "refPrn: " << refPrn.toString() << endl;
     699#endif
     700    vector<t_pppSatObs*> obsVector;
     701    for (unsigned jj = 0; jj < allObs.size(); jj++) {
     702      if (allObs[jj]->prn().system() == sys) {
     703        obsVector.push_back(allObs[jj]);
     704      }
     705    }
     706
     707    vector<t_lc::type> LCs = OPT->LCs(sys);
     708    unsigned usedLCs = LCs.size();
     709    unsigned realUsedLCs = usedLCs;
     710    if (OPT->_pseudoObsIono && !epoch->pseudoObsIono()) {
     711        usedLCs -= 1;  // GIM not used
     712    }
     713    int hlpLCs = 0;
     714    if (OPT->_pseudoObsTropo) {
     715      hlpLCs = -1;
     716      realUsedLCs -= 1;
     717    }
     718    // max Obs
     719    unsigned maxObs = obsVector.size() * (usedLCs + hlpLCs);
     720    if (OPT->_pseudoObsTropo) {
     721      maxObs += 1;
     722    }
     723    if (OPT->_pseudoObsIono && epoch->pseudoObsIono()) {
     724      maxObs -= 1; // pseudo obs iono with respect to refSat
     725    }
     726
     727    Matrix  AA(maxObs, _parlist->nPar());
     728
     729    // Real Observations
     730    // -----------------
     731    int iObs = -1;
     732    for (unsigned ii = 0; ii < obsVector.size(); ii++) {
     733      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          }
     744        }
     745      }
     746    }
     747    _datumTrafo->updateIndices(iObs+1);
     748    _datumTrafo->prepareAA(AA.SubMatrix(1, iObs+1 , 1, _parlist->nPar()), 2);
     749  }
     750
     751  // Datum Transformation
     752  // ====================
     753#ifdef BNC_DEBUG_PPP
     754      LOG << "AA1\n"; _datumTrafo->printMatrix(_datumTrafo->AA1(), _datumTrafo->obsNum(), _datumTrafo->parNum());
     755      LOG << "AA2\n"; _datumTrafo->printMatrix(_datumTrafo->AA2(), _datumTrafo->obsNum(), _datumTrafo->parNum());
     756#endif
    620757  Matrix D21 = _datumTrafo->computeTrafoMatrix();
     758#ifdef BNC_DEBUG_PPP
     759      LOG << "D21" << endl; _datumTrafo->printMatrix(D21, _datumTrafo->parNum(), _datumTrafo->parNum());
     760#endif
     761  ColumnVector    xFltOld = _xFlt;
     762  SymmetricMatrix QFltOld = _QFlt;
     763
    621764  _QFlt << D21 * QFltOld * D21.t();
    622   _xFlt = D21 * xFltOld;
    623 }
    624 
    625 
    626 // Reset Ambiguity Parameter (cycle slip)
    627 ////////////////////////////////////////////////////////////////////////////
    628 t_irc t_pppFilter::addInfiniteNoise(t_pppParam::e_type para) {
    629   t_irc irc = failure;
    630   vector<t_pppParam*>& params = _parlist->params();
    631   for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    632     t_pppParam* par = params[iPar];
    633     if (par->type() == para) {
    634       int ind = par->indexNew();
    635       if (ind < 0) {
    636         return irc;
    637       }
    638       _QFlt(ind+1,ind+1) += par->sigma0() * par->sigma0();
    639       irc = success;
    640     }
    641   }
    642   return irc;
    643 }
     765  _xFlt =  D21 * xFltOld;
     766
     767#ifdef BNC_DEBUG_PPP
     768  LOG << "xFltOld:\n" << xFltOld << endl;
     769  LOG << "xFlt   :\n" << _xFlt   << endl;
     770#endif
     771
     772  // Reset Ambiguities after Datum Transformation
     773  // ============================================
     774  for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
     775    char sys = OPT->systems()[iSys];
     776    t_irc irc = resetAmb(_obsPool->getRefSatMapElementLastEpoch(sys), allObs);
     777    if (OPT->_obsModelType == OPT->DCMcodeBias) {
     778      if (irc == success) {
     779        addNoiseToIono(sys);}
     780    }
     781  }
     782
     783  // switch AA2 to AA1
     784  // =================
     785  _datumTrafo->switchAA();
     786
     787  return success;
     788}
     789
     790// Init datum transformation
     791////////////////////////////////////////////////////////////////////////////
     792void t_pppFilter::initDatumTransformation(const std::vector<t_pppSatObs*>& allObs) {
     793  unsigned realObs = 0;
     794  for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
     795    char system = OPT->systems()[iSys];
     796    int satNum = 0;
     797    for (unsigned jj = 0; jj < allObs.size(); jj++) {
     798      if (allObs[jj]->prn().system() == system) {
     799        satNum++;
     800      }
     801    }
     802    // all LCs
     803    unsigned realUsedLCs = OPT->LCs(system).size();
     804    // exclude pseudo obs GIM
     805    if (OPT->_pseudoObsIono) {
     806      realUsedLCs -= 1;
     807    }
     808    if (OPT->_pseudoObsTropo) {
     809      realUsedLCs -= 1;
     810    }
     811    realObs += satNum * realUsedLCs;
     812  }
     813  _datumTrafo->setObsNum(realObs);
     814  _datumTrafo->setParNum(_parlist->nPar());
     815  _datumTrafo->initAA();
     816}
     817
     818//
     819//////////////////////////////////////////////////////////////////////////////
     820bool t_pppFilter::resetRefSatellitesLastEpoch(std::vector<t_pppSatObs*>& obsVector) {
     821
     822  bool resetRefSat = false;
     823  // reference satellite definition per system
     824  for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
     825    char sys = OPT->systems()[iSys];
     826    t_pppRefSat* refSat = _obsPool->getRefSatMapElement(sys);
     827    t_prn newPrn = refSat->prn();
     828    t_prn oldPrn = _obsPool->getRefSatMapElementLastEpoch(sys);
     829#ifdef BNC_DEBUG_PPP
     830    LOG << "oldPrn: " << oldPrn.toString() << " => newPrn: " << newPrn.toString() << endl;
     831#endif
     832    vector<t_pppSatObs*>::iterator it = obsVector.begin();
     833    while (it != obsVector.end()) {
     834      t_pppSatObs* satObs = *it;
     835      if      (satObs->prn() == newPrn) {
     836        resetRefSat = true;
     837        satObs->setAsReference();
     838      }
     839      else if (satObs->prn() == oldPrn) {
     840        satObs->resetReference();
     841      }
     842     it++;
     843    }
     844  }
     845  return resetRefSat;
     846}
Note: See TracChangeset for help on using the changeset viewer.