Index: trunk/BNC/src/PPP/ambres.cpp =================================================================== --- trunk/BNC/src/PPP/ambres.cpp (revision 10791) +++ trunk/BNC/src/PPP/ambres.cpp (revision 10791) @@ -0,0 +1,260 @@ +#include +#include +#include +#include + +#include "ambres.h" +#include "bncutils.h" +#include "lambda.h" +#include "pppClient.h" + +using namespace std; +using namespace BNC_PPP; + +// +////////////////////////////////////////////////////////////////////////////////////////// +AmbRes::AmbRes() { + reset(); +} + +// +////////////////////////////////////////////////////////////////////////////////////////// +AmbRes::~AmbRes() { +} + +// +////////////////////////////////////////////////////////////////////////////////////////// +bool AmbRes::isResolvable(const t_pppParam* amb) { + if ( amb->ambNumEpo() >= OPT->_ar._minNumEpo && amb->ambEleSat() >= OPT->_ar._minEle) { + return true; + } + else { + return false; + } +} + +// +////////////////////////////////////////////////////////////////////////////////////////// +bool AmbRes::isFixable(double xx, double sigma, bool* checked) { + if (OPT->_ar._maxFrac <= 0.0 || OPT->_ar._maxSig <= 0.0) { + if (checked) { + *checked = false; + } + return true; + } + else { + if (checked) { + *checked = true; + } + return ( fabs(xx - nint(xx)) <= OPT->_ar._maxFrac && sigma <= OPT->_ar._maxSig ); + } +} + +// +////////////////////////////////////////////////////////////////////////////////////////// +void AmbRes::addToGroup(const t_pppParam* amb) { + t_lc LC = amb->LC(); + char sys = amb->prn().system(); + AmbGroup* p_ambGroup = 0; + for (AmbGroup& ambGroup : _ambGroups) { + if (ambGroup._system == sys && ambGroup._LC == LC) { + p_ambGroup = &ambGroup; + break; + } + } + if (!p_ambGroup) { + _ambGroups.emplace_back(AmbGroup(sys, LC)); + p_ambGroup = &_ambGroups.back(); + } + p_ambGroup->_zdAmbs.push_back(ZdAmb(amb)); + + _numZdAmbs += 1; +} + +// Select Reference Ambiguity +//////////////////////////////////////////////////////////////////////////////////////////// +void AmbRes::selectReference(AmbGroup& ambGroup) { + + // Compute sum of variances of all double-difference ambiguities + // ------------------------------------------------------------- + double minVarSum = 0.0; + for (const ZdAmb& zdAmb1 : ambGroup._zdAmbs) { + int i1 = zdAmb1.index(); + double varSum = 0.0; + for (const ZdAmb& zdAmb2 : ambGroup._zdAmbs) { + int i2 = zdAmb2.index(); + varSum += _QQ[i1][i1] - 2.0 * _QQ(i1+1,i2+1) + _QQ[i2][i2]; + } + if (ambGroup._zdAmbRef == 0 || minVarSum > varSum) { + ambGroup._zdAmbRef = &zdAmb1; + minVarSum = varSum; + } + } +} + +// +////////////////////////////////////////////////////////////////////////////////////////// +t_irc AmbRes::run(const bncTime& epoTime, + const std::vector& params, + SymmetricMatrix& QFinal, + ColumnVector& xFinal, + double& fixRatio, + std::ostream& msg) { + reset(); + fixRatio = 0.0; + + msg << string(epoTime) << " Ambiguity Resolution (BIE)" << endl; + + // Resolvable ambiguities + // ---------------------- + map> usedPrns; + for (char system : OPT->_ar._systems) { + for (const auto& par : params) { + if (par->type() == t_pppParam::amb && par->prn().system() == system) { + if (isResolvable(par)) { + addToGroup(par); + } + } + } + } + + // Remove Groups with less than 2 ambiguities + // ------------------------------------------ + auto iGrp = _ambGroups.begin(); + while (iGrp != _ambGroups.end()) { + if (iGrp->_zdAmbs.size() < 2) { + _numZdAmbs -= iGrp->_zdAmbs.size(); + iGrp = _ambGroups.erase(iGrp); + } + else { + ++iGrp; + } + } + + // Check number of zero-difference ambiguities + // ------------------------------------------- + bool hasEnoughAmbs = false; + for (const auto& ambGroup : _ambGroups) { + if (int(ambGroup._zdAmbs.size()) >= OPT->_ar._minNumSat) { + hasEnoughAmbs = true; + break; + } + } + if (!hasEnoughAmbs) { + msg << " not enough resolvable ambiguities" << endl; + return t_irc::failure; + } + + // Construct the variance-covariance matrix of ambiguities + // ------------------------------------------------------- + Matrix ZD(_numZdAmbs, params.size()); ZD = 0.0; + int idx = -1; + for (auto& ambGroup : _ambGroups) { + for (auto& zdAmb : ambGroup._zdAmbs) { + idx += 1; + zdAmb.setIndex(idx); + ZD[idx][zdAmb.indexGlobal()] = 1.0; + } + } + _xx = ZD * xFinal; + _QQ << ZD * QFinal * ZD.t(); + + // Select and Constrain Reference Ambiguities + // ------------------------------------------ + Matrix HH(_ambGroups.size(), _numZdAmbs); HH = 0.0; + ColumnVector hh(_ambGroups.size()); + for (unsigned iGrp = 0; iGrp < _ambGroups.size(); ++iGrp) { + AmbGroup& ambGroup = _ambGroups[iGrp]; + selectReference(ambGroup); + HH[iGrp][ambGroup._zdAmbRef->index()] = 1.0; + hh[iGrp] = nint(_xx[ambGroup._zdAmbRef->index()]); + } + DiagonalMatrix PP(HH.nrows()); PP = 1.0 / (sigCon * sigCon); + kalman(HH, hh, PP, _QQ, _xx); + + // Perform the Lambda (BIE) Search + // ------------------------------- + ColumnVector xBie; + SymmetricMatrix covBie; + Lambda::search(_xx, _QQ, xBie, covBie); + if (xBie.Nrows() == 0) { + msg << " BIE search: failed" << endl; + return t_irc::failure; + } + + // Print BIE Results + // ----------------- + printBieResults(xBie, covBie, msg); + + // Constrain ambiguities + // --------------------- + setConstraints(QFinal, xFinal, xBie, covBie, fixRatio, msg); + + return t_irc::success; +} + +// Print single-difference ambiguities and AR results +//////////////////////////////////////////////////////////////////////////////////////////// +void AmbRes::printBieResults(const ColumnVector& xBie, const SymmetricMatrix& covBie, + ostream& msg) const { + msg.setf(ios::fixed); + for (const AmbGroup& ambGroup : _ambGroups) { + msg << " Group " << ambGroup._LC.toString() << endl; + + for (const ZdAmb& zdAmb: ambGroup._zdAmbs) { + int idx = zdAmb.index(); + msg << " " << zdAmb.prn().toString() << ' ' << setw(7) << setprecision(2) << _xx[idx]; + if (&zdAmb == ambGroup._zdAmbRef) { + msg << " ref"; + } + else { + msg << " +- "; + } + msg << setw(6) << setprecision(2) << sqrtMod(_QQ[idx][idx]) << ' ' + << setw(7) << setprecision(2) << xBie[idx] << ' ' + << setw(9) << setprecision(6) << sqrtMod(covBie[idx][idx]); + if (isFixable(xBie[idx], sqrtMod(covBie[idx][idx]))) { + msg << " fix "; + } + msg << endl; + } + } +} + +// +//////////////////////////////////////////////////////////////////////////////////////////// +void AmbRes::setConstraints(SymmetricMatrix& QFinal, ColumnVector& xFinal, + const ColumnVector& xBie, const SymmetricMatrix& covBie, + double& fixRatio, ostream& /*msg*/) { + + fixRatio = 0.0; + + vector> CC; + vector cc; + vector Sc; + int numSdAmbs = 0; + int numFixSdAll = 0; + for (const AmbGroup& ambGroup : _ambGroups) { + numSdAmbs += ambGroup._zdAmbs.size() - 1; + for (const ZdAmb& zdAmb : ambGroup._zdAmbs) { + int idx = zdAmb.index(); + if (isFixable(xBie[idx], sqrtMod(covBie[idx][idx]))) { + if (ambGroup._zdAmbRef->indexGlobal() != zdAmb.indexGlobal()) { + numFixSdAll += 1; + } + RowVector* pCC; + CC.emplace_back(pCC = new RowVector(xFinal.Nrows())); + cc.push_back(xBie[idx]); + Sc.push_back(sigCon); + (*pCC) = 0.0; + (*pCC)[zdAmb.indexGlobal()] = 1.0; + } + } + } + if (CC.size() > 0) { + kalman(CC, cc, Sc, QFinal, xFinal); + } + if (numSdAmbs > 0) { + fixRatio = double(numFixSdAll) / double(numSdAmbs); + } +} Index: trunk/BNC/src/PPP/ambres.h =================================================================== --- trunk/BNC/src/PPP/ambres.h (revision 10791) +++ trunk/BNC/src/PPP/ambres.h (revision 10791) @@ -0,0 +1,95 @@ +#ifndef AMBRES_H +#define AMBRES_H + +#include +#include +#include +#include "pppInclude.h" +#include "pppParlist.h" + +namespace BNC_PPP { + + class AmbRes { + public: + AmbRes(); + ~AmbRes(); + + t_irc run(const bncTime& epoTime, + const std::vector& params, + SymmetricMatrix& QFinal, + ColumnVector& xFinal, + double& fixRatio, + std::ostream& msg); + + private: + static constexpr double sigCon = 1e-4; + + class ZdAmb { + public: + ZdAmb(const t_pppParam* amb) : _amb(amb), _index(-1) {} + t_prn prn() const {return _amb->prn();} + t_lc LC() const {return _amb->LC();} + double valueApr() const {return _amb->x0();} + int index() const {return _index;} + int indexGlobal() const {return _amb->indexNew();} + void setIndex(int index) {_index = index;} + bncTime firstObsTime() const {return _amb->firstObsTime();} + const t_pppParam* ambGlobal() const {return _amb;} + private: + const t_pppParam* _amb; + int _index; + }; + + class SdAmb { + public: + SdAmb(const ZdAmb& zdAmb1, const ZdAmb& zdAmb2, + const ColumnVector& xBie, const SymmetricMatrix& covBie); + t_lc _LC; + t_prn _prn1; + t_prn _prn2; + double _aprVal; + double _bieValue; + double _bieSigma; + bool _fixed; + bncTime _firstObsTime; + const t_pppParam* _ambGlobal[2]; + }; + + class AmbGroup { + public: + AmbGroup(char sys, const t_lc& LC) : _system(sys), _LC(LC), _zdAmbRef(0) {}; + char _system; + t_lc _LC; + const ZdAmb* _zdAmbRef; + std::vector _zdAmbs; + std::vector _sdAmbs; + }; + + static bool isResolvable(const t_pppParam* amb); + static bool isFixable(double xx, double sigma, bool* checked = 0); + + void reset() { + _numZdAmbs = 0; + _ambGroups.clear(); + _QQ.cleanup(); + _xx.cleanup(); + } + + void addToGroup(const t_pppParam* amb); + void selectReference(AmbGroup& ambGroup); + void printBieResults(const ColumnVector& xBie, const SymmetricMatrix& covBie, + std::ostream& msg) const; + void setConstraints(SymmetricMatrix& QFinal, ColumnVector& xFinal, + const ColumnVector& xBie, const SymmetricMatrix& covBie, + double& fixRatio, std::ostream& msg); + + unsigned _numZdAmbs; + std::vector _ambGroups; + SymmetricMatrix _QQ; + ColumnVector _xx; + + }; + +} + +#endif Index: trunk/BNC/src/PPP/lambda.cpp =================================================================== --- trunk/BNC/src/PPP/lambda.cpp (revision 10791) +++ trunk/BNC/src/PPP/lambda.cpp (revision 10791) @@ -0,0 +1,430 @@ + +#include +#include +#include +#include + +#include "lambda.h" + +using namespace BNC_PPP; +using namespace std; + +// Gauss Error Function +///////////////////////////////////////////////////////////////////////////////////////// +double Lambda::erf(double xx) { + static const double a1 = 0.254829592; + static const double a2 = -0.284496736; + static const double a3 = 1.421413741; + static const double a4 = -1.453152027; + static const double a5 = 1.061405429; + static const double pp = 0.3275911; + + int sign = (xx < 0) ? -1 : 1; + xx = fabs(xx); + + double tt = 1.0/(1.0 + pp*xx); + double yy = 1.0 - (((((a5*tt + a4)*tt) + a3)*tt + a2)*tt + a1)*tt*exp(-xx*xx); + + return sign * yy; +} + +// Cumulative Distribution Function (Normal Distribution) +///////////////////////////////////////////////////////////////////////////////////////// +double Lambda::normcdf(double x, double mu, double sigma) { + static const double root_two = sqrt(2.0); + return 0.5 * ( 1.0 + erf( (x-mu) / (sigma*root_two) ) ); +} + +// Auxiliary functions +///////////////////////////////////////////////////////////////////////////////////////// +void Lambda::swap(double& a, double& b) { + double t(a); a = b; b = t; +} +double Lambda::sign(double a) { + if (a < 0.0) { + return -1.0; + } + else if (a > 0.0) { + return 1.0; + } + else { + return 0.0; + } +} +double Lambda::nint(double val) { + return ((val < 0.0) ? -floor(fabs(val)+0.5) : floor(val+0.5)); +} + +// LAMBDA/BIE Search +///////////////////////////////////////////////////////////////////////////////////////// +void Lambda::search(ColumnVector aFlt, const SymmetricMatrix& QQ, + ColumnVector& aFix, SymmetricMatrix& covBie) { + + int nn = QQ.Nrows(); + + // Remove integer numbers from float solution (for computational convenience only) + // ------------------------------------------------------------------------------- + ColumnVector incr(nn); + for (int ii = 0; ii < nn; ii++) { + incr[ii] = nint(aFlt[ii]); + aFlt[ii] = aFlt[ii] - incr[ii]; + } + + // Compute ZZ matrix based on the decomposition Q=LL^T*D*LL; The transformed + // float solution: zFlt = ZZ^T *aFlt, QzFlt = ZZ^T * QQ * ZZ + // ----------------------------------------------------------------------- + SymmetricMatrix QzFlt; + Matrix ZZ; + LowerTriangularMatrix LL; + DiagonalMatrix DD; + ColumnVector zFlt; + Matrix iZt; + decorrel(QQ, aFlt, QzFlt, ZZ, LL, DD, zFlt, iZt); + + // Perform the search + // ------------------ + Info info; + BIE(zFlt, LL, DD, info); + + // Perform the back-transformation and add the increments + // ------------------------------------------------------ + aFix = iZt * info.zBie + incr; + info.zBie = ZZ.t() * aFix; + info.aFix = iZt * info.zFix; + for (int iCand = 1; iCand <= info.zFix.Ncols(); iCand++) { + info.aFix.column(iCand) += incr; + } + info.zFix = ZZ.t() * info.aFix; + +#ifdef LAMBDA_MAIN_TEST + cout.setf(ios::fixed); + cout << "aFix(3 cand)= \n" << setw(10) << setprecision(2) << info.aFix.columns(1,3) << endl; +#endif + + // Variances of ambiguities + // ------------------------ + static const double sigCon = 1e-6; + DiagonalMatrix covZ(info.zFix.Nrows()); covZ = 0.0; + for (int ia = 0; ia < info.zFix.Nrows(); ia++) { // loop over all ambiguities + for (int ic = 0; ic < info.zFix.Ncols(); ic++) { // loop over all candidates + double dZ = info.zBie[ia] - info.zFix[ia][ic]; + covZ[ia] += info.wgt[ic] * dZ * dZ; + } + if (covZ[ia] < sigCon*sigCon) { // make the matrix positive-definite + covZ[ia] = sigCon*sigCon; + } + } + Matrix invZ = ZZ.i(); + covBie << invZ.t() * covZ * invZ; +} + +// Best Integer Equivariant Estimator +///////////////////////////////////////////////////////////////////////////////////////// +void Lambda::BIE(const ColumnVector& zFlt, // Original ambiguities + const LowerTriangularMatrix& LL, // L matrix from L'DL-decomposition of QzFlt + const DiagonalMatrix& DD, // D matrix from L'DL-decomposition of QzFlt + Info& info) { + + int nn = zFlt.Nrows(); + int ncands = 100; + ColumnVector sqnorm; + info.wgt.ReSize(ncands); info.wgt = 0.0; + info.zBie.ReSize(nn); info.zBie = 0.0; + + ssearch(zFlt, LL, DD, ncands, info.zFix, sqnorm); + + LowerTriangularMatrix Li = LL.i(); + SymmetricMatrix QzFltInv; QzFltInv << Li * DD.i() * Li.t(); + + info.zBie.ReSize(nn); info.zBie = 0.0; + double wgtSum = 0.0; + double norm1 = 0.0; + for (int ic = 1; ic <= ncands; ic++) { + + ColumnVector da = zFlt - info.zFix.column(ic); + double daNorm = DotProduct(da, QzFltInv * da); + + if (ic == 1) { + norm1 = daNorm; + info.wgt(ic) = 1.0; + } + else { + info.wgt(ic) = exp(-0.5 * (daNorm-norm1)); // weights scaled by exp(0.5 * norm1) + } + + wgtSum += info.wgt(ic); + + for (int ia = 1; ia <= info.zFix.Nrows(); ia++) { + info.zBie(ia) += info.wgt(ic) * info.zFix(ia,ic); + } + } + + info.zBie /= wgtSum; + info.wgt /= wgtSum; +} + +// Decomposition Q = L'DL (L is lower triangular) +///////////////////////////////////////////////////////////////////////////////////////// +void Lambda::ldldecom(const SymmetricMatrix& QQ, LowerTriangularMatrix& LL, DiagonalMatrix& DD) { + + const int n = QQ.Nrows(); + + Matrix QC = QQ; + + LL.ReSize(n); LL = 0.0; + DD.resize(n); DD = 0.0; + + for (int i = n-1; i >= 0; i--) { + DD[i] = QC[i][i]; + if ( DD[i] <= 0.0 ) { + throw "ldldecom problem"; + } + double temp = sqrt(DD[i]); + for (int j = 0; j <= i; j++) { + LL[i][j] = QC[i][j]/temp; + } + for (int j = 0; j <= i-1; j++) { + for(int k = 0; k <= j; k++) { + QC[j][k] -= LL[i][k] * LL[i][j]; + } + } + for (int j = 0; j <= i; j++) { + LL[i][j] /= LL[i][i]; + } + } +} + +// +///////////////////////////////////////////////////////////////////////////////////////// +void Lambda::decorrel(const SymmetricMatrix& QQ, // Variance-covariance matrix of ambiguities + const ColumnVector& aFlt, // Original ambiguities + SymmetricMatrix& QzFlt, // Cov. matrix of decorrelated ambiguities + Matrix& ZZ, // ZZ-transformation matrix + LowerTriangularMatrix& LL, // L matrix from L'DL-decomposition of QzFlt + DiagonalMatrix& DD, // D matrix from L'DL-decomposition of QzFlt + ColumnVector& zFlt, // Transformed ambiguities + Matrix& iZt) { // ZZ.t().i() transformation matrix + + // L'DL Decomposition + // ------------------ + ldldecom(QQ, LL, DD); + + // Reduction + // --------- + int n = DD.Nrows(); + + iZt.ReSize(n,n); iZt = 0.0; + for (int i = 0; i < n; i++) { + iZt[i][i] = 1.0; + } + + int i1 = n - 1; + bool sw = true; + while (sw) { + + int i = n; // loop for column from n to 1 + sw = false; + + while ( !sw && i > 1) { + + i = i - 1; // the ith column + if (i <= i1) { + for (int j = i+1; j <= n; j++) { + double mu = nint(LL(j,i)); + if (mu != 0.0) { + for (int k = j; k <= n; k++) { + LL(k,i) = LL(k,i) - mu * LL(k,j); + } + for (int k = 1; k <= n; k++) { + iZt(k,j) = iZt(k,j) + mu * iZt(k,i); + } + } + } + } + + double delta = DD(i) + LL(i+1,i) * LL(i+1,i) * DD(i+1); + if (delta < DD(i+1)) { + double lambda = DD(i+1) * LL(i+1,i) / delta; + double eta = DD(i) / delta; + DD(i) = eta * DD(i+1); + DD(i+1) = delta; + + Matrix hlp(2,2); hlp << -LL(i+1,i) << 1.0 + << eta << lambda; + + LL.submatrix(i,i+1,1,i-1) = hlp * LL.submatrix(i,i+1,1,i-1); + + LL(i+1,i) = lambda; + + for (int k = i+2; k <= n; k++) swap( LL(k,i), LL(k,i+1)); + for (int k = 1; k <= n; k++) swap(iZt(k,i), iZt(k,i+1)); + i1 = i; + sw = true; + } + } + } + + // Transformed Q-matrix, transformation-matrix, and decorrelated ambiguities + // ------------------------------------------------------------------------- + ZZ = iZt.i().t(); + for (int i = 0; i < ZZ.Nrows(); i++) { + for (int j = 0; j < ZZ.Nrows(); j++) { + ZZ[i][j] = nint(ZZ[i][j]); + } + } + + QzFlt << ZZ.t() * QQ * ZZ; // it is also L'DL + zFlt = ZZ.t() * aFlt; +} + +// Integer ambiguity vector search by employing the search-and-shrink technique +///////////////////////////////////////////////////////////////////////////////////////// +void Lambda::ssearch(const ColumnVector& zFlt, // Original ambiguities + const LowerTriangularMatrix& LL, // L matrix from L'DL-decomposition of QzFlt + const DiagonalMatrix& DD, // D matrix from L'DL-decomposition of QzFlt + int ncands, // Number of requested candidates + Matrix& zFix, // estimated integers (n x ncands ) + ColumnVector& sqnorm) { // squared norms (ascendantly sorted) + + // Initialize outputs + // ------------------ + int n = zFlt.Nrows(); + zFix.ReSize(n, ncands); zFix = 0.0; + sqnorm.ReSize(ncands); sqnorm = 0.0; + + // Initializing the variables for searching + // ---------------------------------------- + double Chi2 = 1.0e+18; // start search with an infinite chi^2 + ColumnVector dist(n); dist(n) = 0.0; // dist(k)=sum_{j=k+1}^{n}(a_j-acond_j)^2/d_j + bool endsearch = false; + int count = 0; // the number of candidates + + ColumnVector acond(n); acond(n) = zFlt(n); + ColumnVector zcond(n); zcond(n) = nint(acond(n)); + double left = acond(n) - zcond(n); + ColumnVector step(n); step(n) = sign(left); + + // For a very occasional case when the value of float solution zFlt(n) == 0, we + // compusively give a positive step to continue. + if (step(n) == 0.0) { + step(n) = 1; + } + + int imax = ncands; // initially, the maximum F(z) is at ncands + Matrix SS(n, n); SS = 0.0; // used to compute conditional ambiguities + + int k = n; + + // Start the main search-loop + // -------------------------- + while (!endsearch) { + double newdist = dist(k) + left*left / DD(k); + if (newdist < Chi2) { + if (k != 1) { // Case 1: move down + k = k - 1; + dist(k) = newdist; + for (int j = 1; j <= k; j++) { + SS(k,j) = SS(k+1,j) + (zcond(k+1)-acond(k+1)) * LL(k+1,j); + } + + acond(k) = zFlt(k) + SS(k, k); + zcond(k) = round(acond(k)); + left = acond(k) - zcond(k); + step(k) = sign(left); + + if (step(k) == 0) { + step(k) = 1.0; + } + } + else { // Case 2: store the found candidate and try next valid integer + if (count < ncands - 1) { + count = count + 1; + zFix.column(count) = zcond; + sqnorm(count) = newdist; + } + else { + zFix.column(imax) = zcond; + sqnorm(imax) = newdist; + Chi2 = sqnorm.maximum1(imax); + } + zcond(1) = zcond(1) + step(1); // next valid integer + left = acond(1) - zcond(1); + step(1) = -step(1) - sign(step(1)); + } + } + else { // Case 3: exit or move up + if (k == n) { + endsearch = true; + } + else { + k = k + 1; // move up + zcond(k) = zcond(k) + step(k); // next valid integer + left = acond(k) - zcond(k); + step(k) = -step(k) - sign(step(k)); + } + } + } + + // Sort + // ---- + for (int i = 0; i < ncands-1; i++) { + for (int j = i+1; j < ncands; j++) { + if (sqnorm[i] > sqnorm[j]) { + swap(sqnorm[i],sqnorm[j]); + for (k = 0; k < n; k++) { + swap(zFix[k][i], zFix[k][j]); + } + } + } + } +} + +#ifdef LAMBDA_MAIN_TEST +// Main test program +// compile: g++ -DLAMBDA_MAIN_TEST -I../../newmat ../../newmat/*.cpp arLambda.cpp +///////////////////////////////////////////////////////////////////////////////////////// +int main(int argc, char* argv[]) { + + const int nn = 12; + + SymmetricMatrix QQ(nn); + QQ << 1.90688560e+04 + << -1.57839723e+04 << 5.90277038e+04 + << -1.73342006e+04 << 3.81426928e+04 << 2.81775654e+04 + << 1.44119240e+04 << 5.62717388e+02 << -7.00050220e+03 << 1.56055082e+04 + << 1.00557170e+04 << -1.38300856e+04 << -1.16958674e+04 << 5.03970282e+03 << 6.82077251e+03 + << -1.42592953e+04 << 2.73734263e+04 << 2.18861681e+04 << -9.64896531e+03 << -6.88024051e+03 << 2.32465490e+04 + << 1.48588484e+04 << -1.22991994e+04 << -1.35071695e+04 << 1.12300704e+04 << 7.83562345e+03 << -1.11111394e+04 << 1.15783237e+04 + << -1.22991994e+04 << 4.59956130e+04 << 2.97215786e+04 << 4.38480888e+02 << -1.07766903e+04 << 2.13299424e+04 << -9.58379157e+03 << 3.58407377e+04 + << -1.35071695e+04 << 2.97215786e+04 << 2.19565441e+04 << -5.45493698e+03 << -9.11366311e+03 << 1.70541567e+04 << -1.05250670e+04 << 2.31596718e+04 << 1.71089957e+04 + << 1.12300704e+04 << 4.38480887e+02 << -5.45493698e+03 << 1.21601359e+04 << 3.92704096e+03 << -7.51867446e+03 << 8.75070439e+03 << 3.41673570e+02 << -4.25060009e+03 << 9.47543087e+03 + << 7.83562345e+03 << -1.07766903e+04 << -9.11366311e+03 << 3.92704096e+03 << 5.31488728e+03 << -5.36122657e+03 << 6.10568076e+03 << -8.39742084e+03 << -7.10155552e+03 << 3.06003207e+03 << 4.14147091e+03 + << -1.11111394e+04 << 2.13299424e+04 << 1.70541567e+04 << -7.51867446e+03 << -5.36122657e+03 << 1.81141936e+04 << -8.65803054e+03 << 1.66207345e+04 << 1.32889535e+04 << -5.85870722e+03 << -4.17757899e+03 << 1.41149564e+04; + + ColumnVector aa(nn); + aa << -2.84908567e+04 + << 6.57526299e+04 + << 3.88303667e+04 + << 5.00370834e+03 + << -2.91960699e+04 + << -2.97658932e+02 + << -2.22010284e+04 + << 5.12358375e+04 + << 3.02577810e+04 + << 3.89940332e+03 + << -2.27491854e+04 + << -1.59278780e+02; + + ColumnVector aBie; + SymmetricMatrix covBie; + + Lambda::search(aa, QQ, aBie, covBie); + + cout.setf(ios::fixed); + cout << "aFlt = \n" << setw(10) << setprecision(2) << aa << endl; + cout << "aBie = \n" << setw(10) << setprecision(2) << aBie << endl; + cout << "covBie = \n" << setw( 7) << setprecision(2) << covBie << endl; + + return 0; +} +#endif Index: trunk/BNC/src/PPP/lambda.h =================================================================== --- trunk/BNC/src/PPP/lambda.h (revision 10791) +++ trunk/BNC/src/PPP/lambda.h (revision 10791) @@ -0,0 +1,47 @@ +#ifndef LAMBDA_H +#define LAMBDA_H + +#include + +namespace BNC_PPP { + + class Lambda { + public: + static void search(ColumnVector aFlt, const SymmetricMatrix& QaFlt, + ColumnVector& aFix, SymmetricMatrix& covBie); + + private: + class Info { + public: + Matrix aFix; + Matrix zFix; + ColumnVector zBie; + ColumnVector wgt; + }; + + static double erf(double xx); + + static double normcdf(double xx, double mu = 0.0, double sigma = 1.0); + + static void ldldecom(const SymmetricMatrix& QQ, LowerTriangularMatrix& LL, DiagonalMatrix& DD); + + static void decorrel(const SymmetricMatrix& QaFlt, const ColumnVector& aFlt, + SymmetricMatrix& QzFlt, Matrix& ZZ, LowerTriangularMatrix& LL, + DiagonalMatrix& DD, ColumnVector& zFlt, Matrix& iZt); + + static void ssearch(const ColumnVector& zFlt, const LowerTriangularMatrix& LL, + const DiagonalMatrix& DD, int ncands, Matrix& zFix, ColumnVector& sqnorm); + + static inline double nint(double val); + + static inline void swap(double& a, double& b); + + static inline double sign(double a); + + static void BIE(const ColumnVector& zFlt, const LowerTriangularMatrix& LL, + const DiagonalMatrix& DD, Info& info); + }; + +} + +#endif Index: trunk/BNC/src/PPP/pppClient.cpp =================================================================== --- trunk/BNC/src/PPP/pppClient.cpp (revision 10765) +++ trunk/BNC/src/PPP/pppClient.cpp (revision 10791) @@ -117,4 +117,7 @@ ////////////////////////////////////////////////////////////////////////////// void t_pppClient::putOrbCorrections(const vector& corr) { + if (OPT->_logMode == t_pppOptions::normal && corr.size() > 0) { + LOG << "orbCorrections " << string(corr[0]->_time) << ' ' << corr.size() << endl; + } for (unsigned ii = 0; ii < corr.size(); ii++) { _ephPool->putOrbCorrection(new t_orbCorr(*corr[ii])); @@ -125,4 +128,7 @@ ////////////////////////////////////////////////////////////////////////////// void t_pppClient::putClkCorrections(const vector& corr) { + if (OPT->_logMode == t_pppOptions::normal && corr.size() > 0) { + LOG << "clkCorrections " << string(corr[0]->_time) << ' ' << corr.size() << endl; + } for (unsigned ii = 0; ii < corr.size(); ii++) { _ephPool->putClkCorrection(new t_clkCorr(*corr[ii])); @@ -133,6 +139,16 @@ ////////////////////////////////////////////////////////////////////////////// void t_pppClient::putCodeBiases(const vector& biases) { + if (OPT->_logMode == t_pppOptions::normal && biases.size() > 0) { + LOG << "codeBiases " << string(biases[0]->_time) << ' ' << biases.size() << endl; + } + set systems; for (unsigned ii = 0; ii < biases.size(); ii++) { - _obsPool->putCodeBias(new t_satCodeBias(*biases[ii])); + t_satCodeBias* newBias = new t_satCodeBias(*biases[ii]); + char sys = newBias->_prn.system(); + if (systems.find(sys) == systems.end()) { + _obsPool->clearCodeBiases(sys); + systems.insert(sys); + } + _obsPool->putCodeBias(newBias); } } @@ -141,6 +157,16 @@ ////////////////////////////////////////////////////////////////////////////// void t_pppClient::putPhaseBiases(const vector& biases) { + if (OPT->_logMode == t_pppOptions::normal && biases.size() > 0) { + LOG << "phaseBiases " << string(biases[0]->_time) << ' ' << biases.size() << endl; + } + set systems; for (unsigned ii = 0; ii < biases.size(); ii++) { - _obsPool->putPhaseBias(new t_satPhaseBias(*biases[ii])); + t_satPhaseBias* newBias = new t_satPhaseBias(*biases[ii]); + char sys = newBias->_prn.system(); + if (systems.find(sys) == systems.end()) { + _obsPool->clearPhaseBiases(sys); + systems.insert(sys); + } + _obsPool->putPhaseBias(newBias); } } @@ -211,38 +237,15 @@ } } -/* - vector::iterator it = obsVector.begin(); - while (it != obsVector.end()) { - t_pppSatObs* satObs = *it; - satObs->printObsMinusComputed(); - it++; - } -*/ + + if (OPT->_logMode == t_pppOptions::all) { + vector::iterator it = obsVector.begin(); + while (it != obsVector.end()) { + t_pppSatObs* satObs = *it; + satObs->printObsMinusComputed(); + it++; + } + } return pseudoObsIono; } - -// -////////////////////////////////////////////////////////////////////////////// -void t_pppClient::useObsWithCodeBiasesOnly(std::vector& obsVector) { - - vector::iterator it = obsVector.begin(); - while (it != obsVector.end()) { - t_pppSatObs* pppSatObs = *it; - bool codeBiasesAvailable = false; - if (pppSatObs->getCodeBias(pppSatObs->fType1()) && - pppSatObs->getCodeBias(pppSatObs->fType2())) { - codeBiasesAvailable = true; - } - if (codeBiasesAvailable) { - ++it; - } - else { - it = obsVector.erase(it); - delete pppSatObs; - } - } - return; -} - // Compute the Bancroft position, check for blunders @@ -251,5 +254,5 @@ vector& obsVector, ColumnVector& xyzc, bool print) { - t_lc::type tLC = t_lc::dummy; + int numBancroft = obsVector.size(); @@ -259,16 +262,6 @@ for (unsigned ii = 0; ii < obsVector.size(); ii++) { const t_pppSatObs* satObs = obsVector.at(ii); - if (tLC == t_lc::dummy) { - if (satObs->isValid(t_lc::cIF)) { - tLC = t_lc::cIF; - } - else if (satObs->isValid(t_lc::c1)) { - tLC = t_lc::c1; - } - else if (satObs->isValid(t_lc::c2)) { - tLC = t_lc::c2; - } - } - if ( satObs->isValid(tLC) && + t_lc LC = satObs->rangeLC(); + if ( satObs->isValid(LC) && (!satObs->modelSet() || satObs->eleSat() >= _opt->_minEle) ) { ++iObs; @@ -276,5 +269,5 @@ BB[iObs][1] = satObs->xc()[1]; BB[iObs][2] = satObs->xc()[2]; - BB[iObs][3] = satObs->obsValue(tLC) - satObs->cmpValueForBanc(tLC); + BB[iObs][3] = satObs->obsValue(LC) - satObs->cmpValueForBanc(LC); } } @@ -297,8 +290,9 @@ for (unsigned ii = 0; ii < obsVector.size(); ii++) { const t_pppSatObs* satObs = obsVector.at(ii); - if (satObs->isValid() && + t_lc LC = satObs->rangeLC(); + if (satObs->isValid(LC) && (!satObs->modelSet() || satObs->eleSat() >= _opt->_minEle) ) { ColumnVector rr = satObs->xc().Rows(1,3) - xyzc.Rows(1,3); - double res = rr.NormFrobenius() - satObs->obsValue(tLC) + double res = rr.NormFrobenius() - satObs->obsValue(LC) - (satObs->xc()[3] - xyzc[3]) * t_CST::c; if (fabs(res) > maxRes) { @@ -354,5 +348,5 @@ // ////////////////////////////////////////////////////////////////////////////// -void t_pppClient::finish(t_irc irc, int ind) { +void t_pppClient::finish(t_irc irc, int /*ind*/) { #ifdef BNC_DEBUG_PPP LOG << "t_pppClient::finish(" << ind << "): " << irc << endl; @@ -463,9 +457,10 @@ } } - // use observations only if satellite code biases are available - /* ------------------------------------------------------------ - if (!_opt->_corrMount.empty() { - useObsWithCodeBiasesOnly(_obsRover); - }*/ + + // Use observations only if satellite code biases are available + // ------------------------------------------------------------ + if (_opt->ambRes()) { + useObsWithBiasesOnly(_obsRover); + } if (int(_obsRover.size()) < _opt->_minObs) { @@ -613,2 +608,18 @@ } + +// +////////////////////////////////////////////////////////////////////////////// +void t_pppClient::useObsWithBiasesOnly(vector& obsVector) const { + vector::iterator it = obsVector.begin(); + while (it != obsVector.end()) { + t_pppSatObs* pppSatObs = *it; + if (pppSatObs->hasBiases()) { + ++it; + } + else { + it = obsVector.erase(it); + delete pppSatObs; + } + } +} Index: trunk/BNC/src/PPP/pppClient.h =================================================================== --- trunk/BNC/src/PPP/pppClient.h (revision 10765) +++ trunk/BNC/src/PPP/pppClient.h (revision 10791) @@ -52,5 +52,4 @@ std::vector& obsVector, bncTime& epoTime); bool preparePseudoObs(std::vector& obsVector); - void useObsWithCodeBiasesOnly(std::vector& obsVector); t_irc cmpModel(t_pppStation* station, const ColumnVector& xyzc, std::vector& obsVector); @@ -61,5 +60,5 @@ double cmpOffGal(std::vector& obsVector); double cmpOffBds(std::vector& obsVector); - + void useObsWithBiasesOnly(std::vector& obsVector) const; t_output* _output; @@ -75,5 +74,4 @@ t_tides* _tides; bool _pseudoObsIono; - QMap _usedSystems; bool _running; }; Index: trunk/BNC/src/PPP/pppEphPool.cpp =================================================================== --- trunk/BNC/src/PPP/pppEphPool.cpp (revision 10765) +++ trunk/BNC/src/PPP/pppEphPool.cpp (revision 10791) @@ -16,4 +16,5 @@ #include +#include #include "pppEphPool.h" #include "pppInclude.h" @@ -41,4 +42,7 @@ if (corr) { _satEphPool[corr->_prn.toInt()].putOrbCorrection(corr); + if (OPT->_logMode > t_pppOptions::normal) { + LOG << "orbCorr " << string(corr->_time) << ' ' << corr->_prn.toString() << endl; + } } } @@ -49,4 +53,9 @@ if (corr) { _satEphPool[corr->_prn.toInt()].putClkCorrection(corr); + if (OPT->_logMode > t_pppOptions::normal) { + LOG.setf(ios::fixed); + LOG << "clkCorr " << string(corr->_time) << ' ' << corr->_prn.toString() << ' ' + << setw(7) << setprecision(3) << corr->_dClk * t_CST::c << endl; + } } } Index: trunk/BNC/src/PPP/pppFilter.cpp =================================================================== --- trunk/BNC/src/PPP/pppFilter.cpp (revision 10765) +++ trunk/BNC/src/PPP/pppFilter.cpp (revision 10791) @@ -16,4 +16,5 @@ #include +#include #include #include @@ -28,4 +29,5 @@ #include "pppStation.h" #include "pppClient.h" +#include "ambres.h" using namespace BNC_PPP; @@ -107,9 +109,29 @@ } - // close epoch processing + // Dillution of Precision // ---------------------- cmpDOP(allObs); - _parlist->printResult(_epoTime, _QFlt, _xFlt); + + // Ambiguity Resolution + // -------------------- + if (OPT->ambRes()) { + AmbRes ambRes; + ostringstream msg; + ColumnVector xFix = _xFlt; + SymmetricMatrix QFix = _QFlt; + double fixRatio = 0.0; + ambRes.run(_epoTime, _parlist->params(), QFix, xFix, fixRatio, msg); + LOG << msg.str(); + _parlist->printResult(_epoTime, QFix, xFix, fixRatio); + } + + // Float Solution + // -------------- + else { + _parlist->printResult(_epoTime, _QFlt, _xFlt); + } + _lastEpoTimeOK = _epoTime; // remember time of last successful epoch processing + return success; } @@ -117,5 +139,5 @@ // Process Selected LCs //////////////////////////////////////////////////////////////////////////// -t_irc t_pppFilter::processSystem(const vector &LCs, +t_irc t_pppFilter::processSystem(const vector& LCs, const vector &obsVector, bool pseudoObsIonoAvailable) { @@ -158,5 +180,5 @@ int iObs = -1; vector usedObs; - vector usedTypes; + vector usedTypes; // Real Observations @@ -171,18 +193,21 @@ nSat++; for (unsigned jj = 0; jj < usedLCs; jj++) { - const t_lc::type tLC = LCs[jj]; - if (tLC == t_lc::GIM) { + const t_lc LC = LCs[jj]; + if (LC._type == t_lc::GIM) { + continue; + } + if (LC._frq1 == t_frequency::G5 && !obs->isValid(LC)) { continue; } ++iObs; usedObs.push_back(obs); - usedTypes.push_back(tLC); + usedTypes.push_back(LC); for (unsigned iPar = 0; iPar < nPar; iPar++) { const t_pppParam *par = params[iPar]; - AA[iObs][iPar] = par->partial(_epoTime, obs, tLC); - } - - ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs + 1)); - PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC)); + AA[iObs][iPar] = par->partial(_epoTime, obs, LC); + } + + ll[iObs] = obs->obsValue(LC) - obs->cmpValue(LC) - DotProduct(_x0, AA.Row(iObs + 1)); + PP[iObs] = 1.0 / (obs->sigma(LC) * obs->sigma(LC)); } } @@ -202,6 +227,6 @@ if (!obs->outlier()) { for (unsigned jj = 0; jj < usedLCs; jj++) { - const t_lc::type tLC = LCs[jj]; - if (tLC == t_lc::GIM) { + const t_lc LC = LCs[jj]; + if (LC._type == t_lc::GIM) { ++iObs; } else { @@ -209,11 +234,11 @@ } usedObs.push_back(obs); - usedTypes.push_back(tLC); + usedTypes.push_back(LC); for (unsigned iPar = 0; iPar < nPar; iPar++) { const t_pppParam *par = params[iPar]; - AA[iObs][iPar] = par->partial(_epoTime, obs, tLC); + AA[iObs][iPar] = par->partial(_epoTime, obs, LC); } - ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs + 1)); - PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC)); + ll[iObs] = obs->obsValue(LC) - obs->cmpValue(LC) - DotProduct(_x0, AA.Row(iObs + 1)); + PP[iObs] = 1.0 / (obs->sigma(LC) * obs->sigma(LC)); } } @@ -236,13 +261,13 @@ double maxOutlier = 0.0; int maxOutlierIndex = -1; - t_lc::type maxOutlierLC = t_lc::dummy; + t_lc maxOutlierLC; for (unsigned ii = 0; ii < usedObs.size(); ii++) { - const t_lc::type tLC = usedTypes[ii]; + const t_lc LC = usedTypes[ii]; double res = fabs(vv[ii]); - if (res > usedObs[ii]->maxRes(tLC)) { + if (res > usedObs[ii]->maxRes(LC)) { if (res > fabs(maxOutlier)) { maxOutlier = vv[ii]; maxOutlierIndex = ii; - maxOutlierLC = tLC; + maxOutlierLC = LC; } } @@ -254,12 +279,12 @@ t_pppSatObs *obs = usedObs[maxOutlierIndex]; t_pppParam *par = 0; - LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' ' + LOG << epoTimeStr << " Outlier " << maxOutlierLC.toString() << ' ' << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << maxOutlier << endl; for (unsigned iPar = 0; iPar < nPar; iPar++) { - t_pppParam *hlp = params[iPar]; + t_pppParam* hlp = params[iPar]; if (hlp->type() == t_pppParam::amb && hlp->prn() == obs->prn() && - hlp->tLC() == usedTypes[maxOutlierIndex]) { + hlp->LC() == usedTypes[maxOutlierIndex]) { par = hlp; } @@ -277,10 +302,10 @@ for (unsigned jj = 0; jj < LCs.size(); jj++) { for (unsigned ii = 0; ii < usedObs.size(); ii++) { - const t_lc::type tLC = usedTypes[ii]; - t_pppSatObs *obs = usedObs[ii]; - if (tLC == LCs[jj]) { - obs->setRes(tLC, vv[ii]); + const t_lc LC = usedTypes[ii]; + t_pppSatObs* obs = usedObs[ii]; + if (LC == LCs[jj]) { + obs->setRes(LC, vv[ii]); LOG << epoTimeStr << " RES " << left << setw(3) - << t_lc::toString(tLC) << right << ' ' + << LC.toString() << right << ' ' << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv[ii] << endl; @@ -296,5 +321,5 @@ // Cycle-Slip Detection //////////////////////////////////////////////////////////////////////////// -t_irc t_pppFilter::detectCycleSlips(const vector &LCs, +t_irc t_pppFilter::detectCycleSlips(const vector& LCs, const vector &obsVector) { @@ -309,9 +334,8 @@ SLIP *= fac; string epoTimeStr = string(_epoTime); - const vector ¶ms = _parlist->params(); for (unsigned ii = 0; ii < LCs.size(); ii++) { - const t_lc::type &tLC = LCs[ii]; - if (t_lc::includesPhase(tLC)) { + const t_lc& LC = LCs[ii]; + if (LC.includesPhase()) { for (unsigned iObs = 0; iObs < obsVector.size(); iObs++) { const t_pppSatObs *obs = obsVector[iObs]; @@ -343,24 +367,6 @@ // -------- if (slip) { - resetAmb(obs->prn(), obsVector, tLC); - } - - // Check Pre-Fit Residuals - /* ----------------------- - else { - ColumnVector AA(params.size()); - for (unsigned iPar = 0; iPar < params.size(); iPar++) { - const t_pppParam* par = params[iPar]; - AA[iPar] = par->partial(_epoTime, obs, tLC); - } - double ll = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA); - double vv = DotProduct(AA, _xFlt) - ll; - - if (fabs(vv) > SLIP) { - LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' ' - << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl; - resetAmb(obs->prn(), obsVector, tLC); - } - }*/ + resetAmb(obs->prn(), obsVector, t_lc()); + } } } @@ -371,5 +377,5 @@ // Reset Ambiguity Parameter (cycle slip) //////////////////////////////////////////////////////////////////////////// -t_irc t_pppFilter::resetAmb(const t_prn prn, const vector &obsVector, t_lc::type lc, +t_irc t_pppFilter::resetAmb(const t_prn prn, const vector &obsVector, t_lc lc, SymmetricMatrix *QSav, ColumnVector *xSav) { @@ -386,8 +392,8 @@ (par->firstObsTime().undef()) ? firstObsTime = lastObsTime : firstObsTime = par->firstObsTime(); - t_lc::type tLC = par->tLC(); - if (tLC != lc) {continue;} + t_lc LC = par->LC(); + if (lc.valid() && LC != lc) {continue;} LOG << string(_epoTime) << " RESET " << par->toString() << endl; - delete par; par = new t_pppParam(t_pppParam::amb, prn, tLC, &obsVector); + delete par; par = new t_pppParam(t_pppParam::amb, prn, LC, &obsVector); par->setIndex(ind); par->setFirstObsTime(firstObsTime); @@ -436,5 +442,5 @@ for (unsigned iPar = 0; iPar < numPar; iPar++) { t_pppParam* par = _parlist->params()[iPar]; - AA[_numSat - 1][iPar] = par->partial(_epoTime, obs, t_lc::c1); + AA[_numSat - 1][iPar] = par->partial(_epoTime, obs, obs->rangeLC()); if (par->type() == t_pppParam::crdX) { parX = par; Index: trunk/BNC/src/PPP/pppFilter.h =================================================================== --- trunk/BNC/src/PPP/pppFilter.h (revision 10765) +++ trunk/BNC/src/PPP/pppFilter.h (revision 10791) @@ -73,12 +73,12 @@ }; - t_irc processSystem(const std::vector& LCs, + t_irc processSystem(const std::vector& LCs, const std::vector& obsVector, bool pseudoObsIonoAvailable); - t_irc detectCycleSlips(const std::vector& LCs, + t_irc detectCycleSlips(const std::vector& LCs, const std::vector& obsVector); - t_irc resetAmb(t_prn prn, const std::vector& obsVector, t_lc::type lc, + t_irc resetAmb(t_prn prn, const std::vector& obsVector, t_lc lc, SymmetricMatrix* QSav = 0, ColumnVector* xSav = 0); Index: trunk/BNC/src/PPP/pppObsPool.cpp =================================================================== --- trunk/BNC/src/PPP/pppObsPool.cpp (revision 10765) +++ trunk/BNC/src/PPP/pppObsPool.cpp (revision 10791) @@ -15,5 +15,7 @@ * -----------------------------------------------------------------------*/ +#include #include "pppObsPool.h" +#include "pppClient.h" using namespace BNC_PPP; @@ -74,4 +76,11 @@ delete _satCodeBiases[iPrn]; _satCodeBiases[iPrn] = satCodeBias; + if (OPT->_logMode > t_pppOptions::normal) { + LOG << "codeBias " << string(satCodeBias->_time) << ' ' << satCodeBias->_prn.toString(); + for (const auto& bias : satCodeBias->_bias) { + LOG << " " << bias._rnxType2ch << ' ' << setw(7) << setprecision(3) << bias._value; + } + LOG << endl; + } } @@ -82,4 +91,45 @@ delete _satPhaseBiases[iPrn]; _satPhaseBiases[iPrn] = satPhaseBias; + if (OPT->_logMode > t_pppOptions::normal) { + LOG.setf(ios::fixed); + LOG << "phaseBias " << string(satPhaseBias->_time) << ' ' << satPhaseBias->_prn.toString() + << ' ' << setw(7) << setprecision(2) << satPhaseBias->_yaw * 180.0 / M_PI + << ' ' << setw(7) << setprecision(3) << satPhaseBias->_yawRate * 180.0 / M_PI; + for (const auto& bias : satPhaseBias->_bias) { + LOG << " " << bias._rnxType2ch << ' ' << setw(2) << bias._jumpCounter; + if (bias._fixIndicator) { + LOG << " x "; + } + else { + LOG << " . "; + } + LOG << setw(6) << setprecision(3) << bias._value; + } + LOG << endl; + } +} + +// +///////////////////////////////////////////////////////////////////////////// +void t_pppObsPool::clearCodeBiases(char sys) { + for (unsigned iPrn = 1; iPrn <= t_prn::MAXPRN; ++iPrn) { + t_prn prn; prn.set(iPrn); + if (prn.system() == sys) { + delete _satCodeBiases[iPrn]; + _satCodeBiases[iPrn] = 0; + } + } +} + +// +///////////////////////////////////////////////////////////////////////////// +void t_pppObsPool::clearPhaseBiases(char sys) { + for (unsigned iPrn = 1; iPrn <= t_prn::MAXPRN; ++iPrn) { + t_prn prn; prn.set(iPrn); + if (prn.system() == sys) { + delete _satPhaseBiases[iPrn]; + _satPhaseBiases[iPrn] = 0; + } + } } Index: trunk/BNC/src/PPP/pppObsPool.h =================================================================== --- trunk/BNC/src/PPP/pppObsPool.h (revision 10765) +++ trunk/BNC/src/PPP/pppObsPool.h (revision 10791) @@ -33,4 +33,6 @@ void putPhaseBias(t_satPhaseBias* satPhaseBias); void putTec(t_vTec* _vTec); + void clearCodeBiases(char sys); + void clearPhaseBiases(char sys); void putEpoch(const bncTime& epoTime, std::vector& obsVector, Index: trunk/BNC/src/PPP/pppParlist.cpp =================================================================== --- trunk/BNC/src/PPP/pppParlist.cpp (revision 10765) +++ trunk/BNC/src/PPP/pppParlist.cpp (revision 10791) @@ -34,10 +34,11 @@ // Constructor //////////////////////////////////////////////////////////////////////////// -t_pppParam::t_pppParam(e_type type, const t_prn& prn, t_lc::type tLC, +t_pppParam::t_pppParam(e_type type, const t_prn& prn, t_lc LC, const vector* obsVector) { _type = type; _prn = prn; - _tLC = tLC; + _sys = '?'; + _LC = LC; _x0 = 0.0; _indexOld = -1; @@ -47,70 +48,49 @@ switch (_type) { - case crdX: - _epoSpec = false; - _sigma0 = OPT->_aprSigCrd[0]; - _noise = OPT->_noiseCrd[0]; - break; - case crdY: - _epoSpec = false; - _sigma0 = OPT->_aprSigCrd[1]; - _noise = OPT->_noiseCrd[1]; - break; - case crdZ: - _epoSpec = false; - _sigma0 = OPT->_aprSigCrd[2]; - _noise = OPT->_noiseCrd[2]; - break; - case rClkG: - _epoSpec = true; - _sigma0 = OPT->_aprSigClk; - break; - case rClkR: - _epoSpec = true; - _sigma0 = OPT->_aprSigClk; - break; - case rClkE: - _epoSpec = true; - _sigma0 = OPT->_aprSigClk; - break; - case rClkC: - _epoSpec = true; - _sigma0 = OPT->_aprSigClk; - break; - case amb: - _epoSpec = false; - _sigma0 = OPT->_aprSigAmb; - if (obsVector) { - for (unsigned ii = 0; ii < obsVector->size(); ii++) { - const t_pppSatObs* obs = obsVector->at(ii); - if (obs->prn() == _prn) { - _x0 = floor((obs->obsValue(tLC) - obs->cmpValue(tLC)) / obs->lambda(tLC) + 0.5); - break; - } - } - } - break; - case trp: - _epoSpec = false; - _sigma0 = OPT->_aprSigTrp; - _noise = OPT->_noiseTrp; - break; - case ion: - _epoSpec = false; - _sigma0 = OPT->_aprSigIon; - _noise = OPT->_noiseIon; - break; - case cBiasG1: case cBiasR1: case cBiasE1: case cBiasC1: - case cBiasG2: case cBiasR2: case cBiasE2: case cBiasC2: - _epoSpec = false; - _sigma0 = OPT->_aprSigCodeBias; - _noise = OPT->_noiseCodeBias; - break; - case pBiasG1: case pBiasR1: case pBiasE1: case pBiasC1: - case pBiasG2: case pBiasR2: case pBiasE2: case pBiasC2: - _epoSpec = false; - _sigma0 = OPT->_aprSigPhaseBias; - _noise = OPT->_noisePhaseBias; - break; + case crdX: + _epoSpec = false; + _sigma0 = OPT->_aprSigCrd[0]; + _noise = OPT->_noiseCrd[0]; + break; + case crdY: + _epoSpec = false; + _sigma0 = OPT->_aprSigCrd[1]; + _noise = OPT->_noiseCrd[1]; + break; + case crdZ: + _epoSpec = false; + _sigma0 = OPT->_aprSigCrd[2]; + _noise = OPT->_noiseCrd[2]; + break; + case rClk: + _epoSpec = true; + _sigma0 = OPT->_aprSigClk; + break; + case amb: + _epoSpec = false; + _sigma0 = OPT->_aprSigAmb; + if (obsVector) { + for (unsigned ii = 0; ii < obsVector->size(); ii++) { + const t_pppSatObs* obs = obsVector->at(ii); + if (obs->prn() == _prn) { + _x0 = floor((obs->obsValue(LC) - obs->cmpValue(LC)) / obs->lambda(LC) + 0.5); + break; + } + } + } + break; + case trp: + _epoSpec = false; + _sigma0 = OPT->_aprSigTrp; + _noise = OPT->_noiseTrp; + break; + case ion: + _epoSpec = true; + _sigma0 = OPT->_aprSigIon; + break; + case bias: + _epoSpec = true; + _sigma0 = OPT->_aprSigCodeBias; + break; } } @@ -124,10 +104,21 @@ //////////////////////////////////////////////////////////////////////////// double t_pppParam::partial(const bncTime& /* epoTime */, const t_pppSatObs* obs, - const t_lc::type& tLC) const { + const t_lc& obsLC) const { // Special Case - Melbourne-Wuebbena // --------------------------------- - if (tLC == t_lc::MW && _type != amb) { + if (obsLC._type == t_lc::MW && _type != amb) { return 0.0; + } + + // Special Case - GIM + // ------------------ + if (obsLC._type == t_lc::GIM) { + if (_type == ion) { + return 1.0; + } + else { + return 0.0; + } } @@ -138,126 +129,39 @@ map phaseCoeff; map ionoCoeff; - obs->lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff); - + obs->lcCoeff(obsLC, codeCoeff, phaseCoeff, ionoCoeff); + switch (_type) { case crdX: - if (tLC == t_lc::GIM) {return 0.0;} return (sta->xyzApr()[0] - obs->xc()[0]) / rhoV.NormFrobenius(); case crdY: - if (tLC == t_lc::GIM) {return 0.0;} return (sta->xyzApr()[1] - obs->xc()[1]) / rhoV.NormFrobenius(); case crdZ: - if (tLC == t_lc::GIM) {return 0.0;} return (sta->xyzApr()[2] - obs->xc()[2]) / rhoV.NormFrobenius(); - case rClkG: - if (tLC == t_lc::GIM) {return 0.0;} - return (obs->prn().system() == 'G') ? 1.0 : 0.0; - case rClkR: - if (tLC == t_lc::GIM) {return 0.0;} - return (obs->prn().system() == 'R') ? 1.0 : 0.0; - case rClkE: - if (tLC == t_lc::GIM) {return 0.0;} - return (obs->prn().system() == 'E') ? 1.0 : 0.0; - case rClkC: - if (tLC == t_lc::GIM) {return 0.0;} - return (obs->prn().system() == 'C') ? 1.0 : 0.0; + case rClk: + return (_sys != obsLC.system() || obsLC.isGeometryFree()) ? 0.0 : 1.0; case amb: - if (tLC == t_lc::GIM) { - return 0.0; - } - else { - if (obs->prn() == _prn) { - if (tLC == _tLC) { - return (obs->lambda(tLC)); - } - else if (tLC == t_lc::lIF && _tLC == t_lc::MW) { - return obs->lambda(t_lc::lIF) * obs->lambda(t_lc::MW) / obs->lambda(t_lc::l2); - } - else { - if (_tLC == t_lc::l1) { - return obs->lambda(t_lc::l1) * phaseCoeff[obs->fType1()]; - } - else if (_tLC == t_lc::l2) { - return obs->lambda(t_lc::l2) * phaseCoeff[obs->fType2()]; - } - } + if (obs->prn() == _prn) { + if (obsLC == _LC) { + return (obs->lambda(obsLC)); + } + else if (_LC._type == t_lc::phase) { + return obs->lambda(_LC) * phaseCoeff[_LC._frq1]; + } + else if (obsLC._type == t_lc::phaseIF && _LC._type == t_lc::MW) { + return obs->lambda(obsLC) * obs->lambda(_LC) / obs->lambda(t_lc(t_lc::phase, _LC._frq2)); } } break; case trp: - if (tLC == t_lc::GIM) { - return 0.0; - } - else { - return 1.0 / sin(obs->eleSat()); - } + return 1.0 / sin(obs->eleSat()); case ion: if (obs->prn() == _prn) { - if (tLC == t_lc::c1) { - return ionoCoeff[obs->fType1()]; - } - else if (tLC == t_lc::c2) { - return ionoCoeff[obs->fType2()]; - } - else if (tLC == t_lc::l1) { - return ionoCoeff[obs->fType1()]; - } - else if (tLC == t_lc::l2) { - return ionoCoeff[obs->fType2()]; - } - else if (tLC == t_lc::GIM) { - return 1.0; - } - } - break; - case cBiasG1: - if ((obs->prn().system() == 'G') && (tLC == t_lc::c1)) {return 1.0;} else {return 0.0;} - break; - case cBiasR1: - if ((obs->prn().system() == 'R') && (tLC == t_lc::c1)) {return 1.0;} else {return 0.0;} - break; - case cBiasE1: - if ((obs->prn().system() == 'E') && (tLC == t_lc::c1)) {return 1.0;} else {return 0.0;} - break; - case cBiasC1: - if ((obs->prn().system() == 'C') && (tLC == t_lc::c1)) {return 1.0;} else {return 0.0;} - break; - case cBiasG2: - if ((obs->prn().system() == 'G') && (tLC == t_lc::c2)) {return 1.0;} else {return 0.0;} - break; - case cBiasR2: - if ((obs->prn().system() == 'R') && (tLC == t_lc::c2)) {return 1.0;} else {return 0.0;} - break; - case cBiasE2: - if ((obs->prn().system() == 'E') && (tLC == t_lc::c2)) {return 1.0;} else {return 0.0;} - break; - case cBiasC2: - if ((obs->prn().system() == 'C') && (tLC == t_lc::c2)) {return 1.0;} else {return 0.0;} - break; - case pBiasG1: - if ((obs->prn().system() == 'G') && (tLC == t_lc::l1)) {return 1.0;} else {return 0.0;} - break; - case pBiasR1: - if ((obs->prn().system() == 'R') && (tLC == t_lc::l1)) {return 1.0;} else {return 0.0;} - break; - case pBiasE1: - if ((obs->prn().system() == 'E') && (tLC == t_lc::l1)) {return 1.0;} else {return 0.0;} - break; - case pBiasC1: - if ((obs->prn().system() == 'C') && (tLC == t_lc::l1)) {return 1.0;} else {return 0.0;} - break; - case pBiasG2: - if ((obs->prn().system() == 'G') && (tLC == t_lc::l2)) {return 1.0;} else {return 0.0;} - break; - case pBiasR2: - if ((obs->prn().system() == 'R') && (tLC == t_lc::l2)) {return 1.0;} else {return 0.0;} - break; - case pBiasE2: - if ((obs->prn().system() == 'E') && (tLC == t_lc::l2)) {return 1.0;} else {return 0.0;} - break; - case pBiasC2: - if ((obs->prn().system() == 'C') && (tLC == t_lc::l2)) {return 1.0;} else {return 0.0;} - break; - + if (obsLC._type == t_lc::code || obsLC._type == t_lc::phase) { + return ionoCoeff[obsLC._frq1]; + } + } + break; + case bias: + return (_LC == obsLC ? 1.0 : 0.0); } return 0.0; @@ -278,15 +182,6 @@ ss << "CRD_Z"; break; - case rClkG: - ss << "REC_CLK G "; - break; - case rClkR: - ss << "REC_CLK R "; - break; - case rClkE: - ss << "REC_CLK E "; - break; - case rClkC: - ss << "REC_CLK C "; + case rClk: + ss << "REC_CLK " << _sys << " "; break; case trp: @@ -294,24 +189,12 @@ break; case amb: - ss << "AMB " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString(); + ss << "AMB " << left << setw(3) << _LC.toString() << right << ' ' << _prn.toString(); break; case ion: - ss << "ION " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString(); - break; - case cBiasG1: case pBiasG1: - case cBiasG2: case pBiasG2: - ss << "BIA " << left << setw(3) << t_lc::toString(_tLC) << right << " G "; - break; - case cBiasR1: case pBiasR1: - case cBiasR2: case pBiasR2: - ss << "BIA " << left << setw(3) << t_lc::toString(_tLC) << right << " R "; - break; - case cBiasE1: case pBiasE1: - case cBiasE2: case pBiasE2: - ss << "BIA " << left << setw(3) << t_lc::toString(_tLC) << right << " E "; - break; - case cBiasC1: case pBiasC1: - case cBiasC2: case pBiasC2: - ss << "BIA " << left << setw(3) << t_lc::toString(_tLC) << right << " C "; + ss << "ION " << left << setw(3) << _LC.toString() << right << ' ' << _prn.toString(); + break; + case bias: + char sys = t_frequency::toSystem(_LC._frq1); + ss << "BIA " << left << setw(3) << _LC.toString() << right << ' ' << sys << " "; break; } @@ -353,17 +236,5 @@ par->type() == t_pppParam::crdY || par->type() == t_pppParam::crdZ || - par->type() == t_pppParam::ion || - par->type() == t_pppParam::cBiasC1 || - par->type() == t_pppParam::cBiasC2 || - par->type() == t_pppParam::cBiasE1 || - par->type() == t_pppParam::cBiasE2 || - par->type() == t_pppParam::cBiasR1 || - par->type() == t_pppParam::cBiasR2 || - par->type() == t_pppParam::pBiasC1 || - par->type() == t_pppParam::pBiasC2 || - par->type() == t_pppParam::pBiasE1 || - par->type() == t_pppParam::pBiasE2 || - par->type() == t_pppParam::pBiasR1 || - par->type() == t_pppParam::pBiasR2) { + par->type() == t_pppParam::ion ) { if (par->lastObsTime().valid() && (epoTime - par->lastObsTime() > 60.0)) { remove = true; @@ -383,15 +254,16 @@ } - // check which systems have observations + // Check which systems have observations // ------------------------------------- - _usedSystems['G'] = _usedSystems['R'] = _usedSystems['E'] = _usedSystems['C'] = 0; - for (unsigned jj = 0; jj < obsVector.size(); jj++) { - const t_pppSatObs* satObs = obsVector[jj]; - char sys = satObs->prn().system(); - if (OPT->LCs(sys).size()) { - _usedSystems[sys]++; + vector systems = OPT->systems(); + for (char sys : systems) { + _usedSystems[sys] = 0; + for (unsigned jj = 0; jj < obsVector.size(); jj++) { + const t_pppSatObs* satObs = obsVector[jj]; + if (satObs->prn().system() == sys) { + _usedSystems[sys]++; + } } }; - // Check whether parameters have observations @@ -425,31 +297,22 @@ // Coordinates // ----------- - required.push_back(new t_pppParam(t_pppParam::crdX, t_prn(), t_lc::dummy)); - required.push_back(new t_pppParam(t_pppParam::crdY, t_prn(), t_lc::dummy)); - required.push_back(new t_pppParam(t_pppParam::crdZ, t_prn(), t_lc::dummy)); + required.push_back(new t_pppParam(t_pppParam::crdX, t_prn())); + required.push_back(new t_pppParam(t_pppParam::crdY, t_prn())); + required.push_back(new t_pppParam(t_pppParam::crdZ, t_prn())); // Receiver Clocks // --------------- - if (_usedSystems['G']) { - //if (OPT->useSystem('G')) { - required.push_back(new t_pppParam(t_pppParam::rClkG, t_prn(), t_lc::dummy)); - } - if (_usedSystems['R']) { - //if (OPT->useSystem('R')) { - required.push_back(new t_pppParam(t_pppParam::rClkR, t_prn(), t_lc::dummy)); - } - if (_usedSystems['E']) { - //if (OPT->useSystem('E')) { - required.push_back(new t_pppParam(t_pppParam::rClkE, t_prn(), t_lc::dummy)); - } - if (_usedSystems['C']) { - //if (OPT->useSystem('C')) { - required.push_back(new t_pppParam(t_pppParam::rClkC, t_prn(), t_lc::dummy)); - } + for (const auto& [sys, numObs] : _usedSystems) { + if (numObs > 0) { + t_pppParam* clk = new t_pppParam(t_pppParam::rClk, t_prn()); + clk->setSystem(sys); + required.push_back(clk); + } + } // Troposphere // ----------- if (OPT->estTrp()) { - required.push_back(new t_pppParam(t_pppParam::trp, t_prn(), t_lc::dummy)); + required.push_back(new t_pppParam(t_pppParam::trp, t_prn())); } @@ -459,12 +322,15 @@ for (unsigned jj = 0; jj < obsVector.size(); jj++) { const t_pppSatObs* satObs = obsVector[jj]; - char sys = satObs->prn().system(); - std::vector LCs = OPT->LCs(sys); - if (std::find(LCs.begin(), LCs.end(), t_lc::cIF) == LCs.end() && - std::find(LCs.begin(), LCs.end(), t_lc::lIF) == LCs.end()) { - required.push_back(new t_pppParam(t_pppParam::ion, satObs->prn(), t_lc::dummy)); - } - } - } + std::vector LCs = OPT->LCs(satObs->prn().system()); + for (auto it = LCs.begin(); it != LCs.end(); ++it) { + const t_lc& lc = *it; + if (!lc.isIonoFree()) { + required.push_back(new t_pppParam(t_pppParam::ion, satObs->prn())); + break; + } + } + } + } + // Ambiguities // ----------- @@ -472,116 +338,40 @@ const t_pppSatObs* satObs = obsVector[jj]; char sys = satObs->prn().system(); - const vector& ambLCs = OPT->ambLCs(sys); + const vector& ambLCs = OPT->ambLCs(sys); for (unsigned ii = 0; ii < ambLCs.size(); ii++) { + if (ambLCs[ii]._frq1 == t_frequency::G5 && !satObs->isValid(ambLCs[ii])) { + continue; + } required.push_back(new t_pppParam(t_pppParam::amb, satObs->prn(), ambLCs[ii], &obsVector)); } } - // Receiver Code Biases - // -------------------- - if (OPT->_ionoModelType == OPT->PPP_RTK) { - std::vector lc; - if (OPT->useSystem('G')) { - lc = OPT->LCs('G'); - if (std::find(lc.begin(), lc.end(), t_lc::c1) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::cBiasG1, t_prn(), t_lc::c1)); - } - if (std::find(lc.begin(), lc.end(), t_lc::c2) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::cBiasG2, t_prn(), t_lc::c2)); - } - } - if (OPT->useSystem('R')) { - lc = OPT->LCs('R'); - if (std::find(lc.begin(), lc.end(), t_lc::c1) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::cBiasR1, t_prn(), t_lc::c1)); - } - if (std::find(lc.begin(), lc.end(), t_lc::c2) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::cBiasR2, t_prn(), t_lc::c2)); - } - } - if (OPT->useSystem('E')) { - lc = OPT->LCs('E'); - if (std::find(lc.begin(), lc.end(), t_lc::c1) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::cBiasE1, t_prn(), t_lc::c1)); - } - if (std::find(lc.begin(), lc.end(), t_lc::c2) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::cBiasE2, t_prn(), t_lc::c2)); - } - } - if (OPT->useSystem('C')) { - lc = OPT->LCs('C'); - if (std::find(lc.begin(), lc.end(), t_lc::c1) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::cBiasC1, t_prn(), t_lc::c1)); - } - if (std::find(lc.begin(), lc.end(), t_lc::c2) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::cBiasC2, t_prn(), t_lc::c2)); - } - } - } - - if (OPT->_pseudoObsIono) { - std::vector lc; - if (OPT->useSystem('G')) { - lc = OPT->LCs('G'); - if (std::find(lc.begin(), lc.end(), t_lc::c2) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::cBiasG2, t_prn(), t_lc::c2)); - } - } - if (OPT->useSystem('R')) { - lc = OPT->LCs('R'); - if (std::find(lc.begin(), lc.end(), t_lc::c2) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::cBiasR2, t_prn(), t_lc::c2)); - } - } - if (OPT->useSystem('E')) { - lc = OPT->LCs('E'); - if (std::find(lc.begin(), lc.end(), t_lc::c2) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::cBiasE2, t_prn(), t_lc::c2)); - } - } - if (OPT->useSystem('C')) { - lc = OPT->LCs('C'); - if (std::find(lc.begin(), lc.end(), t_lc::c2) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::cBiasC2, t_prn(), t_lc::c2)); - } - } - } - - if (OPT->_ionoModelType == OPT->PPP_RTK) { - std::vector lc; - if (OPT->useSystem('G')) { - lc = OPT->LCs('G'); - if (std::find(lc.begin(), lc.end(), t_lc::l1) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::pBiasG1, t_prn(), t_lc::l1)); - } - if (std::find(lc.begin(), lc.end(), t_lc::l2) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::pBiasG2, t_prn(), t_lc::l2)); - } - } - if (OPT->useSystem('R')) { - lc = OPT->LCs('R'); - if (std::find(lc.begin(), lc.end(), t_lc::l1) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::pBiasR1, t_prn(), t_lc::l1)); - } - if (std::find(lc.begin(), lc.end(), t_lc::l2) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::pBiasR2, t_prn(), t_lc::l2)); - } - } - if (OPT->useSystem('E')) { - lc = OPT->LCs('E'); - if (std::find(lc.begin(), lc.end(), t_lc::l1) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::pBiasE1, t_prn(), t_lc::l1)); - } - if (std::find(lc.begin(), lc.end(), t_lc::l2) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::pBiasE2, t_prn(), t_lc::l2)); - } - } - if (OPT->useSystem('C')) { - lc = OPT->LCs('C'); - if (std::find(lc.begin(), lc.end(), t_lc::l1) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::pBiasC1, t_prn(), t_lc::l1)); - } - if (std::find(lc.begin(), lc.end(), t_lc::l2) != lc.end()) { - required.push_back(new t_pppParam(t_pppParam::pBiasC2, t_prn(), t_lc::l2)); + // Biases + // ------ + int maxSkip = (OPT->_pseudoObsIono ? 1 : 2); + for (const auto& [sys, numObs] : _usedSystems) { + if (numObs > 0) { + bool ar = OPT->arSystem(sys); + int skip = 0; + vector LCs = OPT->LCs(sys); + for (const t_lc& lc : LCs) { + if (ar) { + if (skip < maxSkip && lc.includesPhase()) { + skip += 1; + } + else { + required.push_back(new t_pppParam(t_pppParam::bias, t_prn(), lc)); + } + } + else { + if (lc.includesCode()) { + if (skip < maxSkip) { + skip += 1; + } + else { + required.push_back(new t_pppParam(t_pppParam::bias, t_prn(), lc)); + } + } + } } } @@ -605,7 +395,4 @@ } else { -#ifdef BNC_DEBUG_PPP - //LOG << "push_back parReq " << parReq->toString() << std::endl; -#endif _params.push_back(parReq); } @@ -634,5 +421,5 @@ //////////////////////////////////////////////////////////////////////////// void t_pppParlist::printResult(const bncTime& epoTime, const SymmetricMatrix& QQ, - const ColumnVector& xx) const { + const ColumnVector& xx, double fixRatio) const { string epoTimeStr = string(epoTime); @@ -703,7 +490,14 @@ << " dU = " << setprecision(4) << neu[2] << " +- " - << setprecision(4) << sqrt(QQneu[2][2]) - - << endl; + << setprecision(4) << sqrt(QQneu[2][2]); + + if (fixRatio > 0.0) { + LOG << " fix " << int(100*fixRatio) << " %"; + } + else { + LOG << " flt "; + } + + LOG << endl; } return; Index: trunk/BNC/src/PPP/pppParlist.h =================================================================== --- trunk/BNC/src/PPP/pppParlist.h (revision 10765) +++ trunk/BNC/src/PPP/pppParlist.h (revision 10791) @@ -14,18 +14,15 @@ class t_pppParam { public: - enum e_type {crdX, crdY, crdZ, rClkG, rClkR, rClkE, rClkC, trp, ion, amb, - cBiasG1, cBiasR1, cBiasE1, cBiasC1, pBiasG1, pBiasR1, pBiasE1, pBiasC1, - cBiasG2, cBiasR2, cBiasE2, cBiasC2, pBiasG2, pBiasR2, pBiasE2, pBiasC2}; + enum e_type {crdX, crdY, crdZ, rClk, trp, ion, amb, bias}; - t_pppParam(e_type type, const t_prn& prn, t_lc::type tLC, const std::vector* obsVector = 0); + t_pppParam(e_type type, const t_prn& prn, t_lc LC = t_lc(), const std::vector* obsVector = 0); ~t_pppParam(); e_type type() const {return _type;} double x0() const {return _x0;} - double partial(const bncTime& epoTime, const t_pppSatObs* obs, - const t_lc::type& tLC) const; + double partial(const bncTime& epoTime, const t_pppSatObs* obs, const t_lc& LC) const; bool epoSpec() const {return _epoSpec;} bool isEqual(const t_pppParam* par2) const { - return (_type == par2->_type && _prn == par2->_prn && _tLC == par2->_tLC); + return (_type == par2->_type && _prn == par2->_prn && _LC == par2->_LC && _sys == par2->_sys); } void setIndex(int indexNew) { @@ -36,10 +33,10 @@ _indexOld = -1; } - int indexOld() const {return _indexOld;} - int indexNew() const {return _indexNew;} - double sigma0() const {return _sigma0;} - double noise() const {return _noise;} - t_lc::type tLC() const {return _tLC;} - t_prn prn() const {return _prn;} + int indexOld() const {return _indexOld;} + int indexNew() const {return _indexNew;} + double sigma0() const {return _sigma0;} + double noise() const {return _noise;} + t_lc LC() const {return _LC;} + t_prn prn() const {return _prn;} std::string toString() const; @@ -55,4 +52,16 @@ unsigned ambNumEpo() const {return _ambInfo ? _ambInfo->_numEpo : 0;} void stepAmbNumEpo() {if (_ambInfo) _ambInfo->_numEpo += 1;} + char system() const { + if (_prn.valid()) { + return _prn.system(); + } + else if (_LC.valid()) { + return _LC.system(); + } + else { + return _sys; + } + } + void setSystem(char sys) {_sys = sys;} static bool sortFunction(const t_pppParam* p1, const t_pppParam* p2) { @@ -60,6 +69,6 @@ return p1->_type < p2->_type; } - else if (p1->_tLC != p2->_tLC) { - return p1->_tLC < p2->_tLC; + else if (p1->_LC != p2->_LC) { + return p1->_LC < p2->_LC; } else if (p1->_prn != p2->_prn) { @@ -84,5 +93,6 @@ e_type _type; t_prn _prn; - t_lc::type _tLC; + char _sys; + t_lc _LC; double _x0; bool _epoSpec; @@ -104,12 +114,11 @@ const std::vector& params() const {return _params;} std::vector& params() {return _params;} - const QMap& usedSystems() const {return _usedSystems;} void printResult(const bncTime& epoTime, const SymmetricMatrix& QQ, - const ColumnVector& xx) const; + const ColumnVector& xx, double fixRatio = 0.0) const; void printParams(const bncTime& epoTime); private: std::vector _params; - QMap _usedSystems; + std::map _usedSystems; }; Index: trunk/BNC/src/PPP/pppSatObs.cpp =================================================================== --- trunk/BNC/src/PPP/pppSatObs.cpp (revision 10765) +++ trunk/BNC/src/PPP/pppSatObs.cpp (revision 10791) @@ -19,4 +19,6 @@ #include #include +#include +#include #include @@ -43,5 +45,4 @@ _reference = false; _stecSat = 0.0; - _signalPriorities = QString::fromStdString(OPT->_signalPriorities); for (unsigned ii = 0; ii < t_frequency::max; ii++) { _obs[ii] = 0; @@ -60,72 +61,77 @@ // //////////////////////////////////////////////////////////////////////////// -void t_pppSatObs::prepareObs(const t_satObs& pppSatObs) { +bool t_pppSatObs::isBetter(const t_frqObs* aa, t_frqObs* bb, const string& trkModes) const { + + if (!trkModes.empty()) { + size_t posA = trkModes.find(aa->trkChar()); + size_t posB = trkModes.find(bb->trkChar()); + if (posA != posB) { + if (posA == string::npos) { + return false; + } + else if (posB == string::npos) { + return true; + } + else { + return posA < posB; + } + } + } + + unsigned numValA = 0; + if (aa->_codeValid) numValA += 1; + if (aa->_phaseValid) numValA += 1; + + unsigned numValB = 0; + if (bb->_codeValid) numValB += 1; + if (bb->_phaseValid) numValB += 1; + + if (numValA != numValB) { + return numValA > numValB; + } + + return false; +} + +// +//////////////////////////////////////////////////////////////////////////// +void t_pppSatObs::prepareObs(const t_satObs& satObs) { _model.reset(); - std::vector bb = OPT->frqBands(_prn.system()); - char frqNum1 = '0'; - if (bb.size() >= 1) { - frqNum1 = bb[0]; - } - char frqNum2 = '0'; - if (bb.size() == 2) { - frqNum2 = bb[1]; - } - - // Select pseudo-ranges and phase observations - // ------------------------------------------- - QStringList priorList = _signalPriorities.split(" ", Qt::SkipEmptyParts); - string preferredAttrib; - for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) { - t_frequency::type frqType = static_cast(iFreq); - char frqSys = t_frequency::toString(frqType)[0]; //cout << "frqSys: " << frqSys << endl; - char frqNum = t_frequency::toString(frqType)[1]; //cout << "frqNum: " << frqNum << endl; - if (frqSys != _prn.system()) { - continue; - } - if (frqNum != frqNum1 && - frqNum != frqNum2 ) { - continue; - } - QStringList hlp; - for (int ii = 0; ii < priorList.size(); ii++) { - if (priorList[ii].indexOf(":") != -1) { - hlp = priorList[ii].split(":", Qt::SkipEmptyParts); - if (hlp.size() == 2 && hlp[0].length() == 1 && hlp[0][0] == frqSys) { - hlp = hlp[1].split("&", Qt::SkipEmptyParts); - } - if (hlp.size() == 2 && hlp[0].indexOf(frqNum) != -1) { - preferredAttrib = hlp[1].toStdString(); //cout << "preferredAttrib: " << preferredAttrib << endl; - } - } - for (unsigned iPref = 0; iPref < preferredAttrib.length(); iPref++) { - QString obsType = QString("%1").arg(frqNum) + preferredAttrib[iPref]; //cout << "obstype: " << obsType.toStdString().c_str() << endl; - if (_obs[iFreq] == 0) { - for (unsigned ii = 0; ii < pppSatObs._obs.size(); ii++) { - const t_frqObs* obs = pppSatObs._obs[ii]; - //cout << "observation2char: " << obs->_rnxType2ch << " vs. " << obsType.toStdString().c_str()<< endl; - if (obs->_rnxType2ch == obsType.toStdString() && - obs->_codeValid && obs->_code && - obs->_phaseValid && obs->_phase) { - _obs[iFreq] = new t_frqObs(*obs); //cout << "================> newObs: " << obs->_rnxType2ch << " obs->_lockTime: " << obs->_lockTime << endl; - } - } - } - } - } - } - - // Used frequency types - // -------------------- - _fType1 = t_frqBand::toFreq(_prn.system(), frqNum1); - _fType2 = t_frqBand::toFreq(_prn.system(), frqNum2); - + const t_pppOptions::SysTrkModes* sysTrkModes = OPT->sysTrkModes(_prn.system()); + using FrqTrkModes = t_pppOptions::SysTrkModes::FrqTrkModes; + + for (const t_frqObs* obs : satObs._obs) { + if ( (obs->_codeValid || obs->_phaseValid)) { + t_frequency::type frq = t_frequency::toFreq(_prn.system(), obs->frqChar()); + if (frq == t_frequency::dummy) { + continue; + } + string trkModes; + if (sysTrkModes) { + const vector& frqTrkModes = sysTrkModes->_frqTrkModes; + auto it = find_if(frqTrkModes.begin(), frqTrkModes.end(), + [frq](const FrqTrkModes& mm){return mm._frq == frq;}); + if (it != frqTrkModes.end()) { + trkModes = it->_trkModes; + } + } + if (_obs[frq] == 0 || isBetter(obs, _obs[frq], trkModes)) { + delete _obs[frq]; + _obs[frq] = new t_frqObs(*obs); + } + } + } + // Check whether all required frequencies available // ------------------------------------------------ - for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) { - t_lc::type tLC = OPT->LCs(_prn.system())[ii]; - if (tLC == t_lc::GIM) {continue;} - if (!isValid(tLC)) { + const std::vector& LCs = OPT->LCs(_prn.system()); + for (unsigned ii = 0; ii < LCs.size(); ii++) { + t_lc LC = LCs[ii]; + if (LC._type == t_lc::GIM) { + continue; + } + if (LC._frq1 != t_frequency::G5 && !isValid(LC)) { _valid = false; return; @@ -148,19 +154,28 @@ bool totOK = false; ColumnVector satPosOld(6); satPosOld = 0.0; - t_lc::type tLC = t_lc::dummy; - if (isValid(t_lc::cIF)) { - tLC = t_lc::cIF; - } - if (tLC == t_lc::dummy && isValid(t_lc::c1)) { - tLC = t_lc::c1; - } - if (tLC == t_lc::dummy && isValid(t_lc::c2)) { - tLC = t_lc::c2; - } - if (tLC == t_lc::dummy) { - _valid = false; + + _valid = false; + t_frequency::type frq1 = t_frequency::dummy; + t_frequency::type frq2 = t_frequency::dummy; + OPT->defaultFrqs(_prn.system(), frq1, frq2); + if (frq1 != t_frequency::dummy) { + t_lc lc1(t_lc::code, frq1); + if (isValid(lc1)) { + _valid = true; + _rangeLC = lc1; + } + if (frq2 != t_frequency::dummy) { + t_lc lcIF(t_lc::codeIF, frq1, frq2); + if (isValid(lcIF)) { + _valid = true; + _rangeLC = lcIF; + } + } + } + if (!_valid) { return; } - double prange = obsValue(tLC); + + double prange = obsValue(_rangeLC); for (int ii = 1; ii <= 10; ii++) { bncTime ToT = _time - prange / t_CST::c - _xcSat[3]; @@ -188,5 +203,5 @@ // //////////////////////////////////////////////////////////////////////////// -void t_pppSatObs::lcCoeff(t_lc::type tLC, +void t_pppSatObs::lcCoeff(t_lc LC, map& codeCoeff, map& phaseCoeff, @@ -197,42 +212,34 @@ ionoCoeff.clear(); - double f1 = t_CST::freq(_fType1, _channel); - double f2 = t_CST::freq(_fType2, _channel); + double f1 = t_CST::freq(LC._frq1, _channel); + double f2 = t_CST::freq(LC._frq2, _channel); double f1GPS = t_CST::freq(t_frequency::G1, 0); - switch (tLC) { - case t_lc::l1: - phaseCoeff[_fType1] = 1.0; - ionoCoeff [_fType1] = -1.0 * pow(f1GPS, 2) / pow(f1, 2); + switch (LC._type) { + case t_lc::phase: + phaseCoeff[LC._frq1] = 1.0; + ionoCoeff [LC._frq1] = -1.0 * pow(f1GPS, 2) / pow(f1, 2); return; - case t_lc::l2: - phaseCoeff[_fType2] = 1.0; - ionoCoeff [_fType2] = -1.0 * pow(f1GPS, 2) / pow(f2, 2); + case t_lc::code: + codeCoeff[LC._frq1] = 1.0; + ionoCoeff[LC._frq1] = pow(f1GPS, 2) / pow(f1, 2); return; - case t_lc::lIF: - phaseCoeff[_fType1] = f1 * f1 / (f1 * f1 - f2 * f2); - phaseCoeff[_fType2] = -f2 * f2 / (f1 * f1 - f2 * f2); + case t_lc::phaseIF: + phaseCoeff[LC._frq1] = f1 * f1 / (f1 * f1 - f2 * f2); + phaseCoeff[LC._frq2] = -f2 * f2 / (f1 * f1 - f2 * f2); + return; + case t_lc::codeIF: + codeCoeff[LC._frq1] = f1 * f1 / (f1 * f1 - f2 * f2); + codeCoeff[LC._frq2] = -f2 * f2 / (f1 * f1 - f2 * f2); return; case t_lc::MW: - phaseCoeff[_fType1] = f1 / (f1 - f2); - phaseCoeff[_fType2] = -f2 / (f1 - f2); - codeCoeff[_fType1] = -f1 / (f1 + f2); - codeCoeff[_fType2] = -f2 / (f1 + f2); + phaseCoeff[LC._frq1] = f1 / (f1 - f2); + phaseCoeff[LC._frq2] = -f2 / (f1 - f2); + codeCoeff [LC._frq1] = -f1 / (f1 + f2); + codeCoeff [LC._frq2] = -f2 / (f1 + f2); return; case t_lc::CL: - phaseCoeff[_fType1] = 0.5; - codeCoeff [_fType1] = 0.5; - return; - case t_lc::c1: - codeCoeff[_fType1] = 1.0; - ionoCoeff[_fType1] = pow(f1GPS, 2) / pow(f1, 2); - return; - case t_lc::c2: - codeCoeff[_fType2] = 1.0; - ionoCoeff[_fType2] = pow(f1GPS, 2) / pow(f2, 2); - return; - case t_lc::cIF: - codeCoeff[_fType1] = f1 * f1 / (f1 * f1 - f2 * f2); - codeCoeff[_fType2] = -f2 * f2 / (f1 * f1 - f2 * f2); + phaseCoeff[LC._frq1] = 0.5; + codeCoeff [LC._frq1] = 0.5; return; case t_lc::GIM: @@ -245,7 +252,7 @@ // //////////////////////////////////////////////////////////////////////////// -bool t_pppSatObs::isValid(t_lc::type tLC) const { +bool t_pppSatObs::isValid(t_lc LC) const { bool valid = true; - obsValue(tLC, &valid); + obsValue(LC, &valid); return valid; @@ -253,5 +260,5 @@ // //////////////////////////////////////////////////////////////////////////// -double t_pppSatObs::obsValue(t_lc::type tLC, bool* valid) const { +double t_pppSatObs::obsValue(t_lc LC, bool* valid) const { double retVal = 0.0; @@ -259,5 +266,5 @@ // Pseudo observations - if (tLC == t_lc::GIM) { + if (LC._type == t_lc::GIM) { if (_stecSat == 0.0) { if (valid) *valid = false; @@ -272,5 +279,5 @@ map phaseCoeff; map ionoCoeff; - lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff); + lcCoeff(LC, codeCoeff, phaseCoeff, ionoCoeff); map::const_iterator it; @@ -303,22 +310,19 @@ // //////////////////////////////////////////////////////////////////////////// -double t_pppSatObs::lambda(t_lc::type tLC) const { - - double f1 = t_CST::freq(_fType1, _channel); - double f2 = t_CST::freq(_fType2, _channel); - - if (tLC == t_lc::l1) { +double t_pppSatObs::lambda(t_lc LC) const { + + double f1 = t_CST::freq(LC._frq1, _channel); + double f2 = t_CST::freq(LC._frq2, _channel); + + if (LC._type == t_lc::phase) { return t_CST::c / f1; } - else if (tLC == t_lc::l2) { - return t_CST::c / f2; - } - else if (tLC == t_lc::lIF) { + else if (LC._type == t_lc::phaseIF) { return t_CST::c / (f1 + f2); } - else if (tLC == t_lc::MW) { + else if (LC._type == t_lc::MW) { return t_CST::c / (f1 - f2); } - else if (tLC == t_lc::CL) { + else if (LC._type == t_lc::CL) { return t_CST::c / f1 / 2.0; } @@ -329,5 +333,5 @@ // //////////////////////////////////////////////////////////////////////////// -double t_pppSatObs::sigma(t_lc::type tLC) const { +double t_pppSatObs::sigma(t_lc LC) const { double retVal = 0.0; @@ -335,7 +339,7 @@ map phaseCoeff; map ionoCoeff; - lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff); - - if (tLC == t_lc::GIM) { + lcCoeff(LC, codeCoeff, phaseCoeff, ionoCoeff); + + if (LC._type == t_lc::GIM) { retVal = OPT->_sigmaGIM * OPT->_sigmaGIM; } @@ -354,5 +358,5 @@ // De-Weight R // ----------- - if (_prn.system() == 'R'&& t_lc::includesCode(tLC)) { + if (_prn.system() == 'R'&& LC.includesCode()) { retVal *= 5.0; } @@ -361,6 +365,6 @@ // ----------------------------- double cEle = 1.0; - if ( (OPT->_eleWgtCode && t_lc::includesCode(tLC)) || - (OPT->_eleWgtPhase && t_lc::includesPhase(tLC)) ) { + if ( (OPT->_eleWgtCode && LC.includesCode()) || + (OPT->_eleWgtPhase && LC.includesPhase()) ) { double eleD = eleSat()*180.0/M_PI; double hlp = fabs(90.0 - eleD); @@ -373,5 +377,5 @@ // //////////////////////////////////////////////////////////////////////////// -double t_pppSatObs::maxRes(t_lc::type tLC) const { +double t_pppSatObs::maxRes(t_lc LC) const { double retVal = 0.0; @@ -379,5 +383,5 @@ map phaseCoeff; map ionoCoeff; - lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff); + lcCoeff(LC, codeCoeff, phaseCoeff, ionoCoeff); map::const_iterator it; @@ -388,5 +392,5 @@ retVal += it->second * it->second * OPT->_maxResL1 * OPT->_maxResL1; } - if (tLC == t_lc::GIM) { + if (LC._type == t_lc::GIM) { retVal = OPT->_maxResGIM * OPT->_maxResGIM + OPT->_maxResGIM * OPT->_maxResGIM; } @@ -517,7 +521,6 @@ } const t_frqObs* obs = _obs[iFreq]; - if (obs && - obs->_rnxType2ch == bias._rnxType2ch) { - _model._codeBias[iFreq] = bias._value; + if (obs && obs->_rnxType2ch == bias._rnxType2ch) { + _model._codeBias[iFreq] = (bias._value != 0.0 ? bias._value : ZEROVALUE); } } @@ -527,25 +530,31 @@ // Phase Biases // ----------- - const t_satPhaseBias* satPhaseBias = PPP_CLIENT->obsPool()->satPhaseBias(_prn); - double yaw = 0.0; - bool ssr = false; - if (satPhaseBias) { - double dt = station->epochTime() - satPhaseBias->_time; - if (satPhaseBias->_updateInt) { - dt -= (0.5 * ssrUpdateInt[satPhaseBias->_updateInt]); - } - yaw = satPhaseBias->_yaw + satPhaseBias->_yawRate * dt; - ssr = true; - for (unsigned ii = 0; ii < satPhaseBias->_bias.size(); ii++) { - const t_frqPhaseBias& bias = satPhaseBias->_bias[ii]; - for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) { - string frqStr = t_frequency::toString(t_frequency::type(iFreq)); - if (frqStr[0] != _prn.system()) { - continue; - } - const t_frqObs* obs = _obs[iFreq]; - if (obs && - obs->_rnxType2ch == bias._rnxType2ch) { - _model._phaseBias[iFreq] = bias._value; + double yaw = 0.0; + bool useYaw = false; + if (OPT->arSystem(_prn.system())) { + const t_satPhaseBias* satPhaseBias = PPP_CLIENT->obsPool()->satPhaseBias(_prn); + if (satPhaseBias) { + if (OPT->_ar._useYaw) { + double dt = station->epochTime() - satPhaseBias->_time; + if (satPhaseBias->_updateInt) { + dt -= (0.5 * ssrUpdateInt[satPhaseBias->_updateInt]); + } + yaw = satPhaseBias->_yaw + satPhaseBias->_yawRate * dt; + useYaw = true; + } + for (unsigned ii = 0; ii < satPhaseBias->_bias.size(); ii++) { + const t_frqPhaseBias& bias = satPhaseBias->_bias[ii]; + if (bias._fixIndicator) { // if AR, biases without fixIndicator not used + for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) { + string frqStr = t_frequency::toString(t_frequency::type(iFreq)); + if (frqStr[0] != _prn.system()) { + continue; + } + t_frqObs* obs = _obs[iFreq]; + if (obs && obs->_rnxType2ch[0] == bias._rnxType2ch[0]) { // allow different tracking mode + _model._phaseBias[iFreq] = (bias._value != 0.0 ? bias._value : ZEROVALUE); + obs->_biasJumpCounter = bias._jumpCounter; + } + } } } @@ -555,5 +564,5 @@ // Phase Wind-Up // ------------- - _model._windUp = station->windUp(_time, _prn, rSat, ssr, yaw, vSat) ; + _model._windUp = station->windUp(_time, _prn, rSat, useYaw, yaw, vSat) ; // Relativistic effect due to earth gravity @@ -574,6 +583,6 @@ bool vTecUsage = true; for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) { - t_lc::type tLC = OPT->LCs(_prn.system())[ii]; - if (tLC == t_lc::cIF || tLC == t_lc::lIF) { + t_lc LC = OPT->LCs(_prn.system())[ii]; + if (LC._type == t_lc::codeIF || LC._type == t_lc::phaseIF) { vTecUsage = false; } @@ -601,5 +610,7 @@ _model._set = true; - //printModel(); + if (OPT->_logMode == t_pppOptions::all) { + printModel(); + } return success; @@ -636,6 +647,6 @@ if (_prn.system() == frqStr[0]) { LOG << "PCO : " << frqStr << setw(12) << setprecision(3) << _model._antPCO[iFreq] << endl - << "BIAS CODE : " << frqStr << setw(12) << setprecision(3) << _model._codeBias[iFreq] << "\t(" << _obs[iFreq]->_rnxType2ch[1] << ") " << endl - << "BIAS PHASE : " << frqStr << setw(12) << setprecision(3) << _model._phaseBias[iFreq] << "\t(" << _obs[iFreq]->_rnxType2ch[1] << ") " << endl + << "BIAS CODE : " << frqStr << setw(12) << setprecision(3) << _model._codeBias[iFreq] << "\t(" << _obs[iFreq]->trkChar() << ") " << endl + << "BIAS PHASE : " << frqStr << setw(12) << setprecision(3) << _model._phaseBias[iFreq] << "\t(" << _obs[iFreq]->trkChar() << ") " << endl << "IONO CODEDELAY: " << frqStr << setw(12) << setprecision(3) << _model._ionoCodeDelay[iFreq]<< endl; } @@ -652,27 +663,27 @@ char sys = _prn.system(); for (unsigned ii = 0; ii < OPT->LCs(sys).size(); ii++) { - t_lc::type tLC = OPT->LCs(sys)[ii]; - LOG << "OBS-CMP " << setw(4) << t_lc::toString(tLC) << ": " << _prn.toString() << " " - << setw(12) << setprecision(3) << obsValue(tLC) << " " - << setw(12) << setprecision(3) << cmpValue(tLC) << " " - << setw(12) << setprecision(3) << obsValue(tLC) - cmpValue(tLC) << endl; - } -} - -// -//////////////////////////////////////////////////////////////////////////// -double t_pppSatObs::cmpValueForBanc(t_lc::type tLC) const { - return cmpValue(tLC) - _model._rho - _model._sagnac - _model._recClkM; -} - -// -//////////////////////////////////////////////////////////////////////////// -double t_pppSatObs::cmpValue(t_lc::type tLC) const { + t_lc LC = OPT->LCs(sys)[ii]; + LOG << "OBS-CMP " << setw(4) << LC.toString() << ": " << _prn.toString() << " " + << setw(12) << setprecision(3) << obsValue(LC) << " " + << setw(12) << setprecision(3) << cmpValue(LC) << " " + << setw(12) << setprecision(3) << obsValue(LC) - cmpValue(LC) << endl; + } +} + +// +//////////////////////////////////////////////////////////////////////////// +double t_pppSatObs::cmpValueForBanc(t_lc LC) const { + return cmpValue(LC) - _model._rho - _model._sagnac - _model._recClkM; +} + +// +//////////////////////////////////////////////////////////////////////////// +double t_pppSatObs::cmpValue(t_lc LC) const { double cmpValue; - if (!isValid(tLC)) { + if (!isValid(LC)) { cmpValue = 0.0; } - else if (tLC == t_lc::GIM) { + else if (LC._type == t_lc::GIM) { cmpValue = _stecSat; } @@ -691,12 +702,9 @@ map phaseCoeff; map ionoCoeff; - lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff); + lcCoeff(LC, codeCoeff, phaseCoeff, ionoCoeff); map::const_iterator it; for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) { t_frequency::type tFreq = it->first; dispPart += it->second * (_model._antPCO[tFreq] - _model._codeBias[tFreq]); - if (OPT->PPP_RTK) { - dispPart += it->second * (_model._ionoCodeDelay[tFreq]); - } } for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) { @@ -704,7 +712,4 @@ dispPart += it->second * (_model._antPCO[tFreq] - _model._phaseBias[tFreq] + _model._windUp * t_CST::lambda(tFreq, _channel)); - if (OPT->PPP_RTK) { - dispPart += it->second * (- _model._ionoCodeDelay[tFreq]); - } } cmpValue = nonDisp + dispPart; @@ -716,12 +721,12 @@ // //////////////////////////////////////////////////////////////////////////// -void t_pppSatObs::setRes(t_lc::type tLC, double res) { - _res[tLC] = res; -} - -// -//////////////////////////////////////////////////////////////////////////// -double t_pppSatObs::getRes(t_lc::type tLC) const { - map::const_iterator it = _res.find(tLC); +void t_pppSatObs::setRes(t_lc LC, double res) { + _res[LC] = res; +} + +// +//////////////////////////////////////////////////////////////////////////// +double t_pppSatObs::getRes(t_lc LC) const { + map::const_iterator it = _res.find(LC); if (it != _res.end()) { return it->second; @@ -742,2 +747,21 @@ return pseudoObsIono; } + +// +//////////////////////////////////////////////////////////////////////////// +bool t_pppSatObs::hasBiases() const { + bool ar = OPT->arSystem(_prn.system()); + set frqs; + for (const auto& lc : OPT->LCs(_prn.system())) { + if (lc._frq1 != t_frequency::dummy) frqs.insert(lc._frq1); + if (lc._frq2 != t_frequency::dummy) frqs.insert(lc._frq2); + } + for (int iFreq : frqs) { + if (_obs[iFreq] != 0) { + if (_model._codeBias[iFreq] == 0 || (ar && _model._phaseBias[iFreq] == 0)) { + return false; + } + } + } + return true; +} Index: trunk/BNC/src/PPP/pppSatObs.h =================================================================== --- trunk/BNC/src/PPP/pppSatObs.h (revision 10765) +++ trunk/BNC/src/PPP/pppSatObs.h (revision 10791) @@ -11,148 +11,148 @@ namespace BNC_PPP { -class t_pppStation; + class t_pppStation; -class t_pppSatObs { - public: - t_pppSatObs(const t_satObs& satObs); - ~t_pppSatObs(); - bool isValid() const {return _valid;}; - bool isValid(t_lc::type tLC) const; - bool isReference() const {return _reference;}; - void setAsReference() {_reference = true;}; - void resetReference() {_reference = false;}; - const t_prn& prn() const {return _prn;} - const ColumnVector& xc() const {return _xcSat;} - const bncTime& time() const {return _time;} - t_irc cmpModel(const t_pppStation* station); - double obsValue(t_lc::type tLC, bool* valid = 0) const; - double cmpValue(t_lc::type tLC) const; - double cmpValueForBanc(t_lc::type tLC) const; - double rho() const {return _model._rho;} - double sagnac() const {return _model._sagnac;} - double eleSat() const {return _model._eleSat;} - bool modelSet() const {return _model._set;} - void printModel() const; - void printObsMinusComputed() const; - void lcCoeff(t_lc::type tLC, - std::map& codeCoeff, - std::map& phaseCoeff, - std::map& ionoCoeff) const; - double lambda(t_lc::type tLC) const; - double sigma(t_lc::type tLC) const; - double maxRes(t_lc::type tLC) const; - bool outlier() const {return _outlier;} - void setOutlier() {_outlier = true;} - void resetOutlier() {_outlier = false;} - void setRes(t_lc::type tLC, double res); - double getRes(t_lc::type tLC) const; - bool setPseudoObsIono(t_frequency::type freq); - double getIonoCodeDelay(t_frequency::type freq) {return _model._ionoCodeDelay[freq];} - double getCodeBias(t_frequency::type freq) {return _model._codeBias[freq];} - t_frequency::type fType1() const {return _fType1;} - t_frequency::type fType2() const {return _fType2;} + class t_pppSatObs { + public: + t_pppSatObs(const t_satObs& satObs); + ~t_pppSatObs(); + bool isValid() const {return _valid;}; + bool isValid(t_lc LC) const; + t_lc rangeLC() const {return _rangeLC;}; + bool isReference() const {return _reference;}; + void setAsReference() {_reference = true;}; + void resetReference() {_reference = false;}; + const t_prn& prn() const {return _prn;} + const ColumnVector& xc() const {return _xcSat;} + const bncTime& time() const {return _time;} + t_irc cmpModel(const t_pppStation* station); + double obsValue(t_lc LC, bool* valid = 0) const; + double cmpValue(t_lc LC) const; + double cmpValueForBanc(t_lc LC) const; + double rho() const {return _model._rho;} + double sagnac() const {return _model._sagnac;} + double eleSat() const {return _model._eleSat;} + bool modelSet() const {return _model._set;} + void printModel() const; + void printObsMinusComputed() const; + void lcCoeff(t_lc LC, + std::map& codeCoeff, + std::map& phaseCoeff, + std::map& ionoCoeff) const; + double lambda(t_lc LC) const; + double sigma(t_lc LC) const; + double maxRes(t_lc LC) const; + bool outlier() const {return _outlier;} + void setOutlier() {_outlier = true;} + void resetOutlier() {_outlier = false;} + void setRes(t_lc LC, double res); + double getRes(t_lc LC) const; + bool setPseudoObsIono(t_frequency::type freq); + double getIonoCodeDelay(t_frequency::type freq) {return _model._ionoCodeDelay[freq];} + double getCodeBias(t_frequency::type freq) {return _model._codeBias[freq];} + bool hasBiases() const; - // RINEX - bool slip() const { - for (unsigned ii = 1; ii < t_frequency::max; ii++) { - if (_obs[ii] && _obs[ii]->_slip) { - return true; + // RINEX + bool slip() const { + for (unsigned ii = 1; ii < t_frequency::max; ii++) { + if (_obs[ii] && _obs[ii]->_slip) { + return true; + } } + return false; } - return false; - } - // RTCM - int slipCounter() const { - int cnt = -1; - for (unsigned ii = 1; ii < t_frequency::max; ii++) { - if (_obs[ii] && _obs[ii]->_slipCounter > cnt) { - cnt = _obs[ii]->_slipCounter; + // RTCM + int slipCounter() const { + int cnt = -1; + for (unsigned ii = 1; ii < t_frequency::max; ii++) { + if (_obs[ii] && _obs[ii]->_slipCounter > cnt) { + cnt = _obs[ii]->_slipCounter; + } } + return cnt; } - return cnt; - } - int biasJumpCounter() const { - int jmp = -1; - for (unsigned ii = 1; ii < t_frequency::max; ii++) { - if (_obs[ii] && _obs[ii]->_biasJumpCounter > jmp) { - jmp = _obs[ii]->_biasJumpCounter; + int biasJumpCounter() const { + int jmp = -1; + for (unsigned ii = 1; ii < t_frequency::max; ii++) { + if (_obs[ii] && _obs[ii]->_biasJumpCounter > jmp) { + jmp = _obs[ii]->_biasJumpCounter; + } } + return jmp; } - return jmp; - } - private: - class t_model { - public: - t_model() {reset();} - ~t_model() {} - void reset() { - _set = false; - _rho = 0.0; - _eleSat = 0.0; - _azSat = 0.0; - _elTx = 0.0; - _azTx = 0.0; - _recClkM = 0.0; - _satClkM = 0.0; - _sagnac = 0.0; - _antEcc = 0.0; - _tropo = 0.0; - _tropo0 = 0.0; - _tideEarth = 0.0; - _tideOcean = 0.0; - _windUp = 0.0; - _rel = 0.0; - for (unsigned ii = 0; ii < t_frequency::max; ii++) { - _antPCO[ii] = 0.0; - _codeBias[ii] = 0.0; - _phaseBias[ii] = 0.0; - _ionoCodeDelay[ii] = 0.0; + private: + class t_model { + public: + t_model() {reset();} + ~t_model() {} + void reset() { + _set = false; + _rho = 0.0; + _eleSat = 0.0; + _azSat = 0.0; + _elTx = 0.0; + _azTx = 0.0; + _recClkM = 0.0; + _satClkM = 0.0; + _sagnac = 0.0; + _antEcc = 0.0; + _tropo = 0.0; + _tropo0 = 0.0; + _tideEarth = 0.0; + _tideOcean = 0.0; + _windUp = 0.0; + _rel = 0.0; + for (unsigned ii = 0; ii < t_frequency::max; ii++) { + _antPCO[ii] = 0.0; + _codeBias[ii] = 0.0; + _phaseBias[ii] = 0.0; + _ionoCodeDelay[ii] = 0.0; + } } - } - bool _set; - double _rho; - double _eleSat; - double _azSat; - double _elTx; - double _azTx; - double _recClkM; - double _satClkM; - double _sagnac; - double _antEcc; - double _tropo; - double _tropo0; - double _tideEarth; - double _tideOcean; - double _windUp; - double _rel; - double _antPCO[t_frequency::max]; - double _codeBias[t_frequency::max]; - double _phaseBias[t_frequency::max]; - double _ionoCodeDelay[t_frequency::max]; + bool _set; + double _rho; + double _eleSat; + double _azSat; + double _elTx; + double _azTx; + double _recClkM; + double _satClkM; + double _sagnac; + double _antEcc; + double _tropo; + double _tropo0; + double _tideEarth; + double _tideOcean; + double _windUp; + double _rel; + double _antPCO[t_frequency::max]; + double _codeBias[t_frequency::max]; + double _phaseBias[t_frequency::max]; + double _ionoCodeDelay[t_frequency::max]; + }; + + void prepareObs(const t_satObs& satObs); + + bool isBetter(const t_frqObs* aa, t_frqObs* bb, const std::string& trkModes) const; + + bool _valid; + bool _reference; + t_lc _rangeLC; + t_prn _prn; + bncTime _time; + int _channel; + t_frqObs* _obs[t_frequency::max]; + ColumnVector _xcSat; + ColumnVector _vvSat; + t_model _model; + bool _outlier; + std::map _res; + double _signalPropagationTime; + double _stecSat; + double _tropo0; }; - - void prepareObs(const t_satObs& satObs); - - bool _valid; - bool _reference; - t_frequency::type _fType1; - t_frequency::type _fType2; - t_prn _prn; - bncTime _time; - int _channel; - t_frqObs* _obs[t_frequency::max]; - ColumnVector _xcSat; - ColumnVector _vvSat; - t_model _model; - bool _outlier; - std::map _res; - double _signalPropagationTime; - double _stecSat; - double _tropo0; - QString _signalPriorities; - }; }