Changeset 8956 in ntrip for trunk


Ignore:
Timestamp:
Jun 23, 2020, 11:58:46 AM (4 years ago)
Author:
stuerze
Message:

update regarding PPP

Location:
trunk/BNC/src
Files:
12 edited

Legend:

Unmodified
Added
Removed

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

    r8912 r8956  
    160160  // -------
    161161  epoTime.reset();
     162  clearObs();
    162163
    163164  // Create vector of valid observations
    … …  
    174175      }
    175176    }
    176   }
    177 
    178   // (re)set reference satellites per system if required
    179   // ---------------------------------------------------
    180   if (_opt->_refSatRequired) {
    181     // reference satellite definition per system
    182     for (unsigned iSys = 0; iSys < _opt->systems().size(); iSys++) {
    183       char system = _opt->systems()[iSys];
    184       bool refSatDefined = false;
    185       t_pppRefSat* refSat = _obsPool->getRefSatMapElement(system);
    186       for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    187         const t_pppSatObs* satObs = obsVector.at(ii);
    188         // reference satellite is unchanged
    189         if      (!_obsPool->refSatChangeRequired() && refSat->prn() == satObs->prn()) {
    190           refSatDefined = true;
    191           obsVector[ii]->setAsReference();
    192           refSat->setStatus(t_pppRefSat::unchanged);
    193         }
    194         // reference satellite has changed
    195         else if ( _obsPool->refSatChangeRequired() && refSat->prn() != satObs->prn()) {
    196           if (satObs->prn().system() == system) {
    197             refSatDefined = true;
    198             obsVector[ii]->setAsReference();
    199             refSat->setStatus(t_pppRefSat::changed);
    200             refSat->setPrn(satObs->prn());
    201           }
    202         }
    203         if (refSatDefined) {
    204           break;
    205         }
    206       }
    207       // reference satellite has to be initialized
    208       if (!refSatDefined) {
    209         for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    210           const t_pppSatObs* satObs = obsVector.at(ii);
    211           if (satObs->prn().system() == system) {
    212             obsVector[ii]->setAsReference();
    213             refSat->setStatus(t_pppRefSat::initialized);
    214             refSat->setPrn(satObs->prn());
    215           }
    216         }
    217       }
    218     }
    219     _obsPool->setRefSatChangeRequired(false);
    220177  }
    221178
    … …  
    497454        satObs->eleSat() >= _opt->_minEle &&
    498455        modelSetup == success) {
    499       if (satObs->isReference() && _opt->_pseudoObsIono) {
    500         char system = satObs->prn().system();
    501         t_pppRefSat* refSat = _obsPool->getRefSatMapElement(system);
    502         refSat->setStecValue(satObs->getIonoCodeDelay(t_frequency::G1));
    503       }
    504456      ++it;
    505457    }
    … …  
    517469void t_pppClient::processEpoch(const vector<t_satObs*>& satObs, t_output* output) {
    518470
     471  _historicalRefSats.clear();
    519472  try {
    520473    initOutput(output);
    … …  
    524477    do {
    525478      num++;
     479      if (num == 1) {
     480        LOG << "\nPPP of Epoch ";
     481        if (!_epoTimeRover.undef()) LOG << string(_epoTimeRover);
     482        LOG << "\n---------------------------------------------------------------\n";
     483      }
    526484
    527485      // Prepare Observations of the Rover
    … …  
    530488        return finish(failure);
    531489      }
    532 
    533       LOG << "\nPPP of Epoch ";
    534       if (!_epoTimeRover.undef()) LOG << string(_epoTimeRover);
    535       LOG << "\n---------------------------------------------------------------\n";
    536490
    537491      for (int iter = 1; iter <= 2; iter++) {
    … …  
    548502      _offGG = cmpOffGG(_obsRover);
    549503
    550       // Prepare Pseudo Observations of the Rover
    551       // ----------------------------------------
    552       _pseudoObsIono = preparePseudoObs(_obsRover);
    553       if (int(_obsRover.size()) < _opt->_minObs) {
    554         LOG << "t_pppClient::processEpoch not enough observations" << endl;
    555         return finish(failure);
    556       }
    557 
    558504      if (_opt->_refSatRequired) {
     505        setRefSatellites(_obsRover);
    559506        LOG.setf(ios::fixed);
    560507        QMapIterator<char, t_pppRefSat*> it(_obsPool->getRefSatMap());
    … …  
    563510          char   sys = it.key();
    564511          string prn = it.value()->prn().toString();
     512          if (num == 1) {
     513            LOG << "set   ";
     514          }
     515          else {
     516            LOG << "reset ";
     517          }
    565518          LOG << string(_epoTimeRover) << " REFSAT  " << sys << ": " << prn << endl;
    566519        }
    567520      }
     521
     522
     523      // Prepare Pseudo Observations of the Rover
     524      // ----------------------------------------
     525      _pseudoObsIono = preparePseudoObs(_obsRover);
     526
     527      if (int(_obsRover.size()) < _opt->_minObs) {
     528        LOG << "t_pppClient::processEpoch not enough observations" << endl;
     529        return finish(failure);
     530      }
     531
    568532
    569533      // Store last epoch of data
    570534      // ------------------------
    571       //_obsRover.resize(2);
    572535      _obsPool->putEpoch(_epoTimeRover, _obsRover, _pseudoObsIono);
     536      _obsPool->setHistoricalRefSatList(_historicalRefSats);
    573537
    574538      // Process Epoch in Filter
    … …  
    582546      if (_obsPool->refSatChangeRequired()) {
    583547        epochReProcessing = true;
     548        _obsPool->deleteLastEpoch();
     549        QMapIterator<char, t_pppRefSat*> it(_obsPool->getRefSatMap());
     550        while (it.hasNext()) {
     551          it.next();
     552          _historicalRefSats.append(it.value()->prn());
     553        }
    584554      }
    585555      else {
    586556        epochReProcessing = false;
    587       }
    588       // Datum transformation required?
    589       // ------------------------------
    590       if (num > 1 && epochReProcessing) {
    591         _filter->datumTransformation();
    592       }
    593 
     557
     558      }
    594559    } while (epochReProcessing);
    595560  }
    … …  
    694659//
    695660//////////////////////////////////////////////////////////////////////////////
     661void t_pppClient::setRefSatellites(std::vector<t_pppSatObs*>& obsVector) {
     662
     663  // reference satellite definition per system
     664  for (unsigned iSys = 0; iSys < _opt->systems().size(); iSys++) {
     665    char system = _opt->systems()[iSys];
     666    bool refSatDefined = false;
     667    t_pppRefSat* refSat = _obsPool->getRefSatMapElement(system);
     668    for (unsigned ii = 0; ii < obsVector.size(); ii++) {
     669      t_pppSatObs* satObs = obsVector.at(ii);
     670      if (satObs->eleSat() < _opt->_minEle) {continue;}
     671      // reference satellite is unchanged
     672      if      (!_obsPool->refSatChangeRequired() &&  refSat->prn() == satObs->prn()) {
     673        refSatDefined = true;
     674        obsVector[ii]->setAsReference();
     675      }
     676      // reference satellite has changed
     677      else if ( _obsPool->refSatChangeRequired() && refSat->prn() != satObs->prn() && !_historicalRefSats.contains(satObs->prn())) {
     678        if (satObs->prn().system() == system) {
     679          refSatDefined = true;
     680          obsVector[ii]->setAsReference();
     681          refSat->setPrn(satObs->prn());
     682        }
     683      }
     684      if (refSatDefined) {
     685        if (OPT->_pseudoObsIono) {
     686          refSat->setStecValue(satObs->getIonoCodeDelay(t_frequency::G1));
     687        }
     688        break;
     689      }
     690    }
     691    // reference satellite has to be initialized
     692    if (!refSatDefined) {
     693      for (unsigned ii = 0; ii < obsVector.size(); ii++) {
     694        t_pppSatObs* satObs = obsVector.at(ii);
     695        if (satObs->eleSat() < _opt->_minEle) {
     696          continue;
     697        }
     698        if (satObs->prn().system() == system) {
     699          obsVector[ii]->setAsReference();
     700          refSat->setPrn(satObs->prn());
     701          if (OPT->_pseudoObsIono) {
     702            refSat->setStecValue(satObs->getIonoCodeDelay(t_frequency::G1));
     703          }
     704          refSatDefined = true;
     705          break;
     706        }
     707      }
     708    }
     709  }
     710  _obsPool->setRefSatChangeRequired(false);
     711
     712}
     713
     714//
     715//////////////////////////////////////////////////////////////////////////////
    696716void t_pppClient::reset() {
    697717

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

    r8911 r8956  
    6060                    ColumnVector& xyzc, bool print);
    6161  double cmpOffGG(std::vector<t_pppSatObs*>& obsVector);
     62  void setRefSatellites(std::vector<t_pppSatObs*>& obsVector);
    6263
    6364  t_output*                 _output;
    … …  
    7475  t_tides*                  _tides;
    7576  bool                      _pseudoObsIono;
     77  QList<t_prn>              _historicalRefSats;
    7678};
    7779

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

    r8915 r8956  
    8383  // Set Parameters
    8484  // --------------
    85   _parlist->set(_epoTime, allObs);
     85  _parlist->set(_epoTime, allObs, _obsPool->getRefSatMap());
    8686  const vector<t_pppParam*>& params = _parlist->params();
    8787
    … …  
    9090  ColumnVector    xFltOld = _xFlt;
    9191  SymmetricMatrix QFltOld = _QFlt;
     92
    9293  _QFlt.ReSize(_parlist->nPar()); _QFlt = 0.0;
    9394  _xFlt.ReSize(_parlist->nPar()); _xFlt = 0.0;
    … …  
    114115  predictCovCrdPart(QFltOld);
    115116
     117  // Init Datum Trafo
     118  // ----------------------------------------
     119  if ((OPT->_obsModelType == OPT->DCMcodeBias ||
     120       OPT->_obsModelType == OPT->DCMphaseBias) &&
     121      (_numEpoProcessing == 1)) {
     122    _numAllUsedLCs = 0;
     123    for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
     124      char system = OPT->systems()[iSys];
     125      _numAllUsedLCs += OPT->LCs(system).size();
     126      if (OPT->_pseudoObsIono && !epoch->pseudoObsIono()) {
     127        _numAllUsedLCs -= 1;  // GIM not used
     128      }
     129    }
     130    int maxObs = allObs.size() * _numAllUsedLCs;
     131    _datumTrafo->initAA(maxObs, _parlist->nPar());
     132  }
    116133
    117134  // Pre-Process Satellite Systems separately
    … …  
    121138      OPT->_obsModelType == OPT->DCMphaseBias) {
    122139    preProcessing = true;
    123     _numAllUsedLCs = 0;
    124140    for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
     141      (iSys) ? _datumTrafo->setFirstSystem(false) : _datumTrafo->setFirstSystem(true);
    125142      char system = OPT->systems()[iSys];
    126       _numAllUsedLCs += OPT->LCs(system).size();
    127       if (OPT->_pseudoObsIono && !epoch->pseudoObsIono()) {
    128         _numAllUsedLCs -= 1;  // GIM not used
    129       }
    130143      if (OPT->_refSatRequired) {
    131144        _refPrn = (_obsPool->getRefSatMapElement(system))->prn();
     145        if (_obsPool->hasHistoricalRefSat(_refPrn)) {
     146          return failure;
     147        }
    132148      }
    133149      vector<t_pppSatObs*> obsVector;
    … …  
    142158      }
    143159    }
    144   }
    145 
    146   if (_numEpoProcessing == 1) {
    147     int maxObs = allObs.size() * _numAllUsedLCs;
    148     _datumTrafo->initAA(maxObs, _parlist->nPar());
     160    // refSat change required?
     161    // -----------------------
     162    if (_obsPool->refSatChangeRequired()) {
     163      _xFlt = xFltOld;
     164      _QFlt = QFltOld;
     165      return success;
     166    }
     167    else if (!_obsPool->refSatChangeRequired() &&
     168             _numEpoProcessing > 1) {
     169      datumTransformation(xFltOld, QFltOld);
     170    }
    149171  }
    150172
    … …  
    154176  for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
    155177    char system = OPT->systems()[iSys];
    156     (iSys) ? _datumTrafo->setFirstSystem(false) :
    157                  _datumTrafo->setFirstSystem(true);
    158178    if (OPT->_refSatRequired) {
    159179      _refPrn = (_obsPool->getRefSatMapElement(system))->prn();
    … …  
    174194  }
    175195
    176   // refSat change required?
    177   // -----------------------
    178   if (_obsPool->refSatChangeRequired()) {
    179     _xFlt = xFltOld;
    180     _QFlt = QFltOld;
    181   }
    182196  // close epoch processing
    183197  // ----------------------
    184   else {
    185     cmpDOP(allObs);
    186     _parlist->printResult(_epoTime, _QFlt, _xFlt);
    187     _lastEpoTimeOK = _epoTime;  // remember time of last successful epoch processing
    188   }
     198  cmpDOP(allObs);
     199  _parlist->printResult(_epoTime, _QFlt, _xFlt);
     200  _lastEpoTimeOK = _epoTime;  // remember time of last successful epoch processing
    189201
    190202  return success;
    … …  
    198210                                 bool pseudoObsIonoAvailable,
    199211                                 bool preProcessing) {
    200   qDebug() << "======t_pppFilter::processSystem=======";
    201212  LOG.setf(ios::fixed);
    202213
    203214  // Detect Cycle Slips
    204215  // ------------------
     216
    205217  if (detectCycleSlips(LCs, obsVector, refPrn, preProcessing) != success) {
    206218    return failure;
    … …  
    221233    maxObs -= 1;
    222234  }
    223   qDebug() << "par.size()      : " << _parlist->nPar() << " LCs.size(): " << usedLCs;
    224   qDebug() << "obsVector.size(): " << obsVector.size() << " maxObs    : " << maxObs;
    225235
    226236  // Outlier Detection Loop
    227237  // ----------------------
    228238  for (unsigned iOutlier = 0; iOutlier < maxObs; iOutlier++) {
    229     qDebug() << "iOutlier: " << iOutlier;
    230239
    231240    if (iOutlier > 0) {
    … …  
    236245    // First-Design Matrix, Terms Observed-Computed, Weight Matrix
    237246    // -----------------------------------------------------------
    238     qDebug() << "A(" << maxObs << "," << _parlist->nPar() << ")";
    239247    Matrix                AA(maxObs, _parlist->nPar());
    240248    ColumnVector          ll(maxObs);
    241249    DiagonalMatrix        PP(maxObs); PP = 0.0;
    242 
    243     // TETSPLOT
    244     for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    245       const t_pppParam* par = params[iPar];
    246       cout << "  " << par->toString() <<endl;
    247     }
    248     cout << endl;
    249     //END TETSPLOT
    250250
    251251    int iObs = -1;
    … …  
    253253    vector<t_lc::type>   usedTypes;
    254254    for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    255       t_pppSatObs* obs = obsVector[ii]; qDebug() << "SATELLITE: " << obs->prn().toString().c_str() << "isRef: " << obs->isReference();
     255      t_pppSatObs* obs = obsVector[ii];
    256256      if (!obs->outlier()) {
    257257        for (unsigned jj = 0; jj < usedLCs; jj++) {
    … …  
    263263          for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    264264            const t_pppParam* par = params[iPar];
    265             AA[iObs][iPar] = par->partial(_epoTime, obs, tLC);
    266             cout << setw(10) << par->partial(_epoTime, obs, tLC);
    267           }
    268           cout << endl;
     265            AA[iObs][iPar] = par->partial(_epoTime, obs, tLC, refPrn);
     266          }
    269267          ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1));
    270268          PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC));
    271           //qDebug() << "ll[iObs]: " << ll[iObs];
    272           //qDebug() << "PP[iObs]: " << PP[iObs];
    273         }
    274       }
    275     }
    276 
    277     if ((!preProcessing) &&
    278        (OPT->_obsModelType == OPT->DCMcodeBias ||
    279         OPT->_obsModelType == OPT->DCMphaseBias)) {
     269        }
     270      }
     271    }
     272
     273    if ((preProcessing && ! iOutlier) &&
     274        (OPT->_obsModelType == OPT->DCMcodeBias || OPT->_obsModelType == OPT->DCMphaseBias)) {
    280275          _datumTrafo->updateIndices(maxObs);
    281276      _datumTrafo->prepareAA(AA, _numEpoProcessing, _parlist->nPar());
     277      if (_obsPool->refSatChangeRequired()) { // from detectCycleSlips()
     278        return success;
     279      }
    282280    }
    283281
    … …  
    331329        }
    332330      }
    333       if      (par && preProcessing) {
    334         if (par->prn() == refPrn) {
     331      if (preProcessing) {
     332        if (par && obs->prn() == refPrn) {
    335333          _obsPool->setRefSatChangeRequired(true);
    336         }
    337       }
    338       else if (par && !preProcessing) {
    339         if (par->ambResetCandidate()) {
    340           resetAmb(par->prn(), obsVector, &QSav, &xSav);
     334          break;
     335        }
     336      }
     337      else {
     338        if (par) {
     339          if (par->ambResetCandidate() || _obsPool->hasHistoricalRefSat(par->prn())) {
     340            resetAmb(par->prn(), obsVector, &QSav, &xSav);
     341          }
     342          else {
     343            par->setAmbResetCandidate();
     344            obs->setOutlier();
     345          }
    341346        }
    342347        else {
    343           par->setAmbResetCandidate();
    344348          obs->setOutlier();
    345349        }
    346350      }
    347       else if (!par && !preProcessing){
    348         obs->setOutlier();
    349       }
    350     }
    351 
     351    }
    352352    // Print Residuals
    353353    // ---------------
    … …  
    381381                                    const t_prn& refPrn,
    382382                                    bool preProcessing) {
    383 
    384   const double            SLIP       = 20.0;  // slip threshold
     383  if (_obsPool->hasHistoricalRefSat(refPrn)) {
     384    return success;
     385  }
     386
     387  double                  SLIP       = 20.0;  // slip threshold
    385388  string                  epoTimeStr = string(_epoTime);
    386389  const vector<t_pppParam*>& params  = _parlist->params();
     390
     391  if (preProcessing) {
     392    SLIP += 10.0;
     393  }
    387394
    388395  for (unsigned ii = 0; ii < LCs.size(); ii++) {
    … …  
    391398      for (unsigned iObs = 0; iObs < obsVector.size(); iObs++) {
    392399        const t_pppSatObs* obs = obsVector[iObs];
     400
     401        if (preProcessing && obs->prn() != refPrn) {continue;}
    393402
    394403        // Check set Slips and Jump Counters
    … …  
    405414          slip = true;
    406415        }
    407         _slips[obs->prn()]._obsSlipCounter = obs->slipCounter();
     416        if (!preProcessing) {
     417          _slips[obs->prn()]._obsSlipCounter = obs->slipCounter();
     418        }
    408419
    409420        if (_slips[obs->prn()]._biasJumpCounter != -1 &&
    … …  
    412423          slip = true;
    413424        }
    414 
     425        if (!preProcessing) {
     426          _slips[obs->prn()]._biasJumpCounter = obs->biasJumpCounter();
     427        }
    415428        // Slip Set
    416429        // --------
    417430        if (slip) {
    418431          if (preProcessing) {
    419             if (obs->prn() == refPrn) {
    420               _obsPool->setRefSatChangeRequired(true);
    421             }
     432            _obsPool->setRefSatChangeRequired(true);
    422433          }
    423434          else {
    … …  
    431442          for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    432443            const t_pppParam* par = params[iPar];
    433             AA[iPar] = par->partial(_epoTime, obs, tLC);
    434           }
     444            AA[iPar] = par->partial(_epoTime, obs, tLC, refPrn);
     445          }
     446
    435447          double ll = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA);
    436448          double vv = DotProduct(AA, _xFlt) - ll;
    437449          if (fabs(vv) > SLIP) {
    438450            if (preProcessing) {
    439               if (obs->prn() == refPrn) {
    440                 _obsPool->setRefSatChangeRequired(true);
    441               }
     451              _obsPool->setRefSatChangeRequired(true);
    442452            }
    443453            else {
    … …  
    451461    }
    452462  }
    453 
    454463  return success;
    455464}
    … …  
    459468t_irc t_pppFilter::resetAmb(t_prn prn, const vector<t_pppSatObs*>& obsVector,
    460469                            SymmetricMatrix* QSav, ColumnVector* xSav) {
     470
    461471  t_irc irc = failure;
    462472  vector<t_pppParam*>& params = _parlist->params();
    … …  
    504514    for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    505515      t_pppSatObs* obs = obsVector[ii];
     516      char system = obs->prn().system();
     517      t_prn refPrn = t_prn();
     518      if (OPT->_refSatRequired) {
     519        refPrn = _obsPool->getRefSatMapElement(system)->prn();
     520      }
    506521      if (obs->isValid() && !obs->outlier()) {
    507522        ++_numSat;
    508523        for (unsigned iPar = 0; iPar < numPar; iPar++) {
    509524          const t_pppParam* par = _parlist->params()[iPar];
    510           AA[_numSat-1][iPar] = par->partial(_epoTime, obs, t_lc::c1);
     525          AA[_numSat-1][iPar] = par->partial(_epoTime, obs, t_lc::c1, refPrn);
    511526        }
    512527      }
    … …  
    571586// Compute datum transformation
    572587////////////////////////////////////////////////////////////////////////////
    573 void t_pppFilter::datumTransformation(void) {
    574   _QFlt = _datumTrafo->varCov(Q());
    575 }
    576 
     588void t_pppFilter::datumTransformation(const ColumnVector& xFltOld, const SymmetricMatrix& QFltOld) {
     589
     590  Matrix D21 = _datumTrafo->computeTrafoMatrix();
     591  _QFlt << D21 * QFltOld * D21.t();
     592  _xFlt = D21 * xFltOld;
     593}
     594
     595
     596// Reset Ambiguity Parameter (cycle slip)
     597////////////////////////////////////////////////////////////////////////////
     598t_irc t_pppFilter::addInfiniteNoise(t_pppParam::e_type para) {
     599  t_irc irc = failure;
     600  vector<t_pppParam*>& params = _parlist->params();
     601  for (unsigned iPar = 0; iPar < params.size(); iPar++) {
     602    t_pppParam* par = params[iPar];
     603    if (par->type() == para) {
     604      int ind = par->indexNew();
     605      if (ind < 0) {
     606        return irc;
     607      }
     608      _QFlt(ind+1,ind+1) += par->sigma0() * par->sigma0();
     609      irc = success;
     610    }
     611  }
     612  return irc;
     613}

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

    r8915 r8956  
    2222  t_irc processEpoch(int num);
    2323
    24   void datumTransformation();
     24  void datumTransformation(const ColumnVector& xFltOld, const SymmetricMatrix& QFltOld);
    2525
    2626  const ColumnVector&    x() const {return _xFlt;}
    … …  
    7575  };
    7676
    77   class t_datumTrafo{
     77  class t_datumTrafo {
    7878  public:
    7979    t_datumTrafo () {initIndices();}
    … …  
    8282    bool firstSystem() {return _firstSys;}
    8383    void updateIndices(int maxObs) {
    84       if (_firstSys) {
     84      if (firstSystem()) {
    8585        initIndices();
    8686      }
    … …  
    9494      _AA2.ReSize(maxObs, numPar); _AA2 = 0.0;
    9595    }
    96     void prepareAA(Matrix& AA, int _numEpoProcessing, int nPar) {
    97       Matrix& Prep = _AA2;
    98       if (_numEpoProcessing == 1) {
    99         Prep = _AA1;
     96    void prepareAA(const Matrix& AA, int numEpoProcessing, int nPar) {
     97      Matrix* Prep = &_AA2;
     98      if (numEpoProcessing == 1) {
     99        Prep = &_AA1;
    100100      }
    101       Prep.SubMatrix(_firstRow, _lastRow, 1, nPar) = AA;
     101      Prep->SubMatrix(_firstRow, _lastRow, 1, nPar) << AA;
    102102    }
    103     Matrix varCov(const SymmetricMatrix& QFlt) {
     103
     104    Matrix computeTrafoMatrix() {
    104105      Matrix D21 = (_AA2.t() * _AA2).i() * _AA2.t() * _AA1;
    105       return D21 * QFlt * D21.t();
     106      return D21;
     107    }
     108
     109    static void printMatrix(const Matrix& X, int nRow, int nCol) {
     110      for (int rr = 0; rr < nRow; rr++) {
     111        for (int cc = 0; cc < nCol; cc++) {
     112          cout << setw(7) << setprecision(4) << X[rr][cc] << " ";
     113        }
     114        cout << endl;      }
     115      cout << endl;
    106116    }
    107117    int     _firstRow;
    … …  
    130140  void predictCovCrdPart(const SymmetricMatrix& QFltOld);
    131141
     142  t_irc addInfiniteNoise(t_pppParam::e_type para);
     143
    132144  bncTime         _epoTime;
    133145  t_pppParlist*   _parlist;

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

    r8911 r8956  
    100100  const unsigned MAXSIZE = 2;
    101101
    102   if (refSatChangeRequired()) {
    103     delete _epochs.back();
    104     _epochs.pop_back();
    105   }
    106   _epochs.push_back(new t_epoch(epoTime, obsVector,pseudoObsIono));
     102  _epochs.push_back(new t_epoch(epoTime, obsVector, pseudoObsIono));
     103
    107104  if (_epochs.size() > MAXSIZE) {
    108105    delete _epochs.front();
    … …  
    110107  }
    111108}
     109
     110//
     111/////////////////////////////////////////////////////////////////////////////
     112void t_pppObsPool::deleteLastEpoch() {
     113
     114  if (!_epochs.empty()) {
     115    delete _epochs.back();
     116    _epochs.pop_back();
     117  }
     118}

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

    r8911 r8956  
    3838                bool pseudoObs);
    3939
     40  void deleteLastEpoch();
     41
    4042  const t_satCodeBias* satCodeBias(const t_prn& prn) const {
    4143    return _satCodeBiases[prn.toInt()];
    … …  
    7678  bool refSatChangeRequired() {return _refSatChangeRequired;}
    7779
     80  void setHistoricalRefSatList(QList<t_prn>& historicalRefSats) {_historicalRefSats = historicalRefSats;}
     81
     82  bool hasHistoricalRefSat(t_prn prn) {return _historicalRefSats.contains(prn);}
     83
    7884 private:
    7985  t_satCodeBias*           _satCodeBiases[t_prn::MAXPRN+1];
    … …  
    8389  QMap<char, t_pppRefSat*> _refSatMap;
    8490  bool                     _refSatChangeRequired;
     91  QList<t_prn>             _historicalRefSats;
    8592};
    8693

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

    r8905 r8956  
    123123////////////////////////////////////////////////////////////////////////////
    124124double t_pppParam::partial(const bncTime& /* epoTime */, const t_pppSatObs* obs,
    125                            const t_lc::type& tLC) const {//qDebug() << "t_pppParam::partial: " << tLC;
     125                           const t_lc::type& tLC, const t_prn refPrn) const {
    126126
    127127  // Special Case - Melbourne-Wuebbena
    … …  
    184184    return  1.0 / sin(obs->eleSat());
    185185  case ion:
    186 //    qDebug() << "refPrn: " << _refPrn.toString().c_str();
     186
    187187    if (obs->prn() == _prn) {
    188188      if      (tLC == t_lc::c1) {
    … …  
    202202      }
    203203    }
    204     if (tLC == t_lc::GIM && _prn == _refPrn) {
     204    if (tLC == t_lc::GIM && _prn == refPrn) {
    205205      return  1.0;
    206206    }
    … …  
    296296////////////////////////////////////////////////////////////////////////////
    297297t_irc t_pppParlist::set(const bncTime& epoTime,
    298     const std::vector<t_pppSatObs*>& obsVector) {
     298    const std::vector<t_pppSatObs*>& obsVector,
     299    const QMap<char, t_pppRefSat*>& refSatMap) {
    299300
    300301  // Remove some Parameters
    … …  
    310311    }
    311312
    312     else if (par->type() == t_pppParam::amb ||
    313              par->type() == t_pppParam::crdX ||
     313    else if (par->type() == t_pppParam::crdX ||
    314314             par->type() == t_pppParam::crdY ||
    315315             par->type() == t_pppParam::crdZ) {
    316316      if (par->lastObsTime().valid() && (epoTime - par->lastObsTime() > 60.0)) {
     317        remove = true;
     318      }
     319    }
     320    else if (par->type() == t_pppParam::amb) {
     321      char system = par->prn().system();
     322      t_prn refPrn = t_prn();
     323      if (OPT->_refSatRequired) {
     324        refPrn = (refSatMap[system])->prn();
     325      }
     326      if ((par->lastObsTime().valid() && (epoTime - par->lastObsTime() > 60.0)) ||
     327          ( refPrn == par->prn())) {
    317328        remove = true;
    318329      }

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

    r8905 r8956  
    77#include "t_prn.h"
    88#include "bnctime.h"
     9#include "pppRefSat.h"
    910
    1011namespace BNC_PPP {
    … …  
    2324  e_type type() const {return _type;}
    2425  double x0()  const {return _x0;}
    25   double partial(const bncTime& epoTime, const t_pppSatObs* obs, const t_lc::type& tLC) const;
     26  double partial(const bncTime& epoTime, const t_pppSatObs* obs,
     27                 const t_lc::type& tLC, const t_prn refPrn) const;
    2628  bool   epoSpec() const {return _epoSpec;}
    2729  bool   isEqual(const t_pppParam* par2) const {
    … …  
    5153  unsigned ambNumEpo() const           {return _ambInfo ? _ambInfo->_numEpo : 0;}
    5254  void     stepAmbNumEpo()             {if (_ambInfo) _ambInfo->_numEpo += 1;}
    53   void     setRefPrn(t_prn prn)        {_refPrn = prn;}
    5455
    5556  static bool sortFunction(const t_pppParam* p1, const t_pppParam* p2) {
    … …  
    8182  e_type       _type;
    8283  t_prn        _prn;
    83   t_prn        _refPrn;
    8484  t_lc::type   _tLC;
    8585  double       _x0;
    … …  
    9999  ~t_pppParlist();
    100100
    101   t_irc set(const bncTime& epoTime, const std::vector<t_pppSatObs*>& obsVector);
     101  t_irc set(const bncTime& epoTime, const std::vector<t_pppSatObs*>& obsVector,
     102            const QMap<char, t_pppRefSat*>& refSatMap);
    102103  unsigned nPar() const {return _params.size();}
    103104  const std::vector<t_pppParam*>& params() const {return _params;}

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

    r8905 r8956  
    1414////////////////////////////////////////////////////////////////////////////
    1515t_pppRefSat::t_pppRefSat() {
    16   _status     = undefined;
    1716  _prn        = t_prn();
    1817  _stecValue  = 0.0;

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

    r8905 r8956  
    1717  t_pppRefSat ();
    1818  ~t_pppRefSat();
    19   enum     e_status {undefined, initialized, changed, unchanged};// ev. wieder löschen
    20   e_status status() const {return _status;}
    21   void     setStatus(e_status status) {_status = status;}
    2219  t_prn    prn() {return _prn;}
    2320  void     setPrn(t_prn prn) {_prn = prn;}
    … …  
    2522  void     setStecValue(double stecValue) {_stecValue = stecValue;}
    2623private:
    27   e_status _status;
    2824  t_prn    _prn;
    2925  double   _stecValue;

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

    r8905 r8956  
    342342  }
    343343  if (tLC == t_lc::GIM) {
    344     retVal = 3 * (OPT->_sigmaGIMdiff * OPT->_sigmaGIMdiff);
     344    retVal = 3.0 * OPT->_sigmaGIMdiff * OPT->_sigmaGIMdiff;
    345345  }
    346346  return sqrt(retVal);
    … …  
    495495  _model._set = true;
    496496
    497   printModel();
     497  //printModel();
    498498
    499499  return success;
    … …  
    503503////////////////////////////////////////////////////////////////////////////
    504504void t_pppSatObs::printModel() const {
    505 // TODO: cout should be LOG
    506   cout.setf(ios::fixed);
    507   cout << "\nMODEL for Satellite " << _prn.toString() << (isReference() ? " (Reference Satellite)" : "") << endl
    508        << "======================= " << endl
    509        << "PPP STRATEGY  : " <<  OPT->_obsmodelTypeStr.at((int)OPT->_obsModelType).toLocal8Bit().constData()
    510        <<  ((OPT->_pseudoObsIono) ? " with pseudo-observations for STEC" : "")  <<  endl
     505
     506  LOG.setf(ios::fixed);
     507  LOG << "\nMODEL for Satellite " << _prn.toString() << (isReference() ? " (Reference Satellite)" : "")
     508
     509      << "======================= " << endl
     510      << "PPP STRATEGY  : " <<  OPT->_obsmodelTypeStr.at((int)OPT->_obsModelType).toLocal8Bit().constData()
     511      <<  ((OPT->_pseudoObsIono) ? " with pseudo-observations for STEC" : "")  <<  endl
    511512      << "RHO           : " << setw(12) << setprecision(3) << _model._rho              << endl
    512513      << "ELE           : " << setw(12) << setprecision(3) << _model._eleSat * RHO_DEG << endl
    … …  
    528529      string frqStr = t_frequency::toString(t_frequency::type(iFreq));
    529530      if (_prn.system() == frqStr[0]) {
    530       cout << "PCO           : " << frqStr << setw(12) << setprecision(3) << _model._antPCO[iFreq]       << endl
     531      LOG << "PCO           : " << frqStr << setw(12) << setprecision(3) << _model._antPCO[iFreq]       << endl
    531532           << "BIAS CODE     : " << frqStr << setw(12) << setprecision(3) << _model._codeBias[iFreq]     << endl
    532533           << "BIAS PHASE    : " << frqStr << setw(12) << setprecision(3) << _model._phaseBias[iFreq]    << endl

  • trunk/BNC/src/pppMain.cpp

    r8905 r8956  
    416416    opt->_aprSigCodeBias  = 1000.0;
    417417    opt->_aprSigPhaseBias = 1000.0;
    418     opt->_noiseIon        = 1.0;
    419     opt->_noiseCodeBias   = 1.0;
    420     opt->_noisePhaseBias  = 0.1;
    421     opt->_sigmaGIMdiff    = 0.05; // pseudo observation GIM: STEC(ref_sat) - STEC(sat)
     418    // TODO: Find realistic values!!!!!!
     419    opt->_noiseIon        = 10.0;
     420    opt->_noiseCodeBias   = 1.00;
     421    opt->_noisePhaseBias  = 5.00;
     422    opt->_sigmaGIMdiff    = 2.00; //pseudo observation GIM: STEC(ref_sat) - STEC(sat)
    422423
    423424    _options << opt;
Note: See TracChangeset for help on using the changeset viewer.