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

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