source: ntrip/trunk/BNC/src/PPP/pppClient.cpp@ 7232

Last change on this file since 7232 was 7231, checked in by stuerze, 10 years ago

some interfaces are added to be able to handle ssr vtec in PPP mode
File size: 13.3 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: t_pppClient
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 "pppClient.h"
46#include "pppUtils.h"
47#include "bncephuser.h"
48#include "bncutils.h"
49
50using namespace BNC_PPP;
51using namespace std;
52
53// Constructor
54////////////////////////////////////////////////////////////////////////////
55t_pppClient::t_pppClient(const t_pppOptions* opt) {
56 _opt = new t_pppOptions(*opt);
57 _filter = new t_pppFilter(this);
58 _epoData = new t_epoData();
59 _log = new ostringstream();
60 _ephUser = new bncEphUser(false);
61 _pppUtils = new t_pppUtils();
62}
63
64// Destructor
65////////////////////////////////////////////////////////////////////////////
66t_pppClient::~t_pppClient() {
67 delete _filter;
68 delete _epoData;
69 delete _opt;
70 delete _ephUser;
71 delete _log;
72}
73
74//
75////////////////////////////////////////////////////////////////////////////
76void t_pppClient::processEpoch(const vector<t_satObs*>& satObs, t_output* output) {
77
78 // Convert and store observations
79 // ------------------------------
80 _epoData->clear();
81 for (unsigned ii = 0; ii < satObs.size(); ii++) {
82 const t_satObs* obs = satObs[ii];
83 t_prn prn = obs->_prn;
84 if (prn.system() == 'E') {prn.setFlags(1);} // force I/NAV usage
85 t_satData* satData = new t_satData();
86
87 if (_epoData->tt.undef()) {
88 _epoData->tt = obs->_time;
89 }
90
91 satData->tt = obs->_time;
92 satData->prn = QString(prn.toInternalString().c_str());
93 satData->slipFlag = false;
94 satData->P1 = 0.0;
95 satData->P2 = 0.0;
96 satData->P5 = 0.0;
97 satData->P7 = 0.0;
98 satData->L1 = 0.0;
99 satData->L2 = 0.0;
100 satData->L5 = 0.0;
101 satData->L7 = 0.0;
102 for (unsigned ifrq = 0; ifrq < obs->_obs.size(); ifrq++) {
103 t_frqObs* frqObs = obs->_obs[ifrq];
104 double cb = 0.0;
105 const t_satCodeBias* satCB = _pppUtils->satCodeBias(prn);
106 if (satCB && satCB->_bias.size()) {
107 for (unsigned ii = 0; ii < satCB->_bias.size(); ii++) {
108
109 const t_frqCodeBias& bias = satCB->_bias[ii];
110 if (frqObs && frqObs->_rnxType2ch == bias._rnxType2ch) {
111 cb = bias._value;
112 }
113 }
114 }
115 if (frqObs->_rnxType2ch[0] == '1') {
116 if (frqObs->_codeValid) satData->P1 = frqObs->_code + cb;
117 if (frqObs->_phaseValid) satData->L1 = frqObs->_phase;
118 if (frqObs->_slip) satData->slipFlag = true;
119 }
120 else if (frqObs->_rnxType2ch[0] == '2') {
121 if (frqObs->_codeValid) satData->P2 = frqObs->_code + cb;
122 if (frqObs->_phaseValid) satData->L2 = frqObs->_phase;
123 if (frqObs->_slip) satData->slipFlag = true;
124 }
125 else if (frqObs->_rnxType2ch[0] == '5') {
126 if (frqObs->_codeValid) satData->P5 = frqObs->_code + cb;
127 if (frqObs->_phaseValid) satData->L5 = frqObs->_phase;
128 if (frqObs->_slip) satData->slipFlag = true;
129 }
130 else if (frqObs->_rnxType2ch[0] == '7') {
131 if (frqObs->_codeValid) satData->P7 = frqObs->_code + cb;
132 if (frqObs->_phaseValid) satData->L7 = frqObs->_phase;
133 if (frqObs->_slip) satData->slipFlag = true;
134 }
135 }
136 putNewObs(satData);
137 }
138
139 // Data Pre-Processing
140 // -------------------
141 QMutableMapIterator<QString, t_satData*> it(_epoData->satData);
142 while (it.hasNext()) {
143 it.next();
144 QString prn = it.key();
145 t_satData* satData = it.value();
146
147 if (cmpToT(satData) != success) {
148 delete satData;
149 it.remove();
150 continue;
151 }
152
153 }
154
155 // Filter Solution
156 // ---------------
157 if (_filter->update(_epoData) == success) {
158 output->_error = false;
159 output->_epoTime = _filter->time();
160 output->_xyzRover[0] = _filter->x();
161 output->_xyzRover[1] = _filter->y();
162 output->_xyzRover[2] = _filter->z();
163 copy(&_filter->Q().data()[0], &_filter->Q().data()[6], output->_covMatrix);
164 output->_neu[0] = _filter->neu()[0];
165 output->_neu[1] = _filter->neu()[1];
166 output->_neu[2] = _filter->neu()[2];
167 output->_numSat = _filter->numSat();
168 output->_pDop = _filter->PDOP();
169 output->_trp0 = _filter->trp0();
170 output->_trp = _filter->trp();
171 }
172 else {
173 output->_error = true;
174 }
175
176 output->_log = _log->str();
177 delete _log; _log = new ostringstream();
178}
179
180//
181////////////////////////////////////////////////////////////////////////////
182void t_pppClient::putNewObs(t_satData* satData) {
183
184 // Set Observations GPS
185 // --------------------
186 if (satData->system() == 'G') {
187 if (satData->P1 != 0.0 && satData->P2 != 0.0 &&
188 satData->L1 != 0.0 && satData->L2 != 0.0 ) {
189 t_frequency::type fType1 = t_lc::toFreq(satData->system(), t_lc::l1);
190 t_frequency::type fType2 = t_lc::toFreq(satData->system(), t_lc::l2);
191 double f1 = t_CST::freq(fType1, 0);
192 double f2 = t_CST::freq(fType2, 0);
193 double a1 = f1 * f1 / (f1 * f1 - f2 * f2);
194 double a2 = - f2 * f2 / (f1 * f1 - f2 * f2);
195 satData->L1 = satData->L1 * t_CST::c / f1;
196 satData->L2 = satData->L2 * t_CST::c / f2;
197 satData->P3 = a1 * satData->P1 + a2 * satData->P2;
198 satData->L3 = a1 * satData->L1 + a2 * satData->L2;
199 satData->lambda3 = a1 * t_CST::c / f1 + a2 * t_CST::c / f2;
200 satData->lkA = a1;
201 satData->lkB = a2;
202 _epoData->satData[satData->prn] = satData;
203 }
204 else {
205 delete satData;
206 }
207 }
208
209 // Set Observations Glonass
210 // ------------------------
211 else if (satData->system() == 'R' && _opt->useSystem('R')) {
212 if (satData->P1 != 0.0 && satData->P2 != 0.0 &&
213 satData->L1 != 0.0 && satData->L2 != 0.0 ) {
214
215 int channel = 0;
216 if (satData->system() == 'R') {
217 const t_eph* eph = _ephUser->ephLast(satData->prn);
218 if (eph) {
219 channel = eph->slotNum();
220 }
221 else {
222 delete satData;
223 return;
224 }
225 }
226
227 t_frequency::type fType1 = t_lc::toFreq(satData->system(), t_lc::l1);
228 t_frequency::type fType2 = t_lc::toFreq(satData->system(), t_lc::l2);
229 double f1 = t_CST::freq(fType1, channel);
230 double f2 = t_CST::freq(fType2, channel);
231 double a1 = f1 * f1 / (f1 * f1 - f2 * f2);
232 double a2 = - f2 * f2 / (f1 * f1 - f2 * f2);
233 satData->L1 = satData->L1 * t_CST::c / f1;
234 satData->L2 = satData->L2 * t_CST::c / f2;
235 satData->P3 = a1 * satData->P1 + a2 * satData->P2;
236 satData->L3 = a1 * satData->L1 + a2 * satData->L2;
237 satData->lambda3 = a1 * t_CST::c / f1 + a2 * t_CST::c / f2;
238 satData->lkA = a1;
239 satData->lkB = a2;
240 _epoData->satData[satData->prn] = satData;
241 }
242 else {
243 delete satData;
244 }
245 }
246
247 // Set Observations Galileo
248 // ------------------------
249 else if (satData->system() == 'E' && _opt->useSystem('E')) {
250 if (satData->P1 != 0.0 && satData->P5 != 0.0 &&
251 satData->L1 != 0.0 && satData->L5 != 0.0 ) {
252 double f1 = t_CST::freq(t_frequency::E1, 0);
253 double f5 = t_CST::freq(t_frequency::E5, 0);
254 double a1 = f1 * f1 / (f1 * f1 - f5 * f5);
255 double a5 = - f5 * f5 / (f1 * f1 - f5 * f5);
256 satData->L1 = satData->L1 * t_CST::c / f1;
257 satData->L5 = satData->L5 * t_CST::c / f5;
258 satData->P3 = a1 * satData->P1 + a5 * satData->P5;
259 satData->L3 = a1 * satData->L1 + a5 * satData->L5;
260 satData->lambda3 = a1 * t_CST::c / f1 + a5 * t_CST::c / f5;
261 satData->lkA = a1;
262 satData->lkB = a5;
263 _epoData->satData[satData->prn] = satData;
264 }
265 else {
266 delete satData;
267 }
268 }
269
270 // Set Observations BDS
271 // ---------------------
272 else if (satData->system() == 'C' && _opt->useSystem('C')) {
273 if (satData->P2 != 0.0 && satData->P7 != 0.0 &&
274 satData->L2 != 0.0 && satData->L7 != 0.0 ) {
275 double f2 = t_CST::freq(t_frequency::C2, 0);
276 double f7 = t_CST::freq(t_frequency::C7, 0);
277 double a2 = f2 * f2 / (f2 * f2 - f7 * f7);
278 double a7 = - f7 * f7 / (f2 * f2 - f7 * f7);
279 satData->L2 = satData->L2 * t_CST::c / f2;
280 satData->L7 = satData->L7 * t_CST::c / f7;
281 satData->P3 = a2 * satData->P2 + a7 * satData->P7;
282 satData->L3 = a2 * satData->L2 + a7 * satData->L7;
283 satData->lambda3 = a2 * t_CST::c / f2 + a7 * t_CST::c / f7;
284 satData->lkA = a2;
285 satData->lkB = a7;
286 _epoData->satData[satData->prn] = satData;
287 }
288 else {
289 delete satData;
290 }
291 }
292}
293
294//
295////////////////////////////////////////////////////////////////////////////
296void t_pppClient::putOrbCorrections(const std::vector<t_orbCorr*>& corr) {
297 for (unsigned ii = 0; ii < corr.size(); ii++) {
298 QString prn = QString(corr[ii]->_prn.toInternalString().c_str());
299 t_eph* eLast = _ephUser->ephLast(prn);
300 t_eph* ePrev = _ephUser->ephPrev(prn);
301 if (eLast && eLast->IOD() == corr[ii]->_iod) {
302 eLast->setOrbCorr(corr[ii]);
303 }
304 else if (ePrev && ePrev->IOD() == corr[ii]->_iod) {
305 ePrev->setOrbCorr(corr[ii]);
306 }
307 }
308}
309
310//
311////////////////////////////////////////////////////////////////////////////
312void t_pppClient::putClkCorrections(const std::vector<t_clkCorr*>& corr) {
313 for (unsigned ii = 0; ii < corr.size(); ii++) {
314 QString prn = QString(corr[ii]->_prn.toInternalString().c_str());
315 t_eph* eLast = _ephUser->ephLast(prn);
316 t_eph* ePrev = _ephUser->ephPrev(prn);
317 if (eLast && eLast->IOD() == corr[ii]->_iod) {
318 eLast->setClkCorr(corr[ii]);
319 }
320 else if (ePrev && ePrev->IOD() == corr[ii]->_iod) {
321 ePrev->setClkCorr(corr[ii]);
322 }
323 }
324}
325
326//
327//////////////////////////////////////////////////////////////////////////////
328void t_pppClient::putCodeBiases(const std::vector<t_satCodeBias*>& satCodeBias) {
329 for (unsigned ii = 0; ii < satCodeBias.size(); ii++) {
330 _pppUtils->putCodeBias(new t_satCodeBias(*satCodeBias[ii]));
331 }
332}
333
334//
335//////////////////////////////////////////////////////////////////////////////
336void t_pppClient::putTec(const t_vTec* vTec) {
337 _pppUtils->putTec(new t_vTec(*vTec));
338}
339
340//
341//////////////////////////////////////////////////////////////////////////////
342void t_pppClient::putEphemeris(const t_eph* eph) {
343 bool check = _opt->_realTime;
344 const t_ephGPS* ephGPS = dynamic_cast<const t_ephGPS*>(eph);
345 const t_ephGlo* ephGlo = dynamic_cast<const t_ephGlo*>(eph);
346 const t_ephGal* ephGal = dynamic_cast<const t_ephGal*>(eph);
347 const t_ephBDS* ephBDS = dynamic_cast<const t_ephBDS*>(eph);
348 if (ephGPS) {
349 _ephUser->putNewEph(new t_ephGPS(*ephGPS), check);
350 }
351 else if (ephGlo) {
352 _ephUser->putNewEph(new t_ephGlo(*ephGlo), check);
353 }
354 else if (ephGal) {
355 _ephUser->putNewEph(new t_ephGal(*ephGal), check);
356 }
357 else if (ephBDS) {
358 _ephUser->putNewEph(new t_ephBDS(*ephBDS), check);
359 }
360}
361
362// Satellite Position
363////////////////////////////////////////////////////////////////////////////
364t_irc t_pppClient::getSatPos(const bncTime& tt, const QString& prn,
365 ColumnVector& xc, ColumnVector& vv) {
366
367 t_eph* eLast = _ephUser->ephLast(prn);
368 t_eph* ePrev = _ephUser->ephPrev(prn);
369 if (eLast && eLast->getCrd(tt, xc, vv, _opt->useOrbClkCorr()) == success) {
370 return success;
371 }
372 else if (ePrev && ePrev->getCrd(tt, xc, vv, _opt->useOrbClkCorr()) == success) {
373 return success;
374 }
375 return failure;
376}
377
378// Correct Time of Transmission
379////////////////////////////////////////////////////////////////////////////
380t_irc t_pppClient::cmpToT(t_satData* satData) {
381
382 double prange = satData->P3;
383 if (prange == 0.0) {
384 return failure;
385 }
386
387 double clkSat = 0.0;
388 for (int ii = 1; ii <= 10; ii++) {
389
390 bncTime ToT = satData->tt - prange / t_CST::c - clkSat;
391
392 ColumnVector xc(4);
393 ColumnVector vv(3);
394 if (getSatPos(ToT, satData->prn, xc, vv) != success) {
395 return failure;
396 }
397
398 double clkSatOld = clkSat;
399 clkSat = xc(4);
400
401 if ( fabs(clkSat-clkSatOld) * t_CST::c < 1.e-4 ) {
402 satData->xx = xc.Rows(1,3);
403 satData->vv = vv;
404 satData->clk = clkSat * t_CST::c;
405 return success;
406 }
407 }
408
409 return failure;
410}
Note: See TracBrowser for help on using the repository browser.