source: ntrip/trunk/BNC/src/PPP_free/pppClient.cpp@ 6663

Last change on this file since 6663 was 6653, checked in by stuerze, 10 years ago

sinex tro file support added,
troposphere results in overall ppp logfile added

File size: 10.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: 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 "bncephuser.h"
47#include "bncutils.h"
48
49using namespace BNC_PPP;
50using namespace std;
51
52// Constructor
53////////////////////////////////////////////////////////////////////////////
54t_pppClient::t_pppClient(const t_pppOptions* opt) {
55
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}
62
63// Destructor
64////////////////////////////////////////////////////////////////////////////
65t_pppClient::~t_pppClient() {
66 delete _filter;
67 delete _epoData;
68 delete _opt;
69 delete _ephUser;
70 delete _log;
71}
72
73//
74////////////////////////////////////////////////////////////////////////////
75void t_pppClient::processEpoch(const vector<t_satObs*>& satObs, t_output* output) {
76
77 // Convert and store observations
78 // ------------------------------
79 _epoData->clear();
80 for (unsigned ii = 0; ii < satObs.size(); ii++) {
81 const t_satObs* obs = satObs[ii];
82 t_satData* satData = new t_satData();
83
84 if (_epoData->tt.undef()) {
85 _epoData->tt = obs->_time;
86 }
87
88 satData->tt = obs->_time;
89 satData->prn = QString(obs->_prn.toString().c_str());
90 satData->slipFlag = false;
91 satData->P1 = 0.0;
92 satData->P2 = 0.0;
93 satData->P5 = 0.0;
94 satData->L1 = 0.0;
95 satData->L2 = 0.0;
96 satData->L5 = 0.0;
97 for (unsigned ifrq = 0; ifrq < obs->_obs.size(); ifrq++) {
98 t_frqObs* frqObs = obs->_obs[ifrq];
99 if (frqObs->_rnxType2ch[0] == '1') {
100 if (frqObs->_codeValid) satData->P1 = frqObs->_code;
101 if (frqObs->_phaseValid) satData->L1 = frqObs->_phase;
102 if (frqObs->_slip) satData->slipFlag = true;
103 }
104 else if (frqObs->_rnxType2ch[0] == '2') {
105 if (frqObs->_codeValid) satData->P2 = frqObs->_code;
106 if (frqObs->_phaseValid) satData->L2 = frqObs->_phase;
107 if (frqObs->_slip) satData->slipFlag = true;
108 }
109 else if (frqObs->_rnxType2ch[0] == '5') {
110 if (frqObs->_codeValid) satData->P5 = frqObs->_code;
111 if (frqObs->_phaseValid) satData->L5 = frqObs->_phase;
112 if (frqObs->_slip) satData->slipFlag = true;
113 }
114 }
115 putNewObs(satData);
116 }
117
118 // Data Pre-Processing
119 // -------------------
120 QMutableMapIterator<QString, t_satData*> it(_epoData->satData);
121 while (it.hasNext()) {
122 it.next();
123 QString prn = it.key();
124 t_satData* satData = it.value();
125
126 if (cmpToT(satData) != success) {
127 delete satData;
128 it.remove();
129 continue;
130 }
131 }
132
133 // Filter Solution
134 // ---------------
135 if (_filter->update(_epoData) == success) {
136 output->_error = false;
137 output->_epoTime = _filter->time();
138 output->_xyzRover[0] = _filter->x();
139 output->_xyzRover[1] = _filter->y();
140 output->_xyzRover[2] = _filter->z();
141 copy(&_filter->Q().data()[0], &_filter->Q().data()[6], output->_covMatrix);
142 output->_neu[0] = _filter->neu()[0];
143 output->_neu[1] = _filter->neu()[1];
144 output->_neu[2] = _filter->neu()[2];
145 output->_numSat = _filter->numSat();
146 output->_pDop = _filter->PDOP();
147 output->_trp0 = _filter->delay_saast(M_PI/2.0);
148 output->_trp = _filter->trp();
149 }
150 else {
151 output->_error = true;
152 }
153
154 output->_log = _log->str();
155 delete _log; _log = new ostringstream();
156}
157
158//
159////////////////////////////////////////////////////////////////////////////
160void t_pppClient::putNewObs(t_satData* satData) {
161
162 // Set Observations GPS and Glonass
163 // --------------------------------
164 if (satData->system() == 'G' || satData->system() == 'R') {
165 if (satData->P1 != 0.0 && satData->P2 != 0.0 &&
166 satData->L1 != 0.0 && satData->L2 != 0.0 ) {
167
168 int channel = 0;
169 if (satData->system() == 'R') {
170 const t_eph* eph = _ephUser->ephLast(satData->prn);
171 if (eph) {
172 channel = eph->slotNum();
173 }
174 else {
175 delete satData;
176 return;
177 }
178 }
179
180 t_frequency::type fType1 = t_lc::toFreq(satData->system(), t_lc::l1);
181 t_frequency::type fType2 = t_lc::toFreq(satData->system(), t_lc::l2);
182 double f1 = t_CST::freq(fType1, channel);
183 double f2 = t_CST::freq(fType2, channel);
184 double a1 = f1 * f1 / (f1 * f1 - f2 * f2);
185 double a2 = - f2 * f2 / (f1 * f1 - f2 * f2);
186 satData->L1 = satData->L1 * t_CST::c / f1;
187 satData->L2 = satData->L2 * t_CST::c / f2;
188 satData->P3 = a1 * satData->P1 + a2 * satData->P2;
189 satData->L3 = a1 * satData->L1 + a2 * satData->L2;
190 satData->lambda3 = a1 * t_CST::c / f1 + a2 * t_CST::c / f2;
191 satData->lkA = a1;
192 satData->lkB = a2;
193 _epoData->satData[satData->prn] = satData;
194 }
195 else {
196 delete satData;
197 }
198 }
199
200 // Set Observations Galileo
201 // ------------------------
202 else if (satData->system() == 'E') {
203 if (satData->P1 != 0.0 && satData->P5 != 0.0 &&
204 satData->L1 != 0.0 && satData->L5 != 0.0 ) {
205 double f1 = t_CST::freq(t_frequency::E1, 0);
206 double f5 = t_CST::freq(t_frequency::E5, 0);
207 double a1 = f1 * f1 / (f1 * f1 - f5 * f5);
208 double a5 = - f5 * f5 / (f1 * f1 - f5 * f5);
209 satData->L1 = satData->L1 * t_CST::c / f1;
210 satData->L5 = satData->L5 * t_CST::c / f5;
211 satData->P3 = a1 * satData->P1 + a5 * satData->P5;
212 satData->L3 = a1 * satData->L1 + a5 * satData->L5;
213 satData->lambda3 = a1 * t_CST::c / f1 + a5 * t_CST::c / f5;
214 satData->lkA = a1;
215 satData->lkB = a5;
216 _epoData->satData[satData->prn] = satData;
217 }
218 else {
219 delete satData;
220 }
221 }
222}
223
224//
225////////////////////////////////////////////////////////////////////////////
226void t_pppClient::putOrbCorrections(const std::vector<t_orbCorr*>& corr) {
227 for (unsigned ii = 0; ii < corr.size(); ii++) {
228 QString prn = QString(corr[ii]->_prn.toString().c_str());
229 t_eph* eLast = _ephUser->ephLast(prn);
230 t_eph* ePrev = _ephUser->ephPrev(prn);
231 if (eLast && eLast->IOD() == corr[ii]->_iod) {
232 eLast->setOrbCorr(corr[ii]);
233 }
234 else if (ePrev && ePrev->IOD() == corr[ii]->_iod) {
235 ePrev->setOrbCorr(corr[ii]);
236 }
237 }
238}
239
240//
241////////////////////////////////////////////////////////////////////////////
242void t_pppClient::putClkCorrections(const std::vector<t_clkCorr*>& corr) {
243 for (unsigned ii = 0; ii < corr.size(); ii++) {
244 QString prn = QString(corr[ii]->_prn.toString().c_str());
245 t_eph* eLast = _ephUser->ephLast(prn);
246 t_eph* ePrev = _ephUser->ephPrev(prn);
247 if (eLast && eLast->IOD() == corr[ii]->_iod) {
248 eLast->setClkCorr(corr[ii]);
249 }
250 else if (ePrev && ePrev->IOD() == corr[ii]->_iod) {
251 ePrev->setClkCorr(corr[ii]);
252 }
253 }
254}
255
256//
257//////////////////////////////////////////////////////////////////////////////
258void t_pppClient::putCodeBiases(const std::vector<t_satCodeBias*>& /* satCodeBias */) {
259}
260
261//
262//////////////////////////////////////////////////////////////////////////////
263void t_pppClient::putEphemeris(const t_eph* eph) {
264 const t_ephGPS* ephGPS = dynamic_cast<const t_ephGPS*>(eph);
265 const t_ephGlo* ephGlo = dynamic_cast<const t_ephGlo*>(eph);
266 const t_ephGal* ephGal = dynamic_cast<const t_ephGal*>(eph);
267 if (ephGPS) {
268 _ephUser->putNewEph(new t_ephGPS(*ephGPS), false);
269 }
270 else if (ephGlo) {
271 _ephUser->putNewEph(new t_ephGlo(*ephGlo), false);
272 }
273 else if (ephGal) {
274 _ephUser->putNewEph(new t_ephGal(*ephGal), false);
275 }
276}
277
278// Satellite Position
279////////////////////////////////////////////////////////////////////////////
280t_irc t_pppClient::getSatPos(const bncTime& tt, const QString& prn,
281 ColumnVector& xc, ColumnVector& vv) {
282
283 t_eph* eLast = _ephUser->ephLast(prn);
284 t_eph* ePrev = _ephUser->ephPrev(prn);
285 if (eLast && eLast->getCrd(tt, xc, vv, _opt->useOrbClkCorr()) == success) {
286 return success;
287 }
288 else if (ePrev && ePrev->getCrd(tt, xc, vv, _opt->useOrbClkCorr()) == success) {
289 return success;
290 }
291 return failure;
292}
293
294// Correct Time of Transmission
295////////////////////////////////////////////////////////////////////////////
296t_irc t_pppClient::cmpToT(t_satData* satData) {
297
298 double prange = satData->P3;
299 if (prange == 0.0) {
300 return failure;
301 }
302
303 double clkSat = 0.0;
304 for (int ii = 1; ii <= 10; ii++) {
305
306 bncTime ToT = satData->tt - prange / t_CST::c - clkSat;
307
308 ColumnVector xc(4);
309 ColumnVector vv(3);
310 if (getSatPos(ToT, satData->prn, xc, vv) != success) {
311 return failure;
312 }
313
314 double clkSatOld = clkSat;
315 clkSat = xc(4);
316
317 if ( fabs(clkSat-clkSatOld) * t_CST::c < 1.e-4 ) {
318 satData->xx = xc.Rows(1,3);
319 satData->vv = vv;
320 satData->clk = clkSat * t_CST::c;
321 return success;
322 }
323 }
324
325 return failure;
326}
327
Note: See TracBrowser for help on using the repository browser.