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

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