Index: /trunk/BNC/src/PPP_SSR_I/pppClient.cpp =================================================================== --- /trunk/BNC/src/PPP_SSR_I/pppClient.cpp (revision 9480) +++ /trunk/BNC/src/PPP_SSR_I/pppClient.cpp (revision 9481) @@ -99,8 +99,10 @@ satData->P2 = 0.0; satData->P5 = 0.0; + satData->P6 = 0.0; satData->P7 = 0.0; satData->L1 = 0.0; satData->L2 = 0.0; satData->L5 = 0.0; + satData->L6 = 0.0; satData->L7 = 0.0; for (unsigned ifrq = 0; ifrq < obs->_obs.size(); ifrq++) { @@ -115,4 +117,6 @@ cb = bias._value; } + // FIXME: use C/Q bias for X observations + // qDebug() << satData->prn << frqObs->_rnxType2ch.c_str(); } } @@ -130,4 +134,9 @@ if (frqObs->_codeValid) satData->P5 = frqObs->_code + cb; if (frqObs->_phaseValid) satData->L5 = frqObs->_phase; + if (frqObs->_slip) satData->slipFlag = true; + } + else if (frqObs->_rnxType2ch[0] == '6') { + if (frqObs->_codeValid) satData->P6 = frqObs->_code + cb; + if (frqObs->_phaseValid) satData->L6 = frqObs->_phase; if (frqObs->_slip) satData->slipFlag = true; } @@ -276,17 +285,17 @@ // --------------------- else if (satData->system() == 'C' && _opt->useSystem('C')) { - if (satData->P2 != 0.0 && satData->P7 != 0.0 && - satData->L2 != 0.0 && satData->L7 != 0.0 ) { + if (satData->P2 != 0.0 && satData->P6 != 0.0 && + satData->L2 != 0.0 && satData->L6 != 0.0 ) { double f2 = t_CST::freq(t_frequency::C2, 0); - double f7 = t_CST::freq(t_frequency::C7, 0); - double a2 = f2 * f2 / (f2 * f2 - f7 * f7); - double a7 = - f7 * f7 / (f2 * f2 - f7 * f7); + double f6 = t_CST::freq(t_frequency::C6, 0); + double a2 = f2 * f2 / (f2 * f2 - f6 * f6); + double a6 = - f6 * f6 / (f2 * f2 - f6 * f6); satData->L2 = satData->L2 * t_CST::c / f2; - satData->L7 = satData->L7 * t_CST::c / f7; - satData->P3 = a2 * satData->P2 + a7 * satData->P7; - satData->L3 = a2 * satData->L2 + a7 * satData->L7; - satData->lambda3 = a2 * t_CST::c / f2 + a7 * t_CST::c / f7; + satData->L6 = satData->L6 * t_CST::c / f6; + satData->P3 = a2 * satData->P2 + a6 * satData->P6; + satData->L3 = a2 * satData->L2 + a6 * satData->L6; + satData->lambda3 = a2 * t_CST::c / f2 + a6 * t_CST::c / f6; satData->lkA = a2; - satData->lkB = a7; + satData->lkB = a6; _epoData->satData[satData->prn] = satData; } Index: /trunk/BNC/src/PPP_SSR_I/pppFilter.cpp =================================================================== --- /trunk/BNC/src/PPP_SSR_I/pppFilter.cpp (revision 9480) +++ /trunk/BNC/src/PPP_SSR_I/pppFilter.cpp (revision 9481) @@ -55,9 +55,6 @@ using namespace std; -const double MAXRES_CODE = 2.98 * 3.0; -const double MAXRES_PHASE_GPS = 0.04; -const double MAXRES_PHASE_GLONASS = 2.98 * 0.03; const double GLONASS_WEIGHT_FACTOR = 5.0; -const double BDS_WEIGHT_FACTOR = 5.0; +const double BDS_WEIGHT_FACTOR = 2.0; // 5.0; #define LOG (_pppClient->log()) @@ -352,6 +349,8 @@ double offset = 0.0; - t_frequency::type frqA = t_frequency::G1; - t_frequency::type frqB = t_frequency::G2; + + t_frequency::type frqA; + t_frequency::type frqB; + if (satData->prn[0] == 'R') { offset = Glonass_offset(); @@ -367,17 +366,82 @@ offset = Bds_offset(); frqA = t_frequency::C2; - frqB = t_frequency::C7; - } + frqB = t_frequency::C6; + } + else { + frqA = t_frequency::G1; + frqB = t_frequency::G2; + } + double phaseCenter = 0.0; + if (_antex) { + + // Satellite correction + // --------------------- + double elTx,azTx; + + // LOS unit vector satellite --> receiver + ColumnVector rho = xRec - satData->xx; + rho /= rho.norm_Frobenius(); + + // Sun unit vector + ColumnVector xSun = t_astro::Sun(satData->tt.mjd()); + xSun /= xSun.norm_Frobenius(); + + // Satellite unit vectors sz, sy, sx + ColumnVector sz = -satData->xx / satData->xx.norm_Frobenius(); + ColumnVector sy = crossproduct(sz, xSun); + ColumnVector sx = crossproduct(sy, sz); + + sx /= sx.norm_Frobenius(); + sy /= sy.norm_Frobenius(); + + // LOS vector in satellite frame + ColumnVector u(3); + u(1) = dotproduct(sx, rho); + u(2) = dotproduct(sy, rho); + u(3) = dotproduct(sz, rho); + + // Azimuth and elevation in satellite antenna frame + elTx = atan2(u(3),sqrt(pow(u(2),2)+pow(u(1),2))); + azTx = atan2(u(2),u(1)); + bool found; - phaseCenter = satData->lkA * _antex->rcvCorr(OPT->_antNameRover, frqA, - satData->eleSat, satData->azSat, - found) - + satData->lkB * _antex->rcvCorr(OPT->_antNameRover, frqB, - satData->eleSat, satData->azSat, - found); + if (OPT->_isAPC) { + phaseCenter += satData->lkB * _antex->satCorr(satData->prn, frqA, elTx, azTx, found); + } + else { + phaseCenter += satData->lkA * _antex->satCorr(satData->prn, frqA, elTx, azTx, found); + } if (!found) { - LOG << "ANTEX: antenna >" << OPT->_antNameRover << "< not found\n"; + LOG << "ANTEX: antenna >" << satData->prn.mid(0,3).toStdString() << " " << frqA << "< not found\n"; + } + + phaseCenter += satData->lkB * _antex->satCorr(satData->prn, frqB, elTx, azTx, found); + if (!found) { + LOG << "ANTEX: antenna >" << satData->prn.mid(0,3).toStdString() << " " << frqB << "< not found\n"; + } + + // Receiver correction + // ------------------- + + phaseCenter += satData->lkA * _antex->rcvCorr(OPT->_antNameRover, frqA, + satData->eleSat, satData->azSat, found); + if (!found) { + phaseCenter += satData->lkA * _antex->rcvCorr(OPT->_antNameRover, t_frequency::G1, + satData->eleSat, satData->azSat, found); + } + if (!found) { + LOG << "ANTEX: antenna >" << OPT->_antNameRover << " " << frqA << "< not found\n"; + } + + phaseCenter += satData->lkB * _antex->rcvCorr(OPT->_antNameRover, frqB, + satData->eleSat, satData->azSat, found); + if (!found) { + phaseCenter += satData->lkB * _antex->rcvCorr(OPT->_antNameRover, t_frequency::G2, + satData->eleSat, satData->azSat, found); + } + if (!found) { + LOG << "ANTEX: antenna >" << OPT->_antNameRover << " " << frqB << "< not found\n"; } } @@ -535,5 +599,8 @@ // -------------- else if (pp->type == t_pppParam::GALILEO_OFFSET) { - _QQ(iPar,iPar) += 0.1 * 0.1; + if (_QQ(iPar,iPar)>pow(1000.0,2)) + _QQ(iPar,iPar) = 1000.0 * 1000.0; + else + _QQ(iPar,iPar) += 0.1 * 0.1; } @@ -541,5 +608,8 @@ // ---------- else if (pp->type == t_pppParam::BDS_OFFSET) { - _QQ(iPar,iPar) += 0.1 * 0.1; //TODO: TEST + if (_QQ(iPar,iPar)>pow(1000.0,2)) + _QQ(iPar,iPar) = 1000.0 * 1000.0; + else + _QQ(iPar,iPar) += 0.1 * 0.1; } } @@ -746,8 +816,18 @@ } +// Iono combi noise factor +//////////////////////////////////////////////////////////////////////////// +double ionFac(const QString prn, QMap& satData) { + if (satData.contains(prn)) + return sqrt(pow(satData.value(prn)->lkA,2) + + pow(satData.value(prn)->lkB,2) ); + else + return 0.0; +}; + // Outlier Detection //////////////////////////////////////////////////////////////////////////// QString t_pppFilter::outlierDetection(int iPhase, const ColumnVector& vv, - QMap& satData) { + QMap& satData) { Tracer tracer("t_pppFilter::outlierDetection"); @@ -755,21 +835,37 @@ QString prnGPS; QString prnGlo; + + double ionFacGPS; + double ionFacGLO; + double maxResGPS = 0.0; // GPS + Galileo double maxResGlo = 0.0; // GLONASS + BDS + findMaxRes(vv, satData, prnGPS, prnGlo, maxResGPS, maxResGlo); + ionFacGLO = ionFac(prnGlo,satData); + if (iPhase == 0) + ionFacGLO *= (prnGlo[0]=='R'? GLONASS_WEIGHT_FACTOR : BDS_WEIGHT_FACTOR); + ionFacGPS = ionFac(prnGPS,satData); + if (iPhase == 1) { - if (maxResGlo > 2.98 * OPT->_maxResL1) { + if (maxResGlo > ionFacGLO * OPT->_maxResL1) { LOG << "Outlier Phase " << prnGlo.mid(0,3).toLatin1().data() << ' ' << maxResGlo << endl; return prnGlo; } - else if (maxResGPS > MAXRES_PHASE_GPS) { + else if (maxResGPS > ionFacGPS * OPT->_maxResL1) { LOG << "Outlier Phase " << prnGPS.mid(0,3).toLatin1().data() << ' ' << maxResGPS << endl; return prnGPS; } } - else if (iPhase == 0 && maxResGPS > 2.98 * OPT->_maxResC1) { - LOG << "Outlier Code " << prnGPS.mid(0,3).toLatin1().data() << ' ' << maxResGPS << endl; - return prnGPS; + else if (iPhase == 0) { + if (maxResGlo > ionFacGLO * OPT->_maxResC1) { + LOG << "Outlier Code " << prnGlo.mid(0,3).toLatin1().data() << ' ' << maxResGlo << endl; + return prnGlo; + } + else if (maxResGPS > ionFacGPS * OPT->_maxResC1) { + LOG << "Outlier Code " << prnGPS.mid(0,3).toLatin1().data() << ' ' << maxResGPS << endl; + return prnGPS; + } } @@ -780,5 +876,5 @@ /////////////////////////////////////////////////////////////////////////// double t_pppFilter::windUp(const QString& prn, const ColumnVector& rSat, - const ColumnVector& rRec) { + const ColumnVector& rRec) { Tracer tracer("t_pppFilter::windUp"); @@ -919,10 +1015,13 @@ satData->obsIndex = iObs; + // Iono-free combination noise factor + // ---------------------------------- + double ionFac = sqrt(pow(satData->lkA,2) + pow(satData->lkB,2)); + // Phase Observations // ------------------ if (iPhase == 1) { - ll(iObs) = satData->L3 - cmpValue(satData, true); - double sigL3 = 2.98 * OPT->_sigmaL1; + double sigL3 = ionFac * ellWgtCoef * OPT->_sigmaL1; if (satData->system() == 'R') { sigL3 *= GLONASS_WEIGHT_FACTOR; @@ -931,5 +1030,7 @@ sigL3 *= BDS_WEIGHT_FACTOR; } - PP(iObs,iObs) = 1.0 / (sigL3 * sigL3) / (ellWgtCoef * ellWgtCoef); + satData->L3sig = sigL3; + ll(iObs) = satData->L3 - cmpValue(satData, true); + PP(iObs,iObs) = 1.0 / (sigL3 * sigL3); for (int iPar = 1; iPar <= _params.size(); iPar++) { if (_params[iPar-1]->type == t_pppParam::AMB_L3 && @@ -944,7 +1045,14 @@ // ----------------- else { - double sigP3 = 2.98 * OPT->_sigmaC1; + double sigP3 = ionFac * ellWgtCoef * OPT->_sigmaC1; + if (satData->system() == 'R') { + sigP3 *= GLONASS_WEIGHT_FACTOR; + } + if (satData->system() == 'C') { + sigP3 *= BDS_WEIGHT_FACTOR; + } + satData->P3sig = sigP3; ll(iObs) = satData->P3 - cmpValue(satData, false); - PP(iObs,iObs) = 1.0 / (sigP3 * sigP3) / (ellWgtCoef * ellWgtCoef); + PP(iObs,iObs) = 1.0 / (sigP3 * sigP3); for (int iPar = 1; iPar <= _params.size(); iPar++) { AA(iObs, iPar) = _params[iPar-1]->partial(satData, false); @@ -973,5 +1081,10 @@ << " RES " << satData->prn.mid(0,3).toLatin1().data() << (iPhase ? " L3 " : " P3 ") - << setw(9) << setprecision(4) << vv(satData->obsIndex) << endl; + << setw(9) << setprecision(4) << vv(satData->obsIndex) << " " + << setprecision(3) + << setw(7) << (iPhase? satData->L3sig : satData->P3sig) << " " + << setprecision(1) + << setw(5) <eleSat * 180 / M_PI + << endl; } } @@ -1166,5 +1279,5 @@ } -// Remeber Original State Vector and Variance-Covariance Matrix +// Remember Original State Vector and Variance-Covariance Matrix //////////////////////////////////////////////////////////////////////////// void t_pppFilter::rememberState(t_epoData* epoData) { Index: /trunk/BNC/src/PPP_SSR_I/pppFilter.h =================================================================== --- /trunk/BNC/src/PPP_SSR_I/pppFilter.h (revision 9480) +++ /trunk/BNC/src/PPP_SSR_I/pppFilter.h (revision 9481) @@ -52,11 +52,15 @@ P2 = 0.0; P5 = 0.0; + P6 = 0.0; P7 = 0.0; P3 = 0.0; + P3sig = 0.0; L1 = 0.0; L2 = 0.0; L5 = 0.0; + L6 = 0.0; L7 = 0.0; L3 = 0.0; + L3sig = 0.0; lkA = 0.0; lkB = 0.0; @@ -74,11 +78,15 @@ double P2; double P5; + double P6; double P7; double P3; + double P3sig; double L1; double L2; double L5; + double L6; double L7; double L3; + double L3sig; ColumnVector xx; ColumnVector vv; Index: /trunk/BNC/src/bncantex.cpp =================================================================== --- /trunk/BNC/src/bncantex.cpp (revision 9480) +++ /trunk/BNC/src/bncantex.cpp (revision 9481) @@ -321,4 +321,19 @@ frqType = t_frequency::R1; } + else if (prn[0] == 'E') { + frqType = t_frequency::E1; + } + else if (prn[0] == 'C') { + frqType = t_frequency::C2; + } + else if (prn[0] == 'S') { + frqType = t_frequency::S1; + } + else if (prn[0] == 'J') { + frqType = t_frequency::J1; + } + else if (prn[0] == 'I') { + frqType = t_frequency::I5; + } QMap::const_iterator it = _maps.find(prn.mid(0,3)); @@ -351,4 +366,45 @@ return failure; +} + +// +//////////////////////////////////////////////////////////////////////////// +double bncAntex::satCorr(const QString& prn, t_frequency::type frqType, + double elTx, double azTx, bool& found) const { + + if (_maps.find(prn.mid(0,3)) == _maps.end()) { + found = false; + return 0.0; + }; + + t_antMap* map = _maps[prn.mid(0,3)]; + + if (map->frqMap.find(frqType) == map->frqMap.end()) { + found = false; + return 0.0; + }; + + t_frqMap* frqMap = map->frqMap[frqType]; + + double var = 0.0; + if (frqMap->pattern.ncols() > 0) { + double zenDiff = 999.999; + double zenTx = 90.0 - elTx * 180.0 / M_PI; + unsigned iZen = 0; + for (double zen = map->zen1; zen <= map->zen2; zen += map->dZen) { + iZen += 1; + double newZenDiff = fabs(zen - zenTx); + if (newZenDiff < zenDiff) { + zenDiff = newZenDiff; + var = frqMap->pattern(iZen); + } + } + } + + found = true; + return var - frqMap->neu[0] * cos(azTx)*cos(elTx) + - frqMap->neu[1] * sin(azTx)*cos(elTx) + - frqMap->neu[2] * sin(elTx); + } Index: /trunk/BNC/src/bncantex.h =================================================================== --- /trunk/BNC/src/bncantex.h (revision 9480) +++ /trunk/BNC/src/bncantex.h (revision 9481) @@ -40,4 +40,6 @@ void print() const; QString pcoSinexString(const std::string& antName, t_frequency::type frqType); + double satCorr(const QString& prn, t_frequency::type frqType, + double eleSat, double azSat, bool& found) const; double rcvCorr(const std::string& antName, t_frequency::type frqType, double eleSat, double azSat, bool& found) const; Index: /trunk/BNC/src/pppMain.cpp =================================================================== --- /trunk/BNC/src/pppMain.cpp (revision 9480) +++ /trunk/BNC/src/pppMain.cpp (revision 9481) @@ -164,4 +164,5 @@ if (_realTime) { opt->_corrMount.assign(settings.value("PPP/corrMount").toString().toStdString()); + opt->_isAPC = (opt->_corrMount.substr(0,4)=="SSRA"); } else { @@ -169,4 +170,6 @@ opt->_rinexNav.assign(settings.value("PPP/rinexNav").toString().toStdString()); opt->_corrFile.assign(settings.value("PPP/corrFile").toString().toStdString()); + QFileInfo tmp = QFileInfo(QString::fromStdString(opt->_corrFile)); + opt->_isAPC = (tmp.baseName().mid(0,4)=="SSRA"); } Index: /trunk/BNC/src/pppOptions.h =================================================================== --- /trunk/BNC/src/pppOptions.h (revision 9480) +++ /trunk/BNC/src/pppOptions.h (revision 9481) @@ -37,4 +37,5 @@ std::string _crdFile; std::string _corrMount; + bool _isAPC; std::string _rinexObs; std::string _rinexNav;