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

Last change on this file since 10092 was 10092, checked in by stuerze, 15 months ago

units are changed in sp3 comparison summary

File size: 15.6 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://www.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("[ ,]"), QString::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("[ ,]"), QString::SkipEmptyParts);
65}
66
67// Destructor
68////////////////////////////////////////////////////////////////////////////
69t_sp3Comp::~t_sp3Comp() {
70 delete _log;
71 delete _logFile;
72}
73
74//
75////////////////////////////////////////////////////////////////////////////
76void t_sp3Comp::run() {
77
78 // Open Log File
79 // -------------
80 _logFile = new QFile(_logFileName);
81 if (_logFile->open(QIODevice::WriteOnly | QIODevice::Text)) {
82 _log = new QTextStream();
83 _log->setDevice(_logFile);
84 }
85 if (!_log) {
86 *_log << "ERROR: SP3Comp requires logfile specification" << "\n";
87 goto exit;
88 }
89
90 for (int ii = 0; ii < _sp3FileNames.size(); ii++) {
91 *_log << "! SP3 File " << ii+1 << ": " << _sp3FileNames[ii] << "\n";
92 }
93 if (_sp3FileNames.size() != 2) {
94 *_log << "ERROR: sp3Comp requires two input SP3 files" << "\n";
95 goto end;
96 }
97
98 try {
99 ostringstream msg;
100 compare(msg);
101 *_log << msg.str().c_str();
102 }
103 catch (const string& error) {
104 *_log << "ERROR: " << error.c_str() << "\n";
105 }
106 catch (const char* error) {
107 *_log << "ERROR: " << error << "\n";
108 }
109 catch (Exception& exc) {
110 *_log << "ERROR: " << exc.what() << "\n";
111 }
112 catch (std::exception& exc) {
113 *_log << "ERROR: " << exc.what() << "\n";
114 }
115 catch (QString error) {
116 *_log << "ERROR: " << error << "\n";
117 }
118 catch (...) {
119 *_log << "ERROR: " << "unknown exception" << "\n";
120 }
121
122 // Exit (thread)
123 // -------------
124 end:
125 _log->flush();
126 exit:
127 // do nothing if no logfile available
128
129 if (BNC_CORE->mode() != t_bncCore::interactive) {
130 qApp->exit(7);
131 msleep(100); //sleep 0.1 sec
132 }
133 else {
134 emit finished();
135 deleteLater();
136 }
137}
138
139// Satellite Index in clkSats set
140////////////////////////////////////////////////////////////////////////////////
141int t_sp3Comp::satIndex(const set<t_prn>& clkSats, const t_prn& prn) const {
142 int ret = 0;
143 for (set<t_prn>::const_iterator it = clkSats.begin(); it != clkSats.end(); it++) {
144 if ( *it == prn) {
145 return ret;
146 }
147 ++ret;
148 }
149 return -1;
150}
151
152// Estimate Clock Offsets
153////////////////////////////////////////////////////////////////////////////////
154void t_sp3Comp::processClocks(const set<t_prn>& clkSats, const vector<t_epoch*>& epochsIn,
155 map<string, t_stat>& stat) const {
156
157 if (clkSats.size() == 0) {
158 return;
159 }
160
161 vector<t_epoch*> epochs;
162 for (unsigned ii = 0; ii < epochsIn.size(); ii++) {
163 unsigned numSatOK = 0;
164 std::map<t_prn, double>::const_iterator it;
165 for (it = epochsIn[ii]->_dc.begin(); it != epochsIn[ii]->_dc.end(); it++) {
166 if (satIndex(clkSats, it->first) != -1) {
167 numSatOK += 1;
168 }
169 }
170 if (numSatOK > 0) {
171 epochs.push_back(epochsIn[ii]);
172 }
173 }
174 if (epochs.size() == 0) {
175 return;
176 }
177
178 int nPar = epochs.size() + clkSats.size();
179 SymmetricMatrix NN(nPar); NN = 0.0;
180 ColumnVector bb(nPar); bb = 0.0;
181
182 // Create Matrix A'A and vector b
183 // ------------------------------
184 for (unsigned ie = 0; ie < epochs.size(); ie++) {
185 const map<t_prn, double>& dc = epochs[ie]->_dc;
186 Matrix AA(dc.size(), nPar); AA = 0.0;
187 ColumnVector ll(dc.size()); ll = 0.0;
188 map<t_prn, double>::const_iterator it;
189 int ii = -1;
190 for (it = dc.begin(); it != dc.end(); it++) {
191 const t_prn& prn = it->first;
192 if (satIndex(clkSats, prn) != -1) {
193 ++ii;
194 int index = epochs.size() + satIndex(clkSats, prn);
195 AA[ii][ie] = 1.0; // epoch-specfic offset (common for all satellites)
196 AA[ii][index] = 1.0; // satellite-specific offset (common for all epochs)
197 ll[ii] = it->second;
198 }
199 }
200 SymmetricMatrix dN; dN << AA.t() * AA;
201 NN += dN;
202 bb += AA.t() * ll;
203 }
204
205 // Regularize NN
206 // -------------
207 RowVector HH(nPar);
208 HH.columns(1, epochs.size()) = 0.0;
209 HH.columns(epochs.size()+1, nPar) = 1.0;
210 SymmetricMatrix dN; dN << HH.t() * HH;
211 NN += dN;
212
213 // Estimate Parameters
214 // -------------------
215 ColumnVector xx = NN.i() * bb;
216
217 // Compute clock residuals
218 // -----------------------
219 const string all = "ZZZ";
220 for (unsigned ie = 0; ie < epochs.size(); ie++) {
221 map<t_prn, double>& dc = epochs[ie]->_dc;
222 map<t_prn, double>& dcRed = epochs[ie]->_dcRed;
223 const map<t_prn, ColumnVector>& dr = epochs[ie]->_dr;
224 for (map<t_prn, double>::iterator it = dc.begin(); it != dc.end(); it++) {
225 const t_prn& prn = it->first;
226 if (satIndex(clkSats, prn) != -1) {
227 int index = epochs.size() + satIndex(clkSats, prn);
228 dc[prn] = it->second - xx[ie] - xx[index];
229 stat[prn.toString()]._offset = xx[index];
230 if (dr.find(prn) != dr.end()){
231 dcRed[prn] = dc[prn] - dr.find(prn)->second[0]; // clock minus radial component
232 stat[prn.toString()]._dcRedMean += dcRed[prn];
233 stat[all ]._dcRedMean += dcRed[prn];
234 stat[prn.toString()]._nc += 1;
235 stat[all ]._nc += 1;
236 }
237 }
238 }
239 }
240
241 // Compute Clock Mean
242 // ------------------
243 for (map<string, t_stat>::iterator it = stat.begin(); it != stat.end(); it++) {
244 t_stat& stat = it->second;
245 if (stat._nc > 0) {
246 stat._dcRedMean = stat._dcRedMean / stat._nc;
247 }
248 }
249}
250
251// Main Routine
252////////////////////////////////////////////////////////////////////////////////
253void t_sp3Comp::compare(ostringstream& out) const {
254
255 // Synchronize reading of two sp3 files
256 // ------------------------------------
257 bncSP3 in1(_sp3FileNames[0]); in1.nextEpoch();
258 bncSP3 in2(_sp3FileNames[1]); in2.nextEpoch();
259
260 vector<t_epoch*> epochs;
261 while (in1.currEpoch() && in2.currEpoch()) {
262 bncTime t1 = in1.currEpoch()->_tt;
263 bncTime t2 = in2.currEpoch()->_tt;
264 if (t1 < t2) {
265 in1.nextEpoch();
266 }
267 else if (t1 > t2) {
268 in2.nextEpoch();
269 }
270 else if (t1 == t2) {
271 t_epoch* epo = new t_epoch; epo->_tt = t1;
272 bool epochOK = false;
273 for (int i1 = 0; i1 < in1.currEpoch()->_sp3Sat.size(); i1++) {
274 bncSP3::t_sp3Sat* sat1 = in1.currEpoch()->_sp3Sat[i1];
275 for (int i2 = 0; i2 < in2.currEpoch()->_sp3Sat.size(); i2++) {
276 bncSP3::t_sp3Sat* sat2 = in2.currEpoch()->_sp3Sat[i2];
277 if (sat1->_prn == sat2->_prn) {
278 epochOK = true;
279 epo->_dr[sat1->_prn] = sat1->_xyz - sat2->_xyz;
280 epo->_xyz[sat1->_prn] = sat1->_xyz;
281 if (sat1->_clkValid && sat2->_clkValid) {
282 epo->_dc[sat1->_prn] = sat1->_clk - sat2->_clk;
283 }
284 }
285 }
286 }
287 if (epochOK) {
288 epochs.push_back(epo);
289 }
290 else {
291 delete epo;
292 }
293 in1.nextEpoch();
294 in2.nextEpoch();
295 }
296 }
297
298 // Transform xyz into radial, along-track, and out-of-plane
299 // --------------------------------------------------------
300 if (epochs.size() < 2) {
301 throw "t_sp3Comp: not enough common epochs";
302 }
303
304 set<t_prn> clkSatsAll;
305
306 for (unsigned ie = 0; ie < epochs.size(); ie++) {
307 t_epoch* epoch = epochs[ie];
308 t_epoch* epoch2 = 0;
309 if (ie == 0) {
310 epoch2 = epochs[ie+1];
311 }
312 else {
313 epoch2 = epochs[ie-1];
314 }
315 double dt = epoch->_tt - epoch2->_tt;
316 map<t_prn, ColumnVector>& dr = epoch->_dr;
317 map<t_prn, ColumnVector>& xyz = epoch->_xyz;
318 map<t_prn, ColumnVector>& xyz2 = epoch2->_xyz;
319 QList<t_prn> satEraseList;
320 for (map<t_prn, ColumnVector>::const_iterator it = dr.begin(); it != dr.end(); it++) {
321 const t_prn& prn = it->first;
322 if (xyz2.find(prn) != xyz2.end()) {
323 const ColumnVector dx = dr[prn];
324 const ColumnVector& x1 = xyz[prn];
325 const ColumnVector& x2 = xyz2[prn];
326 ColumnVector vel = (x1 - x2) / dt;
327 XYZ_to_RSW(x1, vel, dx, dr[prn]);
328 if (epoch->_dc.find(prn) != epoch->_dc.end()) {
329 clkSatsAll.insert(prn);
330 }
331 }
332 else {
333 satEraseList << it->first;
334 }
335 }
336 for (QList<t_prn>::const_iterator it = satEraseList.begin(); it != satEraseList.end(); it++) {
337 epoch->_dc.erase(*it);
338 epoch->_dr.erase(*it);
339 }
340 }
341
342 map<string, t_stat> stat;
343
344 // Estimate Clock Offsets
345 // ----------------------
346 string systems;
347 for (set<t_prn>::const_iterator it = clkSatsAll.begin(); it != clkSatsAll.end(); it++) {
348 if (systems.find(it->system()) == string::npos) {
349 systems += it->system();
350 }
351 }
352 for (unsigned iSys = 0; iSys < systems.size(); iSys++) {
353 char system = systems[iSys];
354 set<t_prn> clkSats;
355 for (set<t_prn>::const_iterator it = clkSatsAll.begin(); it != clkSatsAll.end(); it++) {
356 if (it->system() == system && !excludeSat(*it)) {
357 clkSats.insert(*it);
358 }
359 }
360 processClocks(clkSats, epochs, stat);
361 }
362
363 // Print epoch-wise Clock Residuals
364 // --------------------------------
365 const string all = "ZZZ";
366
367 out.setf(ios::fixed);
368 out << "!\n! Clock residuals and orbit differences in [m]\n"
369 "! ----------------------------------------------------------------------------\n";
370 out << "!\n! Epoch PRN radial along out clk clkRed iPRN"
371 "\n! ----------------------------------------------------------------------------\n";
372 for (unsigned ii = 0; ii < epochs.size(); ii++) {
373 const t_epoch* epo = epochs[ii];
374 const map<t_prn, ColumnVector>& dr = epochs[ii]->_dr;
375 const map<t_prn, double>& dc = epochs[ii]->_dc;
376 const map<t_prn, double>& dcRed = epochs[ii]->_dcRed;
377 for (map<t_prn, ColumnVector>::const_iterator it = dr.begin(); it != dr.end(); it++) {
378 const t_prn& prn = it->first;
379 if (!excludeSat(prn)) {
380 const ColumnVector& rao = it->second;
381 out << setprecision(6) << string(epo->_tt) << ' ' << prn.toString() << ' '
382 << setw(7) << setprecision(4) << rao[0] << ' '
383 << setw(7) << setprecision(4) << rao[1] << ' '
384 << setw(7) << setprecision(4) << rao[2] << " ";
385 stat[prn.toString()]._rao += SP(rao, rao); // Schur product
386 stat[all ]._rao += SP(rao, rao);
387 stat[prn.toString()]._nr += 1;
388 stat[all ]._nr += 1;
389 if (dc.find(prn) != dc.end() && dcRed.find(prn) != dc.end()) {
390 double clkRes = dc.find(prn)->second;
391 double clkResRed = dcRed.find(prn)->second;
392 out << setw(7) << setprecision(4) << clkRes << ' '
393 << setw(7) << setprecision(4) << clkResRed;
394 stat[prn.toString()]._dcRMS += clkRes * clkRes;
395 stat[all ]._dcRMS += clkRes * clkRes;
396 stat[prn.toString()]._dcRedRMS += clkResRed * clkResRed;
397 stat[all ]._dcRedRMS += clkResRed * clkResRed;
398 stat[prn.toString()]._dcRedSig += (clkResRed - stat[prn.toString()]._dcRedMean) *
399 (clkResRed - stat[prn.toString()]._dcRedMean);
400 stat[all ]._dcRedSig += (clkResRed - stat[all ]._dcRedMean) *
401 (clkResRed - stat[all ]._dcRedMean);
402 }
403 else {
404 out << " . . ";
405 }
406 out << " " << setw(2) << int(prn) << "\n";
407 }
408 }
409 delete epo;
410 }
411
412 // Print Summary
413 // -------------
414 out << "!\n! Summary";
415 out << "\n! -----------------------------------------------------------------------------------------------------------------\n";
416 out << "!\n! PRN radialRMS alongRMS outRMS 3DRMS nOrb clkRMS clkRedRMS clkRedSig nClk Offset "
417 "\n! [mm] [mm] [mm] [mm] [-] [ns] [ns] [ns] [-] [ns] "
418 "\n! -----------------------------------------------------------------------------------------------------------------\n";
419 for (map<string, t_stat>::iterator it = stat.begin(); it != stat.end(); it++) {
420 const string& prn = it->first;
421 t_stat& stat = it->second;
422 if (stat._nr > 0) {
423 stat._rao[0] = sqrt(stat._rao[0] / stat._nr);
424 stat._rao[1] = sqrt(stat._rao[1] / stat._nr);
425 stat._rao[2] = sqrt(stat._rao[2] / stat._nr);
426 stat._rao3DRMS = stat._rao.NormFrobenius();
427 if (prn == all) {
428 out << "!\n! Total ";
429 }
430 else {
431 out << "! " << prn << ' ';
432 }
433 // orbit values in millimeters
434 out << setw(10) << setprecision(1) << stat._rao[0] * 1e3 << ' '
435 << setw(10) << setprecision(1) << stat._rao[1] * 1e3 << ' '
436 << setw(10) << setprecision(1) << stat._rao[2] * 1e3 << ' '
437 << setw(10) << setprecision(1) << stat._rao3DRMS * 1e3 << ' '
438 << setw( 7) << stat._nr << " ";
439 if (stat._nc > 0) {
440 stat._dcRMS = sqrt(stat._dcRMS / stat._nc);
441 stat._dcRedRMS = sqrt(stat._dcRedRMS / stat._nc);
442 stat._dcRedSig = sqrt(stat._dcRedSig / stat._nc);
443 // clock values in nano seconds
444 out << setw(10) << setprecision(2) << stat._dcRMS / t_CST::c * 1e9 << ' '
445 << setw(10) << setprecision(2) << stat._dcRedRMS / t_CST::c * 1e9 << ' '
446 << setw(10) << setprecision(2) << stat._dcRedSig / t_CST::c * 1e9 << ' '
447 << setw( 9) << stat._nc << " ";
448 if (prn != all) {
449 out << setw(9) << setprecision(2) << stat._offset / t_CST::c * 1e9;
450 }
451 }
452 else {
453 out << " . . ";
454 }
455 out << "\n";
456 }
457 }
458}
459
460//
461////////////////////////////////////////////////////////////////////////////
462bool t_sp3Comp::excludeSat(const t_prn& prn) const {
463 QStringListIterator it(_excludeSats);
464 while (it.hasNext()) {
465 string prnStr = it.next().toLatin1().data();
466 if (prnStr == prn.toString() || prnStr == prn.toString().substr(0,1)) {
467 return true;
468 }
469 }
470 return false;
471}
472
Note: See TracBrowser for help on using the repository browser.