Changeset 4590 in ntrip for trunk/BNC


Ignore:
Timestamp:
Aug 30, 2012, 11:51:41 AM (12 years ago)
Author:
mervart
Message:
 
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/BNC/src/rinex/reqcanalyze.cpp

    r4584 r4590  
    229229 
    230230        if (obs.satSys == 'R') {
    231           continue; // TODO: set channel number
     231          // TODO: set channel number
    232232        }
    233233 
     
    375375  t_allObs& allObs = _allObsMap[prn];
    376376
     377  bncTime currTime;
     378  bncTime prevTime;
     379
    377380  for (int chunkStart = 0; chunkStart + numEpo < allObs._oneObsVec.size();
    378381       chunkStart += chunkStep) {
    379382
    380     bncTime firstEpoch;
    381383
    382384    // Compute Mean
    383385    // ------------
    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;
    389398
    390399    for (int ii = 0; ii < numEpo; ii++) {
     
    392401      const t_oneObs* oneObs = allObs._oneObsVec[iEpo];
    393402
     403      currTime.set(oneObs->_GPSWeek, oneObs->_GPSWeeks);
     404
     405      // Compute the Azimuth and Zenith Distance
     406      // ---------------------------------------
    394407      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
    409469      // Check Slip
    410470      // ----------
     
    414474        if (fabs(diff1) > SLIPTRESH || fabs(diff2) > SLIPTRESH) {
    415475          slipFlag = true;
    416           break;
    417476        }
    418477      }
    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));
    492513    }
    493514  }
Note: See TracChangeset for help on using the changeset viewer.