Index: /trunk/BNC/src/bncutils.cpp =================================================================== --- /trunk/BNC/src/bncutils.cpp (revision 4337) +++ /trunk/BNC/src/bncutils.cpp (revision 4338) @@ -408,2 +408,33 @@ return ok ? 0 : 1; } + +// Topocentrical Distance and Elevation +//////////////////////////////////////////////////////////////////////////// +void topos(double xRec, double yRec, double zRec, + double xSat, double ySat, double zSat, + double& rho, double& eleSat, double& azSat) { + + double dx[3]; + dx[0] = xSat-xRec; + dx[1] = ySat-yRec; + dx[2] = zSat-zRec; + + rho = sqrt( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] ); + + double xyzRec[3]; + xyzRec[0] = xRec; + xyzRec[1] = yRec; + xyzRec[2] = zRec; + + double Ell[3]; + double neu[3]; + xyz2ell(xyzRec, Ell); + xyz2neu(Ell, dx, neu); + + eleSat = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / rho ); + if (neu[2] < 0) { + eleSat *= -1.0; + } + + azSat = atan2(neu[1], neu[0]); +} Index: /trunk/BNC/src/bncutils.h =================================================================== --- /trunk/BNC/src/bncutils.h (revision 4337) +++ /trunk/BNC/src/bncutils.h (revision 4338) @@ -76,3 +76,7 @@ int readDbl(const QString& str, int pos, int len, double& value); +void topos(double xRec, double yRec, double zRec, + double xSat, double ySat, double zSat, + double& rho, double& eleSat, double& azSat); + #endif Index: /trunk/BNC/src/rinex/reqcanalyze.cpp =================================================================== --- /trunk/BNC/src/rinex/reqcanalyze.cpp (revision 4337) +++ /trunk/BNC/src/rinex/reqcanalyze.cpp (revision 4338) @@ -64,4 +64,6 @@ _currEpo = 0; + _dataMP1 = 0; + _dataMP2 = 0; connect(this, SIGNAL(displayGraph()), this, SLOT(slotDisplayGraph())); @@ -86,45 +88,11 @@ if (((bncApp*) qApp)->mode() == bncApp::interactive) { - QVector* data1 = new QVector; - - //// beg test - { - const QwtInterval zenithInterval(0.0, 90.0); - const QwtInterval azimuthInterval(0.0, 360.0 ); - const int numPoints = 1000; - const double stepA = 4 * azimuthInterval.width() / numPoints; - const double stepR = zenithInterval.width() / numPoints; - for (int ii = 0; ii < numPoints; ii++) { - double aa = azimuthInterval.minValue() + ii * stepA; - double rr = zenithInterval.minValue() + ii * stepR; - double vv = static_cast(ii) / numPoints; - (*data1) << (new t_polarPoint(aa, rr, vv)); - } - } - //// end test - t_polarPlot* plotMP1 = new t_polarPlot(0); - plotMP1->addCurve(data1); - - QVector* data2 = new QVector; - - //// beg test - { - const QwtInterval zenithInterval(0.0, 90.0); - const QwtInterval azimuthInterval(0.0, 360.0 ); - const int numPoints = 1000; - const double stepA = 4 * azimuthInterval.width() / numPoints; - const double stepR = zenithInterval.width() / numPoints; - for (int ii = 0; ii < numPoints; ii++) { - double aa = azimuthInterval.minValue() + ii * stepA; - double rr = zenithInterval.minValue() + ii * stepR; - double vv = static_cast(ii) / numPoints; - (*data2) << (new t_polarPoint(aa, rr, vv)); - } - } - //// end test + plotMP1->addCurve(_dataMP1); + _dataMP1 = 0; t_polarPlot* plotMP2 = new t_polarPlot(0); - plotMP2->addCurve(data2); + plotMP2->addCurve(_dataMP2); + _dataMP2 = 0; QVector plots; @@ -201,6 +169,15 @@ .arg(obs.satNum, 2, 10, QChar('0')); + t_eph* eph = 0; + for (int ie = 0; ie < _ephs.size(); ie++) { + if (_ephs[ie]->prn() == prn) { + eph = _ephs[ie]; + break; + } + } + t_satStat& satStat = _satStat[prn]; - satStat.addObs(obs); + + satStat.addObs(eph, obs); } @@ -211,4 +188,7 @@ _log->setRealNumberNotation(QTextStream::FixedNotation); _log->setRealNumberPrecision(2); + + delete _dataMP1; _dataMP1 = new QVector; + delete _dataMP2; _dataMP2 = new QVector; QMapIterator it(_satStat); @@ -218,23 +198,34 @@ const t_satStat& satStat = it.value(); - for (int im = 1; im <= 2; im++) { - const QVector& MP = (im == 1) ? satStat.MP1 : satStat.MP2; - if (MP.size() > 1) { - double mean = 0.0; - for (int ii = 0; ii < MP.size(); ii++) { - mean += MP[ii]; - } - mean /= MP.size(); - double stddev = 0.0; - for (int ii = 0; ii < MP.size(); ii++) { - double diff = MP[ii] - mean; - stddev += diff * diff; - } - double multipath = sqrt(stddev / (MP.size()-1)); - - *_log << "MP" << im << " " << prn << " " << multipath << endl; - } - } - + int numVal = satStat.anaObs.size(); + if (numVal > 1) { + double mean1 = 0.0; + double mean2 = 0.0; + for (int ii = 0; ii < numVal; ii++) { + const t_anaObs* anaObs = satStat.anaObs[ii]; + mean1 += anaObs->MP1; + mean2 += anaObs->MP2; + } + mean1 /= numVal; + mean2 /= numVal; + double stddev1 = 0.0; + double stddev2 = 0.0; + for (int ii = 0; ii < numVal; ii++) { + const t_anaObs* anaObs = satStat.anaObs[ii]; + double diff1 = anaObs->MP1 - mean1; + double diff2 = anaObs->MP2 - mean2; + stddev1 += diff1 * diff1; + stddev2 += diff2 * diff2; + //// beg test + (*_dataMP1) << (new t_polarPoint(anaObs->az, anaObs->zen, 0.5)); + (*_dataMP2) << (new t_polarPoint(anaObs->az, anaObs->zen, 1.0)); + //// end test + } + double MP1 = sqrt(stddev1 / (numVal-1)); + double MP2 = sqrt(stddev2 / (numVal-1)); + + *_log << "MP1 " << prn << " " << MP1 << endl; + *_log << "MP2 " << prn << " " << MP2 << endl; + } } @@ -246,10 +237,8 @@ // //////////////////////////////////////////////////////////////////////////// -void t_reqcAnalyze::t_satStat::addObs(const t_obs& obs) { - if (currObs) { - delete prevObs; - prevObs = currObs; - } - currObs = new t_anaObs(obs); +void t_reqcAnalyze::t_satStat::addObs(const t_eph* eph, const t_obs& obs) { + + t_anaObs* newObs = new t_anaObs(obs); + anaObs << newObs; // Compute the Multipath @@ -263,11 +252,16 @@ if (obs.p1() != 0.0) { - currObs->M1 = obs.p1() - L1 - 2.0*f2*f2/(f1*f1-f2*f2) * (L1 - L2); - MP1 << currObs->M1; + newObs->MP1 = obs.p1() - L1 - 2.0*f2*f2/(f1*f1-f2*f2) * (L1 - L2); } if (obs.p2() != 0.0) { - currObs->M2 = obs.p2() - L2 - 2.0*f1*f1/(f1*f1-f2*f2) * (L1 - L2); - MP2 << currObs->M2; - } - } -} + newObs->MP2 = obs.p2() - L2 - 2.0*f1*f1/(f1*f1-f2*f2) * (L1 - L2); + } + } + + // Compute the Azimuth and Zenith Distance + // --------------------------------------- + if (eph) { + double xx, yy, zz, cc; + eph->position(obs.GPSWeek, obs.GPSWeeks, xx, yy, zz, cc); + } +} Index: /trunk/BNC/src/rinex/reqcanalyze.h =================================================================== --- /trunk/BNC/src/rinex/reqcanalyze.h (revision 4337) +++ /trunk/BNC/src/rinex/reqcanalyze.h (revision 4338) @@ -32,4 +32,6 @@ #include "RTCM3/ephemeris.h" +class t_polarPoint; + class t_reqcAnalyze : public QThread { Q_OBJECT @@ -54,29 +56,23 @@ class t_anaObs { public: - t_anaObs(const t_obs& obsIn) { - obs = obsIn; - M1 = 0.0; - M2 = 0.0; - } + t_anaObs(const t_obs& obsIn) : + obs(obsIn), az(0.0), zen(0.0), MP1(0.0), MP2(0.0) {} t_obs obs; - double M1; - double M2; + double az; + double zen; + double MP1; + double MP2; }; class t_satStat { public: - t_satStat() { - currObs = 0; - prevObs = 0; + t_satStat() {} + ~t_satStat() { + for (int ii = 0; ii < anaObs.size(); ii++) { + delete anaObs[ii]; + } } - ~t_satStat() { - delete currObs; - delete prevObs; - } - void addObs(const t_obs& obs); - QVector MP1; - QVector MP2; - t_anaObs* currObs; - t_anaObs* prevObs; + void addObs(const t_eph* eph, const t_obs& obs); + QVector anaObs; }; @@ -92,4 +88,6 @@ QMap _satStat; t_rnxObsFile::t_rnxEpo* _currEpo; + QVector* _dataMP1; + QVector* _dataMP2; };