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

Last change on this file since 6630 was 6612, checked in by weber, 10 years ago

Output format modified

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