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

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