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

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