Changeset 4355 in ntrip for trunk/BNC/src/rinex/reqcanalyze.cpp
- Timestamp:
- Jun 24, 2012, 6:25:02 PM (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/BNC/src/rinex/reqcanalyze.cpp
r4354 r4355 175 175 .arg(obs.satNum, 2, 10, QChar('0')); 176 176 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 //////////////////////////////////////////////////////////////////////////// 204 void 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 //////////////////////////////////////////////////////////////////////////// 240 void 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()) { 177 306 t_eph* eph = 0; 178 307 for (int ie = 0; ie < _ephs.size(); ie++) { … … 182 311 } 183 312 } 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; 296 324 } 297 325 } 298 299 if (slipFlag) {300 continue;301 }302 303 mean1 /= LENGTH;304 mean2 /= LENGTH;305 306 // Compute Standard Deviation307 // --------------------------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));320 326 321 327 // Add new Point 322 328 // ------------- 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)); 326 331 327 332 _log->setRealNumberNotation(QTextStream::FixedNotation); 328 333 329 334 _log->setRealNumberPrecision(2); 330 *_log << "MP1 " << prn << " " << a naObs->az << " " << anaObs->zen << " ";335 *_log << "MP1 " << prn << " " << az << " " << zen << " "; 331 336 _log->setRealNumberPrecision(3); 332 337 *_log << MP1 << endl; 333 338 334 339 _log->setRealNumberPrecision(2); 335 *_log << "MP2 " << prn << " " << a naObs->az << " " << anaObs->zen << " ";340 *_log << "MP2 " << prn << " " << az << " " << zen << " "; 336 341 _log->setRealNumberPrecision(3); 337 342 *_log << MP2 << endl;
Note:
See TracChangeset
for help on using the changeset viewer.