Changeset 9481 in ntrip for trunk


Ignore:
Timestamp:
Jul 19, 2021, 10:46:16 PM (3 years ago)
Author:
stuerze
Message:

minor changes

Location:
trunk/BNC/src
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r8542 r9481  
    9999    satData->P2       = 0.0;
    100100    satData->P5       = 0.0;
     101    satData->P6       = 0.0;
    101102    satData->P7       = 0.0;
    102103    satData->L1       = 0.0;
    103104    satData->L2       = 0.0;
    104105    satData->L5       = 0.0;
     106    satData->L6       = 0.0;
    105107    satData->L7       = 0.0;
    106108    for (unsigned ifrq = 0; ifrq < obs->_obs.size(); ifrq++) {
    … …  
    115117            cb  = bias._value;
    116118          }
     119          // FIXME: use C/Q bias for X observations
     120          // qDebug() << satData->prn << frqObs->_rnxType2ch.c_str();
    117121        }
    118122      }
    … …  
    130134        if (frqObs->_codeValid)  satData->P5       = frqObs->_code + cb;
    131135        if (frqObs->_phaseValid) satData->L5       = frqObs->_phase;
     136        if (frqObs->_slip)       satData->slipFlag = true;
     137      }
     138      else if (frqObs->_rnxType2ch[0] == '6') {
     139        if (frqObs->_codeValid)  satData->P6       = frqObs->_code + cb;
     140        if (frqObs->_phaseValid) satData->L6       = frqObs->_phase;
    132141        if (frqObs->_slip)       satData->slipFlag = true;
    133142      }
    … …  
    276285  // ---------------------
    277286  else if (satData->system() == 'C' && _opt->useSystem('C')) {
    278     if (satData->P2 != 0.0 && satData->P7 != 0.0 &&
    279         satData->L2 != 0.0 && satData->L7 != 0.0 ) {
     287    if (satData->P2 != 0.0 && satData->P6 != 0.0 &&
     288        satData->L2 != 0.0 && satData->L6 != 0.0 ) {
    280289      double f2 = t_CST::freq(t_frequency::C2, 0);
    281       double f7 = t_CST::freq(t_frequency::C7, 0);
    282       double a2 =   f2 * f2 / (f2 * f2 - f7 * f7);
    283       double a7 = - f7 * f7 / (f2 * f2 - f7 * f7);
     290      double f6 = t_CST::freq(t_frequency::C6, 0);
     291      double a2 =   f2 * f2 / (f2 * f2 - f6 * f6);
     292      double a6 = - f6 * f6 / (f2 * f2 - f6 * f6);
    284293      satData->L2      = satData->L2 * t_CST::c / f2;
    285       satData->L7      = satData->L7 * t_CST::c / f7;
    286       satData->P3      = a2 * satData->P2 + a7 * satData->P7;
    287       satData->L3      = a2 * satData->L2 + a7 * satData->L7;
    288       satData->lambda3 = a2 * t_CST::c / f2 + a7 * t_CST::c / f7;
     294      satData->L6      = satData->L6 * t_CST::c / f6;
     295      satData->P3      = a2 * satData->P2 + a6 * satData->P6;
     296      satData->L3      = a2 * satData->L2 + a6 * satData->L6;
     297      satData->lambda3 = a2 * t_CST::c / f2 + a6 * t_CST::c / f6;
    289298      satData->lkA     = a2;
    290       satData->lkB     = a7;
     299      satData->lkB     = a6;
    291300      _epoData->satData[satData->prn] = satData;
    292301    }

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

    r9418 r9481  
    5555using namespace std;
    5656
    57 const double   MAXRES_CODE           = 2.98 * 3.0;
    58 const double   MAXRES_PHASE_GPS      = 0.04;
    59 const double   MAXRES_PHASE_GLONASS  = 2.98 * 0.03;
    6057const double   GLONASS_WEIGHT_FACTOR = 5.0;
    61 const double   BDS_WEIGHT_FACTOR     = 5.0;
     58const double   BDS_WEIGHT_FACTOR     = 2.0; // 5.0;
    6259
    6360#define LOG (_pppClient->log())
    … …  
    352349
    353350  double offset = 0.0;
    354   t_frequency::type frqA = t_frequency::G1;
    355   t_frequency::type frqB = t_frequency::G2;
     351
     352  t_frequency::type frqA;
     353  t_frequency::type frqB;
     354
    356355  if      (satData->prn[0] == 'R') {
    357356    offset = Glonass_offset();
    … …  
    367366    offset = Bds_offset();
    368367    frqA = t_frequency::C2;
    369     frqB = t_frequency::C7;
    370   }
     368    frqB = t_frequency::C6;
     369  }
     370  else {
     371    frqA = t_frequency::G1;
     372    frqB = t_frequency::G2;
     373  }
     374
    371375  double phaseCenter = 0.0;
     376
    372377  if (_antex) {
     378
     379    // Satellite correction
     380    // ---------------------
     381    double elTx,azTx;
     382
     383    // LOS unit vector satellite --> receiver
     384    ColumnVector rho = xRec - satData->xx;
     385    rho /= rho.norm_Frobenius();
     386
     387    // Sun unit vector
     388    ColumnVector xSun = t_astro::Sun(satData->tt.mjd());
     389    xSun /= xSun.norm_Frobenius();
     390
     391    // Satellite unit vectors sz, sy, sx
     392    ColumnVector sz = -satData->xx / satData->xx.norm_Frobenius();
     393    ColumnVector sy = crossproduct(sz, xSun);
     394    ColumnVector sx = crossproduct(sy, sz);
     395
     396    sx /= sx.norm_Frobenius();
     397    sy /= sy.norm_Frobenius();
     398
     399    // LOS vector in satellite frame
     400    ColumnVector u(3);
     401    u(1) = dotproduct(sx, rho);
     402    u(2) = dotproduct(sy, rho);
     403    u(3) = dotproduct(sz, rho);
     404
     405    // Azimuth and elevation in satellite antenna frame
     406    elTx = atan2(u(3),sqrt(pow(u(2),2)+pow(u(1),2)));
     407    azTx = atan2(u(2),u(1));
     408
    373409    bool found;
    374     phaseCenter = satData->lkA * _antex->rcvCorr(OPT->_antNameRover, frqA,
    375                                                  satData->eleSat, satData->azSat,
    376                                                  found)
    377                 + satData->lkB * _antex->rcvCorr(OPT->_antNameRover, frqB,
    378                                                  satData->eleSat, satData->azSat,
    379                                                  found);
     410    if (OPT->_isAPC) {
     411      phaseCenter += satData->lkB * _antex->satCorr(satData->prn, frqA, elTx, azTx, found);
     412    }
     413    else {
     414      phaseCenter += satData->lkA * _antex->satCorr(satData->prn, frqA, elTx, azTx, found);
     415    }
    380416    if (!found) {
    381       LOG << "ANTEX: antenna >" << OPT->_antNameRover << "< not found\n";
     417      LOG << "ANTEX: antenna >" << satData->prn.mid(0,3).toStdString() << " " << frqA << "< not found\n";
     418    }
     419
     420    phaseCenter += satData->lkB * _antex->satCorr(satData->prn, frqB, elTx, azTx, found);
     421    if (!found) {
     422      LOG << "ANTEX: antenna >" << satData->prn.mid(0,3).toStdString() << " " << frqB << "< not found\n";
     423    }
     424
     425    // Receiver correction
     426    // -------------------
     427
     428    phaseCenter += satData->lkA * _antex->rcvCorr(OPT->_antNameRover, frqA,
     429                                                  satData->eleSat, satData->azSat, found);
     430    if (!found) {
     431      phaseCenter += satData->lkA * _antex->rcvCorr(OPT->_antNameRover, t_frequency::G1,
     432                                                    satData->eleSat, satData->azSat, found);
     433    }
     434    if (!found) {
     435      LOG << "ANTEX: antenna >" << OPT->_antNameRover << " " << frqA << "< not found\n";
     436    }
     437
     438    phaseCenter += satData->lkB * _antex->rcvCorr(OPT->_antNameRover, frqB,
     439                                                  satData->eleSat, satData->azSat, found);
     440    if (!found) {
     441      phaseCenter += satData->lkB * _antex->rcvCorr(OPT->_antNameRover, t_frequency::G2,
     442                                                    satData->eleSat, satData->azSat, found);
     443    }
     444    if (!found) {
     445      LOG << "ANTEX: antenna >" << OPT->_antNameRover << " " << frqB << "< not found\n";
    382446    }
    383447  }
    … …  
    535599      // --------------
    536600      else if (pp->type == t_pppParam::GALILEO_OFFSET) {
    537         _QQ(iPar,iPar) += 0.1 * 0.1;
     601        if (_QQ(iPar,iPar)>pow(1000.0,2))
     602          _QQ(iPar,iPar) = 1000.0 * 1000.0;
     603        else
     604          _QQ(iPar,iPar) += 0.1 * 0.1;
    538605      }
    539606
    … …  
    541608      // ----------
    542609      else if (pp->type == t_pppParam::BDS_OFFSET) {
    543         _QQ(iPar,iPar) += 0.1 * 0.1;    //TODO: TEST
     610        if (_QQ(iPar,iPar)>pow(1000.0,2))
     611          _QQ(iPar,iPar) = 1000.0 * 1000.0;
     612        else
     613          _QQ(iPar,iPar) += 0.1 * 0.1;
    544614      }
    545615    }
    … …  
    746816}
    747817
     818// Iono combi noise factor
     819////////////////////////////////////////////////////////////////////////////
     820double ionFac(const QString prn, QMap<QString, t_satData*>& satData) {
     821  if (satData.contains(prn))
     822    return sqrt(pow(satData.value(prn)->lkA,2) +
     823                pow(satData.value(prn)->lkB,2)  );
     824  else
     825    return 0.0;
     826};
     827
    748828// Outlier Detection
    749829////////////////////////////////////////////////////////////////////////////
    750830QString t_pppFilter::outlierDetection(int iPhase, const ColumnVector& vv,
    751                                    QMap<QString, t_satData*>& satData) {
     831                                      QMap<QString, t_satData*>& satData) {
    752832
    753833  Tracer tracer("t_pppFilter::outlierDetection");
    … …  
    755835  QString prnGPS;
    756836  QString prnGlo;
     837
     838  double  ionFacGPS;
     839  double  ionFacGLO;
     840
    757841  double  maxResGPS = 0.0; // GPS + Galileo
    758842  double  maxResGlo = 0.0; // GLONASS + BDS
     843
    759844  findMaxRes(vv, satData, prnGPS, prnGlo, maxResGPS, maxResGlo);
    760845
     846  ionFacGLO = ionFac(prnGlo,satData);
     847  if (iPhase == 0)
     848    ionFacGLO *= (prnGlo[0]=='R'? GLONASS_WEIGHT_FACTOR : BDS_WEIGHT_FACTOR);
     849  ionFacGPS = ionFac(prnGPS,satData);
     850
    761851  if      (iPhase == 1) {
    762     if      (maxResGlo > 2.98 * OPT->_maxResL1) {
     852    if      (maxResGlo > ionFacGLO * OPT->_maxResL1) {
    763853      LOG << "Outlier Phase " << prnGlo.mid(0,3).toLatin1().data() << ' ' << maxResGlo << endl;
    764854      return prnGlo;
    765855    }
    766     else if (maxResGPS > MAXRES_PHASE_GPS) {
     856    else if (maxResGPS > ionFacGPS * OPT->_maxResL1) {
    767857      LOG << "Outlier Phase " << prnGPS.mid(0,3).toLatin1().data() << ' ' << maxResGPS << endl;
    768858      return prnGPS;
    769859    }
    770860  }
    771   else if (iPhase == 0 && maxResGPS > 2.98 * OPT->_maxResC1) {
    772     LOG << "Outlier Code  " << prnGPS.mid(0,3).toLatin1().data() << ' ' << maxResGPS << endl;
    773     return prnGPS;
     861  else if (iPhase == 0) {
     862    if (maxResGlo > ionFacGLO * OPT->_maxResC1) {
     863      LOG << "Outlier Code  " << prnGlo.mid(0,3).toLatin1().data() << ' ' << maxResGlo << endl;
     864      return prnGlo;
     865    }
     866    else if (maxResGPS > ionFacGPS * OPT->_maxResC1) {
     867      LOG << "Outlier Code  " << prnGPS.mid(0,3).toLatin1().data() << ' ' << maxResGPS << endl;
     868      return prnGPS;
     869    }
    774870  }
    775871
    … …  
    780876///////////////////////////////////////////////////////////////////////////
    781877double t_pppFilter::windUp(const QString& prn, const ColumnVector& rSat,
    782                         const ColumnVector& rRec) {
     878                           const ColumnVector& rRec) {
    783879
    784880  Tracer tracer("t_pppFilter::windUp");
    … …  
    9191015  satData->obsIndex = iObs;
    9201016
     1017  // Iono-free combination noise factor
     1018  // ----------------------------------
     1019  double ionFac = sqrt(pow(satData->lkA,2) + pow(satData->lkB,2));
     1020
    9211021  // Phase Observations
    9221022  // ------------------
    9231023
    9241024  if (iPhase == 1) {
    925     ll(iObs)      = satData->L3 - cmpValue(satData, true);
    926     double sigL3 = 2.98 * OPT->_sigmaL1;
     1025    double sigL3 = ionFac * ellWgtCoef * OPT->_sigmaL1;
    9271026    if (satData->system() == 'R') {
    9281027      sigL3 *= GLONASS_WEIGHT_FACTOR;
    … …  
    9311030      sigL3 *= BDS_WEIGHT_FACTOR;
    9321031    }
    933     PP(iObs,iObs) = 1.0 / (sigL3 * sigL3) / (ellWgtCoef * ellWgtCoef);
     1032    satData->L3sig = sigL3;
     1033    ll(iObs)      = satData->L3 - cmpValue(satData, true);
     1034    PP(iObs,iObs) = 1.0 / (sigL3 * sigL3);
    9341035    for (int iPar = 1; iPar <= _params.size(); iPar++) {
    9351036      if (_params[iPar-1]->type == t_pppParam::AMB_L3 &&
    … …  
    9441045  // -----------------
    9451046  else {
    946     double sigP3 = 2.98 * OPT->_sigmaC1;
     1047    double sigP3 = ionFac * ellWgtCoef * OPT->_sigmaC1;
     1048    if (satData->system() == 'R') {
     1049      sigP3 *= GLONASS_WEIGHT_FACTOR;
     1050    }
     1051    if  (satData->system() == 'C') {
     1052      sigP3 *= BDS_WEIGHT_FACTOR;
     1053    }
     1054    satData->P3sig = sigP3;
    9471055    ll(iObs)      = satData->P3 - cmpValue(satData, false);
    948     PP(iObs,iObs) = 1.0 / (sigP3 * sigP3) / (ellWgtCoef * ellWgtCoef);
     1056    PP(iObs,iObs) = 1.0 / (sigP3 * sigP3);
    9491057    for (int iPar = 1; iPar <= _params.size(); iPar++) {
    9501058      AA(iObs, iPar) = _params[iPar-1]->partial(satData, false);
    … …  
    9731081          << " RES " << satData->prn.mid(0,3).toLatin1().data()
    9741082          << (iPhase ? "   L3 " : "   P3 ")
    975           << setw(9) << setprecision(4) << vv(satData->obsIndex) << endl;
     1083          << setw(9) << setprecision(4) << vv(satData->obsIndex) << " "
     1084          << setprecision(3)
     1085          << setw(7) << (iPhase? satData->L3sig : satData->P3sig) << " "
     1086          << setprecision(1)
     1087          << setw(5) <<satData->eleSat * 180 / M_PI
     1088          << endl;
    9761089    }
    9771090  }
    … …  
    11661279}
    11671280
    1168 // Remeber Original State Vector and Variance-Covariance Matrix
     1281// Remember Original State Vector and Variance-Covariance Matrix
    11691282////////////////////////////////////////////////////////////////////////////
    11701283void t_pppFilter::rememberState(t_epoData* epoData) {

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

    r8252 r9481  
    5252    P2       = 0.0;
    5353    P5       = 0.0;
     54    P6       = 0.0;
    5455    P7       = 0.0;
    5556    P3       = 0.0;
     57    P3sig    = 0.0;
    5658    L1       = 0.0;
    5759    L2       = 0.0;
    5860    L5       = 0.0;
     61    L6       = 0.0;
    5962    L7       = 0.0;
    6063    L3       = 0.0;
     64    L3sig    = 0.0;
    6165    lkA      = 0.0;
    6266    lkB      = 0.0;
    … …  
    7478  double       P2;
    7579  double       P5;
     80  double       P6;
    7681  double       P7;
    7782  double       P3;
     83  double       P3sig;
    7884  double       L1;
    7985  double       L2;
    8086  double       L5;
     87  double       L6;
    8188  double       L7;
    8289  double       L3;
     90  double       L3sig;
    8391  ColumnVector xx;
    8492  ColumnVector vv;

  • trunk/BNC/src/bncantex.cpp

    r8630 r9481  
    321321    frqType = t_frequency::R1;
    322322  }
     323  else if (prn[0] == 'E') {
     324    frqType = t_frequency::E1;
     325  }
     326  else if (prn[0] == 'C') {
     327    frqType = t_frequency::C2;
     328  }
     329  else if (prn[0] == 'S') {
     330    frqType = t_frequency::S1;
     331  }
     332  else if (prn[0] == 'J') {
     333    frqType = t_frequency::J1;
     334  }
     335  else if (prn[0] == 'I') {
     336    frqType = t_frequency::I5;
     337  }
    323338
    324339  QMap<QString, t_antMap*>::const_iterator it = _maps.find(prn.mid(0,3));
    … …  
    351366
    352367  return failure;
     368}
     369
     370//
     371////////////////////////////////////////////////////////////////////////////
     372double bncAntex::satCorr(const QString& prn, t_frequency::type frqType,
     373                         double elTx, double azTx, bool& found) const {
     374
     375  if (_maps.find(prn.mid(0,3)) == _maps.end()) {
     376    found = false;
     377    return 0.0;
     378  };
     379
     380  t_antMap* map = _maps[prn.mid(0,3)];
     381
     382  if (map->frqMap.find(frqType) == map->frqMap.end()) {
     383    found = false;
     384    return 0.0;
     385  };
     386
     387  t_frqMap* frqMap = map->frqMap[frqType];
     388
     389  double var = 0.0;
     390  if (frqMap->pattern.ncols() > 0) {
     391    double zenDiff = 999.999;
     392    double zenTx  = 90.0 - elTx * 180.0 / M_PI;
     393    unsigned iZen = 0;
     394    for (double zen = map->zen1; zen <= map->zen2; zen += map->dZen) {
     395      iZen += 1;
     396      double newZenDiff = fabs(zen - zenTx);
     397      if (newZenDiff < zenDiff) {
     398        zenDiff = newZenDiff;
     399        var = frqMap->pattern(iZen);
     400      }
     401    }
     402  }
     403
     404  found = true;
     405  return var - frqMap->neu[0] * cos(azTx)*cos(elTx)
     406             - frqMap->neu[1] * sin(azTx)*cos(elTx)
     407             - frqMap->neu[2] * sin(elTx);
     408
    353409}
    354410

  • trunk/BNC/src/bncantex.h

    r8078 r9481  
    4040  void    print() const;
    4141  QString pcoSinexString(const std::string& antName, t_frequency::type frqType);
     42  double  satCorr(const QString& prn, t_frequency::type frqType,
     43                  double eleSat, double azSat, bool& found) const;
    4244  double  rcvCorr(const std::string& antName, t_frequency::type frqType,
    4345                  double eleSat, double azSat, bool& found) const;

  • trunk/BNC/src/pppMain.cpp

    r9442 r9481  
    164164    if (_realTime) {
    165165      opt->_corrMount.assign(settings.value("PPP/corrMount").toString().toStdString());
     166      opt->_isAPC = (opt->_corrMount.substr(0,4)=="SSRA");
    166167    }
    167168    else {
    … …  
    169170      opt->_rinexNav.assign(settings.value("PPP/rinexNav").toString().toStdString());
    170171      opt->_corrFile.assign(settings.value("PPP/corrFile").toString().toStdString());
     172      QFileInfo tmp = QFileInfo(QString::fromStdString(opt->_corrFile));
     173      opt->_isAPC = (tmp.baseName().mid(0,4)=="SSRA");
    171174    }
    172175

  • trunk/BNC/src/pppOptions.h

    r9439 r9481  
    3737  std::string             _crdFile;
    3838  std::string             _corrMount;
     39  bool                    _isAPC;
    3940  std::string             _rinexObs;
    4041  std::string             _rinexNav;
Note: See TracChangeset for help on using the changeset viewer.