source: ntrip/trunk/BNC/src/orbComp/sp3Comp.cpp@ 11076

Last change on this file since 11076 was 11076, checked in by stuerze, 14 hours ago

addded SIS to the SP3 comparison
File size: 18.1 KB
Line 
1// Part of BNC, a utility for retrieving decoding and
2// converting GNSS data streams from NTRIP broadcasters.
3//
4// Copyright (C) 2007
5// German Federal Agency for Cartography and Geodesy (BKG)
6// http://bkg.bund.de
7// Czech Technical University Prague, Department of Geodesy
8// http://www.fsv.cvut.cz
9//
10// Email: euref-ip@bkg.bund.de
11//
12// This program is free software; you can redistribute it and/or
13// modify it under the terms of the GNU General Public License
14// as published by the Free Software Foundation, version 2.
15//
16// This program is distributed in the hope that it will be useful,
17// but WITHOUT ANY WARRANTY; without even the implied warranty of
18// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19// GNU General Public License for more details.
20//
21// You should have received a copy of the GNU General Public License
22// along with this program; if not, write to the Free Software
23// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
24
25/* -------------------------------------------------------------------------
26 * BKG NTRIP Client
27 * -------------------------------------------------------------------------
28 *
29 * Class: t_sp3Comp
30 *
31 * Purpose: Compare SP3 Files
32 *
33 * Author: L. Mervart
34 *
35 * Created: 24-Nov-2014
36 *
37 * Changes:
38 *
39 * -----------------------------------------------------------------------*/
40
41#include <iostream>
42#include <iomanip>
43#include "sp3Comp.h"
44#include "bnccore.h"
45#include "bncsettings.h"
46#include "bncutils.h"
47#include "bncsp3.h"
48
49using namespace std;
50
51// Constructor
52////////////////////////////////////////////////////////////////////////////
53t_sp3Comp::t_sp3Comp(QObject* parent) : QThread(parent) {
54
55 bncSettings settings;
56 _sp3FileNames = settings.value("sp3CompFile").toString().split(QRegExp("[ ,]"), Qt::SkipEmptyParts);
57 for (int ii = 0; ii < _sp3FileNames.size(); ii++) {
58 expandEnvVar(_sp3FileNames[ii]);
59 }
60 _logFileName = settings.value("sp3CompOutLogFile").toString(); expandEnvVar(_logFileName);
61 _logFile = 0;
62 _log = 0;
63
64 _excludeSats = settings.value("sp3CompExclude").toString().split(QRegExp("[ ,]"), Qt::SkipEmptyParts);
65
66 _summaryOnly = (Qt::CheckState(settings.value("sp3CompSummaryOnly").toInt()) == Qt::Checked);
67}
68
69// Destructor
70////////////////////////////////////////////////////////////////////////////
71t_sp3Comp::~t_sp3Comp() {
72 delete _log;
73 delete _logFile;
74}
75
76//
77////////////////////////////////////////////////////////////////////////////
78void t_sp3Comp::run() {
79
80 // Open Log File
81 // -------------
82 _logFile = new QFile(_logFileName);
83 if (_logFile->open(QIODevice::WriteOnly | QIODevice::Text)) {
84 _log = new QTextStream();
85 _log->setDevice(_logFile);
86 }
87 if (!_log) {
88 *_log << "ERROR: SP3Comp requires logfile specification" << "\n";
89 goto exit;
90 }
91
92 for (int ii = 0; ii < _sp3FileNames.size(); ii++) {
93 *_log << "! SP3 File " << ii+1 << ": " << _sp3FileNames[ii] << "\n";
94 }
95 if (_sp3FileNames.size() != 2) {
96 *_log << "ERROR: sp3Comp requires two input SP3 files" << "\n";
97 goto end;
98 }
99
100 try {
101 ostringstream msg;
102 compare(msg);
103 *_log << msg.str().c_str();
104 }
105 catch (const string& error) {
106 *_log << "ERROR: " << error.c_str() << "\n";
107 }
108 catch (const char* error) {
109 *_log << "ERROR: " << error << "\n";
110 }
111 catch (Exception& exc) {
112 *_log << "ERROR: " << exc.what() << "\n";
113 }
114 catch (std::exception& exc) {
115 *_log << "ERROR: " << exc.what() << "\n";
116 }
117 catch (QString error) {
118 *_log << "ERROR: " << error << "\n";
119 }
120 catch (...) {
121 *_log << "ERROR: " << "unknown exception" << "\n";
122 }
123
124 // Exit (thread)
125 // -------------
126 end:
127 _log->flush();
128 exit:
129 // do nothing if no logfile available
130
131 if (BNC_CORE->mode() != t_bncCore::interactive) {
132 qApp->exit(7);
133 msleep(100); //sleep 0.1 sec
134 }
135 else {
136 emit finished();
137 deleteLater();
138 }
139}
140
141// Satellite Index in clkSats set
142////////////////////////////////////////////////////////////////////////////////
143int t_sp3Comp::satIndex(const set<t_prn>& clkSats, const t_prn& prn) const {
144 int ret = 0;
145 for (set<t_prn>::const_iterator it = clkSats.begin(); it != clkSats.end(); it++) {
146 if ( *it == prn) {
147 return ret;
148 }
149 ++ret;
150 }
151 return -1;
152}
153
154// Estimate Clock Offsets
155////////////////////////////////////////////////////////////////////////////////
156void t_sp3Comp::processClocks(const set<t_prn>& clkSats, const vector<t_epoch*>& epochsIn,
157 map<string, t_stat>& stat) const {
158
159 if (clkSats.size() == 0) {
160 return;
161 }
162
163 vector<t_epoch*> epochs;
164 for (unsigned ii = 0; ii < epochsIn.size(); ii++) {
165 unsigned numSatOK = 0;
166 std::map<t_prn, double>::const_iterator it;
167 for (it = epochsIn[ii]->_dc.begin(); it != epochsIn[ii]->_dc.end(); it++) {
168 if (satIndex(clkSats, it->first) != -1) {
169 numSatOK += 1;
170 }
171 }
172 if (numSatOK > 0) {
173 epochs.push_back(epochsIn[ii]);
174 }
175 }
176 if (epochs.size() == 0) {
177 return;
178 }
179
180 int nPar = epochs.size() + clkSats.size();
181 SymmetricMatrix NN(nPar); NN = 0.0;
182 ColumnVector bb(nPar); bb = 0.0;
183
184 // Create Matrix A'A and vector b
185 // ------------------------------
186 for (unsigned ie = 0; ie < epochs.size(); ie++) {
187 const map<t_prn, double>& dc = epochs[ie]->_dc;
188 Matrix AA(dc.size(), nPar); AA = 0.0;
189 ColumnVector ll(dc.size()); ll = 0.0;
190 map<t_prn, double>::const_iterator it;
191 int ii = -1;
192 for (it = dc.begin(); it != dc.end(); it++) {
193 const t_prn& prn = it->first;
194 if (satIndex(clkSats, prn) != -1) {
195 ++ii;
196 int index = epochs.size() + satIndex(clkSats, prn);
197 AA[ii][ie] = 1.0; // epoch-specfic offset (common for all satellites)
198 AA[ii][index] = 1.0; // satellite-specific offset (common for all epochs)
199 ll[ii] = it->second;
200 }
201 }
202 SymmetricMatrix dN; dN << AA.t() * AA;
203 NN += dN;
204 bb += AA.t() * ll;
205 }
206
207 // Regularize NN
208 // -------------
209 RowVector HH(nPar);
210 HH.columns(1, epochs.size()) = 0.0;
211 HH.columns(epochs.size()+1, nPar) = 1.0;
212 SymmetricMatrix dN; dN << HH.t() * HH;
213 NN += dN;
214
215 // Estimate Parameters
216 // -------------------
217 ColumnVector xx = NN.i() * bb;
218
219 // Compute clock residuals
220 // -----------------------
221 for (unsigned ie = 0; ie < epochs.size(); ie++) {
222 map<t_prn, double>& dc = epochs[ie]->_dc;
223 map<t_prn, double>& dcRed = epochs[ie]->_dcRed;
224 const map<t_prn, ColumnVector>& dr = epochs[ie]->_dr;
225 for (map<t_prn, double>::iterator it = dc.begin(); it != dc.end(); it++) {
226 const t_prn& prn = it->first;
227 stringstream all; all << prn.system() << 99;
228 if (satIndex(clkSats, prn) != -1) {
229 int index = epochs.size() + satIndex(clkSats, prn);
230 dc[prn] = it->second - xx[ie] - xx[index];
231 stat[prn.toString()]._offset = xx[index];
232 if (dr.find(prn) != dr.end()){
233 // signal-in-space difference: radial minus clock, only the
234 // epoch-wise clock offset (datum) removed, i.e. including the
235 // satellite-specific clock offset
236 epochs[ie]->_sis[prn] = dr.find(prn)->second[0] - (dc[prn] + xx[index]);
237 dcRed[prn] = dc[prn] - dr.find(prn)->second[0]; // clock minus radial component
238 stat[prn.toString()]._dcRedMean += dcRed[prn];
239 stat[all.str() ]._dcRedMean += dcRed[prn];
240 stat[prn.toString()]._nc += 1;
241 stat[all.str() ]._nc += 1;
242 }
243 }
244 }
245 }
246
247 // Compute Clock Mean
248 // ------------------
249 for (map<string, t_stat>::iterator it = stat.begin(); it != stat.end(); it++) {
250 t_stat& stat = it->second;
251 if (stat._nc > 0) {
252 stat._dcRedMean = stat._dcRedMean / stat._nc;
253 }
254 }
255}
256
257// Main Routine
258////////////////////////////////////////////////////////////////////////////////
259void t_sp3Comp::compare(ostringstream& out) const {
260
261 // Synchronize reading of two sp3 files
262 // ------------------------------------
263 bncSP3 in1(_sp3FileNames[0]); in1.nextEpoch();
264 bncSP3 in2(_sp3FileNames[1]); in2.nextEpoch();
265
266 vector<t_epoch*> epochs;
267 while (in1.currEpoch() && in2.currEpoch()) {
268 bncTime t1 = in1.currEpoch()->_tt;
269 bncTime t2 = in2.currEpoch()->_tt;
270 if (t1 < t2) {
271 in1.nextEpoch();
272 }
273 else if (t1 > t2) {
274 in2.nextEpoch();
275 }
276 else if (t1 == t2) {
277 t_epoch* epo = new t_epoch; epo->_tt = t1;
278 bool epochOK = false;
279 for (int i1 = 0; i1 < in1.currEpoch()->_sp3Sat.size(); i1++) {
280 bncSP3::t_sp3Sat* sat1 = in1.currEpoch()->_sp3Sat[i1];
281 for (int i2 = 0; i2 < in2.currEpoch()->_sp3Sat.size(); i2++) {
282 bncSP3::t_sp3Sat* sat2 = in2.currEpoch()->_sp3Sat[i2];
283 if (sat1->_prn == sat2->_prn &&
284 sat1->_clkValid && sat2->_clkValid) {
285 epochOK = true;
286 epo->_dr[sat1->_prn] = sat1->_xyz - sat2->_xyz;
287 epo->_xyz[sat1->_prn] = sat1->_xyz;
288 epo->_dc[sat1->_prn] = sat1->_clk - sat2->_clk;
289 }
290 }
291 }
292 if (epochOK) {
293 epochs.push_back(epo);
294 }
295 else {
296 delete epo;
297 }
298 in1.nextEpoch();
299 in2.nextEpoch();
300 }
301 }
302
303 // Transform xyz into radial, along-track, and out-of-plane
304 // --------------------------------------------------------
305 if (epochs.size() < 2) {
306 throw "t_sp3Comp: not enough common epochs";
307 }
308
309 set<t_prn> clkSatsAll;
310
311 for (unsigned ie = 0; ie < epochs.size(); ie++) {
312 t_epoch* epoch = epochs[ie];
313 t_epoch* epoch2 = 0;
314 if (ie == 0) {
315 epoch2 = epochs[ie+1];
316 }
317 else {
318 epoch2 = epochs[ie-1];
319 }
320 double dt = epoch->_tt - epoch2->_tt;
321 map<t_prn, ColumnVector>& dr = epoch->_dr;
322 map<t_prn, ColumnVector>& xyz = epoch->_xyz;
323 map<t_prn, ColumnVector>& xyz2 = epoch2->_xyz;
324 QList<t_prn> satEraseList;
325 for (map<t_prn, ColumnVector>::const_iterator it = dr.begin(); it != dr.end(); it++) {
326 const t_prn& prn = it->first;
327 if (xyz2.find(prn) != xyz2.end()) {
328 const ColumnVector dx = dr[prn];
329 const ColumnVector& x1 = xyz[prn];
330 const ColumnVector& x2 = xyz2[prn];
331 ColumnVector vel = (x1 - x2) / dt;
332 XYZ_to_RSW(x1, vel, dx, dr[prn]);
333 if (epoch->_dc.find(prn) != epoch->_dc.end()) {
334 clkSatsAll.insert(prn);
335 }
336 }
337 else {
338 satEraseList << it->first;
339 }
340 }
341 for (QList<t_prn>::const_iterator it = satEraseList.begin(); it != satEraseList.end(); it++) {
342 epoch->_dc.erase(*it);
343 epoch->_dr.erase(*it);
344 }
345 }
346
347 map<string, t_stat> stat;
348
349 // Estimate Clock Offsets
350 // ----------------------
351 string systems;
352 for (set<t_prn>::const_iterator it = clkSatsAll.begin(); it != clkSatsAll.end(); it++) {
353 if (systems.find(it->system()) == string::npos) {
354 systems += it->system();
355 }
356 }
357 for (unsigned iSys = 0; iSys < systems.size(); iSys++) {
358 char system = systems[iSys];
359 set<t_prn> clkSats;
360 for (set<t_prn>::const_iterator it = clkSatsAll.begin(); it != clkSatsAll.end(); it++) {
361 if (it->system() == system && !excludeSat(*it)) {
362 clkSats.insert(*it);
363 }
364 }
365 processClocks(clkSats, epochs, stat);
366 }
367
368 // Print epoch-wise Clock Residuals
369 // --------------------------------
370 out.setf(ios::fixed);
371 if (!_summaryOnly) {
372 out << "!\n! Clock residuals and orbit differences in [m]\n"
373 "! ----------------------------------------------------------------------------\n";
374 out << "!\n! Epoch PRN radial along out clk clkRed iPRN SIS"
375 "\n! -------------------------------------------------------------------------------------\n";
376 }
377 for (unsigned ii = 0; ii < epochs.size(); ii++) {
378 const t_epoch* epo = epochs[ii];
379 const map<t_prn, ColumnVector>& dr = epochs[ii]->_dr;
380 const map<t_prn, double>& dc = epochs[ii]->_dc;
381 const map<t_prn, double>& dcRed = epochs[ii]->_dcRed;
382 const map<t_prn, double>& sis = epochs[ii]->_sis;
383 for (map<t_prn, ColumnVector>::const_iterator it = dr.begin(); it != dr.end(); it++) {
384 const t_prn& prn = it->first;
385 stringstream all; all << prn.system() << 99;
386 if (!excludeSat(prn)) {
387 const ColumnVector& rao = it->second;
388 if (!_summaryOnly) {
389 out << setprecision(6) << string(epo->_tt) << ' ' << prn.toString() << ' '
390 << setw(7) << setprecision(4) << rao[0] << ' '
391 << setw(7) << setprecision(4) << rao[1] << ' '
392 << setw(7) << setprecision(4) << rao[2] << " ";
393 }
394 stat[prn.toString()]._rao += SP(rao, rao); // Schur product
395 stat[all.str() ]._rao += SP(rao, rao);
396 stat[prn.toString()]._nr += 1;
397 stat[all.str() ]._nr += 1;
398 if (dc.find(prn) != dc.end() && dcRed.find(prn) != dc.end()) {
399 double clkRes = dc.find(prn)->second;
400 double clkResRed = dcRed.find(prn)->second;
401 if (!_summaryOnly) {
402 out << setw(7) << setprecision(4) << clkRes << ' '
403 << setw(7) << setprecision(4) << clkResRed;
404 }
405 stat[prn.toString()]._dcRMS += clkRes * clkRes;
406 stat[all.str() ]._dcRMS += clkRes * clkRes;
407 stat[prn.toString()]._dcRedRMS += clkResRed * clkResRed;
408 stat[all.str() ]._dcRedRMS += clkResRed * clkResRed;
409 stat[prn.toString()]._dcRedSig += (clkResRed - stat[prn.toString()]._dcRedMean) *
410 (clkResRed - stat[prn.toString()]._dcRedMean);
411 stat[all.str() ]._dcRedSig += (clkResRed - stat[all.str() ]._dcRedMean) *
412 (clkResRed - stat[all.str() ]._dcRedMean);
413 }
414 else {
415 if (!_summaryOnly) {
416 out << " . . ";
417 }
418 }
419 if (!_summaryOnly) {
420 out << " " << setw(2) << int(prn);
421 }
422 if (sis.find(prn) != sis.end()) {
423 double sisVal = sis.find(prn)->second;
424 stat[prn.toString()]._sisSum += sisVal;
425 stat[all.str() ]._sisSum += sisVal;
426 stat[prn.toString()]._sisSqr += sisVal * sisVal;
427 stat[all.str() ]._sisSqr += sisVal * sisVal;
428 stat[prn.toString()]._ns += 1;
429 stat[all.str() ]._ns += 1;
430 if (!_summaryOnly) {
431 out << " " << setw(7) << setprecision(4) << sisVal;
432 }
433 }
434 if (!_summaryOnly) {
435 out << "\n";
436 }
437 }
438 }
439 delete epo;
440 }
441
442 // Print Summary
443 // -------------
444 out << "!\n! Summary";
445 out << "\n! -----------------------------------------------------------------------------------------------------------------\n";
446 out << "!\n! PRN radialRMS alongRMS outRMS 3DRMS nOrb clkRMS clkRedRMS clkRedSig nClk Offset SISRMS SISSig"
447 "\n! [mm] [mm] [mm] [mm] [-] [ns] [ns] [ns] [-] [ns] [mm] [mm]"
448 "\n! -------------------------------------------------------------------------------------------------------------------------------------\n";
449 // SIS standard deviation of a system: pooled about the satellite means
450 map<char, double> sisVarSys;
451 for (map<string, t_stat>::const_iterator it = stat.begin(); it != stat.end(); it++) {
452 const t_stat& st = it->second;
453 stringstream all; all << it->first[0] << 99;
454 if (it->first != all.str() && st._ns > 0) {
455 sisVarSys[it->first[0]] += st._sisSqr - st._sisSum * st._sisSum / st._ns;
456 }
457 }
458 for (map<string, t_stat>::iterator it = stat.begin(); it != stat.end(); it++) {
459 const string& prn = it->first;
460 t_stat& stat = it->second;
461 char sys = prn[0];
462 stringstream all; all << sys << 99;
463 if (stat._nr > 0) {
464 stat._rao[0] = sqrt(stat._rao[0] / stat._nr);
465 stat._rao[1] = sqrt(stat._rao[1] / stat._nr);
466 stat._rao[2] = sqrt(stat._rao[2] / stat._nr);
467 stat._rao3DRMS = stat._rao.NormFrobenius();
468 // orbit values in millimeter
469 if (prn != all.str()) {
470 (_summaryOnly) ? out << " " << prn << ' ':
471 out << "! " << prn << ' ';
472 }
473 else {
474 QString sys; sys = all.str()[0];
475 (_summaryOnly) ? out << " " << QString("%1").arg(sys,3,' ').toStdString().c_str() << " ":
476 out << "! " << QString("%1").arg(sys,3,' ').toStdString().c_str() << " ";
477 }
478 out << setw(10) << setprecision(1) << stat._rao[0] * 1e3 << ' '
479 << setw(10) << setprecision(1) << stat._rao[1] * 1e3 << ' '
480 << setw(10) << setprecision(1) << stat._rao[2] * 1e3 << ' '
481 << setw(10) << setprecision(1) << stat._rao3DRMS * 1e3 << ' '
482 << setw( 7) << stat._nr << " ";
483 // clock values in nano seconds
484 if (stat._nc > 0) {
485 stat._dcRMS = sqrt(stat._dcRMS / stat._nc);
486 stat._dcRedRMS = sqrt(stat._dcRedRMS / stat._nc);
487 stat._dcRedSig = sqrt(stat._dcRedSig / stat._nc);
488 out << setw(10) << setprecision(2) << stat._dcRMS / t_CST::c * 1e9 << ' '
489 << setw(10) << setprecision(2) << stat._dcRedRMS / t_CST::c * 1e9 << ' '
490 << setw(10) << setprecision(2) << stat._dcRedSig / t_CST::c * 1e9 << ' '
491 << setw( 9) << stat._nc << " ";
492 if (prn != "G99" && prn != "R99" && prn != "E99" && prn != "C98" && prn != "C99") {
493 out << setw( 9) << setprecision(2) << stat._offset / t_CST::c * 1e9;
494 }
495 else {
496 out << setw( 9) << " ";
497 }
498 }
499 if (stat._ns > 0) {
500 double var = (prn == all.str()) ? sisVarSys[sys] / stat._ns
501 : stat._sisSqr / stat._ns - pow(stat._sisSum / stat._ns, 2);
502 out << setw(11) << setprecision(1) << sqrt(stat._sisSqr / stat._ns) * 1e3 << ' '
503 << setw( 9) << setprecision(1) << sqrt(max(0.0, var)) * 1e3;
504 }
505 out << "\n";
506 }
507 }
508}
509
510//
511////////////////////////////////////////////////////////////////////////////
512bool t_sp3Comp::excludeSat(const t_prn& prn) const {
513 QStringListIterator it(_excludeSats);
514 while (it.hasNext()) {
515 string prnStr = it.next().toLatin1().data();
516 if (prnStr == prn.toString() || // prn
517 prnStr == prn.toString().substr(0,1)) { // sys
518 return true;
519 }
520 }
521 return false;
522}
523
Note: See TracBrowser for help on using the repository browser.