- Timestamp:
- Aug 30, 2012, 11:51:41 AM (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/BNC/src/rinex/reqcanalyze.cpp
r4584 r4590 229 229 230 230 if (obs.satSys == 'R') { 231 continue;// TODO: set channel number231 // TODO: set channel number 232 232 } 233 233 … … 375 375 t_allObs& allObs = _allObsMap[prn]; 376 376 377 bncTime currTime; 378 bncTime prevTime; 379 377 380 for (int chunkStart = 0; chunkStart + numEpo < allObs._oneObsVec.size(); 378 381 chunkStart += chunkStep) { 379 382 380 bncTime firstEpoch;381 383 382 384 // Compute Mean 383 385 // ------------ 384 bool slipFlag = false; 385 double mean1 = 0.0; 386 double mean2 = 0.0; 387 double SNR1 = 0.0; 388 double SNR2 = 0.0; 386 bncTime chunkStartTime; 387 bool availL1 = false; 388 bool availL2 = false; 389 bool missL1 = false; 390 bool missL2 = false; 391 bool slipFlag = false; 392 double meanMP1 = 0.0; 393 double meanMP2 = 0.0; 394 double minSNR1 = 0.0; 395 double minSNR2 = 0.0; 396 double aziDeg = 0.0; 397 double zenDeg = 0.0; 389 398 390 399 for (int ii = 0; ii < numEpo; ii++) { … … 392 401 const t_oneObs* oneObs = allObs._oneObsVec[iEpo]; 393 402 403 currTime.set(oneObs->_GPSWeek, oneObs->_GPSWeeks); 404 405 // Compute the Azimuth and Zenith Distance 406 // --------------------------------------- 394 407 if (ii == 0) { 395 bncTime epoTime(oneObs->_GPSWeek, oneObs->_GPSWeeks); 396 _availDataMap[prn]._epoL1 << epoTime.mjddec(); 397 } 398 399 mean1 += oneObs->_MP1; 400 mean2 += oneObs->_MP2; 401 402 if ( oneObs->_SNR1 > 0 && (SNR1 == 0 || SNR1 > oneObs->_SNR1) ) { 403 SNR1 = oneObs->_SNR1; 404 } 405 if ( oneObs->_SNR2 > 0 && (SNR2 == 0 || SNR2 > oneObs->_SNR2) ) { 406 SNR2 = oneObs->_SNR2; 407 } 408 408 chunkStartTime = currTime; 409 410 if (xyz.size()) { 411 t_eph* eph = 0; 412 for (int ie = 0; ie < _ephs.size(); ie++) { 413 if (_ephs[ie]->prn() == prn) { 414 eph = _ephs[ie]; 415 break; 416 } 417 } 418 419 if (eph) { 420 double xSat, ySat, zSat, clkSat; 421 eph->position(oneObs->_GPSWeek, oneObs->_GPSWeeks, 422 xSat, ySat, zSat, clkSat); 423 424 double rho, eleSat, azSat; 425 topos(xyz(1), xyz(2), xyz(3), xSat, ySat, zSat, rho, eleSat, azSat); 426 427 aziDeg = azSat * 180.0/M_PI; 428 zenDeg = 90.0 - eleSat * 180.0/M_PI; 429 } 430 } 431 } 432 433 // Check Interval 434 // -------------- 435 double dt = 0.0; 436 if (prevTime.valid()) { 437 dt = currTime - prevTime; 438 } 439 if (dt != obsInterval) { 440 missL1 = true; 441 missL2 = true; 442 } 443 prevTime = currTime; 444 445 // Check L1 and L2 availability 446 // ---------------------------- 447 if (oneObs->_hasL1) { 448 availL1 = true; 449 } 450 else { 451 missL1 = true; 452 } 453 if (oneObs->_hasL2) { 454 availL2 = true; 455 } 456 else { 457 missL2 = true; 458 } 459 460 // Check Minimal Signal-to-Noise Ratio 461 // ----------------------------------- 462 if ( oneObs->_SNR1 > 0 && (minSNR1 == 0 || minSNR1 > oneObs->_SNR1) ) { 463 minSNR1 = oneObs->_SNR1; 464 } 465 if ( oneObs->_SNR2 > 0 && (minSNR2 == 0 || minSNR2 > oneObs->_SNR2) ) { 466 minSNR2 = oneObs->_SNR2; 467 } 468 409 469 // Check Slip 410 470 // ---------- … … 414 474 if (fabs(diff1) > SLIPTRESH || fabs(diff2) > SLIPTRESH) { 415 475 slipFlag = true; 416 break;417 476 } 418 477 } 419 } 420 421 if (slipFlag) { 422 continue; 423 } 424 425 mean1 /= numEpo; 426 mean2 /= numEpo; 427 428 // Compute Standard Deviation 429 // -------------------------- 430 double stddev1 = 0.0; 431 double stddev2 = 0.0; 432 for (int ii = 0; ii < numEpo; ii++) { 433 int iEpo = chunkStart + ii; 434 const t_oneObs* oneObs = allObs._oneObsVec[iEpo]; 435 double diff1 = oneObs->_MP1 - mean1; 436 double diff2 = oneObs->_MP2 - mean2; 437 stddev1 += diff1 * diff1; 438 stddev2 += diff2 * diff2; 439 } 440 double MP1 = sqrt(stddev1 / (numEpo-1)); 441 double MP2 = sqrt(stddev2 / (numEpo-1)); 442 443 const t_oneObs* oneObs0 = allObs._oneObsVec[chunkStart]; 444 445 // Compute the Azimuth and Zenith Distance 446 // --------------------------------------- 447 double az = 0.0; 448 double zen = 0.0; 449 if (xyz.size()) { 450 t_eph* eph = 0; 451 for (int ie = 0; ie < _ephs.size(); ie++) { 452 if (_ephs[ie]->prn() == prn) { 453 eph = _ephs[ie]; 454 break; 455 } 456 } 457 458 if (eph) { 459 double xSat, ySat, zSat, clkSat; 460 eph->position(oneObs0->_GPSWeek, oneObs0->_GPSWeeks, 461 xSat, ySat, zSat, clkSat); 462 463 double rho, eleSat, azSat; 464 topos(xyz(1), xyz(2), xyz(3), xSat, ySat, zSat, rho, eleSat, azSat); 465 466 az = azSat * 180.0/M_PI; 467 zen = 90.0 - eleSat * 180.0/M_PI; 468 } 469 } 470 471 // Add new Point 472 // ------------- 473 (*dataMP1) << (new t_polarPoint(az, zen, MP1)); 474 (*dataMP2) << (new t_polarPoint(az, zen, MP2)); 475 (*dataSNR1) << (new t_polarPoint(az, zen, SNR1)); 476 (*dataSNR2) << (new t_polarPoint(az, zen, SNR2)); 477 478 if (_log) { 479 _log->setRealNumberNotation(QTextStream::FixedNotation); 480 481 _log->setRealNumberPrecision(2); 482 *_log << "MP1 " << prn << " " << az << " " << zen << " "; 483 _log->setRealNumberPrecision(3); 484 *_log << MP1 << endl; 485 486 _log->setRealNumberPrecision(2); 487 *_log << "MP2 " << prn << " " << az << " " << zen << " "; 488 _log->setRealNumberPrecision(3); 489 *_log << MP2 << endl; 490 491 _log->flush(); 478 479 meanMP1 += oneObs->_MP1; 480 meanMP2 += oneObs->_MP2; 481 } 482 483 // Availability Plot Data 484 // ---------------------- 485 if (availL1) { 486 _availDataMap[prn]._epoL1 << chunkStartTime.mjddec(); 487 } 488 489 // Signal-to-Noise Ration Plot Data 490 // -------------------------------- 491 (*dataSNR1) << (new t_polarPoint(aziDeg, zenDeg, minSNR1)); 492 (*dataSNR2) << (new t_polarPoint(aziDeg, zenDeg, minSNR2)); 493 494 // Compute the Multipath 495 // --------------------- 496 if (!slipFlag) { 497 meanMP1 /= numEpo; 498 meanMP2 /= numEpo; 499 double MP1 = 0.0; 500 double MP2 = 0.0; 501 for (int ii = 0; ii < numEpo; ii++) { 502 int iEpo = chunkStart + ii; 503 const t_oneObs* oneObs = allObs._oneObsVec[iEpo]; 504 double diff1 = oneObs->_MP1 - meanMP1; 505 double diff2 = oneObs->_MP2 - meanMP2; 506 MP1 += diff1 * diff1; 507 MP2 += diff2 * diff2; 508 } 509 MP1 = sqrt(MP1 / (numEpo-1)); 510 MP2 = sqrt(MP2 / (numEpo-1)); 511 (*dataMP1) << (new t_polarPoint(aziDeg, zenDeg, MP1)); 512 (*dataMP2) << (new t_polarPoint(aziDeg, zenDeg, MP2)); 492 513 } 493 514 }
Note:
See TracChangeset
for help on using the changeset viewer.