source: ntrip/trunk/BNC/bncpppclient.cpp@ 2710

Last change on this file since 2710 was 2710, checked in by mervart, 13 years ago
File size: 14.9 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: bncPPPclient
30 *
31 * Purpose: Precise Point Positioning
32 *
33 * Author: L. Mervart
34 *
35 * Created: 21-Nov-2009
36 *
37 * Changes:
38 *
39 * -----------------------------------------------------------------------*/
40
41#include <newmatio.h>
42#include <iomanip>
43#include <sstream>
44
45#include "bncpppclient.h"
46#include "bncapp.h"
47#include "bncutils.h"
48#include "bncconst.h"
49#include "bncmodel.h"
50#include "bncsettings.h"
51
52extern "C" {
53#include "clock_orbit_rtcm.h"
54}
55
56using namespace std;
57
58// Constructor
59////////////////////////////////////////////////////////////////////////////
60bncPPPclient::bncPPPclient(QByteArray staID) {
61
62 bncSettings settings;
63
64 if ( Qt::CheckState(settings.value("pppGLONASS").toInt()) == Qt::Checked) {
65 _useGlonass = true;
66 }
67 else {
68 _useGlonass = false;
69 }
70
71 if (settings.value("pppSPP").toString() == "PPP") {
72 _pppMode = true;
73 }
74 else {
75 _pppMode = false;
76 }
77
78 connect(this, SIGNAL(newMessage(QByteArray,bool)),
79 ((bncApp*)qApp), SLOT(slotMessage(const QByteArray,bool)));
80
81 connect(((bncApp*)qApp), SIGNAL(newEphGPS(gpsephemeris)),
82 this, SLOT(slotNewEphGPS(gpsephemeris)));
83
84 connect(((bncApp*)qApp), SIGNAL(newEphGlonass(glonassephemeris)),
85 this, SLOT(slotNewEphGlonass(glonassephemeris)));
86
87 connect(((bncApp*)qApp), SIGNAL(newCorrections(QList<QString>)),
88 this, SLOT(slotNewCorrections(QList<QString>)));
89
90 _staID = staID;
91 _epoData = 0;
92 _model = new bncModel(staID);
93 connect(_model, SIGNAL(newNMEAstr(QByteArray)),
94 this, SIGNAL(newNMEAstr(QByteArray)));
95}
96
97// Destructor
98////////////////////////////////////////////////////////////////////////////
99bncPPPclient::~bncPPPclient() {
100 delete _model;
101 delete _epoData;
102 QMapIterator<QString, t_ephPair*> it(_eph);
103 while (it.hasNext()) {
104 it.next();
105 delete it.value();
106 }
107 QMapIterator<QString, t_corr*> ic(_corr);
108 while (ic.hasNext()) {
109 ic.next();
110 delete ic.value();
111 }
112 QMapIterator<QString, t_bias*> ib(_bias);
113 while (ib.hasNext()) {
114 ib.next();
115 delete ib.value();
116 }
117}
118
119//
120////////////////////////////////////////////////////////////////////////////
121void bncPPPclient::putNewObs(t_obs* obs) {
122 QMutexLocker locker(&_mutex);
123
124 if (obs->satSys != 'G' && !_useGlonass) {
125 return;
126 }
127
128 t_satData* satData = new t_satData();
129
130 // Satellite Number
131 // ----------------
132 if (obs->satSys == 'G') {
133 QString prn = QString("G%1").arg(obs->satNum, 2, 10, QChar('0'));
134 satData->prn = prn;
135 }
136 else if (obs->satSys == 'R') {
137 QString prn = QString("R%1").arg(obs->satNum, 2, 10, QChar('0'));
138 satData->prn = prn;
139 }
140
141 // Check Slips
142 // -----------
143 slipInfo& sInfo = _slips[satData->prn];
144 if ( sInfo.slipCntL1 == obs->slip_cnt_L1 &&
145 sInfo.slipCntL2 == obs->slip_cnt_L2 ) {
146 satData->slipFlag = false;
147 }
148 else {
149 satData->slipFlag = true;
150 }
151 sInfo.slipCntL1 = obs->slip_cnt_L1;
152 sInfo.slipCntL2 = obs->slip_cnt_L2;
153
154 // Handle Code Biases
155 // ------------------
156 t_bias* bb = 0;
157 if (_bias.contains(satData->prn)) {
158 bb = _bias.value(satData->prn);
159 }
160
161 // Set Code Observations
162 // ---------------------
163 if (obs->P1) {
164 satData->P1 = obs->P1 + (bb ? bb->p1 : 0.0);
165 satData->codeTypeF1 = t_satData::P_CODE;
166 }
167 else if (obs->C1) {
168 satData->P1 = obs->C1 + (bb ? bb->c1 : 0.0);
169 satData->codeTypeF1 = t_satData::C_CODE;
170 }
171 else {
172 delete satData;
173 return;
174 }
175
176 if (obs->P2) {
177 satData->P2 = obs->P2 + (bb ? bb->p2 : 0.0);
178 satData->codeTypeF2 = t_satData::P_CODE;
179 }
180 else if (obs->C2) {
181 satData->P2 = obs->C2;
182 satData->codeTypeF2 = t_satData::C_CODE;
183 }
184 else {
185 delete satData;
186 return;
187 }
188
189 double f1 = t_CST::freq1;
190 double f2 = t_CST::freq2;
191
192 if (obs->satSys == 'R') {
193 f1 = 1602000000.0 + 562500.0 * obs->slotNum;
194 f2 = 1246000000.0 + 437500.0 * obs->slotNum;
195 }
196
197 // Ionosphere-Free Combination
198 // ---------------------------
199 double c1 = f1 * f1 / (f1 * f1 - f2 * f2);
200 double c2 = - f2 * f2 / (f1 * f1 - f2 * f2);
201
202 satData->P3 = c1 * satData->P1 + c2 * satData->P2;
203
204 // Set Phase Observations
205 // ----------------------
206 if (obs->L1() && obs->L2()) {
207 satData->L1 = obs->L1() * t_CST::c / f1;
208 satData->L2 = obs->L2() * t_CST::c / f2;
209 }
210 else {
211 delete satData;
212 return;
213 }
214 satData->L3 = c1 * satData->L1 + c2 * satData->L2;
215
216 // Set Ionosphere-Free Wavelength
217 // ------------------------------
218 satData->lambda3 = c1 * t_CST::c / f1 + c2 * t_CST::c / f2;
219
220 // Add new Satellite to the epoch
221 // ------------------------------
222 bncTime tt(obs->GPSWeek, obs->GPSWeeks);
223
224 if (!_epoData) {
225 _epoData = new t_epoData();
226 _epoData->tt = tt;
227 }
228 else if (tt != _epoData->tt) {
229 processEpoch();
230 delete _epoData;
231 _epoData = new t_epoData();
232 _epoData->tt = tt;
233 }
234
235 if (obs->satSys == 'G') {
236 _epoData->satDataGPS[satData->prn] = satData;
237 }
238 else if (obs->satSys == 'R') {
239 _epoData->satDataGlo[satData->prn] = satData;
240 }
241}
242
243//
244////////////////////////////////////////////////////////////////////////////
245void bncPPPclient::slotNewEphGPS(gpsephemeris gpseph) {
246 QMutexLocker locker(&_mutex);
247
248 QString prn = QString("G%1").arg(gpseph.satellite, 2, 10, QChar('0'));
249
250 if (_eph.contains(prn)) {
251 t_ephGPS* eLast = static_cast<t_ephGPS*>(_eph.value(prn)->last);
252 if ( (eLast->GPSweek() < gpseph.GPSweek) ||
253 (eLast->GPSweek() == gpseph.GPSweek &&
254 eLast->TOC() < gpseph.TOC) ) {
255 delete static_cast<t_ephGPS*>(_eph.value(prn)->prev);
256 _eph.value(prn)->prev = _eph.value(prn)->last;
257 _eph.value(prn)->last = new t_ephGPS();
258 static_cast<t_ephGPS*>(_eph.value(prn)->last)->set(&gpseph);
259 }
260 }
261 else {
262 t_ephGPS* eLast = new t_ephGPS();
263 eLast->set(&gpseph);
264 _eph.insert(prn, new t_ephPair());
265 _eph[prn]->last = eLast;
266 }
267}
268
269//
270////////////////////////////////////////////////////////////////////////////
271void bncPPPclient::slotNewEphGlonass(glonassephemeris gloeph) {
272 QMutexLocker locker(&_mutex);
273
274 QString prn = QString("R%1").arg(gloeph.almanac_number, 2, 10, QChar('0'));
275
276 if (_eph.contains(prn)) {
277 int ww = gloeph.GPSWeek;
278 int tow = gloeph.GPSTOW;
279 updatetime(&ww, &tow, gloeph.tb*1000, 0); // Moscow -> GPS
280 t_ephGlo* eLast = static_cast<t_ephGlo*>(_eph.value(prn)->last);
281 if (eLast->GPSweek() < ww ||
282 (eLast->GPSweek() == ww && eLast->GPSweeks() < tow)) {
283 delete static_cast<t_ephGlo*>(_eph.value(prn)->prev);
284 _eph.value(prn)->prev = _eph.value(prn)->last;
285 _eph.value(prn)->last = new t_ephGlo();
286 static_cast<t_ephGlo*>(_eph.value(prn)->last)->set(&gloeph);
287 }
288 }
289 else {
290 t_ephGlo* eLast = new t_ephGlo();
291 eLast->set(&gloeph);
292 _eph.insert(prn, new t_ephPair());
293 _eph[prn]->last = eLast;
294 }
295}
296
297//
298////////////////////////////////////////////////////////////////////////////
299void bncPPPclient::slotNewCorrections(QList<QString> corrList) {
300 QMutexLocker locker(&_mutex);
301
302 if (corrList.size() == 0) {
303 return;
304 }
305
306 // Remove All Corrections
307 // ----------------------
308 QMapIterator<QString, t_corr*> ic(_corr);
309 while (ic.hasNext()) {
310 ic.next();
311 delete ic.value();
312 }
313 _corr.clear();
314
315 QListIterator<QString> it(corrList);
316 while (it.hasNext()) {
317 QTextStream in(it.next().toAscii());
318 int messageType;
319 int updateInterval;
320 int GPSweek;
321 double GPSweeks;
322 QString prn;
323 in >> messageType >> updateInterval >> GPSweek >> GPSweeks >> prn;
324 if ( messageType == COTYPE_GPSCOMBINED ||
325 messageType == COTYPE_GLONASSCOMBINED ||
326 messageType == COTYPE_GPSORBIT ||
327 messageType == COTYPE_GPSCLOCK ||
328 messageType == COTYPE_GLONASSORBIT ||
329 messageType == COTYPE_GLONASSCLOCK ) {
330
331 t_corr* cc = 0;
332 if (_corr.contains(prn)) {
333 cc = _corr.value(prn);
334 }
335 else {
336 cc = new t_corr();
337 _corr[prn] = cc;
338 }
339
340 cc->tt.set(GPSweek, GPSweeks);
341
342 if ( messageType == COTYPE_GPSCOMBINED ||
343 messageType == COTYPE_GLONASSCOMBINED ) {
344 cc->rao.ReSize(3); cc->rao = 0.0;
345 cc->dotRao.ReSize(3); cc->dotRao = 0.0;
346 cc->dotDotRao.ReSize(3); cc->dotDotRao = 0.0;
347 cc->dClk = 0.0;
348 cc->dotDClk = 0.0;
349 cc->dotDotDClk = 0.0;
350 in >> cc->iod
351 >> cc->dClk >> cc->rao[0] >> cc->rao[1] >> cc->rao[2]
352 >> cc->dotDClk >> cc->dotRao[0] >> cc->dotRao[1] >> cc->dotRao[2]
353 >> cc->dotDotDClk >> cc->dotDotRao[0] >> cc->dotDotRao[1] >> cc->dotDotRao[2];
354 cc->dClk /= t_CST::c;
355 cc->dotDClk /= t_CST::c;
356 cc->dotDotDClk /= t_CST::c;
357 cc->raoSet = true;
358 cc->dClkSet = true;
359 }
360 else if ( messageType == COTYPE_GPSORBIT ||
361 messageType == COTYPE_GLONASSORBIT ) {
362 cc->rao.ReSize(3); cc->rao = 0.0;
363 cc->dotRao.ReSize(3); cc->dotRao = 0.0;
364 cc->dotDotRao.ReSize(3); cc->dotDotRao = 0.0;
365 in >> cc->iod
366 >> cc->rao[0] >> cc->rao[1] >> cc->rao[2]
367 >> cc->dotRao[0] >> cc->dotRao[1] >> cc->dotRao[2]
368 >> cc->dotDotRao[0] >> cc->dotDotRao[1] >> cc->dotDotRao[2];
369 cc->raoSet = true;
370 }
371 else if ( messageType == COTYPE_GPSCLOCK ||
372 messageType == COTYPE_GLONASSCLOCK ) {
373 int dummyIOD;
374 cc->dClk = 0.0;
375 cc->dotDClk = 0.0;
376 cc->dotDotDClk = 0.0;
377 in >> dummyIOD >> cc->dClk >> cc->dotDClk >> cc->dotDotDClk;
378 cc->dClk /= t_CST::c;
379 cc->dotDClk /= t_CST::c;
380 cc->dotDotDClk /= t_CST::c;
381 cc->dClkSet = true;
382 }
383 }
384 else if ( messageType == BTYPE_GPS ) {
385
386 t_bias* bb = 0;
387 if (_bias.contains(prn)) {
388 bb = _bias.value(prn);
389 }
390 else {
391 bb = new t_bias();
392 _bias[prn] = bb;
393 }
394
395 bb->tt.set(GPSweek, GPSweeks);
396
397 int numBiases;
398 in >> numBiases;
399 for (int ii = 0; ii < numBiases; ++ii) {
400 int bType;
401 double bValue;
402 in >> bType >> bValue;
403 if (bType == CODETYPEGPS_L1_Z) {
404 bb->p1 = bValue;
405 }
406 else if (bType == CODETYPEGPS_L1_CA) {
407 bb->c1 = bValue;
408 }
409 else if (bType == CODETYPEGPS_L2_Z) {
410 bb->p2 = bValue;
411 }
412 }
413 }
414 }
415
416 QMutableMapIterator<QString, t_corr*> im(_corr);
417 while (im.hasNext()) {
418 im.next();
419 t_corr* cc = im.value();
420 if (!cc->ready()) {
421 delete cc;
422 im.remove();
423 }
424 }
425}
426
427// Satellite Position
428////////////////////////////////////////////////////////////////////////////
429t_irc bncPPPclient::getSatPos(const bncTime& tt, const QString& prn,
430 ColumnVector& xc, ColumnVector& vv) {
431
432 const double MAXAGE = 120.0;
433
434 if (_eph.contains(prn)) {
435
436 if (_pppMode) {
437 if (_corr.contains(prn)) {
438 t_corr* cc = _corr.value(prn);
439 if (tt - cc->tt < MAXAGE) {
440 t_eph* eLast = _eph.value(prn)->last;
441 t_eph* ePrev = _eph.value(prn)->prev;
442 if (eLast && eLast->IOD() == cc->iod) {
443 eLast->position(tt.gpsw(), tt.gpssec(), xc.data(), vv.data());
444 applyCorr(tt, cc, xc, vv);
445 return success;
446 }
447 else if (ePrev && ePrev->IOD() == cc->iod) {
448 ePrev->position(tt.gpsw(), tt.gpssec(), xc.data(), vv.data());
449 applyCorr(tt, cc, xc, vv);
450 return success;
451 }
452 }
453 }
454 }
455
456 else {
457 t_eph* ee = _eph.value(prn)->last;
458 ee->position(tt.gpsw(), tt.gpssec(), xc.data(), vv.data());
459 return success;
460 }
461 }
462
463 return failure;
464}
465
466//
467////////////////////////////////////////////////////////////////////////////
468void bncPPPclient::applyCorr(const bncTime& tt, const t_corr* cc,
469 ColumnVector& xc, ColumnVector& vv) {
470 ColumnVector dx(3);
471
472 double dt = tt - cc->tt;
473 ColumnVector raoHlp = cc->rao + cc->dotRao * dt + cc->dotDotRao * dt * dt;
474
475 RSW_to_XYZ(xc.Rows(1,3), vv, raoHlp, dx);
476
477 xc[0] -= dx[0];
478 xc[1] -= dx[1];
479 xc[2] -= dx[2];
480 xc[3] += cc->dClk + cc->dotDClk * dt + cc->dotDotDClk * dt * dt;
481}
482
483// Correct Time of Transmission
484////////////////////////////////////////////////////////////////////////////
485t_irc bncPPPclient::cmpToT(t_satData* satData) {
486
487 double prange = satData->P3;
488 if (prange == 0.0) {
489 return failure;
490 }
491
492 double clkSat = 0.0;
493 for (int ii = 1; ii <= 10; ii++) {
494
495 bncTime ToT = _epoData->tt - prange / t_CST::c - clkSat;
496
497 ColumnVector xc(4);
498 ColumnVector vv(3);
499 if (getSatPos(ToT, satData->prn, xc, vv) != success) {
500 return failure;
501 }
502
503 double clkSatOld = clkSat;
504 clkSat = xc(4);
505
506 if ( fabs(clkSat-clkSatOld) * t_CST::c < 1.e-4 ) {
507 satData->xx = xc.Rows(1,3);
508 satData->vv = vv;
509 satData->clk = clkSat * t_CST::c;
510 return success;
511 }
512 }
513
514 return failure;
515}
516
517//
518////////////////////////////////////////////////////////////////////////////
519void bncPPPclient::processEpoch() {
520
521 // Data Pre-Processing
522 // -------------------
523 QMutableMapIterator<QString, t_satData*> iGPS(_epoData->satDataGPS);
524 while (iGPS.hasNext()) {
525 iGPS.next();
526 QString prn = iGPS.key();
527 t_satData* satData = iGPS.value();
528
529 if (cmpToT(satData) != success) {
530 delete satData;
531 iGPS.remove();
532 continue;
533 }
534 }
535
536 QMutableMapIterator<QString, t_satData*> iGlo(_epoData->satDataGlo);
537 while (iGlo.hasNext()) {
538 iGlo.next();
539 QString prn = iGlo.key();
540 t_satData* satData = iGlo.value();
541
542 if (cmpToT(satData) != success) {
543 delete satData;
544 iGlo.remove();
545 continue;
546 }
547 }
548
549 // Filter Solution
550 // ---------------
551 if (_model->update(_epoData) == success) {
552 emit newPosition(_model->time(), _model->x(), _model->y(), _model->z());
553 }
554}
555
Note: See TracBrowser for help on using the repository browser.