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

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