Changeset 8905 in ntrip for trunk/BNC/src/PPP/pppSatObs.cpp
- Timestamp:
- Mar 18, 2020, 11:13:50 AM (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/BNC/src/PPP/pppSatObs.cpp
r8619 r8905 35 35 //////////////////////////////////////////////////////////////////////////// 36 36 t_pppSatObs::t_pppSatObs(const t_satObs& pppSatObs) { 37 _prn = pppSatObs._prn; 38 _time = pppSatObs._time; 39 _outlier = false; 40 _valid = true; 37 _prn = pppSatObs._prn; 38 _time = pppSatObs._time; 39 _outlier = false; 40 _valid = true; 41 _reference = false; 42 _stecRefSat = 0.0; 43 _stecSat = 0.0; 41 44 for (unsigned ii = 0; ii < t_frequency::max; ii++) { 42 45 _obs[ii] = 0; … … 61 64 // Select pseudoranges and phase observations 62 65 // ------------------------------------------ 63 const string preferredAttrib = "CWPXI_"; 64 //const string preferredAttrib = "G:12&PWCSLXYN G:5&IQX R:12&PC R:3&IQX E:16&BCX E:578&IQX J:1&SLXCZ J:26&SLX J:5&IQX C:IQX I:ABCX S:1&C S:5&IQX"; 66 const string preferredAttrib = "G:12&PWCSLXYN G:5&IQX R:12&PC R:3&IQX E:16&BCX E:578&IQX J:1&SLXCZ J:26&SLX J:5&IQX C:IQX I:ABCX S:1&C S:5&IQX"; 65 67 66 68 for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) { … … 90 92 for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) { 91 93 t_lc::type tLC = OPT->LCs(_prn.system())[ii]; 94 if (tLC == t_lc::GIM) {continue;} 92 95 if (!isValid(tLC)) { 93 96 _valid = false; … … 121 124 ColumnVector dx = _xcSat - satPosOld; 122 125 dx[3] *= t_CST::c; 123 if (dx. norm_Frobenius() < 1.e-4) {126 if (dx.NormFrobenius() < 1.e-4) { 124 127 totOK = true; 125 128 break; … … 140 143 void t_pppSatObs::lcCoeff(t_lc::type tLC, 141 144 map<t_frequency::type, double>& codeCoeff, 142 map<t_frequency::type, double>& phaseCoeff) const { 145 map<t_frequency::type, double>& phaseCoeff, 146 map<t_frequency::type, double>& ionoCoeff) const { 143 147 144 148 codeCoeff.clear(); 145 149 phaseCoeff.clear(); 150 ionoCoeff.clear(); 146 151 147 152 double f1 = t_CST::freq(_fType1, _channel); 148 153 double f2 = t_CST::freq(_fType2, _channel); 154 double f1GPS = t_CST::freq(t_frequency::G1, 0); 149 155 150 156 switch (tLC) { 151 157 case t_lc::l1: 152 phaseCoeff[_fType1] = 1.0; 158 phaseCoeff[_fType1] = 1.0; 159 ionoCoeff [_fType1] = -1.0 * pow(f1GPS, 2) / pow(f1, 2); 153 160 return; 154 161 case t_lc::l2: 155 phaseCoeff[_fType2] = 1.0; 162 phaseCoeff[_fType2] = 1.0; 163 ionoCoeff [_fType2] = -1.0 * pow(f1GPS, 2) / pow(f2, 2); 156 164 return; 157 165 case t_lc::lIF: … … 167 175 case t_lc::CL: 168 176 phaseCoeff[_fType1] = 0.5; 169 codeCoeff [_fType1]= 0.5;177 codeCoeff [_fType1] = 0.5; 170 178 return; 171 179 case t_lc::c1: 172 180 codeCoeff[_fType1] = 1.0; 181 ionoCoeff[_fType1] = pow(f1GPS, 2) / pow(f1, 2); 173 182 return; 174 183 case t_lc::c2: 175 184 codeCoeff[_fType2] = 1.0; 185 ionoCoeff[_fType2] = pow(f1GPS, 2) / pow(f2, 2); 176 186 return; 177 187 case t_lc::cIF: … … 179 189 codeCoeff[_fType2] = -f2 * f2 / (f1 * f1 - f2 * f2); 180 190 return; 191 case t_lc::GIM: 181 192 case t_lc::dummy: 182 193 case t_lc::maxLc: … … 190 201 bool valid = true; 191 202 obsValue(tLC, &valid); 203 //qDebug() << "tLC: " << tLC << " valid: " << valid; 192 204 return valid; 193 205 } … … 195 207 //////////////////////////////////////////////////////////////////////////// 196 208 double t_pppSatObs::obsValue(t_lc::type tLC, bool* valid) const { 209 210 double retVal = 0.0; 211 if (valid) *valid = true; 212 213 // Pseudo observations 214 if (tLC == t_lc::GIM) { 215 if (_stecRefSat == 0.0 || _stecSat == 0.0) { 216 if (valid) *valid = false; 217 return 0.0; 218 } 219 else { 220 return _stecRefSat - _stecSat; 221 } 222 } 197 223 198 224 map<t_frequency::type, double> codeCoeff; 199 225 map<t_frequency::type, double> phaseCoeff; 200 lcCoeff(tLC, codeCoeff, phaseCoeff); 201 202 double retVal = 0.0; 203 if (valid) *valid = true; 226 map<t_frequency::type, double> ionoCoeff; 227 lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff); 204 228 205 229 map<t_frequency::type, double>::const_iterator it; 230 231 // Code observations 206 232 for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) { 207 233 t_frequency::type tFreq = it->first; … … 214 240 } 215 241 } 242 // Phase observations 216 243 for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) { 217 244 t_frequency::type tFreq = it->first; … … 224 251 } 225 252 } 226 227 253 return retVal; 228 254 } … … 258 284 double t_pppSatObs::sigma(t_lc::type tLC) const { 259 285 286 double retVal = 0.0; 260 287 map<t_frequency::type, double> codeCoeff; 261 288 map<t_frequency::type, double> phaseCoeff; 262 lcCoeff(tLC, codeCoeff, phaseCoeff); 263 264 double retVal = 0.0; 289 map<t_frequency::type, double> ionoCoeff; 290 lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff); 291 292 if (tLC == t_lc::GIM) { 293 retVal = OPT->_sigmaGIMdiff * OPT->_sigmaGIMdiff; 294 } 265 295 266 296 map<t_frequency::type, double>::const_iterator it; 267 for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {297 for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {//qDebug() << "codeCoeff : " << t_frequency::toString(it->first).c_str() << ": " << it->second; 268 298 retVal += it->second * it->second * OPT->_sigmaC1 * OPT->_sigmaC1; 269 299 } 270 for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) { 300 301 for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {//qDebug() << "phaseCoeff: " << t_frequency::toString(it->first).c_str() << ": " << it->second; 271 302 retVal += it->second * it->second * OPT->_sigmaL1 * OPT->_sigmaL1; 272 303 } … … 295 326 // 296 327 //////////////////////////////////////////////////////////////////////////// 297 double t_pppSatObs::maxRes(t_lc::type tLC) const { 328 double t_pppSatObs::maxRes(t_lc::type tLC) const {//qDebug() << "t_pppSatObs::maxRes(t_lc::type tLC)"; 329 double retVal = 0.0; 298 330 299 331 map<t_frequency::type, double> codeCoeff; 300 332 map<t_frequency::type, double> phaseCoeff; 301 lcCoeff(tLC, codeCoeff, phaseCoeff); 302 303 double retVal = 0.0; 333 map<t_frequency::type, double> ionoCoeff; 334 lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff); 304 335 305 336 map<t_frequency::type, double>::const_iterator it; 306 for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {337 for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {//qDebug() << "codeCoeff: " << it->first << ": " << it->second; 307 338 retVal += it->second * it->second * OPT->_maxResC1 * OPT->_maxResC1; 308 339 } 309 for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) { 340 for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {//qDebug() << "phaseCoeff: " << it->first << ": " << it->second; 310 341 retVal += it->second * it->second * OPT->_maxResL1 * OPT->_maxResL1; 311 342 } 312 343 if (tLC == t_lc::GIM) { 344 retVal = 3 * (OPT->_sigmaGIMdiff * OPT->_sigmaGIMdiff); 345 } 313 346 return sqrt(retVal); 314 347 } … … 326 359 // ------------------------------ 327 360 ColumnVector rSat = _xcSat.Rows(1,3); 328 ColumnVector rhoV = rSat - station->xyzApr(); 329 _model._rho = rhoV.norm_Frobenius(); 361 ColumnVector rRec = station->xyzApr(); 362 ColumnVector rhoV = rSat - rRec; 363 _model._rho = rhoV.NormFrobenius(); 330 364 331 365 ColumnVector vSat = _vvSat; … … 334 368 xyz2neu(station->ellApr().data(), rhoV.data(), neu.data()); 335 369 336 _model._eleSat = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / _model._rho);370 _model._eleSat = acos(sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / _model._rho); 337 371 if (neu[2] < 0) { 338 372 _model._eleSat *= -1.0; … … 354 388 Omega[1] = 0.0; 355 389 Omega[2] = t_CST::omega / t_CST::c; 356 _model._sagnac = DotProduct(Omega, crossproduct(rSat, station->xyzApr()));390 _model._sagnac = DotProduct(Omega, crossproduct(rSat, rRec)); 357 391 358 392 // Antenna Eccentricity … … 373 407 // Tropospheric Delay 374 408 // ------------------ 375 _model._tropo = t_tropo::delay_saast( station->xyzApr(), _model._eleSat);409 _model._tropo = t_tropo::delay_saast(rRec, _model._eleSat); 376 410 377 411 // Code Biases … … 392 426 // Phase Biases 393 427 // ----------- 394 // TODO: consideration of fix indicators and jump counter395 428 const t_satPhaseBias* satPhaseBias = PPP_CLIENT->obsPool()->satPhaseBias(_prn); 396 429 double yaw = 0.0; 397 430 bool ssr = false; 398 431 if (satPhaseBias) { 399 yaw = satPhaseBias->_yaw; 432 double dt = station->epochTime() - satPhaseBias->_time; 433 if (satPhaseBias->_updateInt) { 434 dt -= (0.5 * ssrUpdateInt[satPhaseBias->_updateInt]); 435 } 436 yaw = satPhaseBias->_yaw + satPhaseBias->_yawRate * dt; 400 437 ssr = true; 401 438 for (unsigned ii = 0; ii < satPhaseBias->_bias.size(); ii++) { … … 414 451 _model._windUp = station->windUp(_time, _prn, rSat, ssr, yaw, vSat) ; 415 452 453 // Relativistic effect due to earth gravity 454 // ---------------------------------------- 455 double a = rSat.NormFrobenius() + rRec.NormFrobenius(); 456 double b = (rSat - rRec).NormFrobenius(); 457 double gm = 3.986004418e14; // m3/s2 458 _model._rel = 2 * gm / t_CST::c / t_CST::c * log((a + b) / (a - b)); 416 459 417 460 // Tidal Correction 418 461 // ---------------- 419 _model._tide = -DotProduct(station->tideDspl(), rhoV) / _model._rho; 462 _model._tideEarth = -DotProduct(station->tideDsplEarth(), rhoV) / _model._rho; 463 _model._tideOcean = -DotProduct(station->tideDsplOcean(), rhoV) / _model._rho; 420 464 421 465 // Ionospheric Delay … … 429 473 } 430 474 } 475 431 476 if (vTecUsage && vTec) { 432 double stec = station->stec(vTec, _signalPropagationTime, rSat); 477 double stec = station->stec(vTec, _signalPropagationTime, rSat); 478 double f1GPS = t_CST::freq(t_frequency::G1, 0); 433 479 for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) { 434 t_frequency::type frqType = static_cast<t_frequency::type>(iFreq);435 _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(t_CST::freq(frqType, _channel), 2) * stec;436 }437 }438 439 // Relativistic effect due to earth gravity440 // ----------------------------------------441 // TODO442 443 // Ocean Loading444 // -------------445 // TODO480 if (OPT->_pseudoObsIono) { // DCMcodeBias, DCMphaseBias 481 // For scaling the slant ionospheric delays the trick is to be consistent with units! 482 // The conversion of TECU into meters requires the frequency of the signal. 483 // Hence, GPS L1 frequency is used for all systems. The same is true for mu_i in lcCoeff(). 484 _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(f1GPS, 2) * stec; 485 } 486 else { // PPP-RTK 487 t_frequency::type frqType = static_cast<t_frequency::type>(iFreq); 488 _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(t_CST::freq(frqType, _channel), 2) * stec; 489 } 490 } 491 } 446 492 447 493 // Set Model Set Flag … … 449 495 _model._set = true; 450 496 451 //printModel();497 printModel(); 452 498 453 499 return success; … … 457 503 //////////////////////////////////////////////////////////////////////////// 458 504 void t_pppSatObs::printModel() const { 459 460 LOG.setf(ios::fixed); 461 LOG << "MODEL for Satellite " << _prn.toString() << endl 462 << "RHO: " << setw(12) << setprecision(3) << _model._rho << endl 463 << "ELE: " << setw(12) << setprecision(3) << _model._eleSat * 180.0 / M_PI << endl 464 << "AZI: " << setw(12) << setprecision(3) << _model._azSat * 180.0 / M_PI << endl 465 << "SATCLK: " << setw(12) << setprecision(3) << _model._satClkM << endl 466 << "RECCLK: " << setw(12) << setprecision(3) << _model._recClkM << endl 467 << "SAGNAC: " << setw(12) << setprecision(3) << _model._sagnac << endl 468 << "ANTECC: " << setw(12) << setprecision(3) << _model._antEcc << endl 469 << "TROPO: " << setw(12) << setprecision(3) << _model._tropo << endl 470 << "WINDUP: " << setw(12) << setprecision(3) << _model._windUp << endl 471 << "TIDES: " << setw(12) << setprecision(3) << _model._tide << endl; 505 // TODO: cout should be LOG 506 cout.setf(ios::fixed); 507 cout << "\nMODEL for Satellite " << _prn.toString() << (isReference() ? " (Reference Satellite)" : "") << endl 508 << "======================= " << endl 509 << "PPP STRATEGY : " << OPT->_obsmodelTypeStr.at((int)OPT->_obsModelType).toLocal8Bit().constData() 510 << ((OPT->_pseudoObsIono) ? " with pseudo-observations for STEC" : "") << endl 511 << "RHO : " << setw(12) << setprecision(3) << _model._rho << endl 512 << "ELE : " << setw(12) << setprecision(3) << _model._eleSat * RHO_DEG << endl 513 << "AZI : " << setw(12) << setprecision(3) << _model._azSat * RHO_DEG << endl 514 << "SATCLK : " << setw(12) << setprecision(3) << _model._satClkM << endl 515 << "RECCLK : " << setw(12) << setprecision(3) << _model._recClkM << endl 516 << "SAGNAC : " << setw(12) << setprecision(3) << _model._sagnac << endl 517 << "ANTECC : " << setw(12) << setprecision(3) << _model._antEcc << endl 518 << "TROPO : " << setw(12) << setprecision(3) << _model._tropo << endl 519 << "WINDUP : " << setw(12) << setprecision(3) << _model._windUp << endl 520 << "REL : " << setw(12) << setprecision(3) << _model._rel << endl 521 << "EARTH TIDES : " << setw(12) << setprecision(3) << _model._tideEarth << endl 522 << "OCEAN TIDES : " << setw(12) << setprecision(3) << _model._tideOcean << endl 523 << endl 524 << "FREQUENCY DEPENDENT CORRECTIONS:" << endl 525 << "-------------------------------" << endl; 472 526 for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) { 473 527 if (_obs[iFreq]) { 474 528 string frqStr = t_frequency::toString(t_frequency::type(iFreq)); 475 529 if (_prn.system() == frqStr[0]) { 476 LOG << "PCO : " << frqStr << setw(12) << setprecision(3) << _model._antPCO[iFreq] << endl 477 << "BIAS CODE : " << frqStr << setw(12) << setprecision(3) << _model._codeBias[iFreq] << endl 478 << "BIAS PHASE : " << frqStr << setw(12) << setprecision(3) << _model._phaseBias[iFreq] << endl 479 << "IONO CODEDELAY: " << frqStr << setw(12) << setprecision(3) << _model._ionoCodeDelay[iFreq] << endl; 480 } 481 } 482 } 530 cout << "PCO : " << frqStr << setw(12) << setprecision(3) << _model._antPCO[iFreq] << endl 531 << "BIAS CODE : " << frqStr << setw(12) << setprecision(3) << _model._codeBias[iFreq] << endl 532 << "BIAS PHASE : " << frqStr << setw(12) << setprecision(3) << _model._phaseBias[iFreq] << endl 533 << "IONO CODEDELAY: " << frqStr << setw(12) << setprecision(3) << _model._ionoCodeDelay[iFreq]<< endl; 534 } 535 } 536 } 537 } 538 539 // 540 //////////////////////////////////////////////////////////////////////////// 541 void t_pppSatObs::printObsMinusComputed() const { 542 // TODO: cout should be LOG 543 cout.setf(ios::fixed); 544 cout << "\nOBS-COMP for Satellite " << _prn.toString() << (isReference() ? " (Reference Satellite)" : "") << endl 545 << "========================== " << endl; 483 546 for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) { 484 547 t_lc::type tLC = OPT->LCs(_prn.system())[ii]; 485 LOG << "OBS-CMP "<< t_lc::toString(tLC) << ": " << _prn.toString() << " "548 cout << "OBS-CMP " << setw(4) << t_lc::toString(tLC) << ": " << _prn.toString() << " " 486 549 << setw(12) << setprecision(3) << obsValue(tLC) << " " 487 550 << setw(12) << setprecision(3) << cmpValue(tLC) << " " 488 551 << setw(12) << setprecision(3) << obsValue(tLC) - cmpValue(tLC) << endl; 489 490 } 491 LOG << "OBS-CMP MW: " << _prn.toString() << " " 492 << setw(12) << setprecision(3) << obsValue(t_lc::MW) << " " 493 << setw(12) << setprecision(3) << cmpValue(t_lc::MW) << " " 494 << setw(12) << setprecision(3) << obsValue(t_lc::MW) - cmpValue(t_lc::MW) << endl; 495 } 552 } 553 } 554 496 555 497 556 // … … 504 563 //////////////////////////////////////////////////////////////////////////// 505 564 double t_pppSatObs::cmpValue(t_lc::type tLC) const { 506 507 if (!isValid(tLC)) { 508 return 0.0; 509 } 510 511 // Non-Dispersive Part 512 // ------------------- 513 double nonDisp = _model._rho + _model._recClkM - _model._satClkM 514 + _model._sagnac + _model._antEcc + _model._tropo 515 + _model._tide; 516 517 // Add Dispersive Part 518 // ------------------- 519 map<t_frequency::type, double> codeCoeff; 520 map<t_frequency::type, double> phaseCoeff; 521 lcCoeff(tLC, codeCoeff, phaseCoeff); 522 523 double dispPart = 0.0; 524 525 map<t_frequency::type, double>::const_iterator it; 526 for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) { 527 t_frequency::type tFreq = it->first; 528 dispPart += it->second * (_model._antPCO[tFreq] - _model._codeBias[tFreq] + 529 _model._ionoCodeDelay[tFreq]); 530 } 531 for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) { 532 t_frequency::type tFreq = it->first; 533 dispPart += it->second * (_model._antPCO[tFreq] - _model._phaseBias[tFreq] + 534 _model._windUp * t_CST::lambda(tFreq, _channel) - 535 _model._ionoCodeDelay[tFreq]); 536 } 537 538 return nonDisp + dispPart; 565 double cmpValue; 566 567 if (!isValid(tLC)) { 568 cmpValue = 0.0; 569 } 570 else if (tLC == t_lc::GIM) { 571 cmpValue = 0.0; 572 } 573 else { 574 // Non-Dispersive Part 575 // ------------------- 576 double nonDisp = _model._rho 577 + _model._recClkM - _model._satClkM 578 + _model._sagnac + _model._antEcc + _model._tropo 579 + _model._tideEarth + _model._tideOcean + _model._rel; 580 581 // Add Dispersive Part 582 // ------------------- 583 double dispPart = 0.0; 584 map<t_frequency::type, double> codeCoeff; 585 map<t_frequency::type, double> phaseCoeff; 586 map<t_frequency::type, double> ionoCoeff; 587 lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff); 588 map<t_frequency::type, double>::const_iterator it; 589 for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) { 590 t_frequency::type tFreq = it->first; 591 dispPart += it->second * (_model._antPCO[tFreq] - _model._codeBias[tFreq]); 592 if (OPT->PPPRTK) { 593 dispPart += it->second * (_model._ionoCodeDelay[tFreq]); 594 } 595 } 596 for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) { 597 t_frequency::type tFreq = it->first; 598 dispPart += it->second * (_model._antPCO[tFreq] - _model._phaseBias[tFreq] + 599 _model._windUp * t_CST::lambda(tFreq, _channel)); 600 if (OPT->PPPRTK) { 601 dispPart += it->second * (- _model._ionoCodeDelay[tFreq]); 602 } 603 } 604 cmpValue = nonDisp + dispPart; 605 } 606 607 return cmpValue; 539 608 } 540 609 … … 556 625 } 557 626 } 627 628 // 629 //////////////////////////////////////////////////////////////////////////// 630 void t_pppSatObs::setPseudoObsIono(t_frequency::type freq, double stecRefSat) { 631 _stecSat = _model._ionoCodeDelay[freq]; 632 _stecRefSat = stecRefSat; 633 }
Note:
See TracChangeset
for help on using the changeset viewer.