Changeset 8905 in ntrip


Ignore:
Timestamp:
Mar 18, 2020, 11:13:50 AM (5 years ago)
Author:
stuerze
Message:

some developments regarding PPP, not completed!

Location:
trunk/BNC/src
Files:
2 added
21 edited

Legend:

Unmodified
Added
Removed

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

    r8453 r8905  
    5555  _obsPool  = new t_pppObsPool();
    5656  _staRover = new t_pppStation();
    57   _filter   = new t_pppFilter();
     57  _filter   = new t_pppFilter(_obsPool);
    5858  _tides    = new t_tides();
    59 
     59  _antex    = 0;
    6060  if (!_opt->_antexFileName.empty()) {
    6161    _antex = new bncAntex(_opt->_antexFileName.c_str());
    6262  }
    63   else {
    64     _antex = 0;
    65   }
    66 
    6763  if (!_opt->_blqFileName.empty()) {
    68     _loading = new t_loading(_opt->_blqFileName.c_str());
    69   }
    70   else {
    71     _loading = 0;
    72   }
    73 
     64    if (_tides->readBlqFile(_opt->_blqFileName.c_str()) == success) {
     65      //_tides->printAllBlqSets();
     66    }
     67  }
     68
     69  if (OPT->_refSatRequired) {
     70    for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
     71      char system = OPT->systems()[iSys];
     72      _obsPool->initRefSatMapElement(system);
     73    }
     74  }
     75  _offGG = 0.0;
    7476  CLIENTS.setLocalData(this);  // CLIENTS takes ownership over "this"
    7577}
    … …  
    8082  delete _log;
    8183  delete _opt;
     84  delete _filter;
    8285  delete _ephPool;
    8386  delete _obsPool;
    … …  
    8689    delete _antex;
    8790  }
    88   delete _filter;
    8991  delete _tides;
    90   if (_loading) {
    91     delete _loading;
    92   }
    9392  clearObs();
    9493}
    … …  
    157156t_irc t_pppClient::prepareObs(const vector<t_satObs*>& satObs,
    158157                              vector<t_pppSatObs*>& obsVector, bncTime& epoTime) {
     158
    159159  // Default
    160160  // -------
    … …  
    176176  }
    177177
     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->epoReProcessing() && 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->epoReProcessing() && 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->setEpoReProcessing(false); //TODO: später erst nach Trafo false setzen
     220  }
     221
    178222  // Check whether data are synchronized, compute epoTime
    179223  // ----------------------------------------------------
    … …  
    202246}
    203247
     248//
     249//////////////////////////////////////////////////////////////////////////////
     250bool t_pppClient::preparePseudoObs(std::vector<t_pppSatObs*>& obsVector) {
     251
     252  bool pseudoObsIono = false;
     253
     254  if (OPT->_pseudoObsIono) {
     255    vector<t_pppSatObs*>::iterator it = obsVector.begin();
     256    while (it != obsVector.end()) {
     257      t_pppSatObs* satObs = *it;
     258      char system = satObs->prn().system();
     259      t_pppRefSat* refSat = _obsPool->getRefSatMapElement(system);
     260      double stecRef = refSat->stecValue();
     261      if (stecRef && !satObs->isReference()) {
     262        pseudoObsIono = true;
     263        satObs->setPseudoObsIono(t_frequency::G1, stecRef);
     264      }
     265      satObs->printObsMinusComputed();
     266      it++;
     267    }
     268  }
     269
     270  return pseudoObsIono;
     271}
     272
    204273// Compute the Bancroft position, check for blunders
    205274//////////////////////////////////////////////////////////////////////////////
    … …  
    250319    for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    251320      const t_pppSatObs* satObs = obsVector.at(ii);
    252       if ( satObs->isValid() && satObs->prn().system() == 'G' &&
    253            (!satObs->modelSet() || satObs->eleSat() >= OPT->_minEle) ) {
     321      if (satObs->isValid() &&
     322          satObs->prn().system() == 'G' &&
     323          (!satObs->modelSet() || satObs->eleSat() >= OPT->_minEle) ) {
    254324        ColumnVector rr = satObs->xc().Rows(1,3) - xyzc.Rows(1,3);
    255         double res = rr.norm_Frobenius() - satObs->obsValue(tLC)
     325        double res = rr.NormFrobenius() - satObs->obsValue(tLC)
    256326                   - (satObs->xc()[3] - xyzc[3]) * t_CST::c;
    257327        if (std::isnan(res) || fabs(res) > maxRes) {
    … …  
    393463t_irc t_pppClient::cmpModel(t_pppStation* station, const ColumnVector& xyzc,
    394464                               vector<t_pppSatObs*>& obsVector) {
    395 
    396465  bncTime time;
    397466  time = _epoTimeRover;
    398467  station->setName(OPT->_roverName);
    399468  station->setAntName(OPT->_antNameRover);
     469  station->setEpochTime(time);
    400470  if (OPT->xyzAprRoverSet()) {
    401471    station->setXyzApr(OPT->_xyzAprRover);
    … …  
    412482  // Tides
    413483  // -----
    414   station->setTideDspl( _tides->displacement(time, station->xyzApr()) );
    415 
    416   // Ionosphere
    417   // ----------
    418   station->setIonoEpochTime(time);
     484  station->setTideDsplEarth(_tides->earth(time, station->xyzApr()));
     485  station->setTideDsplOcean(_tides->ocean(time, station->xyzApr(), station->name()));
    419486
    420487  // Observation model
    … …  
    430497        satObs->eleSat() >= OPT->_minEle &&
    431498        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      }
    432504      ++it;
    433505    }
    … …  
    447519  try {
    448520    initOutput(output);
    449 
    450     // Prepare Observations of the Rover
    451     // ---------------------------------
    452     if (prepareObs(satObs, _obsRover, _epoTimeRover) != success) {
    453       return finish(failure);
    454     }
    455 
    456     LOG << "\nPPP of Epoch ";
    457     if (!_epoTimeRover.undef()) LOG << string(_epoTimeRover);
    458     LOG << "\n---------------------------------------------------------------\n";
    459 
    460     for (int iter = 1; iter <= 2; iter++) {
    461       ColumnVector xyzc(4); xyzc = 0.0;
    462       bool print = (iter == 2);
    463       if (cmpBancroft(_epoTimeRover, _obsRover, xyzc, print) != success) {
     521    _num = 0;
     522    _obsPool->setEpoReProcessing(false); // initialize for epoch
     523
     524    do {
     525      _num++;
     526
     527      // Prepare Observations of the Rover
     528      // ---------------------------------
     529      if (prepareObs(satObs, _obsRover, _epoTimeRover) != success) {
    464530        return finish(failure);
    465531      }
    466       if (cmpModel(_staRover, xyzc, _obsRover) != success) {
     532
     533      LOG << "\nPPP of Epoch ";
     534      if (!_epoTimeRover.undef()) LOG << string(_epoTimeRover);
     535      LOG << "\n---------------------------------------------------------------\n";
     536
     537      for (int iter = 1; iter <= 2; iter++) {
     538        ColumnVector xyzc(4); xyzc = 0.0;
     539        bool print = (iter == 2);
     540        if (cmpBancroft(_epoTimeRover, _obsRover, xyzc, print) != success) {
     541          return finish(failure);
     542        }
     543        if (cmpModel(_staRover, xyzc, _obsRover) != success) {
     544          return finish(failure);
     545        }
     546      }
     547
     548      _offGG = cmpOffGG(_obsRover);
     549
     550      // Prepare Pseudo Observations of the Rover
     551      // ----------------------------------------
     552      _pseudoObsIono = preparePseudoObs(_obsRover);//qDebug() << "_pseudoObsIonoAvailable: " << _pseudoObsIono;
     553
     554      if (int(_obsRover.size()) < OPT->_minObs) {
     555        LOG << "t_pppClient::processEpoch not enough observations" << endl;
    467556        return finish(failure);
    468557      }
    469     }
    470 
    471     _offGG = cmpOffGG(_obsRover);
    472 
    473     if (int(_obsRover.size()) < OPT->_minObs) {
    474       LOG << "t_pppClient::processEpoch not enough observations" << endl;
    475       return finish(failure);
    476     }
    477 
    478     // Store last epoch of data
    479     // ------------------------
    480     _obsPool->putEpoch(_epoTimeRover, _obsRover);
    481 
    482     // Process Epoch in Filter
    483     // -----------------------
    484     if (_filter->processEpoch(_obsPool) != success) {
    485       return finish(failure);
    486     }
     558
     559      if (OPT->_refSatRequired) {
     560        LOG.setf(ios::fixed);
     561        QMapIterator<char, t_pppRefSat*> it(_obsPool->getRefSatMap());
     562        while (it.hasNext()) {
     563          it.next();
     564          char   sys = it.key();
     565          string prn = it.value()->prn().toString();
     566          LOG << string(_epoTimeRover) << " REFSAT  " << sys << ": " << prn << endl;
     567        }
     568      }
     569
     570      // Store last epoch of data
     571      // ------------------------
     572      //_obsRover.resize(2);
     573      _obsPool->putEpoch(_epoTimeRover, _obsRover, _pseudoObsIono);
     574
     575      // Process Epoch in Filter
     576      // -----------------------
     577      if (_filter->processEpoch(_num) != success) {
     578        return finish(failure);
     579      }
     580      // if num > 1 und !obsPool->epoReProcessing() => filter ->datumTransformation
     581    } while (_obsPool->epoReProcessing());
    487582  }
    488583  catch (Exception& exc) {
    … …  
    588683void t_pppClient::reset() {
    589684
    590   // to delete all parameters
    591   delete _filter;
    592   _filter   = new t_pppFilter();
    593 
    594685  // to delete old orbit and clock corrections
    595686  delete _ephPool;
    … …  
    599690  delete _obsPool;
    600691  _obsPool  = new t_pppObsPool();
    601 }
     692
     693  // to delete all parameters
     694  delete _filter;
     695  _filter   = new t_pppFilter(_obsPool);
     696
     697}

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

    r7996 r8905  
    1010
    1111class bncAntex;
     12class t_pppRefSat;
    1213
    1314namespace BNC_PPP {
    … …  
    3435  const t_pppEphPool* ephPool() const {return _ephPool;}
    3536  const t_pppObsPool* obsPool() const {return _obsPool;}
    36   const bncAntex*  antex() const {return _antex;}
     37  const bncAntex*     antex() const {return _antex;}
    3738  const t_pppStation* staRover() const {return _staRover;}
    38   double           offGG() const {return _offGG;}
     39  double              offGG() const {return _offGG;}
    3940
    4041  std::ostringstream& log() {return *_log;}
    … …  
    5354  t_irc prepareObs(const std::vector<t_satObs*>& satObs,
    5455                   std::vector<t_pppSatObs*>& obsVector, bncTime& epoTime);
     56  bool  preparePseudoObs(std::vector<t_pppSatObs*>& obsVector);
    5557  t_irc cmpModel(t_pppStation* station, const ColumnVector& xyzc,
    5658                 std::vector<t_pppSatObs*>& obsVector);
    … …  
    7173  t_pppOptions*             _opt;
    7274  t_tides*                  _tides;
    73   t_loading*                _loading;
     75  bool                      _pseudoObsIono;
     76  int                       _num;
    7477};
    7578

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

    r8329 r8905  
    3434// Constructor
    3535////////////////////////////////////////////////////////////////////////////
    36 t_pppFilter::t_pppFilter() {
     36t_pppFilter::t_pppFilter(t_pppObsPool* obsPool) {
    3737  _parlist = 0;
     38  _numSat = 0;
     39  _obsPool = obsPool;
    3840}
    3941
    … …  
    4648// Process Single Epoch
    4749////////////////////////////////////////////////////////////////////////////
    48 t_irc t_pppFilter::processEpoch(t_pppObsPool* obsPool) {
    49 
     50t_irc t_pppFilter::processEpoch(int num) {
    5051  _numSat     = 0;
    5152  const double maxSolGap = 60.0;
    … …  
    5758  // Vector of all Observations
    5859  // --------------------------
    59   t_pppObsPool::t_epoch* epoch = obsPool->lastEpoch();
     60  t_pppObsPool::t_epoch* epoch = _obsPool->lastEpoch();
    6061  if (!epoch) {
    6162    return failure;
    … …  
    7576  string epoTimeStr = string(_epoTime);
    7677
     78  //--
    7779  // Set Parameters
    7880  // --------------
    … …  
    9193  for (unsigned ii = 0; ii < params.size(); ii++) {
    9294    const t_pppParam* par1 = params[ii];
    93 
    9495    _x0[ii] = par1->x0();
    95 
    9696    int iOld = par1->indexOld();
    9797    if (iOld < 0) {
    … …  
    113113  predictCovCrdPart(QFltOld);
    114114
     115
     116  // Pre-Process Satellite Systems separately
     117  // ----------------------------------------
     118  bool preProcessing = false;
     119  if (OPT->_obsModelType == OPT->DCMcodeBias ||
     120      OPT->_obsModelType == OPT->DCMphaseBias) {
     121    preProcessing = true;
     122    for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
     123      char system = OPT->systems()[iSys];
     124      t_prn refPrn = (_obsPool->getRefSatMapElement(system))->prn();
     125      vector<t_pppSatObs*> obsVector;
     126      for (unsigned jj = 0; jj < allObs.size(); jj++) {
     127        if (allObs[jj]->prn().system() == system) {
     128          obsVector.push_back(allObs[jj]);
     129        }
     130      }
     131      if (processSystem(OPT->LCs(system), obsVector, refPrn,
     132                        epoch->pseudoObsIono(), preProcessing) != success) {
     133        return failure;
     134      }
     135    }
     136  }
     137
    115138  // Process Satellite Systems separately
    116139  // ------------------------------------
    117140  for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
    118141    char system = OPT->systems()[iSys];
     142    t_prn refPrn = (_obsPool->getRefSatMapElement(system))->prn();
    119143    unsigned int num = 0;
    120144    vector<t_pppSatObs*> obsVector;
    … …  
    126150    }
    127151    LOG << epoTimeStr << " SATNUM " << system << ' ' << right << setw(2) << num << endl;
    128     if ( processSystem(OPT->LCs(system), obsVector) != success ) {
     152    if (processSystem(OPT->LCs(system), obsVector, refPrn,
     153                      epoch->pseudoObsIono(), preProcessing) != success) {
    129154      return failure;
    130155    }
    131156  }
    132 
    133   cmpDOP(allObs);
    134 
    135   _parlist->printResult(_epoTime, _QFlt, _xFlt);
    136   _lastEpoTimeOK = _epoTime;  // remember time of last successful epoch processing
     157  if (_obsPool->epoReProcessing()) {
     158    // set A1 und A2 abhängig von num
     159    // if num == 1 => A1
     160    // if num >  1 => A2
     161  }
     162  else {
     163    cmpDOP(allObs);
     164    _parlist->printResult(_epoTime, _QFlt, _xFlt);
     165    _lastEpoTimeOK = _epoTime;  // remember time of last successful epoch processing
     166  }
     167
    137168  return success;
    138169}
    … …  
    141172////////////////////////////////////////////////////////////////////////////
    142173t_irc t_pppFilter::processSystem(const vector<t_lc::type>& LCs,
    143                                  const vector<t_pppSatObs*>& obsVector) {
    144 
     174                                 const vector<t_pppSatObs*>& obsVector,
     175                                 const t_prn& refPrn,
     176                                 bool pseudoObsIonoAvailable,
     177                                 bool preProcessing) {
     178  qDebug() << "======t_pppFilter::processSystem=======";
    145179  LOG.setf(ios::fixed);
    146180
    147181  // Detect Cycle Slips
    148182  // ------------------
    149   if (detectCycleSlips(LCs, obsVector) != success) {
     183  if (detectCycleSlips(LCs, obsVector, refPrn, preProcessing) != success) {
    150184    return failure;
    151185  }
    152186
    153   ColumnVector            xSav       = _xFlt;
    154   SymmetricMatrix         QSav       = _QFlt;
    155   string                  epoTimeStr = string(_epoTime);
     187  unsigned usedLCs = LCs.size(); //qDebug() << "usedLCs: " << usedLCs;
     188  if (OPT->_pseudoObsIono && !pseudoObsIonoAvailable) {
     189    usedLCs -= 1;  // GIM not used
     190  }
     191  ColumnVector               xSav       = _xFlt;
     192  SymmetricMatrix            QSav       = _QFlt;
     193  string                     epoTimeStr = string(_epoTime);
    156194  const vector<t_pppParam*>& params     = _parlist->params();
    157   unsigned                maxObs     = obsVector.size() * LCs.size();
     195  unsigned                   maxObs     = obsVector.size() * usedLCs;
     196  //unsigned                   maxObs     = 2 * usedLCs;
     197
     198  if (OPT->_pseudoObsIono && pseudoObsIonoAvailable) { // stecDiff w.r.t refSat
     199    maxObs -= 1;
     200  }
     201  qDebug() << "par.size()      : " << _parlist->nPar() << " LCs.size(): " << usedLCs;
     202  qDebug() << "obsVector.size(): " << obsVector.size() << " maxObs    : " << maxObs;
    158203
    159204  // Outlier Detection Loop
    160205  // ----------------------
    161206  for (unsigned iOutlier = 0; iOutlier < maxObs; iOutlier++) {
     207    qDebug() << "iOutlier: " << iOutlier;
    162208
    163209    if (iOutlier > 0) {
    … …  
    168214    // First-Design Matrix, Terms Observed-Computed, Weight Matrix
    169215    // -----------------------------------------------------------
     216    qDebug() << "A(" << maxObs << "," << _parlist->nPar() << ")";
    170217    Matrix                AA(maxObs, _parlist->nPar());
    171218    ColumnVector          ll(maxObs);
    172219    DiagonalMatrix        PP(maxObs); PP = 0.0;
     220    //TETSPLOT
     221    for (unsigned iPar = 0; iPar < params.size(); iPar++) {
     222      const t_pppParam* par = params[iPar];
     223      cout << "  " << par->toString() <<endl;
     224    }
     225    cout << endl;
    173226
    174227    int iObs = -1;
    … …  
    176229    vector<t_lc::type>   usedTypes;
    177230    for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    178       t_pppSatObs* obs = obsVector[ii];
     231      t_pppSatObs* obs = obsVector[ii]; qDebug() << "SATELLITE: " << obs->prn().toString().c_str() << "isRef: " << obs->isReference();
    179232      if (!obs->outlier()) {
    180         for (unsigned jj = 0; jj < LCs.size(); jj++) {
     233        for (unsigned jj = 0; jj < usedLCs; jj++) {
    181234          const t_lc::type tLC = LCs[jj];
     235          if (tLC == t_lc::GIM &&  obs->isReference()) {continue;}
    182236          ++iObs;
    183237          usedObs.push_back(obs);
    … …  
    186240            const t_pppParam* par = params[iPar];
    187241            AA[iObs][iPar] = par->partial(_epoTime, obs, tLC);
     242            cout << setw(10) << par->partial(_epoTime, obs, tLC);
    188243          }
     244          cout << endl;
    189245          ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1));
    190246          PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC));
     247          //qDebug() << "ll[iObs]: " << ll[iObs];
     248          //qDebug() << "PP[iObs]: " << PP[iObs];
    191249        }
    192250      }
    … …  
    229287      t_pppSatObs* obs = usedObs[maxOutlierIndex];
    230288      t_pppParam* par = 0;
    231       LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' '
    232           << obs->prn().toString()                        << ' '
    233           << setw(8) << setprecision(4) << maxOutlier << endl;
     289      if (!preProcessing) {
     290        LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' '
     291            << obs->prn().toString()                        << ' '
     292            << setw(8) << setprecision(4) << maxOutlier << endl;
     293      }
    234294      for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    235295        t_pppParam* hlp = params[iPar];
    236         if (hlp->type() == t_pppParam::amb && hlp->prn()  == obs->prn() &&
     296        if (hlp->type() == t_pppParam::amb &&
     297            hlp->prn()  == obs->prn() &&
    237298            hlp->tLC()  == usedTypes[maxOutlierIndex]) {
    238299          par = hlp;
    239300        }
    240301      }
    241       if (par) {
     302      if      (par && preProcessing) {
     303        if (par->prn() == refPrn) {
     304          _obsPool->setEpoReProcessing(true);
     305        }
     306      }
     307      else if (par && !preProcessing) {
    242308        if (par->ambResetCandidate()) {
    243309          resetAmb(par->prn(), obsVector, &QSav, &xSav);
    … …  
    248314        }
    249315      }
    250       else {
     316      else if (!par && !preProcessing){
    251317        obs->setOutlier();
    252318      }
    … …  
    256322    // ---------------
    257323    else {
    258       for (unsigned jj = 0; jj < LCs.size(); jj++) {
    259         for (unsigned ii = 0; ii < usedObs.size(); ii++) {
    260           const t_lc::type tLC = usedTypes[ii];
    261           t_pppSatObs*        obs = usedObs[ii];
    262           if (tLC == LCs[jj]) {
    263             obs->setRes(tLC, vv[ii]);
    264             LOG << epoTimeStr << " RES "
    265                 << left << setw(3) << t_lc::toString(tLC) << right << ' '
    266                 << obs->prn().toString().substr(0,3) << ' '
    267                 << setw(8) << setprecision(4) << vv[ii] << endl;
     324      if (!preProcessing) {
     325        for (unsigned jj = 0; jj < LCs.size(); jj++) {
     326          for (unsigned ii = 0; ii < usedObs.size(); ii++) {
     327            const t_lc::type tLC = usedTypes[ii];
     328            t_pppSatObs*        obs = usedObs[ii];
     329            if (tLC == LCs[jj]) {
     330              obs->setRes(tLC, vv[ii]);
     331              LOG << epoTimeStr << " RES "
     332                  << left << setw(3) << t_lc::toString(tLC) << right << ' '
     333                  << obs->prn().toString().substr(0,3) << ' '
     334                  << setw(8) << setprecision(4) << vv[ii] << endl;
     335            }
    268336          }
    269337        }
    … …  
    279347////////////////////////////////////////////////////////////////////////////
    280348t_irc t_pppFilter::detectCycleSlips(const vector<t_lc::type>& LCs,
    281                                     const vector<t_pppSatObs*>& obsVector) {
     349                                    const vector<t_pppSatObs*>& obsVector,
     350                                    const t_prn& refPrn,
     351                                    bool preProcessing) {
    282352
    283353  const double            SLIP       = 20.0;  // slip threshold
    284354  string                  epoTimeStr = string(_epoTime);
    285   const vector<t_pppParam*>& params     = _parlist->params();
     355  const vector<t_pppParam*>& params  = _parlist->params();
    286356
    287357  for (unsigned ii = 0; ii < LCs.size(); ii++) {
    … …  
    294364        // ---------------------------------
    295365        bool slip = false;
    296 
    297366        if (obs->slip()) {
    298367          LOG << epoTimeStr  << "cycle slip set (obs) " << obs->prn().toString()  << endl;
    … …  
    312381          slip = true;
    313382        }
    314         _slips[obs->prn()]._biasJumpCounter = obs->biasJumpCounter();
    315383
    316384        // Slip Set
    317385        // --------
    318386        if (slip) {
    319           resetAmb(obs->prn(), obsVector);
    320         }
    321 
     387          if (preProcessing) {
     388            if (obs->prn() == refPrn) {
     389              _obsPool->setEpoReProcessing(true);
     390            }
     391          }
     392          else {
     393            resetAmb(obs->prn(), obsVector);
     394          }
     395        }
    322396        // Check Pre-Fit Residuals
    323397        // -----------------------
    … …  
    328402            AA[iPar] = par->partial(_epoTime, obs, tLC);
    329403          }
    330 
    331404          double ll = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA);
    332405          double vv = DotProduct(AA, _xFlt) - ll;
    333 
    334406          if (fabs(vv) > SLIP) {
    335             LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' '
    336                 << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl;
    337             resetAmb(obs->prn(), obsVector);
     407            if (preProcessing) {
     408              if (obs->prn() == refPrn) {
     409                _obsPool->setEpoReProcessing(true);
     410              }
     411            }
     412            else {
     413              LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' '
     414                  << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl;
     415              resetAmb(obs->prn(), obsVector);
     416            }
    338417          }
    339418        }

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

    r7928 r8905  
    1717class t_pppFilter {
    1818 public:
    19   t_pppFilter();
     19  t_pppFilter(t_pppObsPool* obsPool);
    2020  ~t_pppFilter();
    2121
    22   t_irc processEpoch(t_pppObsPool* obsPool);
     22  t_irc processEpoch(int num);
    2323
    2424  const ColumnVector&    x() const {return _xFlt;}
    … …  
    7474
    7575  t_irc processSystem(const std::vector<t_lc::type>& LCs,
    76                       const std::vector<t_pppSatObs*>& obsVector);
     76                      const std::vector<t_pppSatObs*>& obsVector,
     77                      const t_prn& refPrn,
     78                      bool pseudoObsIonoAvailable,
     79                      bool preProcessing);
    7780
    7881  t_irc detectCycleSlips(const std::vector<t_lc::type>& LCs,
    79                          const std::vector<t_pppSatObs*>& obsVector);
     82                         const std::vector<t_pppSatObs*>& obsVector,
     83                         const t_prn& refPrn,
     84                         bool preProcessing);
    8085
    8186  t_irc resetAmb(t_prn prn, const std::vector<t_pppSatObs*>& obsVector,
    … …  
    8893  bncTime         _epoTime;
    8994  t_pppParlist*   _parlist;
     95  t_pppObsPool*   _obsPool;
    9096  SymmetricMatrix _QFlt;
    9197  ColumnVector    _xFlt;

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

    r7643 r8905  
    2222// Constructor
    2323/////////////////////////////////////////////////////////////////////////////
    24 t_pppObsPool::t_epoch::t_epoch(const bncTime& epoTime, vector<t_pppSatObs*>& obsVector) {
    25   _epoTime   = epoTime;
     24t_pppObsPool::t_epoch::t_epoch(const bncTime& epoTime, vector<t_pppSatObs*>& obsVector,
     25                               bool pseudoObsIono) {
     26  _epoTime        = epoTime;
     27  _pseudoObsIono  = pseudoObsIono;
    2628  for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    2729    _obsVector.push_back(obsVector[ii]);
    … …  
    4850  }
    4951  _vTec = 0;
     52  _epoReProcessing = false;
    5053}
    5154
    … …  
    6467    _epochs.pop_front();
    6568  }
     69  clearRefSatMap();
    6670}
    6771
    … …  
    9296//
    9397/////////////////////////////////////////////////////////////////////////////
    94 void t_pppObsPool::putEpoch(const bncTime& epoTime, vector<t_pppSatObs*>& obsVector) {
     98void t_pppObsPool::putEpoch(const bncTime& epoTime, vector<t_pppSatObs*>& obsVector,
     99    bool pseudoObsIono) {
    95100  const unsigned MAXSIZE = 2;
    96   _epochs.push_back(new t_epoch(epoTime, obsVector));
     101
     102  if (epoReProcessing()) {
     103    delete _epochs.back();
     104    _epochs.pop_back();
     105  }
     106  _epochs.push_back(new t_epoch(epoTime, obsVector,pseudoObsIono));
    97107  if (_epochs.size() > MAXSIZE) {
    98108    delete _epochs.front();

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

    r7288 r8905  
    44#include <vector>
    55#include <deque>
     6#include <QMap>
    67#include "pppSatObs.h"
    78#include "bnctime.h"
     9#include "pppRefSat.h"
    810
    911namespace BNC_PPP {
    … …  
    1416  class t_epoch {
    1517   public:
    16     t_epoch(const bncTime& epoTime, std::vector<t_pppSatObs*>& obsVector);
     18    t_epoch(const bncTime& epoTime, std::vector<t_pppSatObs*>& obsVector,
     19            bool pseudoObsIono);
    1720    ~t_epoch();
    18     std::vector<t_pppSatObs*>& obsVector() {return _obsVector;}
     21          std::vector<t_pppSatObs*>& obsVector() {return _obsVector;}
    1922    const std::vector<t_pppSatObs*>& obsVector() const {return _obsVector;}
    2023    const bncTime& epoTime() const {return _epoTime;}
     24    bool pseudoObsIono() const {return _pseudoObsIono;}
    2125   private:
    22     bncTime                _epoTime;
     26    bncTime                   _epoTime;
    2327    std::vector<t_pppSatObs*> _obsVector;
     28    bool                      _pseudoObsIono;
    2429  };
    2530
    … …  
    3035  void putTec(t_vTec* _vTec);
    3136
    32   void putEpoch(const bncTime& epoTime, std::vector<t_pppSatObs*>& obsVector);
     37  void putEpoch(const bncTime& epoTime, std::vector<t_pppSatObs*>& obsVector,
     38                bool pseudoObs);
    3339
    3440  const t_satCodeBias* satCodeBias(const t_prn& prn) const {
    … …  
    4955  }
    5056
     57  void initRefSatMapElement(char system) {
     58    _refSatMap[system] = new t_pppRefSat();
     59  }
     60  void clearRefSatMap() {
     61    QMapIterator<char, t_pppRefSat*> it(_refSatMap);
     62    while (it.hasNext()) {
     63      it.next();
     64      delete it.value();
     65    }
     66    _refSatMap.clear();
     67  }
     68  t_pppRefSat* getRefSatMapElement(char system) {
     69    return _refSatMap[system];
     70  }
     71  QMap<char, t_pppRefSat*> getRefSatMap() {return _refSatMap;}
     72
     73  void setEpoReProcessing(bool epoReProcessing) {
     74    _epoReProcessing = epoReProcessing;
     75  }
     76  bool epoReProcessing() {return _epoReProcessing;}
     77
    5178 private:
    52   t_satCodeBias*       _satCodeBiases[t_prn::MAXPRN+1];
    53   t_satPhaseBias*      _satPhaseBiases[t_prn::MAXPRN+1];
    54   t_vTec*              _vTec;
    55   std::deque<t_epoch*> _epochs;
     79  t_satCodeBias*           _satCodeBiases[t_prn::MAXPRN+1];
     80  t_satPhaseBias*          _satPhaseBiases[t_prn::MAXPRN+1];
     81  t_vTec*                  _vTec;
     82  std::deque<t_epoch*>     _epochs;
     83  QMap<char, t_pppRefSat*> _refSatMap;
     84  bool                     _epoReProcessing;
    5685};
    5786

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

    r7988 r8905  
    9494     _noise   = OPT->_noiseTrp;
    9595     break;
     96    case ion:
     97     _epoSpec = false;
     98     _sigma0  = OPT->_aprSigIon;
     99     _noise   = OPT->_noiseIon;
     100     break;
     101   case cBias1:
     102   case cBias2:
     103     _epoSpec = false;
     104     _sigma0  = OPT->_aprSigCodeBias;
     105     _noise   = OPT->_noiseCodeBias;
     106     break;
     107   case pBias1:
     108   case pBias2:
     109     _epoSpec = false;
     110     _sigma0  = OPT->_aprSigPhaseBias;
     111     _noise   = OPT->_noisePhaseBias;
     112     break;
    96113  }
    97114}
    … …  
    106123////////////////////////////////////////////////////////////////////////////
    107124double t_pppParam::partial(const bncTime& /* epoTime */, const t_pppSatObs* obs,
    108                         const t_lc::type& tLC) const {
     125                           const t_lc::type& tLC) const {//qDebug() << "t_pppParam::partial: " << tLC;
    109126
    110127  // Special Case - Melbourne-Wuebbena
    … …  
    115132
    116133  const t_pppStation* sta  = PPP_CLIENT->staRover();
    117   ColumnVector     rhoV = sta->xyzApr() - obs->xc().Rows(1,3);
     134  ColumnVector        rhoV = sta->xyzApr() - obs->xc().Rows(1,3);
     135
     136  map<t_frequency::type, double> codeCoeff;
     137  map<t_frequency::type, double> phaseCoeff;
     138  map<t_frequency::type, double> ionoCoeff;
     139  obs->lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
    118140
    119141  switch (_type) {
    120142  case crdX:
    121     return (sta->xyzApr()[0] - obs->xc()[0]) / rhoV.norm_Frobenius();
     143    if (tLC == t_lc::GIM) {return 0.0;}
     144    return (sta->xyzApr()[0] - obs->xc()[0]) / rhoV.NormFrobenius();
    122145  case crdY:
    123     return (sta->xyzApr()[1] - obs->xc()[1]) / rhoV.norm_Frobenius();
     146    if (tLC == t_lc::GIM) {return 0.0;}
     147    return (sta->xyzApr()[1] - obs->xc()[1]) / rhoV.NormFrobenius();
    124148  case crdZ:
    125     return (sta->xyzApr()[2] - obs->xc()[2]) / rhoV.norm_Frobenius();
     149    if (tLC == t_lc::GIM) {return 0.0;}
     150    return (sta->xyzApr()[2] - obs->xc()[2]) / rhoV.NormFrobenius();
    126151  case clkR:
     152    if (tLC == t_lc::GIM) {return 0.0;}
    127153    return 1.0;
    128154  case offGG:
     155    if (tLC == t_lc::GIM) {return 0.0;}
    129156    return (obs->prn().system() == 'R') ? 1.0 : 0.0;
    130157  case amb:
     158    if      (tLC == t_lc::GIM) {return 0.0;}
     159    else if ((OPT->_obsModelType == OPT->IF)     ||
     160             (OPT->_obsModelType == OPT->PPPRTK) ||
     161             (OPT->_obsModelType == OPT->UncombPPP) ||
     162             (OPT->_obsModelType == OPT->DCMcodeBias  && !obs->isReference()) ||
     163             (OPT->_obsModelType == OPT->DCMphaseBias && !obs->isReference())   ) {
     164      if (obs->prn() == _prn) {
     165        if      (tLC == _tLC) {
     166          return (obs->lambda(tLC));
     167        }
     168        else if (tLC == t_lc::lIF && _tLC == t_lc::MW) {
     169          return obs->lambda(t_lc::lIF) * obs->lambda(t_lc::MW) / obs->lambda(t_lc::l2);
     170        }
     171        else {
     172          if      (_tLC == t_lc::l1) {
     173            return obs->lambda(t_lc::l1) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
     174          }
     175          else if (_tLC == t_lc::l2) {
     176            return obs->lambda(t_lc::l2) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
     177          }
     178        }
     179      }
     180    }
     181    break;
     182  case trp:
     183    if (tLC == t_lc::GIM) {return 0.0;}
     184    return  1.0 / sin(obs->eleSat());
     185  case ion:
     186//    qDebug() << "refPrn: " << _refPrn.toString().c_str();
    131187    if (obs->prn() == _prn) {
    132       if      (tLC == _tLC) {
    133         return (obs->lambda(tLC));
    134       }
    135       else if (tLC == t_lc::lIF && _tLC == t_lc::MW) {
    136         return obs->lambda(t_lc::lIF) * obs->lambda(t_lc::MW) / obs->lambda(t_lc::l2);
    137       }
    138       else {
    139         map<t_frequency::type, double> codeCoeff;
    140         map<t_frequency::type, double> phaseCoeff;
    141         obs->lcCoeff(tLC, codeCoeff, phaseCoeff);
    142         if      (_tLC == t_lc::l1) {
    143           return obs->lambda(t_lc::l1) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
    144         }
    145         else if (_tLC == t_lc::l2) {
    146           return obs->lambda(t_lc::l2) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
    147         }
    148       }
    149     }
    150     return 0.0;
    151   case trp:
    152     return 1.0 / sin(obs->eleSat());
    153   }
    154 
     188      if      (tLC == t_lc::c1) {
     189        return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::c1)];
     190      }
     191      else if (tLC == t_lc::c2) {
     192        return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::c2)];
     193      }
     194      else if (tLC == t_lc::l1) {
     195        return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
     196      }
     197      else if (tLC == t_lc::l2) {
     198        return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
     199      }
     200      else if (tLC == t_lc::GIM) {
     201        return -1.0;
     202      }
     203    }
     204    if (tLC == t_lc::GIM && _prn == _refPrn) {
     205      return  1.0;
     206    }
     207    break;
     208  case cBias1:
     209    if  (tLC == t_lc::c1) {
     210      return  1.0;
     211    }
     212    else {
     213      return 0.0;
     214    }
     215    break;
     216  case cBias2:
     217     if (tLC == t_lc::c2) {
     218      return  1.0;
     219    }
     220    else {
     221      return 0.0;
     222    }
     223    break;
     224  case pBias1:
     225    if  (tLC == t_lc::l1) {
     226      return  1.0;
     227    }
     228    else {
     229      return 0.0;
     230    }
     231    break;
     232  case pBias2:
     233    if  (tLC == t_lc::l2) {
     234      return  1.0;
     235    }
     236    else {
     237      return 0.0;
     238    }
     239  }
    155240  return 0.0;
    156241}
    … …  
    171256    break;
    172257  case clkR:
    173     ss << "CLK        ";
    174     break;
    175   case amb:
    176     ss << "AMB " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
     258    ss << "REC_CLK    ";
    177259    break;
    178260  case offGG:
    … …  
    181263  case trp:
    182264    ss << "TRP        ";
     265    break;
     266  case amb:
     267    ss << "AMB  " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
     268    break;
     269  case ion:
     270    ss << "ION  " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
     271    break;
     272  case cBias1:
     273  case cBias2:
     274  case pBias1:
     275  case pBias2:
     276    ss << "BIAS " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << "REC";
    183277    break;
    184278  }
    … …  
    201295//
    202296////////////////////////////////////////////////////////////////////////////
    203 t_irc t_pppParlist::set(const bncTime& epoTime, const std::vector<t_pppSatObs*>& obsVector) {
     297t_irc t_pppParlist::set(const bncTime& epoTime,
     298    const std::vector<t_pppSatObs*>& obsVector) {
    204299
    205300  // Remove some Parameters
    … …  
    271366  required.push_back(new t_pppParam(t_pppParam::clkR, t_prn(), t_lc::dummy));
    272367
    273   // GPS-Glonass Clock Offset
     368  // GPS-GLONASS Clock Offset
    274369  // ------------------------
    275370  if (OPT->useSystem('R')) {
    … …  
    283378  }
    284379
     380  // Ionosphere
     381  // ----------
     382  if (OPT->_obsModelType == OPT->UncombPPP    ||
     383      OPT->_obsModelType == OPT->DCMcodeBias  ||
     384      OPT->_obsModelType == OPT->DCMphaseBias   ) {
     385    for (unsigned jj = 0; jj < obsVector.size(); jj++) {
     386      const t_pppSatObs*        satObs = obsVector[jj];
     387      required.push_back(new t_pppParam(t_pppParam::ion, satObs->prn(), t_lc::dummy));
     388    }
     389  }
     390
    285391  // Ambiguities
    286392  // -----------
    287393  for (unsigned jj = 0; jj < obsVector.size(); jj++) {
    288394    const t_pppSatObs*        satObs = obsVector[jj];
    289     const vector<t_lc::type>& ambLCs = OPT->ambLCs(satObs->prn().system());
    290     for (unsigned ii = 0; ii < ambLCs.size(); ii++) {
    291       required.push_back(new t_pppParam(t_pppParam::amb, satObs->prn(), ambLCs[ii], &obsVector));
    292     }
     395    if ((OPT->_obsModelType == OPT->IF)        ||
     396        (OPT->_obsModelType == OPT->PPPRTK)    ||
     397        (OPT->_obsModelType == OPT->UncombPPP) ||
     398        (OPT->_obsModelType == OPT->DCMcodeBias  && !satObs->isReference()) ||
     399        (OPT->_obsModelType == OPT->DCMphaseBias && !satObs->isReference())   ) {
     400      const vector<t_lc::type>& ambLCs = OPT->ambLCs(satObs->prn().system());
     401      for (unsigned ii = 0; ii < ambLCs.size(); ii++) {
     402        required.push_back(new t_pppParam(t_pppParam::amb, satObs->prn(), ambLCs[ii], &obsVector));
     403      }
     404    }
     405  }
     406
     407  // Receiver Code Biases
     408  // --------------------
     409  if (OPT->_obsModelType == OPT->DCMcodeBias) {
     410    required.push_back(new t_pppParam(t_pppParam::cBias1, t_prn(), t_lc::c1));
     411    required.push_back(new t_pppParam(t_pppParam::cBias2, t_prn(), t_lc::c2));
     412  }
     413
     414  // Receiver Phase Biases
     415  // ---------------------
     416  if ((OPT->_obsModelType == OPT->DCMphaseBias) ||
     417      (OPT->_obsModelType == OPT->PPPRTK)     ) {
     418    required.push_back(new t_pppParam(t_pppParam::pBias1, t_prn(), t_lc::l1));
     419    required.push_back(new t_pppParam(t_pppParam::pBias2, t_prn(), t_lc::l2));
    293420  }
    294421

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

    r7237 r8905  
    1414class t_pppParam {
    1515 public:
    16   enum e_type {crdX, crdY, crdZ, clkR, amb, offGG, trp};
     16  enum e_type {crdX, crdY, crdZ, clkR, amb, offGG, trp, ion, cBias1, cBias2, pBias1, pBias2};
    1717
    1818  t_pppParam(e_type type, const t_prn& prn, t_lc::type tLC,
    … …  
    2020
    2121  ~t_pppParam();
     22
    2223  e_type type() const {return _type;}
    2324  double x0()  const {return _x0;}
    24   double partial(const bncTime& epoTime, const t_pppSatObs* obs,
    25                  const t_lc::type& tLC) const;
     25  double partial(const bncTime& epoTime, const t_pppSatObs* obs, const t_lc::type& tLC) const;
    2626  bool   epoSpec() const {return _epoSpec;}
    2727  bool   isEqual(const t_pppParam* par2) const {
    … …  
    5151  unsigned ambNumEpo() const           {return _ambInfo ? _ambInfo->_numEpo : 0;}
    5252  void     stepAmbNumEpo()             {if (_ambInfo) _ambInfo->_numEpo += 1;}
     53  void     setRefPrn(t_prn prn)        {_refPrn = prn;}
    5354
    5455  static bool sortFunction(const t_pppParam* p1, const t_pppParam* p2) {
    … …  
    7879    unsigned _numEpo;
    7980  };
    80   e_type     _type;
    81   t_prn      _prn;
    82   t_lc::type _tLC;
    83   double     _x0;
    84   bool       _epoSpec;
    85   int        _indexOld;
    86   int        _indexNew;
    87   double     _sigma0;
    88   double     _noise;
    89   t_ambInfo* _ambInfo;
    90   bncTime    _lastObsTime;
    91   bncTime    _firstObsTime;
     81  e_type       _type;
     82  t_prn        _prn;
     83  t_prn        _refPrn;
     84  t_lc::type   _tLC;
     85  double       _x0;
     86  bool         _epoSpec;
     87  int          _indexOld;
     88  int          _indexNew;
     89  double       _sigma0;
     90  double       _noise;
     91  t_ambInfo*   _ambInfo;
     92  bncTime      _lastObsTime;
     93  bncTime      _firstObsTime;
    9294};
    9395
    … …  
    100102  unsigned nPar() const {return _params.size();}
    101103  const std::vector<t_pppParam*>& params() const {return _params;}
    102   std::vector<t_pppParam*>& params() {return _params;}
    103   void printResult(const bncTime& epoTime, const SymmetricMatrix& QQ, 
     104        std::vector<t_pppParam*>& params()      {return _params;}
     105  void printResult(const bncTime& epoTime, const SymmetricMatrix& QQ,
    104106                   const ColumnVector& xx) const;
     107
    105108 private:
    106109  std::vector<t_pppParam*> _params;

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

    r8619 r8905  
    3535////////////////////////////////////////////////////////////////////////////
    3636t_pppSatObs::t_pppSatObs(const t_satObs& pppSatObs) {
    37   _prn     = pppSatObs._prn;
    38   _time    = pppSatObs._time;
    39   _outlier = false;
    40   _valid   = true;
     37  _prn        = pppSatObs._prn;
     38  _time       = pppSatObs._time;
     39  _outlier    = false;
     40  _valid      = true;
     41  _reference  = false;
     42  _stecRefSat = 0.0;
     43  _stecSat    = 0.0;
    4144  for (unsigned ii = 0; ii < t_frequency::max; ii++) {
    4245    _obs[ii] = 0;
    … …  
    6164  // Select pseudoranges and phase observations
    6265  // ------------------------------------------
    63   const string preferredAttrib = "CWPXI_";
    64   //const string preferredAttrib = "G:12&PWCSLXYN G:5&IQX R:12&PC R:3&IQX E:16&BCX E:578&IQX J:1&SLXCZ J:26&SLX J:5&IQX C:IQX I:ABCX S:1&C S:5&IQX";
     66  const string preferredAttrib = "G:12&PWCSLXYN G:5&IQX R:12&PC R:3&IQX E:16&BCX E:578&IQX J:1&SLXCZ J:26&SLX J:5&IQX C:IQX I:ABCX S:1&C S:5&IQX";
    6567
    6668  for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
    … …  
    9092  for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) {
    9193    t_lc::type tLC = OPT->LCs(_prn.system())[ii];
     94    if (tLC == t_lc::GIM) {continue;}
    9295    if (!isValid(tLC)) {
    9396      _valid = false;
    … …  
    121124    ColumnVector dx = _xcSat - satPosOld;
    122125    dx[3] *= t_CST::c;
    123     if (dx.norm_Frobenius() < 1.e-4) {
     126    if (dx.NormFrobenius() < 1.e-4) {
    124127      totOK = true;
    125128      break;
    … …  
    140143void t_pppSatObs::lcCoeff(t_lc::type tLC,
    141144                          map<t_frequency::type, double>& codeCoeff,
    142                           map<t_frequency::type, double>& phaseCoeff) const {
     145                          map<t_frequency::type, double>& phaseCoeff,
     146                          map<t_frequency::type, double>& ionoCoeff) const {
    143147
    144148  codeCoeff.clear();
    145149  phaseCoeff.clear();
     150  ionoCoeff.clear();
    146151
    147152  double f1 = t_CST::freq(_fType1, _channel);
    148153  double f2 = t_CST::freq(_fType2, _channel);
     154  double f1GPS = t_CST::freq(t_frequency::G1, 0);
    149155
    150156  switch (tLC) {
    151157  case t_lc::l1:
    152     phaseCoeff[_fType1] = 1.0;
     158    phaseCoeff[_fType1] =  1.0;
     159    ionoCoeff [_fType1] = -1.0 * pow(f1GPS, 2) / pow(f1, 2);
    153160    return;
    154161  case t_lc::l2:
    155     phaseCoeff[_fType2] = 1.0;
     162    phaseCoeff[_fType2] =  1.0;
     163    ionoCoeff [_fType2] = -1.0 * pow(f1GPS, 2) / pow(f2, 2);
    156164    return;
    157165  case t_lc::lIF:
    … …  
    167175  case t_lc::CL:
    168176    phaseCoeff[_fType1] =  0.5;
    169     codeCoeff[_fType1] =  0.5;
     177    codeCoeff [_fType1] =  0.5;
    170178    return;
    171179  case t_lc::c1:
    172180    codeCoeff[_fType1] = 1.0;
     181    ionoCoeff[_fType1] = pow(f1GPS, 2) / pow(f1, 2);
    173182    return;
    174183  case t_lc::c2:
    175184    codeCoeff[_fType2] = 1.0;
     185    ionoCoeff[_fType2] = pow(f1GPS, 2) / pow(f2, 2);
    176186    return;
    177187  case t_lc::cIF:
    … …  
    179189    codeCoeff[_fType2] = -f2 * f2 / (f1 * f1 - f2 * f2);
    180190    return;
     191  case t_lc::GIM:
    181192  case t_lc::dummy:
    182193  case t_lc::maxLc:
    … …  
    190201  bool valid = true;
    191202  obsValue(tLC, &valid);
     203  //qDebug() << "tLC: " << tLC << "  valid: " << valid;
    192204  return valid;
    193205}
    … …  
    195207////////////////////////////////////////////////////////////////////////////
    196208double t_pppSatObs::obsValue(t_lc::type tLC, bool* valid) const {
     209
     210  double retVal = 0.0;
     211  if (valid) *valid = true;
     212
     213  // Pseudo observations
     214  if (tLC == t_lc::GIM) {
     215    if (_stecRefSat == 0.0 || _stecSat == 0.0) {
     216      if (valid) *valid = false;
     217      return 0.0;
     218    }
     219    else {
     220      return _stecRefSat - _stecSat;
     221    }
     222  }
    197223
    198224  map<t_frequency::type, double> codeCoeff;
    199225  map<t_frequency::type, double> phaseCoeff;
    200   lcCoeff(tLC, codeCoeff, phaseCoeff);
    201 
    202   double retVal = 0.0;
    203   if (valid) *valid = true;
     226  map<t_frequency::type, double> ionoCoeff;
     227  lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
    204228
    205229  map<t_frequency::type, double>::const_iterator it;
     230
     231  // Code observations
    206232  for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
    207233    t_frequency::type tFreq = it->first;
    … …  
    214240    }
    215241  }
     242  // Phase observations
    216243  for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
    217244    t_frequency::type tFreq = it->first;
    … …  
    224251    }
    225252  }
    226 
    227253  return retVal;
    228254}
    … …  
    258284double t_pppSatObs::sigma(t_lc::type tLC) const {
    259285
     286  double retVal = 0.0;
    260287  map<t_frequency::type, double> codeCoeff;
    261288  map<t_frequency::type, double> phaseCoeff;
    262   lcCoeff(tLC, codeCoeff, phaseCoeff);
    263 
    264   double retVal = 0.0;
     289  map<t_frequency::type, double> ionoCoeff;
     290  lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
     291
     292  if (tLC == t_lc::GIM) {
     293    retVal = OPT->_sigmaGIMdiff * OPT->_sigmaGIMdiff;
     294  }
    265295
    266296  map<t_frequency::type, double>::const_iterator it;
    267   for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
     297  for (it = codeCoeff.begin(); it != codeCoeff.end(); it++)   {//qDebug() << "codeCoeff : " << t_frequency::toString(it->first).c_str() << ": " << it->second;
    268298    retVal += it->second * it->second * OPT->_sigmaC1 * OPT->_sigmaC1;
    269299  }
    270   for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
     300
     301  for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {//qDebug() << "phaseCoeff: " << t_frequency::toString(it->first).c_str()  << ": " << it->second;
    271302    retVal += it->second * it->second * OPT->_sigmaL1 * OPT->_sigmaL1;
    272303  }
    … …  
    295326//
    296327////////////////////////////////////////////////////////////////////////////
    297 double t_pppSatObs::maxRes(t_lc::type tLC) const {
     328double t_pppSatObs::maxRes(t_lc::type tLC) const {//qDebug() << "t_pppSatObs::maxRes(t_lc::type tLC)";
     329  double retVal = 0.0;
    298330
    299331  map<t_frequency::type, double> codeCoeff;
    300332  map<t_frequency::type, double> phaseCoeff;
    301   lcCoeff(tLC, codeCoeff, phaseCoeff);
    302 
    303   double retVal = 0.0;
     333  map<t_frequency::type, double> ionoCoeff;
     334  lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
    304335
    305336  map<t_frequency::type, double>::const_iterator it;
    306   for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
     337  for (it = codeCoeff.begin(); it != codeCoeff.end(); it++)   {//qDebug() << "codeCoeff: " << it->first << ": " << it->second;
    307338    retVal += it->second * it->second * OPT->_maxResC1 * OPT->_maxResC1;
    308339  }
    309   for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
     340  for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {//qDebug() << "phaseCoeff: " << it->first << ": " << it->second;
    310341    retVal += it->second * it->second * OPT->_maxResL1 * OPT->_maxResL1;
    311342  }
    312 
     343  if (tLC == t_lc::GIM) {
     344    retVal = 3 * (OPT->_sigmaGIMdiff * OPT->_sigmaGIMdiff);
     345  }
    313346  return sqrt(retVal);
    314347}
    … …  
    326359  // ------------------------------
    327360  ColumnVector rSat = _xcSat.Rows(1,3);
    328   ColumnVector rhoV = rSat - station->xyzApr();
    329   _model._rho = rhoV.norm_Frobenius();
     361  ColumnVector rRec = station->xyzApr();
     362  ColumnVector rhoV = rSat - rRec;
     363  _model._rho = rhoV.NormFrobenius();
    330364
    331365  ColumnVector vSat = _vvSat;
    … …  
    334368  xyz2neu(station->ellApr().data(), rhoV.data(), neu.data());
    335369
    336   _model._eleSat = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / _model._rho );
     370  _model._eleSat = acos(sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / _model._rho);
    337371  if (neu[2] < 0) {
    338372    _model._eleSat *= -1.0;
    … …  
    354388  Omega[1] = 0.0;
    355389  Omega[2] = t_CST::omega / t_CST::c;
    356   _model._sagnac = DotProduct(Omega, crossproduct(rSat, station->xyzApr()));
     390  _model._sagnac = DotProduct(Omega, crossproduct(rSat, rRec));
    357391
    358392  // Antenna Eccentricity
    … …  
    373407  // Tropospheric Delay
    374408  // ------------------
    375   _model._tropo = t_tropo::delay_saast(station->xyzApr(), _model._eleSat);
     409  _model._tropo = t_tropo::delay_saast(rRec, _model._eleSat);
    376410
    377411  // Code Biases
    … …  
    392426  // Phase Biases
    393427  // -----------
    394   // TODO: consideration of fix indicators and jump counter
    395428  const t_satPhaseBias* satPhaseBias = PPP_CLIENT->obsPool()->satPhaseBias(_prn);
    396429  double yaw = 0.0;
    397430  bool ssr = false;
    398431  if (satPhaseBias) {
    399     yaw = satPhaseBias->_yaw;
     432    double dt = station->epochTime() - satPhaseBias->_time;
     433    if (satPhaseBias->_updateInt) {
     434      dt -= (0.5 * ssrUpdateInt[satPhaseBias->_updateInt]);
     435    }
     436    yaw = satPhaseBias->_yaw + satPhaseBias->_yawRate * dt;
    400437    ssr = true;
    401438    for (unsigned ii = 0; ii < satPhaseBias->_bias.size(); ii++) {
    … …  
    414451  _model._windUp = station->windUp(_time, _prn, rSat, ssr, yaw, vSat) ;
    415452
     453  // Relativistic effect due to earth gravity
     454  // ----------------------------------------
     455  double a = rSat.NormFrobenius() + rRec.NormFrobenius();
     456  double b = (rSat - rRec).NormFrobenius();
     457  double gm = 3.986004418e14; // m3/s2
     458  _model._rel = 2 * gm / t_CST::c / t_CST::c * log((a + b) / (a - b));
    416459
    417460  // Tidal Correction
    418461  // ----------------
    419   _model._tide = -DotProduct(station->tideDspl(), rhoV) / _model._rho;
     462  _model._tideEarth = -DotProduct(station->tideDsplEarth(), rhoV) / _model._rho;
     463  _model._tideOcean = -DotProduct(station->tideDsplOcean(), rhoV) / _model._rho;
    420464
    421465  // Ionospheric Delay
    … …  
    429473    }
    430474  }
     475
    431476  if (vTecUsage && vTec) {
    432     double stec = station->stec(vTec, _signalPropagationTime, rSat);
     477    double stec  = station->stec(vTec, _signalPropagationTime, rSat);
     478    double f1GPS = t_CST::freq(t_frequency::G1, 0);
    433479    for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
    434       t_frequency::type frqType = static_cast<t_frequency::type>(iFreq);
    435       _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(t_CST::freq(frqType, _channel), 2) * stec;
    436     }
    437   }
    438 
    439   // Relativistic effect due to earth gravity
    440   // ----------------------------------------
    441   // TODO
    442 
    443   // Ocean Loading
    444   // -------------
    445   // TODO
     480      if (OPT->_pseudoObsIono) { // DCMcodeBias, DCMphaseBias
     481        // For scaling the slant ionospheric delays the trick is to be consistent with units!
     482        // The conversion of TECU into meters requires the frequency of the signal.
     483        // Hence, GPS L1 frequency is used for all systems. The same is true for mu_i in lcCoeff().
     484        _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(f1GPS, 2) * stec;
     485      }
     486      else { // PPP-RTK
     487        t_frequency::type frqType = static_cast<t_frequency::type>(iFreq);
     488        _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(t_CST::freq(frqType, _channel), 2) * stec;
     489      }
     490    }
     491  }
    446492
    447493  // Set Model Set Flag
    … …  
    449495  _model._set = true;
    450496
    451   //printModel();
     497  printModel();
    452498
    453499  return success;
    … …  
    457503////////////////////////////////////////////////////////////////////////////
    458504void t_pppSatObs::printModel() const {
    459 
    460   LOG.setf(ios::fixed);
    461   LOG << "MODEL for Satellite " << _prn.toString() << endl
    462       << "RHO:        " << setw(12) << setprecision(3) << _model._rho     << endl
    463       << "ELE:        " << setw(12) << setprecision(3) << _model._eleSat * 180.0 / M_PI << endl
    464       << "AZI:        " << setw(12) << setprecision(3) << _model._azSat  * 180.0 / M_PI << endl
    465       << "SATCLK:     " << setw(12) << setprecision(3) << _model._satClkM << endl
    466       << "RECCLK:     " << setw(12) << setprecision(3) << _model._recClkM << endl
    467       << "SAGNAC:     " << setw(12) << setprecision(3) << _model._sagnac  << endl
    468       << "ANTECC:     " << setw(12) << setprecision(3) << _model._antEcc  << endl
    469       << "TROPO:      " << setw(12) << setprecision(3) << _model._tropo   << endl
    470       << "WINDUP:     " << setw(12) << setprecision(3) << _model._windUp  << endl
    471       << "TIDES:      " << setw(12) << setprecision(3) << _model._tide    << endl;
     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
     511      << "RHO           : " << setw(12) << setprecision(3) << _model._rho              << endl
     512      << "ELE           : " << setw(12) << setprecision(3) << _model._eleSat * RHO_DEG << endl
     513      << "AZI           : " << setw(12) << setprecision(3) << _model._azSat  * RHO_DEG << endl
     514      << "SATCLK        : " << setw(12) << setprecision(3) << _model._satClkM          << endl
     515      << "RECCLK        : " << setw(12) << setprecision(3) << _model._recClkM          << endl
     516      << "SAGNAC        : " << setw(12) << setprecision(3) << _model._sagnac           << endl
     517      << "ANTECC        : " << setw(12) << setprecision(3) << _model._antEcc           << endl
     518      << "TROPO         : " << setw(12) << setprecision(3) << _model._tropo            << endl
     519      << "WINDUP        : " << setw(12) << setprecision(3) << _model._windUp           << endl
     520      << "REL           : " << setw(12) << setprecision(3) << _model._rel              << endl
     521      << "EARTH TIDES   : " << setw(12) << setprecision(3) << _model._tideEarth        << endl
     522      << "OCEAN TIDES   : " << setw(12) << setprecision(3) << _model._tideOcean        << endl
     523      << endl
     524      << "FREQUENCY DEPENDENT CORRECTIONS:" << endl
     525      << "-------------------------------" << endl;
    472526  for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
    473527    if (_obs[iFreq]) {
    474528      string frqStr = t_frequency::toString(t_frequency::type(iFreq));
    475529      if (_prn.system() == frqStr[0]) {
    476       LOG << "PCO           : " << frqStr << setw(12) << setprecision(3) << _model._antPCO[iFreq]    << endl
    477           << "BIAS CODE     : " << frqStr << setw(12) << setprecision(3) << _model._codeBias[iFreq]  << endl
    478           << "BIAS PHASE    : " << frqStr << setw(12) << setprecision(3) << _model._phaseBias[iFreq]  << endl
    479           << "IONO CODEDELAY: " << frqStr << setw(12) << setprecision(3) << _model._ionoCodeDelay[iFreq] << endl;
    480       }
    481     }
    482   }
     530      cout << "PCO           : " << frqStr << setw(12) << setprecision(3) << _model._antPCO[iFreq]       << endl
     531           << "BIAS CODE     : " << frqStr << setw(12) << setprecision(3) << _model._codeBias[iFreq]     << endl
     532           << "BIAS PHASE    : " << frqStr << setw(12) << setprecision(3) << _model._phaseBias[iFreq]    << endl
     533           << "IONO CODEDELAY: " << frqStr << setw(12) << setprecision(3) << _model._ionoCodeDelay[iFreq]<< endl;
     534      }
     535    }
     536  }
     537}
     538
     539//
     540////////////////////////////////////////////////////////////////////////////
     541void t_pppSatObs::printObsMinusComputed() const {
     542// TODO: cout should be LOG
     543  cout.setf(ios::fixed);
     544  cout << "\nOBS-COMP for Satellite " << _prn.toString() << (isReference() ? " (Reference Satellite)" : "") << endl
     545       << "========================== " << endl;
    483546  for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) {
    484547    t_lc::type tLC = OPT->LCs(_prn.system())[ii];
    485     LOG << "OBS-CMP " << t_lc::toString(tLC) << ": " << _prn.toString() << " "
     548    cout << "OBS-CMP " << setw(4) << t_lc::toString(tLC) << ": " << _prn.toString() << " "
    486549        << setw(12) << setprecision(3) << obsValue(tLC) << " "
    487550        << setw(12) << setprecision(3) << cmpValue(tLC) << " "
    488551        << setw(12) << setprecision(3) << obsValue(tLC) - cmpValue(tLC) << endl;
    489 
    490   }
    491   LOG << "OBS-CMP MW: " << _prn.toString() << " "
    492       << setw(12) << setprecision(3) << obsValue(t_lc::MW) << " "
    493       << setw(12) << setprecision(3) << cmpValue(t_lc::MW) << " "
    494       << setw(12) << setprecision(3) << obsValue(t_lc::MW) - cmpValue(t_lc::MW) << endl;
    495 }
     552  }
     553}
     554
    496555
    497556//
    … …  
    504563////////////////////////////////////////////////////////////////////////////
    505564double t_pppSatObs::cmpValue(t_lc::type tLC) const {
    506 
    507   if (!isValid(tLC)) {
    508     return 0.0;
    509   }
    510 
    511   // Non-Dispersive Part
    512   // -------------------
    513   double nonDisp = _model._rho    + _model._recClkM - _model._satClkM
    514                  + _model._sagnac + _model._antEcc  + _model._tropo
    515                  + _model._tide;
    516 
    517   // Add Dispersive Part
    518   // -------------------
    519   map<t_frequency::type, double> codeCoeff;
    520   map<t_frequency::type, double> phaseCoeff;
    521   lcCoeff(tLC, codeCoeff, phaseCoeff);
    522 
    523   double dispPart = 0.0;
    524 
    525   map<t_frequency::type, double>::const_iterator it;
    526   for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
    527     t_frequency::type tFreq = it->first;
    528     dispPart += it->second * (_model._antPCO[tFreq] - _model._codeBias[tFreq] +
    529                               _model._ionoCodeDelay[tFreq]);
    530   }
    531   for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
    532     t_frequency::type tFreq = it->first;
    533     dispPart += it->second * (_model._antPCO[tFreq] - _model._phaseBias[tFreq] +
    534                               _model._windUp * t_CST::lambda(tFreq, _channel)  -
    535                               _model._ionoCodeDelay[tFreq]);
    536   }
    537 
    538     return nonDisp + dispPart;
     565  double cmpValue;
     566
     567  if      (!isValid(tLC)) {
     568    cmpValue =  0.0;
     569  }
     570  else if (tLC == t_lc::GIM) {
     571    cmpValue = 0.0;
     572  }
     573  else {
     574    // Non-Dispersive Part
     575    // -------------------
     576    double nonDisp = _model._rho
     577                   + _model._recClkM   - _model._satClkM
     578                   + _model._sagnac    + _model._antEcc    + _model._tropo
     579                   + _model._tideEarth + _model._tideOcean + _model._rel;
     580
     581    // Add Dispersive Part
     582    // -------------------
     583    double dispPart = 0.0;
     584    map<t_frequency::type, double> codeCoeff;
     585    map<t_frequency::type, double> phaseCoeff;
     586    map<t_frequency::type, double> ionoCoeff;
     587    lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
     588    map<t_frequency::type, double>::const_iterator it;
     589    for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
     590      t_frequency::type tFreq = it->first;
     591      dispPart += it->second * (_model._antPCO[tFreq] - _model._codeBias[tFreq]);
     592      if (OPT->PPPRTK) {
     593        dispPart += it->second * (_model._ionoCodeDelay[tFreq]);
     594      }
     595    }
     596    for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
     597      t_frequency::type tFreq = it->first;
     598      dispPart += it->second * (_model._antPCO[tFreq] - _model._phaseBias[tFreq] +
     599                                _model._windUp * t_CST::lambda(tFreq, _channel));
     600      if (OPT->PPPRTK) {
     601        dispPart += it->second * (- _model._ionoCodeDelay[tFreq]);
     602      }
     603    }
     604    cmpValue = nonDisp + dispPart;
     605  }
     606
     607  return cmpValue;
    539608}
    540609
    … …  
    556625  }
    557626}
     627
     628//
     629////////////////////////////////////////////////////////////////////////////
     630void  t_pppSatObs::setPseudoObsIono(t_frequency::type freq, double stecRefSat) {
     631  _stecSat    = _model._ionoCodeDelay[freq];
     632  _stecRefSat = stecRefSat;
     633}

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

    r8619 r8905  
    1919  bool                isValid() const {return _valid;};
    2020  bool                isValid(t_lc::type tLC) const;
     21  bool                isReference() const {return _reference;};
     22  void                setAsReference() {_reference = true;};
    2123  const t_prn&        prn() const {return _prn;}
    2224  const ColumnVector& xc() const {return _xcSat;}
    … …  
    3133  bool                modelSet() const {return _model._set;}
    3234  void                printModel() const;
     35  void                printObsMinusComputed() const;
    3336  void                lcCoeff(t_lc::type tLC,
    3437                              std::map<t_frequency::type, double>& codeCoeff,
    35                               std::map<t_frequency::type, double>& phaseCoeff) const;
     38                              std::map<t_frequency::type, double>& phaseCoeff,
     39                              std::map<t_frequency::type, double>& ionoCoeff) const;
    3640  double              lambda(t_lc::type tLC) const;
    3741  double              sigma(t_lc::type tLC) const;
    … …  
    4145  void                setRes(t_lc::type tLC, double res);
    4246  double              getRes(t_lc::type tLC) const;
     47  void                setPseudoObsIono(t_frequency::type freq, double stecRefSat);
     48  double              getIonoCodeDelay(t_frequency::type freq) {return _model._ionoCodeDelay[freq];}
    4349
    4450  // RINEX
    … …  
    7985    ~t_model() {}
    8086    void reset() {
    81       _set     = false;
    82       _rho     = 0.0;
    83       _eleSat  = 0.0;
    84       _azSat   = 0.0;
    85       _recClkM = 0.0;
    86       _satClkM = 0.0;
    87       _sagnac  = 0.0;
    88       _antEcc  = 0.0;
    89       _tropo   = 0.0;
    90       _tide    = 0.0;
    91       _windUp  = 0.0;
    92       _rel    = 0.0;
     87      _set       = false;
     88      _rho       = 0.0;
     89      _eleSat    = 0.0;
     90      _azSat     = 0.0;
     91      _recClkM   = 0.0;
     92      _satClkM   = 0.0;
     93      _sagnac    = 0.0;
     94      _antEcc    = 0.0;
     95      _tropo     = 0.0;
     96      _tideEarth = 0.0;
     97      _tideOcean = 0.0;
     98      _windUp    = 0.0;
     99      _rel       = 0.0;
    93100      for (unsigned ii = 0; ii < t_frequency::max; ii++) {
    94101        _antPCO[ii]        = 0.0;
    … …  
    107114    double _antEcc;
    108115    double _tropo;
    109     double _tide;
     116    double _tideEarth;
     117    double _tideOcean;
    110118    double _windUp;
    111119    double _rel;
    … …  
    119127
    120128  bool                         _valid;
     129  bool                         _reference;
    121130  t_frequency::type            _fType1;
    122131  t_frequency::type            _fType2;
    … …  
    131140  std::map<t_lc::type, double> _res;
    132141  double                       _signalPropagationTime;
     142  double                       _stecRefSat;
     143  double                       _stecSat;
    133144};
    134145

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

    r8619 r8905  
    2121  void setNeuEcc(const ColumnVector& neuEcc);
    2222  void setDClk(double dClk) {_dClk = dClk;}
    23   void setTideDspl(const ColumnVector& tideDspl) {_tideDspl = tideDspl;}
    24   void setIonoEpochTime(const bncTime& epochTime) {_epochTime = epochTime;}
     23  void setTideDsplEarth(const ColumnVector& tideDsplEarth) {_tideDsplEarth = tideDsplEarth;}
     24  void setTideDsplOcean(const ColumnVector& tideDsplOcean) {_tideDsplOcean = tideDsplOcean;}
     25  void setEpochTime(const bncTime& epochTime) {_epochTime = epochTime;}
    2526  const std::string&  name()      const {return _name;}
    2627  const std::string&  antName()   const {return _antName;}
    … …  
    2930  const ColumnVector& neuEcc()    const {return _neuEcc;}
    3031  const ColumnVector& xyzEcc()    const {return _xyzEcc;}
    31   const ColumnVector& tideDspl()  const {return _tideDspl;}
     32  const ColumnVector& tideDsplEarth()  const {return _tideDsplEarth;}
     33  const ColumnVector& tideDsplOcean()  const {return _tideDsplOcean;}
     34  const bncTime       epochTime() const {return _epochTime;}
    3235  double dClk() const {return _dClk;}
    3336  double windUp(const bncTime& time, t_prn prn, const ColumnVector& rSat, bool ssr,
    … …  
    4245  ColumnVector      _neuEcc;
    4346  ColumnVector      _xyzEcc;
    44   ColumnVector      _tideDspl;
     47  ColumnVector      _tideDsplEarth;
     48  ColumnVector      _tideDsplOcean;
    4549  double            _dClk;
    4650  mutable t_windUp* _windUp;

  • trunk/BNC/src/bncwindow.cpp

    r8733 r8905  
    10231023  pppLayout3->addWidget(_pppWidgets._minEle,                ir, 7);
    10241024  ++ir;
    1025   pppLayout3->addWidget(new QLabel("Wait for clock corr."), ir, 0, Qt::AlignLeft);
    1026   pppLayout3->addWidget(_pppWidgets._corrWaitTime,          ir, 1);
    1027   pppLayout3->addWidget(new QLabel("Seeding (sec)"),        ir, 3, Qt::AlignLeft);
    1028   pppLayout3->addWidget(_pppWidgets._seedingTime,           ir, 4);
     1025  pppLayout3->addWidget(new QLabel("Model Obs"),            ir, 0, Qt::AlignLeft);
     1026  pppLayout3->addWidget(_pppWidgets._modelObs,              ir, 1);
     1027  pppLayout3->addWidget(new QLabel("Wait for clock corr."), ir, 3, Qt::AlignLeft);
     1028  pppLayout3->addWidget(_pppWidgets._corrWaitTime,          ir, 4);
     1029  pppLayout3->addWidget(new QLabel("Seeding (sec)"),        ir, 6, Qt::AlignLeft);
     1030  pppLayout3->addWidget(_pppWidgets._seedingTime,           ir, 7);
     1031  ++ir;
     1032  pppLayout3->addWidget(new QLabel("Pseudo Obs"),           ir, 0, Qt::AlignLeft);
     1033  pppLayout3->addWidget(_pppWidgets._pseudoObs,             ir, 1);
    10291034  ++ir;
    10301035  pppLayout3->addWidget(new QLabel(""),                     ir, 8);
    10311036  pppLayout3->setColumnStretch(8, 999);
    10321037  ++ir;
    1033   pppLayout3->addWidget(new QLabel(""),                      ir, 1);
     1038  pppLayout3->addWidget(new QLabel(""),                     ir, 1);
    10341039  pppLayout3->setRowStretch(ir, 999);
    10351040
    … …  
    14021407  _pppWidgets._lcGLONASS->setWhatsThis(tr("<p>Specify which kind of GLONASS observations you want to use and on which kind of ionosphere-free linear combination of GLONASS observations you want to base ambiguity resolutions.</p><p><ul><li>Specifying 'P3' means that you request BNC to use code data and so-called P3 ionosphere-free linear combination of code observations.</li><li>'L3' means that you request BNC to use phase data and so-called L3 ionosphere-free linear combination of phase observations.</li> <li>'P3&L3' means that you request BNC to use both, code and phase data and so-called P3 and L3 ionosphere-free linear combination of code and phase observations.</li></ul></p><p>Specifying 'no' means that you don't want BNC to use GLONASS data. <i>[key: PPP/lcGLONASS]</i></p>"));
    14031408  _pppWidgets._lcGalileo->setWhatsThis(tr("<p>Specify which kind of Galileo observations you want to use and on which kind of of ionosphere-free linear combination of Galileo observations you want to base ambiguity resolutions.</p><p><ul><li>Specifying 'P3' means that you request BNC to use code data and so-called P3 ionosphere-free linear combination of code observations.</li><li>'L3' means that you request BNC to use phase data and so-called L3 ionosphere-free linear combination of phase observations.</li> <li>'P3&L3' means that you request BNC to use both, code and phase data and so-called P3 and L3 ionosphere-free linear combination of code and phase observations.</li></ul></p><p>Specifying on of these options makes only sense if Galileo data are part of the processed observation stream.</p><p>Specifying 'no' means that you don't want BNC to use Galileo data. <i>[key: PPP/lcGalileo]</i></p>"));
    1404   _pppWidgets._lcBDS->setWhatsThis(tr("<p>Specify which kind of BDS observations you want to use and on which kind of ionosphere-free linear combination of BDS observations you want to base ambiguity resolutions.</p><p><ul><li>Specifying 'P3' means that you request BNC to use code data and so-called P3 ionosphere-free linear combination of code observations.</li><li>'L3' means that you request BNC to use phase data and so-called L3 ionosphere-free linear combination of phase observations.</li> <li>'P3&L3' means that you request BNC to use both, code and phase data and so-called P3 and L3 ionosphere-free linear combination of code and phase observations.</li></ul></p><p>Specifying on of these options makes only sense if BDS data are part of the processed observation stream.</p><p>Specifying 'no' means that you don't want to use BDS data. <i>[key: PPP/lcBDS]</i></p>"));
     1409  _pppWidgets._lcBDS->setWhatsThis(tr("<p>Specify which kind of model for observations you want to use and on which kind of ionosphere-free linear combination of BDS observations you want to base ambiguity resolutions.</p><p><ul><li>Specifying 'P3' means that you request BNC to use code data and so-called P3 ionosphere-free linear combination of code observations.</li><li>'L3' means that you request BNC to use phase data and so-called L3 ionosphere-free linear combination of phase observations.</li> <li>'P3&L3' means that you request BNC to use both, code and phase data and so-called P3 and L3 ionosphere-free linear combination of code and phase observations.</li></ul></p><p>Specifying on of these options makes only sense if BDS data are part of the processed observation stream.</p><p>Specifying 'no' means that you don't want to use BDS data. <i>[key: PPP/lcBDS]</i></p>"));
     1410  _pppWidgets._modelObs->setWhatsThis(tr("<p>Specify which kind of PPP model you want to use. <i>[key: PPP/modelObs]</i></p>"));
    14051411  _pppWidgets._sigmaC1->setWhatsThis(tr("<p>Enter a Sigma for GNSS C1 code observations in meters.</p><p>The higher the sigma you enter, the less the contribution of C1 code observations to a PPP solution from combined code and phase data. 2.0 is likely to be an appropriate choice.</p><p>Default is an empty option field, meaning<br>'Sigma C1 = 2.0' <i>[key: PPP/sigmaC1]</i></p>"));
    14061412  _pppWidgets._sigmaL1->setWhatsThis(tr("<p>Enter a Sigma for GNSS L1 phase observations in meters.</p><p>The higher the sigma you enter, the less the contribution of L1 phase observations to a PPP solutions from combined code and phase data. 0.01 is likely to be an appropriate choice.</p><p>Default is an empty option field, meaning<br>'Sigma L1 = 0.01' <i>[key: PPP/sigmaL1]</i></p>"));
    … …  
    28552861      QComboBox* system = new QComboBox();
    28562862      system->setEditable(false);
    2857       system->addItems(QString("ALL,GPS,GLONASS,Galileo,BDS,QZSS,SBAS").split(","));
     2863      system->addItems(QString("ALL,GPS,GLONASS,Galileo,BDS,QZSS,SBAS,IRNSS").split(","));
    28582864      system->setFrame(false);
    28592865      _uploadEphTable->setCellWidget(iRow, iCol, system);
    … …  
    29412947        QComboBox* system = new QComboBox();
    29422948        system->setEditable(false);
    2943         system->addItems(QString("ALL,GPS,GLONASS,Galileo,BDS,QZSS,SBAS").split(","));
     2949        system->addItems(QString("ALL,GPS,GLONASS,Galileo,BDS,QZSS,SBAS,IRNSS").split(","));
    29442950        system->setFrame(false);
    29452951        system->setCurrentIndex(system->findText(hlp[iCol]));

  • trunk/BNC/src/pppInclude.h

    r7926 r8905  
    4242class t_lc {
    4343 public:
    44   enum type {dummy = 0, l1, l2, c1, c2, lIF, cIF, MW, CL, maxLc};
     44  enum type {dummy = 0, l1, l2, c1, c2, lIF, cIF, MW, CL, GIM, maxLc};
    4545
    4646  static bool includesPhase(type tt) {
    … …  
    5656    case cIF:
    5757      return false;
    58     case dummy: case maxLc: return false;
     58    case dummy:
     59    case maxLc:
     60    case GIM:
     61      return false;
    5962    }
    6063    return false;
    … …  
    7376    case lIF:
    7477      return false;
    75     case dummy: case maxLc: return false;
     78    case dummy:
     79    case maxLc:
     80    case GIM:
     81      return false;
    7682    }
    7783    return false;
    … …  
    94100    case lIF: case cIF: case MW: case CL:
    95101      return t_frequency::dummy;
    96     case dummy: case maxLc: return t_frequency::dummy;
     102    case dummy:
     103    case maxLc:
     104    case GIM:
     105      return t_frequency::dummy;
    97106    }
    98107    return t_frequency::dummy;
    … …  
    109118    case c2:  return "c2";
    110119    case cIF: return "cIF";
    111     case dummy: case maxLc: return "";
     120    case GIM: return "GIM";
     121    case dummy:
     122    case maxLc:
     123      return "";
    112124    }
    113125    return "";

  • trunk/BNC/src/pppMain.cpp

    r8773 r8905  
    109109        }
    110110        delete pppThread;
    111       }
     111     }
    112112#endif
    113113    }
    114114    _pppThreads.clear();
    115   }
     115 }
    116116
    117117  _running = false;
    … …  
    183183      opt->_corrWaitTime = 0;
    184184    }
     185    opt->_obsModelType   = t_pppOptions::IF;
     186    opt->_pseudoObsIono  = false;
     187    opt->_refSatRequired = false;
     188#ifdef USE_PPP
     189    // Pseudo Observations
     190    if      (settings.value("PPP/pseudoObs").toString() == "Ionosphere") {
     191      opt->_pseudoObsIono = true;
     192    }
     193    else if (settings.value("PPP/pseudoObs").toString() == "no") {
     194      opt->_pseudoObsIono = false;
     195    }
     196    // Observation Model
     197    if      (settings.value("PPP/modelObs").toString() == "Ionosphere-free PPP") {
     198      opt->_obsModelType = t_pppOptions::IF;
     199      opt->_pseudoObsIono = false;
     200    }
     201    else if (settings.value("PPP/modelObs").toString() == "PPP-RTK") {
     202      opt->_obsModelType = t_pppOptions::PPPRTK;
     203      opt->_pseudoObsIono = false;
     204    }
     205    else if (settings.value("PPP/modelObs").toString() == "Uncombined PPP") {
     206      opt->_obsModelType = t_pppOptions::UncombPPP;
     207      if (opt->_pseudoObsIono) {
     208        opt->_refSatRequired = true;
     209      }
     210    }
     211    else if (settings.value("PPP/modelObs").toString() == "DCM with Code Biases") {
     212      opt->_obsModelType = t_pppOptions::DCMcodeBias;
     213      opt->_refSatRequired = true;
     214    }
     215    else if (settings.value("PPP/modelObs").toString() == "DCM with Phase Biases") {
     216      opt->_obsModelType = t_pppOptions::DCMphaseBias;
     217      opt->_refSatRequired = true;
     218    }
     219#endif
    185220    // GPS
    186     if      (settings.value("PPP/lcGPS").toString() == "Pi") {
    187       opt->_LCsGPS.push_back(t_lc::c1);
    188       opt->_LCsGPS.push_back(t_lc::c2);
     221    if (settings.value("PPP/lcGPS").toString() == "Pi") {
     222      if (opt->_obsModelType == t_pppOptions::IF) {
     223        opt->_LCsGPS.push_back(t_lc::cIF);
     224      }
     225      else {
     226        opt->_LCsGPS.push_back(t_lc::c1);
     227        opt->_LCsGPS.push_back(t_lc::c2);
     228        if (opt->_pseudoObsIono) {
     229          opt->_LCsGPS.push_back(t_lc::GIM);
     230        }
     231      }
    189232    }
    190233    else if (settings.value("PPP/lcGPS").toString() == "Li") {
    191       opt->_LCsGPS.push_back(t_lc::l1);
    192       opt->_LCsGPS.push_back(t_lc::l2);
     234      if (opt->_obsModelType == t_pppOptions::IF) {
     235        opt->_LCsGPS.push_back(t_lc::lIF);
     236      }
     237      else {
     238        opt->_LCsGPS.push_back(t_lc::l1);
     239        opt->_LCsGPS.push_back(t_lc::l2);
     240        if (opt->_pseudoObsIono) {
     241          opt->_LCsGPS.push_back(t_lc::GIM);
     242        }
     243      }
    193244    }
    194245    else if (settings.value("PPP/lcGPS").toString() == "Pi&Li") {
    195       opt->_LCsGPS.push_back(t_lc::c1);
    196       opt->_LCsGPS.push_back(t_lc::c2);
    197       opt->_LCsGPS.push_back(t_lc::l1);
    198       opt->_LCsGPS.push_back(t_lc::l2);
    199     }
    200     if      (settings.value("PPP/lcGPS").toString() == "P3") {
    201       opt->_LCsGPS.push_back(t_lc::cIF);
    202     }
    203     else if (settings.value("PPP/lcGPS").toString() == "L3") {
    204       opt->_LCsGPS.push_back(t_lc::lIF);
    205     }
    206     else if (settings.value("PPP/lcGPS").toString() == "P3&L3") {
    207       opt->_LCsGPS.push_back(t_lc::cIF);
    208       opt->_LCsGPS.push_back(t_lc::lIF);
     246      if (opt->_obsModelType == t_pppOptions::IF) {
     247        opt->_LCsGPS.push_back(t_lc::cIF);
     248        opt->_LCsGPS.push_back(t_lc::lIF);
     249      }
     250      else {
     251        opt->_LCsGPS.push_back(t_lc::c1);
     252        opt->_LCsGPS.push_back(t_lc::c2);
     253        opt->_LCsGPS.push_back(t_lc::l1);
     254        opt->_LCsGPS.push_back(t_lc::l2);
     255        if (opt->_pseudoObsIono) {
     256          opt->_LCsGPS.push_back(t_lc::GIM);
     257        }
     258      }
    209259    }
    210260    // GLONASS
    211     if      (settings.value("PPP/lcGLONASS").toString() == "Pi") {
    212       opt->_LCsGLONASS.push_back(t_lc::c1);
    213       opt->_LCsGLONASS.push_back(t_lc::c2);
     261    if (settings.value("PPP/lcGLONASS").toString() == "Pi") {
     262      if (opt->_obsModelType == t_pppOptions::IF) {
     263        opt->_LCsGLONASS.push_back(t_lc::cIF);
     264      }
     265      else {
     266        opt->_LCsGLONASS.push_back(t_lc::c1);
     267        opt->_LCsGLONASS.push_back(t_lc::c2);
     268        if (opt->_pseudoObsIono) {
     269          opt->_LCsGLONASS.push_back(t_lc::GIM);
     270        }
     271      }
    214272    }
    215273    else if (settings.value("PPP/lcGLONASS").toString() == "Li") {
    216       opt->_LCsGLONASS.push_back(t_lc::l1);
    217       opt->_LCsGLONASS.push_back(t_lc::l2);
     274      if (opt->_obsModelType == t_pppOptions::IF) {
     275        opt->_LCsGLONASS.push_back(t_lc::lIF);
     276      }
     277      else {
     278        opt->_LCsGLONASS.push_back(t_lc::l1);
     279        opt->_LCsGLONASS.push_back(t_lc::l2);
     280        if (opt->_obsModelType == t_pppOptions::IF) {
     281          opt->_LCsGLONASS.push_back(t_lc::GIM);
     282        }
     283      }
    218284    }
    219285    else if (settings.value("PPP/lcGLONASS").toString() == "Pi&Li") {
    220       opt->_LCsGLONASS.push_back(t_lc::c1);
    221       opt->_LCsGLONASS.push_back(t_lc::c2);
    222       opt->_LCsGLONASS.push_back(t_lc::l1);
    223       opt->_LCsGLONASS.push_back(t_lc::l2);
    224     }
    225     if      (settings.value("PPP/lcGLONASS").toString() == "P3") {
    226       opt->_LCsGLONASS.push_back(t_lc::cIF);
    227     }
    228     else if (settings.value("PPP/lcGLONASS").toString() == "L3") {
    229       opt->_LCsGLONASS.push_back(t_lc::lIF);
    230     }
    231     else if (settings.value("PPP/lcGLONASS").toString() == "P3&L3") {
    232       opt->_LCsGLONASS.push_back(t_lc::cIF);
    233       opt->_LCsGLONASS.push_back(t_lc::lIF);
     286      if (opt->_obsModelType == t_pppOptions::IF) {
     287        opt->_LCsGLONASS.push_back(t_lc::cIF);
     288        opt->_LCsGLONASS.push_back(t_lc::lIF);
     289      }
     290      else {
     291        opt->_LCsGLONASS.push_back(t_lc::c1);
     292        opt->_LCsGLONASS.push_back(t_lc::c2);
     293        opt->_LCsGLONASS.push_back(t_lc::l1);
     294        opt->_LCsGLONASS.push_back(t_lc::l2);
     295        if (opt->_pseudoObsIono) {
     296          opt->_LCsGLONASS.push_back(t_lc::GIM);
     297        }
     298      }
    234299    }
    235300    // Galileo
    236     if      (settings.value("PPP/lcGalileo").toString() == "Pi") {
    237       opt->_LCsGalileo.push_back(t_lc::c1);
    238       opt->_LCsGalileo.push_back(t_lc::c2);
     301    if (settings.value("PPP/lcGalileo").toString() == "Pi") {
     302      if (opt->_obsModelType == t_pppOptions::IF) {
     303        opt->_LCsGalileo.push_back(t_lc::cIF);
     304      }
     305      else {
     306        opt->_LCsGalileo.push_back(t_lc::c1);
     307        opt->_LCsGalileo.push_back(t_lc::c2);
     308        if (opt->_pseudoObsIono) {
     309          opt->_LCsGalileo.push_back(t_lc::GIM);
     310        }
     311      }
    239312    }
    240313    else if (settings.value("PPP/lcGalileo").toString() == "Li") {
    241       opt->_LCsGalileo.push_back(t_lc::l1);
    242       opt->_LCsGalileo.push_back(t_lc::l2);
     314      if (opt->_obsModelType == t_pppOptions::IF) {
     315        opt->_LCsGalileo.push_back(t_lc::lIF);
     316      }
     317      else {
     318        opt->_LCsGalileo.push_back(t_lc::l1);
     319        opt->_LCsGalileo.push_back(t_lc::l2);
     320        if (opt->_pseudoObsIono) {
     321          opt->_LCsGalileo.push_back(t_lc::GIM);
     322        }
     323      }
    243324    }
    244325    else if (settings.value("PPP/lcGalileo").toString() == "Pi&Li") {
    245       opt->_LCsGalileo.push_back(t_lc::c1);
    246       opt->_LCsGalileo.push_back(t_lc::c2);
    247       opt->_LCsGalileo.push_back(t_lc::l1);
    248       opt->_LCsGalileo.push_back(t_lc::l2);
    249     }
    250     if      (settings.value("PPP/lcGalileo").toString() == "P3") {
    251       opt->_LCsGalileo.push_back(t_lc::cIF);
    252     }
    253     else if (settings.value("PPP/lcGalileo").toString() == "L3") {
    254       opt->_LCsGalileo.push_back(t_lc::lIF);
    255     }
    256     else if (settings.value("PPP/lcGalileo").toString() == "P3&L3") {
    257       opt->_LCsGalileo.push_back(t_lc::cIF);
    258       opt->_LCsGalileo.push_back(t_lc::lIF);
     326      if (opt->_obsModelType == t_pppOptions::IF) {
     327        opt->_LCsGalileo.push_back(t_lc::cIF);
     328        opt->_LCsGalileo.push_back(t_lc::lIF);
     329      }
     330      else {
     331        opt->_LCsGalileo.push_back(t_lc::c1);
     332        opt->_LCsGalileo.push_back(t_lc::c2);
     333        opt->_LCsGalileo.push_back(t_lc::l1);
     334        opt->_LCsGalileo.push_back(t_lc::l2);
     335        if (opt->_pseudoObsIono) {
     336          opt->_LCsGalileo.push_back(t_lc::GIM);
     337        }
     338      }
    259339    }
    260340    // BDS
    261     if      (settings.value("PPP/lcBDS").toString() == "Pi") {
    262       opt->_LCsBDS.push_back(t_lc::c1);
    263       opt->_LCsBDS.push_back(t_lc::c2);
     341    if (settings.value("PPP/lcBDS").toString() == "Pi") {
     342      if (opt->_obsModelType == t_pppOptions::IF) {
     343        opt->_LCsBDS.push_back(t_lc::cIF);
     344      }
     345      else {
     346        opt->_LCsBDS.push_back(t_lc::c1);
     347        opt->_LCsBDS.push_back(t_lc::c2);
     348        if (opt->_pseudoObsIono) {
     349          opt->_LCsBDS.push_back(t_lc::GIM);
     350        }
     351      }
    264352    }
    265353    else if (settings.value("PPP/lcBDS").toString() == "Li") {
    266       opt->_LCsBDS.push_back(t_lc::l1);
    267       opt->_LCsBDS.push_back(t_lc::l2);
     354      if (opt->_obsModelType == t_pppOptions::IF) {
     355        opt->_LCsBDS.push_back(t_lc::lIF);
     356      }
     357      else {
     358        opt->_LCsBDS.push_back(t_lc::l1);
     359        opt->_LCsBDS.push_back(t_lc::l2);
     360        if (opt->_pseudoObsIono) {
     361          opt->_LCsBDS.push_back(t_lc::GIM);
     362        }
     363      }
    268364    }
    269365    else if (settings.value("PPP/lcBDS").toString() == "Pi&Li") {
    270       opt->_LCsBDS.push_back(t_lc::c1);
    271       opt->_LCsBDS.push_back(t_lc::c2);
    272       opt->_LCsBDS.push_back(t_lc::l1);
    273       opt->_LCsBDS.push_back(t_lc::l2);
    274     }
    275     if      (settings.value("PPP/lcBDS").toString() == "P3") {
    276       opt->_LCsBDS.push_back(t_lc::cIF);
    277     }
    278     else if (settings.value("PPP/lcBDS").toString() == "L3") {
    279       opt->_LCsBDS.push_back(t_lc::lIF);
    280     }
    281     else if (settings.value("PPP/lcBDS").toString() == "P3&L3") {
    282       opt->_LCsBDS.push_back(t_lc::cIF);
    283       opt->_LCsBDS.push_back(t_lc::lIF);
     366      if (opt->_obsModelType == t_pppOptions::IF) {
     367        opt->_LCsBDS.push_back(t_lc::cIF);
     368        opt->_LCsBDS.push_back(t_lc::lIF);
     369      }
     370      else {
     371        opt->_LCsBDS.push_back(t_lc::c1);
     372        opt->_LCsBDS.push_back(t_lc::c2);
     373        opt->_LCsBDS.push_back(t_lc::l1);
     374        opt->_LCsBDS.push_back(t_lc::l2);
     375        if (opt->_pseudoObsIono) {
     376          opt->_LCsBDS.push_back(t_lc::GIM);
     377        }
     378      }
    284379    }
    285380
    … …  
    316411    // Some default values
    317412    // -------------------
    318     opt->_aprSigAmb   = 1000.0;
    319     opt->_noiseClk    = 1000.0;
     413    opt->_aprSigAmb       = 1000.0;
     414    opt->_aprSigIon       = 1000.0;
     415    opt->_noiseClk        = 1000.0;
     416    opt->_aprSigCodeBias  = 1000.0;
     417    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)
    320422
    321423    _options << opt;

  • trunk/BNC/src/pppModel.cpp

    r8774 r8905  
    1 
    21// Part of BNC, a utility for retrieving decoding and
    32// converting GNSS data streams from NTRIP broadcasters.
    … …  
    4039 * -----------------------------------------------------------------------*/
    4140
    42 
    4341#include <cmath>
    4442
    … …  
    4846using namespace std;
    4947
    50 const double t_astro::RHO_DEG   = 180.0 / M_PI;
    51 const double t_astro::RHO_SEC   = 3600.0 * 180.0 / M_PI;
    52 const double t_astro::MJD_J2000 = 51544.5;
     48
    5349
    5450Matrix t_astro::rotX(double Angle) {
    5551  const double C = cos(Angle);
    5652  const double S = sin(Angle);
    57   Matrix UU(3,3);
    58   UU[0][0] = 1.0;  UU[0][1] = 0.0;  UU[0][2] = 0.0;
    59   UU[1][0] = 0.0;  UU[1][1] =  +C;  UU[1][2] =  +S;
    60   UU[2][0] = 0.0;  UU[2][1] =  -S;  UU[2][2] =  +C;
     53  Matrix UU(3, 3);
     54  UU[0][0] = 1.0;
     55  UU[0][1] = 0.0;
     56  UU[0][2] = 0.0;
     57  UU[1][0] = 0.0;
     58  UU[1][1] = +C;
     59  UU[1][2] = +S;
     60  UU[2][0] = 0.0;
     61  UU[2][1] = -S;
     62  UU[2][2] = +C;
    6163  return UU;
    6264}
    … …  
    6567  const double C = cos(Angle);
    6668  const double S = sin(Angle);
    67   Matrix UU(3,3);
    68   UU[0][0] =  +C;  UU[0][1] = 0.0;  UU[0][2] =  -S;
    69   UU[1][0] = 0.0;  UU[1][1] = 1.0;  UU[1][2] = 0.0;
    70   UU[2][0] =  +S;  UU[2][1] = 0.0;  UU[2][2] =  +C;
     69  Matrix UU(3, 3);
     70  UU[0][0] = +C;
     71  UU[0][1] = 0.0;
     72  UU[0][2] = -S;
     73  UU[1][0] = 0.0;
     74  UU[1][1] = 1.0;
     75  UU[1][2] = 0.0;
     76  UU[2][0] = +S;
     77  UU[2][1] = 0.0;
     78  UU[2][2] = +C;
    7179  return UU;
    7280}
    … …  
    7583  const double C = cos(Angle);
    7684  const double S = sin(Angle);
    77   Matrix UU(3,3);
    78   UU[0][0] =  +C;  UU[0][1] =  +S;  UU[0][2] = 0.0;
    79   UU[1][0] =  -S;  UU[1][1] =  +C;  UU[1][2] = 0.0;
    80   UU[2][0] = 0.0;  UU[2][1] = 0.0;  UU[2][2] = 1.0;
     85  Matrix UU(3, 3);
     86  UU[0][0] = +C;
     87  UU[0][1] = +S;
     88  UU[0][2] = 0.0;
     89  UU[1][0] = -S;
     90  UU[1][1] = +C;
     91  UU[1][2] = 0.0;
     92  UU[2][0] = 0.0;
     93  UU[2][1] = 0.0;
     94  UU[2][2] = 1.0;
    8195  return UU;
    8296}
    … …  
    89103
    90104  double Mjd_0 = floor(Mjd_UT1);
    91   double UT1   = Secs*(Mjd_UT1-Mjd_0);
    92   double T_0   = (Mjd_0  -MJD_J2000)/36525.0;
    93   double T     = (Mjd_UT1-MJD_J2000)/36525.0;
    94 
    95   double gmst  = 24110.54841 + 8640184.812866*T_0 + 1.002737909350795*UT1
    96                  + (0.093104-6.2e-6*T)*T*T;
    97 
    98   return  2.0*M_PI*Frac(gmst/Secs);
     105  double UT1 = Secs * (Mjd_UT1 - Mjd_0);
     106  double T_0 = (Mjd_0 - MJD_J2000) / 36525.0;
     107  double T = (Mjd_UT1 - MJD_J2000) / 36525.0;
     108
     109  double gmst = 24110.54841 + 8640184.812866 * T_0 + 1.002737909350795 * UT1
     110      + (0.093104 - 6.2e-6 * T) * T * T;
     111
     112  return 2.0 * M_PI * Frac(gmst / Secs);
    99113}
    100114
    … …  
    103117Matrix t_astro::NutMatrix(double Mjd_TT) {
    104118
    105   const double T  = (Mjd_TT-MJD_J2000)/36525.0;
    106 
    107   double ls = 2.0*M_PI*Frac(0.993133+  99.997306*T);
    108   double D  = 2.0*M_PI*Frac(0.827362+1236.853087*T);
    109   double F  = 2.0*M_PI*Frac(0.259089+1342.227826*T);
    110   double N  = 2.0*M_PI*Frac(0.347346-   5.372447*T);
    111 
    112   double dpsi = ( -17.200*sin(N)   - 1.319*sin(2*(F-D+N)) - 0.227*sin(2*(F+N))
    113                 + 0.206*sin(2*N) + 0.143*sin(ls) ) / RHO_SEC;
    114   double deps = ( + 9.203*cos(N)   + 0.574*cos(2*(F-D+N)) + 0.098*cos(2*(F+N))
    115                 - 0.090*cos(2*N)                 ) / RHO_SEC;
    116 
    117   double eps  = 0.4090928-2.2696E-4*T;
    118 
    119   return  rotX(-eps-deps)*rotZ(-dpsi)*rotX(+eps);
     119  const double T = (Mjd_TT - MJD_J2000) / 36525.0;
     120
     121  double ls = 2.0 * M_PI * Frac(0.993133 + 99.997306 * T);
     122  double D = 2.0 * M_PI * Frac(0.827362 + 1236.853087 * T);
     123  double F = 2.0 * M_PI * Frac(0.259089 + 1342.227826 * T);
     124  double N = 2.0 * M_PI * Frac(0.347346 - 5.372447 * T);
     125
     126  double dpsi = (-17.200 * sin(N) - 1.319 * sin(2 * (F - D + N))
     127      - 0.227 * sin(2 * (F + N))
     128      + 0.206 * sin(2 * N) + 0.143 * sin(ls)) / RHO_SEC;
     129  double deps = (+9.203 * cos(N) + 0.574 * cos(2 * (F - D + N))
     130      + 0.098 * cos(2 * (F + N))
     131      - 0.090 * cos(2 * N)) / RHO_SEC;
     132
     133  double eps = 0.4090928 - 2.2696E-4 * T;
     134
     135  return rotX(-eps - deps) * rotZ(-dpsi) * rotX(+eps);
    120136}
    121137
    … …  
    124140Matrix t_astro::PrecMatrix(double Mjd_1, double Mjd_2) {
    125141
    126   const double T  = (Mjd_1-MJD_J2000)/36525.0;
    127   const double dT = (Mjd_2-Mjd_1)/36525.0;
    128 
    129   double zeta  =  ( (2306.2181+(1.39656-0.000139*T)*T)+
    130                         ((0.30188-0.000344*T)+0.017998*dT)*dT )*dT/RHO_SEC;
    131   double z     =  zeta + ( (0.79280+0.000411*T)+0.000205*dT)*dT*dT/RHO_SEC;
    132   double theta =  ( (2004.3109-(0.85330+0.000217*T)*T)-
    133                         ((0.42665+0.000217*T)+0.041833*dT)*dT )*dT/RHO_SEC;
     142  const double T = (Mjd_1 - MJD_J2000) / 36525.0;
     143  const double dT = (Mjd_2 - Mjd_1) / 36525.0;
     144
     145  double zeta = ((2306.2181 + (1.39656 - 0.000139 * T) * T) +
     146      ((0.30188 - 0.000344 * T) + 0.017998 * dT) * dT) * dT / RHO_SEC;
     147  double z = zeta
     148      + ((0.79280 + 0.000411 * T) + 0.000205 * dT) * dT * dT / RHO_SEC;
     149  double theta = ((2004.3109 - (0.85330 + 0.000217 * T) * T) -
     150      ((0.42665 + 0.000217 * T) + 0.041833 * dT) * dT) * dT / RHO_SEC;
    134151
    135152  return rotZ(-z) * rotY(theta) * rotZ(-zeta);
    … …  
    140157ColumnVector t_astro::Sun(double Mjd_TT) {
    141158
    142   const double eps = 23.43929111/RHO_DEG;
    143   const double T   = (Mjd_TT-MJD_J2000)/36525.0;
    144 
    145   double M = 2.0*M_PI * Frac ( 0.9931267 + 99.9973583*T);
    146   double L = 2.0*M_PI * Frac ( 0.7859444 + M/2.0/M_PI +
    147                         (6892.0*sin(M)+72.0*sin(2.0*M)) / 1296.0e3);
    148   double r = 149.619e9 - 2.499e9*cos(M) - 0.021e9*cos(2*M);
     159  const double eps = 23.43929111 / RHO_DEG;
     160  const double T = (Mjd_TT - MJD_J2000) / 36525.0;
     161
     162  double M = 2.0 * M_PI * Frac(0.9931267 + 99.9973583 * T);
     163  double L = 2.0 * M_PI * Frac(0.7859444 + M / 2.0 / M_PI +
     164      (6892.0 * sin(M) + 72.0 * sin(2.0 * M)) / 1296.0e3);
     165  double r = 149.619e9 - 2.499e9 * cos(M) - 0.021e9 * cos(2 * M);
    149166
    150167  ColumnVector r_Sun(3);
    151   r_Sun << r*cos(L) << r*sin(L) << 0.0; r_Sun = rotX(-eps) * r_Sun;
    152 
    153   return    rotZ(GMST(Mjd_TT))
    154           * NutMatrix(Mjd_TT)
    155           * PrecMatrix(MJD_J2000, Mjd_TT)
    156           * r_Sun;
     168  r_Sun << r * cos(L) << r * sin(L) << 0.0;
     169  r_Sun = rotX(-eps) * r_Sun;
     170
     171  return rotZ(GMST(Mjd_TT))
     172      * NutMatrix(Mjd_TT)
     173      * PrecMatrix(MJD_J2000, Mjd_TT)
     174      * r_Sun;
    157175}
    158176
    … …  
    161179ColumnVector t_astro::Moon(double Mjd_TT) {
    162180
    163   const double eps = 23.43929111/RHO_DEG;
    164   const double T   = (Mjd_TT-MJD_J2000)/36525.0;
    165 
    166   double L_0 = Frac ( 0.606433 + 1336.851344*T );
    167   double l   = 2.0*M_PI*Frac ( 0.374897 + 1325.552410*T );
    168   double lp  = 2.0*M_PI*Frac ( 0.993133 +   99.997361*T );
    169   double D   = 2.0*M_PI*Frac ( 0.827361 + 1236.853086*T );
    170   double F   = 2.0*M_PI*Frac ( 0.259086 + 1342.227825*T );
    171 
    172   double dL = +22640*sin(l) - 4586*sin(l-2*D) + 2370*sin(2*D) +  769*sin(2*l)
    173               -668*sin(lp) - 412*sin(2*F) - 212*sin(2*l-2*D)- 206*sin(l+lp-2*D)
    174               +192*sin(l+2*D) - 165*sin(lp-2*D) - 125*sin(D) - 110*sin(l+lp)
    175               +148*sin(l-lp) - 55*sin(2*F-2*D);
    176 
    177   double L = 2.0*M_PI * Frac( L_0 + dL/1296.0e3 );
    178 
    179   double S  = F + (dL+412*sin(2*F)+541*sin(lp)) / RHO_SEC;
    180   double h  = F-2*D;
    181   double N  = -526*sin(h) + 44*sin(l+h) - 31*sin(-l+h) - 23*sin(lp+h)
    182               +11*sin(-lp+h) - 25*sin(-2*l+F) + 21*sin(-l+F);
    183 
    184   double B = ( 18520.0*sin(S) + N ) / RHO_SEC;
     181  const double eps = 23.43929111 / RHO_DEG;
     182  const double T = (Mjd_TT - MJD_J2000) / 36525.0;
     183
     184  double L_0 = Frac(0.606433 + 1336.851344 * T);
     185  double l = 2.0 * M_PI * Frac(0.374897 + 1325.552410 * T);
     186  double lp = 2.0 * M_PI * Frac(0.993133 + 99.997361 * T);
     187  double D = 2.0 * M_PI * Frac(0.827361 + 1236.853086 * T);
     188  double F = 2.0 * M_PI * Frac(0.259086 + 1342.227825 * T);
     189
     190  double dL = +22640 * sin(l) - 4586 * sin(l - 2 * D) + 2370 * sin(2 * D)
     191      + 769 * sin(2 * l)
     192      - 668 * sin(lp) - 412 * sin(2 * F) - 212 * sin(2 * l - 2 * D)
     193      - 206 * sin(l + lp - 2 * D)
     194      + 192 * sin(l + 2 * D) - 165 * sin(lp - 2 * D) - 125 * sin(D)
     195      - 110 * sin(l + lp)
     196      + 148 * sin(l - lp) - 55 * sin(2 * F - 2 * D);
     197
     198  double L = 2.0 * M_PI * Frac(L_0 + dL / 1296.0e3);
     199
     200  double S = F + (dL + 412 * sin(2 * F) + 541 * sin(lp)) / RHO_SEC;
     201  double h = F - 2 * D;
     202  double N = -526 * sin(h) + 44 * sin(l + h) - 31 * sin(-l + h)
     203      - 23 * sin(lp + h)
     204      + 11 * sin(-lp + h) - 25 * sin(-2 * l + F) + 21 * sin(-l + F);
     205
     206  double B = (18520.0 * sin(S) + N) / RHO_SEC;
    185207
    186208  double cosB = cos(B);
    187209
    188   double R = 385000e3 - 20905e3*cos(l) - 3699e3*cos(2*D-l) - 2956e3*cos(2*D)
    189       -570e3*cos(2*l) + 246e3*cos(2*l-2*D) - 205e3*cos(lp-2*D)
    190       -171e3*cos(l+2*D) - 152e3*cos(l+lp-2*D);
     210  double R = 385000e3 - 20905e3 * cos(l) - 3699e3 * cos(2 * D - l)
     211      - 2956e3 * cos(2 * D)
     212      - 570e3 * cos(2 * l) + 246e3 * cos(2 * l - 2 * D)
     213      - 205e3 * cos(lp - 2 * D)
     214      - 171e3 * cos(l + 2 * D) - 152e3 * cos(l + lp - 2 * D);
    191215
    192216  ColumnVector r_Moon(3);
    193   r_Moon << R*cos(L)*cosB << R*sin(L)*cosB << R*sin(B);
     217  r_Moon << R * cos(L) * cosB << R * sin(L) * cosB << R * sin(B);
    194218  r_Moon = rotX(-eps) * r_Moon;
    195219
    196   return    rotZ(GMST(Mjd_TT))
    197           * NutMatrix(Mjd_TT)
    198           * PrecMatrix(MJD_J2000, Mjd_TT)
    199           * r_Moon;
     220  return rotZ(GMST(Mjd_TT))
     221      * NutMatrix(Mjd_TT)
     222      * PrecMatrix(MJD_J2000, Mjd_TT)
     223      * r_Moon;
    200224}
    201225
    202226// Tidal Correction
    203227////////////////////////////////////////////////////////////////////////////
    204 ColumnVector t_tides::displacement(const bncTime& time, const ColumnVector& xyz) {
     228ColumnVector t_tides::earth(const bncTime& time, const ColumnVector& xyz) {
    205229
    206230  if (time.undef()) {
    207     ColumnVector dX(3); dX = 0.0;
     231    ColumnVector dX(3);
     232    dX = 0.0;
    208233    return dX;
    209234  }
    … …  
    214239    _lastMjd = Mjd;
    215240    _xSun = t_astro::Sun(Mjd);
    216     _rSun = sqrt(DotProduct(_xSun,_xSun));
     241    _rSun = sqrt(DotProduct(_xSun, _xSun));
    217242    _xSun /= _rSun;
    218243    _xMoon = t_astro::Moon(Mjd);
    219     _rMoon = sqrt(DotProduct(_xMoon,_xMoon));
     244    _rMoon = sqrt(DotProduct(_xMoon, _xMoon));
    220245    _xMoon /= _rMoon;
    221246  }
    222247
    223   double       rRec    = sqrt(DotProduct(xyz, xyz));
     248  double rRec = sqrt(DotProduct(xyz, xyz));
    224249  ColumnVector xyzUnit = xyz / rRec;
    225250
    … …  
    231256  // Tidal Displacement
    232257  // ------------------
    233   double scSun  = DotProduct(xyzUnit, _xSun);
     258  double scSun = DotProduct(xyzUnit, _xSun);
    234259  double scMoon = DotProduct(xyzUnit, _xMoon);
    235260
    236   double p2Sun  = 3.0 * (H2/2.0-L2) * scSun  * scSun  - H2/2.0;
    237   double p2Moon = 3.0 * (H2/2.0-L2) * scMoon * scMoon - H2/2.0;
    238 
    239   double x2Sun  = 3.0 * L2 * scSun;
     261  double p2Sun = 3.0 * (H2 / 2.0 - L2) * scSun * scSun - H2 / 2.0;
     262  double p2Moon = 3.0 * (H2 / 2.0 - L2) * scMoon * scMoon - H2 / 2.0;
     263
     264  double x2Sun = 3.0 * L2 * scSun;
    240265  double x2Moon = 3.0 * L2 * scMoon;
    241266
    242267  const double gmWGS = 398.6005e12;
    243   const double gms   = 1.3271250e20;
    244   const double gmm   = 4.9027890e12;
    245 
    246   double facSun  = gms / gmWGS *
    247                    (rRec * rRec * rRec * rRec) / (_rSun * _rSun * _rSun);
     268  const double gms = 1.3271250e20;
     269  const double gmm = 4.9027890e12;
     270
     271  double facSun = gms / gmWGS *
     272      (rRec * rRec * rRec * rRec) / (_rSun * _rSun * _rSun);
    248273
    249274  double facMoon = gmm / gmWGS *
    250                    (rRec * rRec * rRec * rRec) / (_rMoon * _rMoon * _rMoon);
    251 
    252   ColumnVector dX = facSun  * (x2Sun  * _xSun  + p2Sun  * xyzUnit) +
    253                     facMoon * (x2Moon * _xMoon + p2Moon * xyzUnit);
     275      (rRec * rRec * rRec * rRec) / (_rMoon * _rMoon * _rMoon);
     276
     277  ColumnVector dX = facSun * (x2Sun * _xSun + p2Sun * xyzUnit)
     278                  + facMoon * (x2Moon * _xMoon + p2Moon * xyzUnit);
    254279
    255280  return dX;
    … …  
    258283// Constructor
    259284///////////////////////////////////////////////////////////////////////////
    260 t_loading::t_loading(const QString& fileName) {
    261 
     285t_tides::t_tides() {
     286  _lastMjd = 0.0;
     287  _rSun = 0.0;
     288  _rMoon = 0.0;
    262289  newBlqData = 0;
    263   readFile(fileName);
    264   printAll();
    265 }
    266 
    267 t_loading::~t_loading() {
     290}
     291
     292t_tides::~t_tides() {
    268293
    269294  if (newBlqData) {
    … …  
    278303}
    279304
    280 t_irc t_loading::readFile(const QString& fileName) {
     305t_irc t_tides::readBlqFile(const char* fileName) {
    281306  QFile inFile(fileName);
    282307  inFile.open(QIODevice::ReadOnly | QIODevice::Text);
    … …  
    285310  QString site = QString();
    286311
    287   while ( !in.atEnd() ) {
     312  while (!in.atEnd()) {
    288313
    289314    QString line = in.readLine();
    290315
    291316    // skip empty lines and comments
    292     if (line.indexOf("$$") != -1 || line.size() == 0) {
     317    if (line.indexOf("$$") != -1) {
    293318      continue;
    294319    }
    295 
    296320    line = line.trimmed();
    297321    QTextStream inLine(line.toLatin1(), QIODevice::ReadOnly);
    298 
    299322    switch (row) {
    300323      case 0:
    301324        site = line;
    302325        site = site.toUpper();
    303         if (!newBlqData) {
    304           newBlqData = new t_blqData;
    305           newBlqData->amplitudes.ReSize(3,11);
    306           newBlqData->phases.ReSize(3,11);
     326        newBlqData = new t_blqData;
     327        newBlqData->amplitudes.ReSize(3, 11);
     328        newBlqData->phases.ReSize(3, 11);
     329        break;
     330      case 1:
     331      case 2:
     332      case 3:
     333        for (int ii = 0; ii < 11; ii++) {
     334          inLine >> newBlqData->amplitudes[row - 1][ii];
    307335        }
    308336        break;
    309       case 1: case 2: case 3:
     337      case 4:
     338      case 5:
    310339        for (int ii = 0; ii < 11; ii++) {
    311           inLine >> newBlqData->amplitudes[row-1][ii];
     340          inLine >> newBlqData->phases[row - 4][ii];
    312341        }
    313342        break;
    314       case 4: case 5: case 6:
     343      case 6:
    315344        for (int ii = 0; ii < 11; ii++) {
    316           inLine >> newBlqData->phases[row-4][ii];
     345          inLine >> newBlqData->phases[row - 4][ii];
    317346        }
    318         break;
    319       case 7:
    320347        if (newBlqData && !site.isEmpty()) {
    321348          blqMap[site] = newBlqData;
    322349          site = QString();
    323           row = -1;
    324350          newBlqData = 0;
    325351        }
     352        row = -1;
    326353        break;
    327354    }
    … …  
    332359}
    333360
     361ColumnVector t_tides::ocean(const bncTime& time,  const ColumnVector& xyz,
     362    const std::string& station) {
     363  ColumnVector dX(3); dX = 0.0;
     364  if (time.undef()) {
     365    return dX;
     366  }
     367  QString stationQ = station.c_str();
     368  if (blqMap.find(stationQ) == blqMap.end()) {
     369    return dX;
     370  }
     371  t_blqData* blqSet = blqMap[stationQ];  //printBlqSet(station, blqSet);
     372
     373  // angular argument: see arg2.f from IERS Conventions software collection
     374  double speed[11] = {1.40519e-4, 1.45444e-4, 1.3788e-4, 1.45842e-4, 7.2921e-5,
     375                      6.7598e-5,  7.2523e-5,  6.4959e-5, 5.3234e-6,  2.6392e-6, 3.982e-7};
     376
     377  double angfac[4][11];
     378  angfac[0][0] = 2.0;
     379  angfac[1][0] =-2.0;
     380  angfac[2][0] = 0.0;
     381  angfac[3][0] = 0.0;
     382
     383  angfac[0][1] = 0.0;
     384  angfac[1][1] = 0.0;
     385  angfac[2][1] = 0.0;
     386  angfac[3][1] = 0.0;
     387
     388  angfac[0][2] = 2.0;
     389  angfac[1][2] =-3.0;
     390  angfac[2][2] = 1.0;
     391  angfac[3][2] = 0.0;
     392
     393  angfac[0][3] = 2.0;
     394  angfac[1][3] = 0.0;
     395  angfac[2][3] = 0.0;
     396  angfac[3][3] = 0.0;
     397
     398  angfac[0][4] = 1.0;
     399  angfac[1][4] = 0.0;
     400  angfac[2][4] = 0.0;
     401  angfac[3][4] = .25;
     402
     403  angfac[0][5] = 1.0;
     404  angfac[1][5] =-2.0;
     405  angfac[2][5] = 0.0;
     406  angfac[3][5] =-.25;
     407
     408  angfac[0][6] =-1.0;
     409  angfac[1][6] = 0.0;
     410  angfac[2][6] = 0.0;
     411  angfac[3][6] =-.25;
     412
     413  angfac[0][7] = 1.0;
     414  angfac[1][7] =-3.0;
     415  angfac[2][7] = 1.0;
     416  angfac[3][7] =-.25;
     417
     418  angfac[0][8] = 0.0;
     419  angfac[1][8] = 2.0;
     420  angfac[2][8] = 0.0;
     421  angfac[3][8] = 0.0;
     422
     423  angfac[0][9] = 0.0;
     424  angfac[1][9] = 1.0;
     425  angfac[2][9] =-1.0;
     426  angfac[3][9] = 0.0;
     427
     428  angfac[0][10] = 2.0;
     429  angfac[1][10] = 0.0;
     430  angfac[2][10] = 0.0;
     431  angfac[3][10] = 0.0;
     432
     433  double twopi = 6.283185307179586476925287e0;
     434  double dtr = 0.0174532925199;
     435
     436  //  fractional part of the day in seconds
     437  unsigned int year, month, day;
     438  time.civil_date(year, month, day);
     439  int iyear = year - 2000;
     440  QDateTime datTim = QDateTime::fromString(QString::fromStdString(time.datestr()), Qt::ISODate);
     441  int doy = datTim.date().dayOfYear();
     442  double fday = time.daysec();
     443  int   icapd = doy + 365 * (iyear - 75) + ((iyear - 73) / 4);
     444  double capt = (icapd * 1.000000035 + 27392.500528) / 36525.0;
     445
     446  // mean longitude of the sun at the beginning of the day
     447  double h0 = (279.69668e0 + (36000.768930485e0 + 3.03e-4 * capt) * capt) * dtr;
     448
     449  // mean longitude of moon at the beginning of the day
     450  double s0 = (((1.9e-6 * capt - .001133e0) * capt + 481267.88314137e0) * capt + 270.434358e0) * dtr;
     451
     452  // mean longitude of lunar perigee at the beginning of the day
     453  double p0 =  (((-1.2e-5 * capt - .010325e0) * capt + 4069.0340329577e0) * capt + 334.329653e0) * dtr;
     454
     455  // tidal angle arguments
     456  double angle[11];
     457  for (int k = 0; k < 11; ++k) {
     458    angle[k] = speed[k] * fday
     459             + angfac[0][k] * h0
     460             + angfac[1][k] * s0
     461             + angfac[2][k] * p0
     462             + angfac[3][k] * twopi;
     463    angle[k] = fmod(angle[k], twopi);
     464    if (angle[k] < 0.0) {
     465      angle[k] += twopi;
     466    }
     467  }
     468
     469  // displacement by 11 constituents
     470  ColumnVector rwsSta(3); rwsSta = 0.0; // radial, west, south
     471  for (int rr = 0; rr < 3; rr++) {
     472    for (int cc = 0; cc < 11; cc++) {
     473      rwsSta[rr] += blqSet->amplitudes[rr][cc] * cos((angle[cc] - (blqSet->phases[rr][cc]/RHO_DEG)));
     474    }
     475  }
     476
     477  // neu2xyz
     478  ColumnVector dneu(3); // neu
     479  dneu[0] = -rwsSta[2];
     480  dneu[1] = -rwsSta[1];
     481  dneu[2] =  rwsSta[0];
     482  double recEll[3]; xyz2ell(xyz.data(), recEll) ;
     483  neu2xyz(recEll, dneu.data(), dX.data());
     484
     485  return dX;
     486}
     487
    334488// Print
    335489////////////////////////////////////////////////////////////////////////////
    336 void t_loading::printAll() const {
     490void t_tides::printAllBlqSets() const {
    337491
    338492  QMapIterator<QString, t_blqData*> it(blqMap);
    … …  
    341495    t_blqData* blq = it.value();
    342496    QString site = it.key();
    343     cout << site.toStdString().c_str() << "\n";
     497    cout << site.toStdString().c_str() << "\n===============\n";
    344498    for (int rr = 0; rr < 3; rr++) {
    345499      for (int cc = 0; cc < 11; cc++) {
    … …  
    357511}
    358512
     513// Print
     514////////////////////////////////////////////////////////////////////////////
     515void t_tides::printBlqSet(const std::string& station, t_blqData* blq) {
     516  cout << station << endl;
     517  for (int rr = 0; rr < 3; rr++) {
     518    for (int cc = 0; cc < 11; cc++) {
     519      cout << blq->amplitudes[rr][cc] << " ";
     520    }
     521    cout << endl;
     522  }
     523  for (int rr = 0; rr < 3; rr++) {
     524    for (int cc = 0; cc < 11; cc++) {
     525      cout << blq->phases[rr][cc] << " ";
     526    }
     527    cout << endl;
     528  }
     529}
     530
    359531// Constructor
    360532///////////////////////////////////////////////////////////////////////////
    361533t_windUp::t_windUp() {
    362534  for (unsigned ii = 0; ii <= t_prn::MAXPRN; ii++) {
    363     sumWind[ii]   = 0.0;
     535    sumWind[ii] = 0.0;
    364536    lastEtime[ii] = 0.0;
    365537  }
    … …  
    369541///////////////////////////////////////////////////////////////////////////
    370542double t_windUp::value(const bncTime& etime, const ColumnVector& rRec,
    371                        t_prn prn, const ColumnVector& rSat, bool ssr,
    372                        double yaw, const ColumnVector& vSat) {
     543    t_prn prn, const ColumnVector& rSat, bool ssr,
     544    double yaw, const ColumnVector& vSat) {
    373545
    374546  if (etime.mjddec() != lastEtime[prn.toInt()]) {
    … …  
    377549    // --------------------------------------
    378550    ColumnVector rho = rRec - rSat;
    379     rho /= rho.norm_Frobenius();
     551    rho /= rho.NormFrobenius();
    380552
    381553    // GPS Satellite unit Vectors sz, sy, sx
    … …  
    392564      sHlp = vSat + crossproduct(Omega, rSat);
    393565    }
    394     sHlp /= sHlp.norm_Frobenius();
    395 
    396     ColumnVector sz = -rSat / rSat.norm_Frobenius();
     566    sHlp /= sHlp.NormFrobenius();
     567
     568    ColumnVector sz = -rSat / rSat.NormFrobenius();
    397569    ColumnVector sy = crossproduct(sz, sHlp);
    398570    ColumnVector sx = crossproduct(sy, sz);
    399571
    400     // Yaw consideration
    401     sx = t_astro::rotZ(yaw) * sx;
    402 
     572    if (ssr) {
     573      // Yaw angle consideration
     574      Matrix SXYZ(3, 3);
     575      SXYZ.Column(1) = sx;
     576      SXYZ.Column(2) = sy;
     577      SXYZ.Column(3) = sz;
     578      SXYZ = DotProduct(t_astro::rotZ(yaw), SXYZ);
     579      sx = SXYZ.Column(1);
     580      sy = SXYZ.Column(2);
     581      sz = SXYZ.Column(3);
     582    }
    403583    // Effective Dipole of the GPS Satellite Antenna
    404584    // ---------------------------------------------
    405     ColumnVector dipSat = sx - rho * DotProduct(rho,sx) - crossproduct(rho, sy);
     585    ColumnVector dipSat = sx - rho * DotProduct(rho, sx)
     586        - crossproduct(rho, sy);
    406587
    407588    // Receiver unit Vectors rx, ry
    … …  
    409590    ColumnVector rx(3);
    410591    ColumnVector ry(3);
    411     double recEll[3]; xyz2ell(rRec.data(), recEll) ;
     592    double recEll[3];
     593    xyz2ell(rRec.data(), recEll);
    412594    double neu[3];
    413595
    … …  
    417599    neu2xyz(recEll, neu, rx.data());
    418600
    419     neu[0] =  0.0;
     601    neu[0] = 0.0;
    420602    neu[1] = -1.0;
    421     neu[2] =  0.0;
     603    neu[2] = 0.0;
    422604    neu2xyz(recEll, neu, ry.data());
    423605
    424606    // Effective Dipole of the Receiver Antenna
    425607    // ----------------------------------------
    426     ColumnVector dipRec = rx - rho * DotProduct(rho,rx) + crossproduct(rho, ry);
     608    ColumnVector dipRec = rx - rho * DotProduct(rho, rx)
     609        + crossproduct(rho, ry);
    427610
    428611    // Resulting Effect
    429612    // ----------------
    430     double alpha = DotProduct(dipSat,dipRec) /
    431                       (dipSat.norm_Frobenius() * dipRec.norm_Frobenius());
    432 /*
    433     if (alpha >  1.0) alpha =  1.0;
    434     if (alpha < -1.0) alpha = -1.0;
    435 */
     613    double alpha = DotProduct(dipSat, dipRec)
     614        / (dipSat.NormFrobenius() * dipRec.NormFrobenius());
     615
     616    if (alpha > 1.0)
     617      alpha = 1.0;
     618    if (alpha < -1.0)
     619      alpha = -1.0;
     620
    436621    double dphi = acos(alpha) / 2.0 / M_PI;  // in cycles
    437622
    438     if ( DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0 ) {
     623    if (DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0) {
    439624      dphi = -dphi;
    440625    }
    … …  
    467652  double height = ell[2];
    468653
    469   double pp =  1013.25 * pow(1.0 - 2.26e-5 * height, 5.225);
    470   double TT =  18.0 - height * 0.0065 + 273.15;
    471   double hh =  50.0 * exp(-6.396e-4 * height);
    472   double ee =  hh / 100.0 * exp(-37.2465 + 0.213166*TT - 0.000256908*TT*TT);
     654  double pp = 1013.25 * pow(1.0 - 2.26e-5 * height, 5.225);
     655  double TT = 18.0 - height * 0.0065 + 273.15;
     656  double hh = 50.0 * exp(-6.396e-4 * height);
     657  double ee = hh / 100.0
     658      * exp(-37.2465 + 0.213166 * TT - 0.000256908 * TT * TT);
    473659
    474660  double h_km = height / 1000.0;
    475661
    476   if (h_km < 0.0) h_km = 0.0;
    477   if (h_km > 5.0) h_km = 5.0;
    478   int    ii   = int(h_km + 1); if (ii > 5) ii = 5;
     662  if (h_km < 0.0)
     663    h_km = 0.0;
     664  if (h_km > 5.0)
     665    h_km = 5.0;
     666  int ii = int(h_km + 1);
     667  if (ii > 5)
     668    ii = 5;
    479669  double href = ii - 1;
    480670
    … …  
    487677  bCor[5] = 0.563;
    488678
    489   double BB = bCor[ii-1] + (bCor[ii]-bCor[ii-1]) * (h_km - href);
    490 
    491   double zen  = M_PI/2.0 - Ele;
    492 
    493   return (0.002277/cos(zen)) * (pp + ((1255.0/TT)+0.05)*ee - BB*(tan(zen)*tan(zen)));
     679  double BB = bCor[ii - 1] + (bCor[ii] - bCor[ii - 1]) * (h_km - href);
     680
     681  double zen = M_PI / 2.0 - Ele;
     682
     683  return (0.002277 / cos(zen))
     684      * (pp + ((1255.0 / TT) + 0.05) * ee - BB * (tan(zen) * tan(zen)));
    494685}
    495686
    … …  
    500691}
    501692
    502 t_iono::~t_iono() {}
     693t_iono::~t_iono() {
     694}
    503695
    504696double t_iono::stec(const t_vTec* vTec, double signalPropagationTime,
    505       const ColumnVector& rSat, const bncTime& epochTime,
    506       const ColumnVector& xyzSta) {
     697    const ColumnVector& rSat, const bncTime& epochTime,
     698    const ColumnVector& xyzSta) {
    507699
    508700  // Latitude, longitude, height are defined with respect to a spherical earth model
    … …  
    524716  // -------------------------------------------------------------
    525717  ColumnVector rhoV = xyzSat - xyzSta;
    526   double rho = rhoV.norm_Frobenius();
     718  double rho = rhoV.NormFrobenius();
    527719  ColumnVector neu(3);
    528720  xyz2neu(geocSta.data(), rhoV.data(), neu.data());
    529   double sphEle = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / rho );
     721  double sphEle = acos(sqrt(neu[0] * neu[0] + neu[1] * neu[1]) / rho);
    530722  if (neu[2] < 0) {
    531723    sphEle *= -1.0;
    … …  
    537729  double stec = 0.0;
    538730  for (unsigned ii = 0; ii < vTec->_layers.size(); ii++) {
    539     piercePoint(vTec->_layers[ii]._height, epoch, geocSta.data(), sphEle, sphAzi);
     731    piercePoint(vTec->_layers[ii]._height, epoch, geocSta.data(), sphEle,
     732        sphAzi);
    540733    double vtec = vtecSingleLayerContribution(vTec->_layers[ii]);
    541734    stec += vtec * sin(sphEle + _psiPP);
    … …  
    547740
    548741  double vtec = 0.0;
    549   int N = vTecLayer._C.Nrows()-1;
    550   int M = vTecLayer._C.Ncols()-1;
     742  int N = vTecLayer._C.Nrows() - 1;
     743  int M = vTecLayer._C.Ncols() - 1;
    551744  double fac;
    552745
    … …  
    576769}
    577770
    578 void t_iono::piercePoint(double layerHeight, double epoch, const double* geocSta,
     771void t_iono::piercePoint(double layerHeight, double epoch,
     772    const double* geocSta,
    579773    double sphEle, double sphAzi) {
    580774
    581775  double q = (t_CST::rgeoc + geocSta[2]) / (t_CST::rgeoc + layerHeight);
    582776
    583   _psiPP = M_PI/2 - sphEle - asin(q * cos(sphEle));
    584 
    585   _phiPP = asin(sin(geocSta[0]) * cos(_psiPP) + cos(geocSta[0]) * sin(_psiPP) * cos(sphAzi));
    586 
    587   if (( (geocSta[0]*180.0/M_PI > 0) && (  tan(_psiPP) * cos(sphAzi)  > tan(M_PI/2 - geocSta[0])) )  ||
    588       ( (geocSta[0]*180.0/M_PI < 0) && (-(tan(_psiPP) * cos(sphAzi)) > tan(M_PI/2 + geocSta[0])) ))  {
    589     _lambdaPP = geocSta[1] + M_PI - asin((sin(_psiPP) * sin(sphAzi) / cos(_phiPP)));
    590   } else {
    591     _lambdaPP = geocSta[1]        + asin((sin(_psiPP) * sin(sphAzi) / cos(_phiPP)));
    592   }
    593 
    594   _lonS = fmod((_lambdaPP + (epoch - 50400) * M_PI / 43200), 2*M_PI);
     777  _psiPP = M_PI / 2 - sphEle - asin(q * cos(sphEle));
     778
     779  _phiPP = asin(
     780      sin(geocSta[0]) * cos(_psiPP)
     781          + cos(geocSta[0]) * sin(_psiPP) * cos(sphAzi));
     782
     783  if (((geocSta[0] * 180.0 / M_PI > 0)
     784      && (tan(_psiPP) * cos(sphAzi) > tan(M_PI / 2 - geocSta[0])))
     785      ||
     786      ((geocSta[0] * 180.0 / M_PI < 0)
     787          && (-(tan(_psiPP) * cos(sphAzi)) > tan(M_PI / 2 + geocSta[0])))) {
     788    _lambdaPP = geocSta[1] + M_PI
     789        - asin((sin(_psiPP) * sin(sphAzi) / cos(_phiPP)));
     790  }
     791  else {
     792    _lambdaPP = geocSta[1] + asin((sin(_psiPP) * sin(sphAzi) / cos(_phiPP)));
     793  }
     794
     795  _lonS = fmod((_lambdaPP + (epoch - 50400) * M_PI / 43200), 2 * M_PI);
    595796
    596797  return;

  • trunk/BNC/src/pppModel.h

    r8619 r8905  
    55#include <newmat.h>
    66#include <iostream>
     7#include <string>
    78#include "bnctime.h"
    89#include "t_prn.h"
    … …  
    2122
    2223 private:
    23   static const double RHO_DEG;
    24   static const double RHO_SEC;
    25   static const double MJD_J2000;
    26 
    2724  static double GMST(double Mjd_UT1);
    2825  static Matrix NutMatrix(double Mjd_TT);
    … …  
    3229class t_tides {
    3330 public:
    34   t_tides() {
    35     _lastMjd = 0.0;
    36     _rSun    = 0.0;
    37     _rMoon   = 0.0;
    38   }
    39   ~t_tides() {}
    40   ColumnVector displacement(const bncTime& time, const ColumnVector& xyz);
     31  t_tides();
     32  ~t_tides();
     33  ColumnVector earth(const bncTime& time, const ColumnVector& xyz);
     34  ColumnVector ocean(const bncTime& time, const ColumnVector& xyz, const std::string& station);
     35  t_irc        readBlqFile(const char* fileName);
     36  void         printAllBlqSets() const;
    4137 private:
    4238  double       _lastMjd;
    … …  
    4541  double       _rSun;
    4642  double       _rMoon;
    47 };
    4843
    49 class t_loading {
    50  public:
    51   t_loading(const QString& fileName);
    52   ~t_loading();
    53   t_irc   readFile(const QString& fileName);
    54   void    printAll() const;
    55 
    56  private:
    5744  class t_blqData {
    5845   public:
    … …  
    6350  t_blqData*                 newBlqData;
    6451  QMap <QString, t_blqData*> blqMap;
    65 
     52  void         printBlqSet(const std::string& station, t_blqData* blq);
    6653};
    6754
    … …  
    9885  double _lambdaPP;
    9986  double _lonS;
    100 
    101 
    10287};
    10388

  • trunk/BNC/src/pppOptions.cpp

    r7048 r8905  
    3535 * Created:    29-Jul-2014
    3636 *
    37  * Changes:   
     37 * Changes:
    3838 *
    3939 * -----------------------------------------------------------------------*/
    … …  
    5959}
    6060
    61 // 
     61//
    6262//////////////////////////////////////////////////////////////////////////////
    6363const std::vector<t_lc::type>& t_pppOptions::LCs(char system) const {
    … …  
    7777}
    7878
    79 // 
     79//
    8080//////////////////////////////////////////////////////////////////////////////
    8181bool t_pppOptions::useOrbClkCorr() const {
    … …  
    9999}
    100100
    101 // 
     101//
    102102/////////////////////////////////////////////////////////////////////////////
    103103vector<t_lc::type> t_pppOptions::ambLCs(char system) const {
    104  
     104
    105105  const vector<t_lc::type>& allLCs = LCs(system);
    106106  vector<t_lc::type>        phaseLCs;
    … …  
    146146  return answ;
    147147}
     148
     149

  • trunk/BNC/src/pppOptions.h

    r7961 r8905  
    1111class t_pppOptions {
    1212 public:
     13  enum e_type {IF, UncombPPP, PPPRTK, DCMcodeBias, DCMphaseBias};
    1314  t_pppOptions();
    1415  ~t_pppOptions();
    … …  
    1819  std::vector<t_lc::type>        ambLCs(char system) const;
    1920  std::vector<t_lc::type>        codeLCs(char system) const;
     21  std::vector<t_lc::type>        ionoLCs(char system) const;
    2022  bool useSystem(char system) const {return LCs(system).size() > 0;}
    2123  bool useOrbClkCorr() const;
    … …  
    2527  }
    2628
     29  e_type                  _obsModelType;
     30  QStringList             _obsmodelTypeStr = QStringList()
     31      << "IF PPP"
     32      << "Uncombined PPP"
     33      << "PPP-RTK"
     34      << "DCM with Code Biases"
     35      << "DCM with Phase Biases";
    2736  bool                    _realTime;
    2837  std::string             _crdFile;
    … …  
    4352  double                  _maxResC1;
    4453  double                  _maxResL1;
     54  double                  _sigmaGIMdiff;
     55  double                  _maxResGIMdiff;
    4556  bool                    _eleWgtCode;
    4657  bool                    _eleWgtPhase;
    … …  
    5263  double                  _aprSigTrp;
    5364  double                  _noiseTrp;
     65  double                  _aprSigIon;
     66  double                  _noiseIon;
     67  double                  _aprSigCodeBias;
     68  double                  _noiseCodeBias;
     69  double                  _aprSigPhaseBias;
     70  double                  _noisePhaseBias;
    5471  int                     _nmeaPort;
    5572  double                  _aprSigAmb;
    … …  
    5976  std::vector<t_lc::type> _LCsGalileo;
    6077  std::vector<t_lc::type> _LCsBDS;
     78  bool                    _pseudoObsIono;
     79  bool                    _refSatRequired;
    6180};
    6281

  • trunk/BNC/src/pppWidgets.cpp

    r8773 r8905  
    8282  _lcGalileo    = new QComboBox();     _lcGalileo   ->setObjectName("PPP/lcGalileo");    _widgets << _lcGalileo;
    8383  _lcBDS        = new QComboBox();     _lcBDS       ->setObjectName("PPP/lcBDS");        _widgets << _lcBDS;
     84  _modelObs     = new QComboBox();     _modelObs    ->setObjectName("PPP/modelObs");     _widgets << _modelObs;
     85  _pseudoObs    = new QComboBox();     _pseudoObs   ->setObjectName("PPP/pseudoObs");    _widgets << _pseudoObs;
    8486  _sigmaC1      = new QLineEdit();     _sigmaC1     ->setObjectName("PPP/sigmaC1");      _widgets << _sigmaC1;
    8587  _sigmaL1      = new QLineEdit();     _sigmaL1     ->setObjectName("PPP/sigmaL1");      _widgets << _sigmaL1;
    … …  
    111113  _dataSource->addItems(QString(",Real-Time Streams,RINEX Files").split(","));
    112114  connect(_dataSource, SIGNAL(currentIndexChanged(const QString&)), this, SLOT(slotEnableWidgets()));
    113 
    114   connect(_snxtroPath, SIGNAL(textChanged(const QString &)),
    115          this, SLOT(slotPPPTextChanged()));
    116 
    117   connect(_snxtroAc, SIGNAL(textChanged(const QString &)),
    118          this, SLOT(slotPPPTextChanged()));
    119 
    120   connect(_snxtroSol, SIGNAL(textChanged(const QString &)),
    121          this, SLOT(slotPPPTextChanged()));
     115  connect(_modelObs, SIGNAL(currentIndexChanged(const QString&)), this, SLOT(slotEnableWidgets()));
     116  connect(_snxtroPath, SIGNAL(textChanged(const QString &)), this, SLOT(slotPPPTextChanged()));
     117  connect(_snxtroAc, SIGNAL(textChanged(const QString &)), this, SLOT(slotPPPTextChanged()));
     118  connect(_snxtroSol, SIGNAL(textChanged(const QString &)), this, SLOT(slotPPPTextChanged()));
    122119
    123120  slotEnableWidgets();
    … …  
    127124  _lcGPS->addItems(QString("P3,P3&L3").split(","));
    128125#else
    129   _lcGPS->addItems(QString("no,Pi,Li,Pi&Li,P3,L3,P3&L3").split(","));
     126  _lcGPS->addItems(QString("no,Pi,Li,Pi&Li").split(","));
    130127#endif
    131128
    … …  
    134131   _lcGLONASS->addItems(QString("no,P3,L3,P3&L3").split(","));
    135132#else
    136   _lcGLONASS->addItems(QString("no,Pi,Li,Pi&Li,P3,L3,P3&L3").split(","));
     133  _lcGLONASS->addItems(QString("no,Pi,Li,Pi&Li").split(","));
    137134#endif
    138135
    … …  
    141138  _lcGalileo->addItems(QString("no,P3,L3,P3&L3").split(","));
    142139#else
    143   _lcGalileo->addItems(QString("no,Pi,Li,Pi&Li,P3,L3,P3&L3").split(","));
     140  _lcGalileo->addItems(QString("no,Pi,Li,Pi&Li").split(","));
    144141#endif
    145142
    … …  
    148145  _lcBDS->addItems(QString("no,P3,L3,P3&L3").split(","));
    149146#else
    150   _lcBDS->addItems(QString("no,Pi,Li,Pi&Li,P3,L3,P3&L3").split(","));
     147  _lcBDS->addItems(QString("no,Pi,Li,Pi&Li").split(","));
     148#endif
     149
     150  _modelObs->setEditable(false);
     151  _pseudoObs->setEditable(false);
     152#ifdef USE_PPP_SSR_I
     153  _modelObs->addItems(QString("Ionosphere-free PPP").split(","));
     154  _pseudoObs->addItems(QString("no").split(","));
     155#else
     156  _modelObs->addItems(QString("Ionosphere-free PPP,Uncombined PPP,PPP-RTK,DCM with Code Biases,DCM with Phase Biases").split(","));
     157  _pseudoObs->addItems(QString("no,Ionosphere").split(","));
    151158#endif
    152159
    … …  
    247254  delete _lcGalileo;
    248255  delete _lcBDS;
     256  delete _modelObs;
     257  delete _pseudoObs;
    249258  delete _sigmaC1;
    250259  delete _sigmaL1;
    … …  
    297306    _lcBDS->setCurrentIndex(ii);
    298307  }
     308  ii = _modelObs->findText(settings.value(_modelObs->objectName()).toString());
     309  if (ii != -1) {
     310    _modelObs->setCurrentIndex(ii);
     311  }
     312  ii = _pseudoObs->findText(settings.value(_pseudoObs->objectName()).toString());
     313  if (ii != -1) {
     314    _pseudoObs->setCurrentIndex(ii);
     315  }
    299316  ii = _snxtroIntr->findText(settings.value(_snxtroIntr->objectName()).toString());
    300317  if (ii != -1) {
    … …  
    433450  settings.setValue(_lcGalileo   ->objectName(), _lcGalileo   ->currentText());
    434451  settings.setValue(_lcBDS       ->objectName(), _lcBDS       ->currentText());
     452  settings.setValue(_modelObs    ->objectName(), _modelObs    ->currentText());
     453  settings.setValue(_pseudoObs   ->objectName(), _pseudoObs   ->currentText());
    435454  settings.setValue(_sigmaC1     ->objectName(), _sigmaC1     ->text());
    436455  settings.setValue(_sigmaL1     ->objectName(), _sigmaL1     ->text());
    … …  
    477496  bool realTime    = _dataSource->currentText() == "Real-Time Streams";
    478497  bool rinexFiles  = _dataSource->currentText() == "RINEX Files";
     498  bool enablePseudoObs;
     499  if (_modelObs->currentText() == "PPP-RTK" ||
     500    _modelObs->currentText() == "Ionosphere-free PPP") {
     501    enablePseudoObs = false;
     502  }
     503  else {
     504    enablePseudoObs = true;
     505  }
    479506
    480507  QListIterator<QWidget*> it(_widgets);
    … …  
    490517  }
    491518  else if (rinexFiles) {
    492     _corrMount->setEnabled(false);
    493 //  _plotCoordinates->setEnabled(false);
    494 //  _audioResponse->setEnabled(false);
     519    _corrMount    ->setEnabled(false);
     520    _audioResponse->setEnabled(false);
    495521  }
    496522
    … …  
    506532    _snxtroAc   ->setEnabled(false);
    507533    _snxtroSol  ->setEnabled(false);
     534  }
     535
     536  if (enablePseudoObs) {
     537    _pseudoObs->setEnabled(true);
     538  } else {
     539    _pseudoObs->setEnabled(false);
    508540  }
    509541
    … …  
    592624    }
    593625  }
    594 }
     626
     627
     628}

  • trunk/BNC/src/pppWidgets.h

    r8403 r8905  
    6666  QComboBox*     _lcGalileo;
    6767  QComboBox*     _lcBDS;
     68  QComboBox*     _modelObs;
     69  QComboBox*     _pseudoObs;
    6870  QLineEdit*     _sigmaC1;
    6971  QLineEdit*     _sigmaL1;

  • trunk/BNC/src/src.pri

    r8267 r8905  
    122122  HEADERS += PPP/pppClient.h    PPP/pppObsPool.h   PPP/pppEphPool.h   \
    123123             PPP/pppStation.h   PPP/pppFilter.h    PPP/pppParlist.h   \
    124              PPP/pppSatObs.h
     124             PPP/pppSatObs.h    PPP/pppRefSat.h
    125125  SOURCES += PPP/pppClient.cpp  PPP/pppObsPool.cpp PPP/pppEphPool.cpp \
    126126             PPP/pppStation.cpp PPP/pppFilter.cpp  PPP/pppParlist.cpp \
    127              PPP/pppSatObs.cpp
     127             PPP/pppSatObs.cpp  PPP/pppRefSat.cpp
    128128}
    129129else {
Note: See TracChangeset for help on using the changeset viewer.