Changeset 8905 in ntrip for trunk/BNC/src/PPP/pppSatObs.cpp


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

some developments regarding PPP, not completed!

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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}
Note: See TracChangeset for help on using the changeset viewer.