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

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