Changeset 4355 in ntrip for trunk/BNC/src/rinex/reqcanalyze.cpp


Ignore:
Timestamp:
Jun 24, 2012, 6:25:02 PM (12 years ago)
Author:
mervart
Message:
 
File:
1 edited

Legend:

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

    r4354 r4355  
    175175                                   .arg(obs.satNum, 2, 10, QChar('0'));
    176176
     177      t_satStat& satStat = _satStat[prn];
     178
     179      satStat.addObs(obs);
     180    }
     181
     182  } // while (_currEpo)
     183
     184  // Analyze the Multipath
     185  // ---------------------
     186  QVector<t_polarPoint*>* dataMP1 = new QVector<t_polarPoint*>;
     187  QVector<t_polarPoint*>* dataMP2 = new QVector<t_polarPoint*>;
     188
     189  QMapIterator<QString, t_satStat> it(_satStat);
     190  while (it.hasNext()) {
     191    it.next();
     192    QString          prn     = it.key();
     193    const t_satStat& satStat = it.value();
     194    analyzeMultipath(prn, satStat, xyz, dataMP1, dataMP2);
     195  }
     196
     197  emit displayGraph(dataMP1, dataMP2);
     198
     199  _log->flush();
     200}
     201
     202// 
     203////////////////////////////////////////////////////////////////////////////
     204void t_reqcAnalyze::t_satStat::addObs(const t_obs& obs) {
     205
     206  t_anaObs* newObs = new t_anaObs(obs.GPSWeek, obs.GPSWeeks);
     207  bool      okFlag = false;
     208
     209  // Compute the Multipath
     210  // ----------------------
     211  if (obs.l1() != 0.0 && obs.l2() != 0.0) {
     212    double f1 = t_CST::f1(obs.satSys, obs.slotNum);
     213    double f2 = t_CST::f2(obs.satSys, obs.slotNum);
     214
     215    double L1 = obs.l1() * t_CST::c / f1;
     216    double L2 = obs.l2() * t_CST::c / f2;
     217
     218    if (obs.p1() != 0.0) {
     219      newObs->_MP1 = obs.p1() - L1 - 2.0*f2*f2/(f1*f1-f2*f2) * (L1 - L2);
     220      okFlag = true;
     221    }
     222    if (obs.p2() != 0.0) {
     223      newObs->_MP2 = obs.p2() - L2 - 2.0*f1*f1/(f1*f1-f2*f2) * (L1 - L2);
     224      okFlag = true;
     225    }
     226  }
     227
     228  // Remember the Observation
     229  // ------------------------
     230  if (okFlag) {
     231    anaObs << newObs;
     232  }
     233  else {
     234    delete newObs;
     235  }
     236}
     237
     238// 
     239////////////////////////////////////////////////////////////////////////////
     240void t_reqcAnalyze::analyzeMultipath(const QString& prn,
     241                                     const t_satStat& satStat,
     242                                     const ColumnVector& xyz,
     243                                     QVector<t_polarPoint*>* dataMP1,
     244                                     QVector<t_polarPoint*>* dataMP2) {
     245
     246  const int    LENGTH = 10;  // number of epochs in one chunk
     247  const double SLIP   = 5.0; // cycle-slip threshold
     248
     249  int numEpo = satStat.anaObs.size();
     250
     251  for (int chunkStart = 0; chunkStart + LENGTH < numEpo; chunkStart += LENGTH) {
     252
     253    // Compute Mean
     254    // ------------
     255    bool   slipFlag = false;
     256    double mean1    = 0.0;
     257    double mean2    = 0.0;
     258
     259    for (int ii = 0; ii < LENGTH; ii++) {
     260      int iEpo = chunkStart + ii;
     261      const t_anaObs* anaObs = satStat.anaObs[iEpo];
     262      mean1 += anaObs->_MP1;
     263      mean2 += anaObs->_MP2;
     264 
     265      // Check Slip
     266      // ----------
     267      if (ii > 0) {
     268        double diff1 = anaObs->_MP1 - satStat.anaObs[iEpo-1]->_MP1;
     269        double diff2 = anaObs->_MP2 - satStat.anaObs[iEpo-1]->_MP2;
     270        if (fabs(diff1) > SLIP || fabs(diff2) > SLIP) {
     271          slipFlag = true;
     272          break;
     273        }
     274      }
     275    }
     276
     277    if (slipFlag) {
     278      continue;
     279    }
     280
     281    mean1 /= LENGTH;
     282    mean2 /= LENGTH;
     283
     284    // Compute Standard Deviation
     285    // --------------------------
     286    double stddev1 = 0.0;
     287    double stddev2 = 0.0;
     288    for (int ii = 0; ii < LENGTH; ii++) {
     289      int iEpo = chunkStart + ii;
     290      const t_anaObs* anaObs = satStat.anaObs[iEpo];
     291      double diff1 = anaObs->_MP1 - mean1;
     292      double diff2 = anaObs->_MP2 - mean2;
     293      stddev1 += diff1 * diff1;
     294      stddev2 += diff2 * diff2;
     295    }
     296    double MP1 = sqrt(stddev1 / (LENGTH-1));
     297    double MP2 = sqrt(stddev2 / (LENGTH-1));
     298
     299    const t_anaObs* anaObs0 = satStat.anaObs[chunkStart];
     300
     301    // Compute the Azimuth and Zenith Distance
     302    // ---------------------------------------
     303    double az  = 0.0;
     304    double zen = 0.0;
     305    if (xyz.size()) {
    177306      t_eph* eph = 0;
    178307      for (int ie = 0; ie < _ephs.size(); ie++) {
     
    182311        }
    183312      }
    184 
    185       t_satStat& satStat = _satStat[prn];
    186 
    187       satStat.addObs(obs, eph, xyz);
    188     }
    189 
    190   } // while (_currEpo)
    191 
    192   // Analyze the Multipath
    193   // ---------------------
    194   QVector<t_polarPoint*>* dataMP1 = new QVector<t_polarPoint*>;
    195   QVector<t_polarPoint*>* dataMP2 = new QVector<t_polarPoint*>;
    196 
    197   QMapIterator<QString, t_satStat> it(_satStat);
    198   while (it.hasNext()) {
    199     it.next();
    200     QString          prn     = it.key();
    201     const t_satStat& satStat = it.value();
    202     analyzeMultipath(prn, satStat, dataMP1, dataMP2);
    203   }
    204 
    205   emit displayGraph(dataMP1, dataMP2);
    206 
    207   _log->flush();
    208 }
    209 
    210 // 
    211 ////////////////////////////////////////////////////////////////////////////
    212 void t_reqcAnalyze::t_satStat::addObs(const t_obs& obs, const t_eph* eph,
    213                                       const ColumnVector& xyz) {
    214 
    215   t_anaObs* newObs = new t_anaObs(obs);
    216   bool      okFlag = false;
    217 
    218   // Compute the Multipath
    219   // ----------------------
    220   if (obs.l1() != 0.0 && obs.l2() != 0.0) {
    221     double f1 = t_CST::f1(obs.satSys, obs.slotNum);
    222     double f2 = t_CST::f2(obs.satSys, obs.slotNum);
    223 
    224     double L1 = obs.l1() * t_CST::c / f1;
    225     double L2 = obs.l2() * t_CST::c / f2;
    226 
    227     if (obs.p1() != 0.0) {
    228       newObs->MP1 = obs.p1() - L1 - 2.0*f2*f2/(f1*f1-f2*f2) * (L1 - L2);
    229       okFlag = true;
    230     }
    231     if (obs.p2() != 0.0) {
    232       newObs->MP2 = obs.p2() - L2 - 2.0*f1*f1/(f1*f1-f2*f2) * (L1 - L2);
    233       okFlag = true;
    234     }
    235   }
    236 
    237   // Remember the Observation
    238   // ------------------------
    239   if (okFlag) {
    240     anaObs << newObs;
    241   }
    242   else {
    243     delete newObs;
    244     return;
    245   }
    246 
    247   // Compute the Azimuth and Zenith Distance
    248   // ---------------------------------------
    249   if (eph && xyz.size()) {
    250     double xSat, ySat, zSat, clkSat;
    251     eph->position(obs.GPSWeek, obs.GPSWeeks, xSat, ySat, zSat, clkSat);
    252 
    253     double rho, eleSat, azSat;
    254     topos(xyz(1), xyz(2), xyz(3), xSat, ySat, zSat, rho, eleSat, azSat);
    255 
    256     newObs->az  = azSat * 180.0/M_PI;
    257     newObs->zen = 90.0 - eleSat * 180.0/M_PI;
    258   }
    259 }
    260 
    261 // 
    262 ////////////////////////////////////////////////////////////////////////////
    263 void t_reqcAnalyze::analyzeMultipath(const QString& prn,
    264                                      const t_satStat& satStat,
    265                                      QVector<t_polarPoint*>* dataMP1,
    266                                      QVector<t_polarPoint*>* dataMP2) {
    267 
    268   const int    LENGTH = 10;  // number of epochs in one chunk
    269   const double SLIP   = 5.0; // cycle-slip threshold
    270 
    271   int numEpo = satStat.anaObs.size();
    272 
    273   for (int chunkStart = 0; chunkStart + LENGTH < numEpo; chunkStart += LENGTH) {
    274 
    275     // Compute Mean
    276     // ------------
    277     bool   slipFlag = false;
    278     double mean1    = 0.0;
    279     double mean2    = 0.0;
    280 
    281     for (int ii = 0; ii < LENGTH; ii++) {
    282       int iEpo = chunkStart + ii;
    283       const t_anaObs* anaObs = satStat.anaObs[iEpo];
    284       mean1 += anaObs->MP1;
    285       mean2 += anaObs->MP2;
    286  
    287       // Check Slip
    288       // ----------
    289       if (ii > 0) {
    290         double diff1 = anaObs->MP1 - satStat.anaObs[iEpo-1]->MP1;
    291         double diff2 = anaObs->MP2 - satStat.anaObs[iEpo-1]->MP2;
    292         if (fabs(diff1) > SLIP || fabs(diff2) > SLIP) {
    293           slipFlag = true;
    294           break;
    295         }
     313     
     314      if (eph) {
     315        double xSat, ySat, zSat, clkSat;
     316        eph->position(anaObs0->_GPSWeek, anaObs0->_GPSWeeks,
     317                      xSat, ySat, zSat, clkSat);
     318     
     319        double rho, eleSat, azSat;
     320        topos(xyz(1), xyz(2), xyz(3), xSat, ySat, zSat, rho, eleSat, azSat);
     321     
     322        az  = azSat * 180.0/M_PI;
     323        zen = 90.0 - eleSat * 180.0/M_PI;
    296324      }
    297325    }
    298 
    299     if (slipFlag) {
    300       continue;
    301     }
    302 
    303     mean1 /= LENGTH;
    304     mean2 /= LENGTH;
    305 
    306     // Compute Standard Deviation
    307     // --------------------------
    308     double stddev1 = 0.0;
    309     double stddev2 = 0.0;
    310     for (int ii = 0; ii < LENGTH; ii++) {
    311       int iEpo = chunkStart + ii;
    312       const t_anaObs* anaObs = satStat.anaObs[iEpo];
    313       double diff1 = anaObs->MP1 - mean1;
    314       double diff2 = anaObs->MP2 - mean2;
    315       stddev1 += diff1 * diff1;
    316       stddev2 += diff2 * diff2;
    317     }
    318     double MP1 = sqrt(stddev1 / (LENGTH-1));
    319     double MP2 = sqrt(stddev2 / (LENGTH-1));
    320326
    321327    // Add new Point
    322328    // -------------
    323     const t_anaObs* anaObs = satStat.anaObs[chunkStart];
    324     (*dataMP1) << (new t_polarPoint(anaObs->az, anaObs->zen, MP1));
    325     (*dataMP2) << (new t_polarPoint(anaObs->az, anaObs->zen, MP2));
     329    (*dataMP1) << (new t_polarPoint(az, zen, MP1));
     330    (*dataMP2) << (new t_polarPoint(az, zen, MP2));
    326331
    327332    _log->setRealNumberNotation(QTextStream::FixedNotation);
    328333
    329334    _log->setRealNumberPrecision(2);
    330     *_log << "MP1 " << prn << " " << anaObs->az << " " << anaObs->zen << " ";
     335    *_log << "MP1 " << prn << " " << az << " " << zen << " ";
    331336    _log->setRealNumberPrecision(3);
    332337    *_log << MP1 << endl;
    333338
    334339    _log->setRealNumberPrecision(2);
    335     *_log << "MP2 " << prn << " " << anaObs->az << " " << anaObs->zen << " ";
     340    *_log << "MP2 " << prn << " " << az << " " << zen << " ";
    336341    _log->setRealNumberPrecision(3);
    337342    *_log << MP2 << endl;
Note: See TracChangeset for help on using the changeset viewer.