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

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