source: ntrip/trunk/BNC/src/PPP/pppParlist.cpp@ 6027

Last change on this file since 6027 was 6021, checked in by mervart, 10 years ago
File size: 12.0 KB
Line 
1
2// Part of BNC, a utility for retrieving decoding and
3// converting GNSS data streams from NTRIP broadcasters.
4//
5// Copyright (C) 2007
6// German Federal Agency for Cartography and Geodesy (BKG)
7// http://www.bkg.bund.de
8// Czech Technical University Prague, Department of Geodesy
9// http://www.fsv.cvut.cz
10//
11// Email: euref-ip@bkg.bund.de
12//
13// This program is free software; you can redistribute it and/or
14// modify it under the terms of the GNU General Public License
15// as published by the Free Software Foundation, version 2.
16//
17// This program is distributed in the hope that it will be useful,
18// but WITHOUT ANY WARRANTY; without even the implied warranty of
19// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20// GNU General Public License for more details.
21//
22// You should have received a copy of the GNU General Public License
23// along with this program; if not, write to the Free Software
24// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
25
26/* -------------------------------------------------------------------------
27 * BKG NTRIP Client
28 * -------------------------------------------------------------------------
29 *
30 * Class: t_pppParlist
31 *
32 * Purpose: List of estimated parameters
33 *
34 * Author: L. Mervart
35 *
36 * Created: 29-Jul-2014
37 *
38 * Changes:
39 *
40 * -----------------------------------------------------------------------*/
41
42#include <cmath>
43#include <iostream>
44#include <sstream>
45#include <iomanip>
46#include <algorithm>
47#include <newmatio.h>
48
49#include "pppParlist.h"
50#include "pppSatObs.h"
51#include "pppStation.h"
52#include "bncutils.h"
53#include "bncconst.h"
54#include "pppClient.h"
55
56using namespace BNC_PPP;
57using namespace std;
58
59// Constructor
60////////////////////////////////////////////////////////////////////////////
61t_pppParam::t_pppParam(e_type type, const t_prn& prn, t_lc::type tLC,
62 const vector<t_pppSatObs*>* obsVector) {
63
64 _type = type;
65 _prn = prn;
66 _tLC = tLC;
67 _x0 = 0.0;
68 _indexOld = -1;
69 _indexNew = -1;
70 _noise = 0.0;
71 _ambInfo = 0;
72
73 switch (_type) {
74 case crdX:
75 _epoSpec = false;
76 _sigma0 = OPT->_aprSigCrd[0];
77 _noise = OPT->_noiseCrd[0];
78 break;
79 case crdY:
80 _epoSpec = false;
81 _sigma0 = OPT->_aprSigCrd[1];
82 _noise = OPT->_noiseCrd[1];
83 break;
84 case crdZ:
85 _epoSpec = false;
86 _sigma0 = OPT->_aprSigCrd[2];
87 _noise = OPT->_noiseCrd[2];
88 break;
89 case clkR:
90 _epoSpec = true;
91 _sigma0 = OPT->_noiseClk;
92 break;
93 case amb:
94 _epoSpec = false;
95 _sigma0 = OPT->_aprSigAmb;
96 _ambInfo = new t_ambInfo();
97 if (obsVector) {
98 for (unsigned ii = 0; ii < obsVector->size(); ii++) {
99 const t_pppSatObs* obs = obsVector->at(ii);
100 if (obs->prn() == _prn) {
101 double offGG = 0;
102 if (_prn.system() == 'R' && tLC != t_lc::MW) {
103 offGG = PPP_CLIENT->offGG();
104 }
105 _x0 = floor((obs->obsValue(tLC) - offGG - obs->cmpValue(tLC)) / obs->lambda(tLC) + 0.5);
106 break;
107 }
108 }
109 }
110 break;
111 case offGG:
112 _epoSpec = true;
113 _sigma0 = 1000.0;
114 _x0 = PPP_CLIENT->offGG();
115 break;
116 case trp:
117 _epoSpec = false;
118 _sigma0 = OPT->_aprSigTrp;
119 _noise = OPT->_noiseTrp;
120 break;
121 }
122}
123
124// Destructor
125////////////////////////////////////////////////////////////////////////////
126t_pppParam::~t_pppParam() {
127 delete _ambInfo;
128}
129
130//
131////////////////////////////////////////////////////////////////////////////
132double t_pppParam::partial(const bncTime& /* epoTime */, const t_pppSatObs* obs,
133 const t_lc::type& tLC) const {
134
135 // Special Case - Melbourne-Wuebbena
136 // ---------------------------------
137 if (tLC == t_lc::MW && _type != amb) {
138 return 0.0;
139 }
140
141 const t_pppStation* sta = PPP_CLIENT->staRover();
142 ColumnVector rhoV = sta->xyzApr() - obs->xc().Rows(1,3);
143
144 switch (_type) {
145 case crdX:
146 return (sta->xyzApr()[0] - obs->xc()[0]) / rhoV.norm_Frobenius();
147 case crdY:
148 return (sta->xyzApr()[1] - obs->xc()[1]) / rhoV.norm_Frobenius();
149 case crdZ:
150 return (sta->xyzApr()[2] - obs->xc()[2]) / rhoV.norm_Frobenius();
151 case clkR:
152 return 1.0;
153 case offGG:
154 return (obs->prn().system() == 'R') ? 1.0 : 0.0;
155 case amb:
156 if (obs->prn() == _prn) {
157 if (tLC == _tLC) {
158 return (obs->lambda(tLC));
159 }
160 else if (tLC == t_lc::lIF && _tLC == t_lc::MW) {
161 return obs->lambda(t_lc::lIF) * obs->lambda(t_lc::MW) / obs->lambda(t_lc::l2);
162 }
163 else {
164 map<t_frequency::type, double> codeCoeff;
165 map<t_frequency::type, double> phaseCoeff;
166 obs->lcCoeff(tLC, codeCoeff, phaseCoeff);
167 if (_tLC == t_lc::l1) {
168 return obs->lambda(t_lc::l1) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
169 }
170 else if (_tLC == t_lc::l2) {
171 return obs->lambda(t_lc::l2) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
172 }
173 }
174 }
175 return 0.0;
176 case trp:
177 return 1.0 / sin(obs->eleSat());
178 }
179
180 return 0.0;
181}
182
183//
184////////////////////////////////////////////////////////////////////////////
185string t_pppParam::toString() const {
186 stringstream ss;
187 switch (_type) {
188 case crdX:
189 ss << "CRD_X";
190 break;
191 case crdY:
192 ss << "CRD_Y";
193 break;
194 case crdZ:
195 ss << "CRD_Z";
196 break;
197 case clkR:
198 ss << "CLK ";
199 break;
200 case amb:
201 ss << "AMB " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
202 break;
203 case offGG:
204 ss << "OGG ";
205 break;
206 case trp:
207 ss << "TRP ";
208 break;
209 }
210 return ss.str();
211}
212
213// Constructor
214////////////////////////////////////////////////////////////////////////////
215t_pppParlist::t_pppParlist() {
216}
217
218// Destructor
219////////////////////////////////////////////////////////////////////////////
220t_pppParlist::~t_pppParlist() {
221 for (unsigned ii = 0; ii < _params.size(); ii++) {
222 delete _params[ii];
223 }
224}
225
226//
227////////////////////////////////////////////////////////////////////////////
228t_irc t_pppParlist::set(const bncTime& epoTime, const std::vector<t_pppSatObs*>& obsVector) {
229
230 // Remove some Parameters
231 // ----------------------
232 vector<t_pppParam*>::iterator it = _params.begin();
233 while (it != _params.end()) {
234 t_pppParam* par = *it;
235
236 bool remove = false;
237
238 if (par->epoSpec()) {
239 remove = true;
240 }
241
242 else if (par->type() == t_pppParam::amb) {
243 if (par->lastObsTime().valid() && (epoTime - par->lastObsTime() > 120.0)) {
244 remove = true;
245 }
246 }
247
248 else if (par->type() == t_pppParam::amb) {
249 if (par->lastObsTime().valid() && (epoTime - par->lastObsTime() > 3600.0)) {
250 remove = true;
251 }
252 }
253
254 if (remove) {
255 delete par;
256 it = _params.erase(it);
257 }
258 else {
259 ++it;
260 }
261 }
262
263 // Check whether parameters have observations
264 // ------------------------------------------
265 for (unsigned ii = 0; ii < _params.size(); ii++) {
266 t_pppParam* par = _params[ii];
267 if (par->prn() == 0) {
268 par->setLastObsTime(epoTime);
269 if (par->firstObsTime().undef()) {
270 par->setFirstObsTime(epoTime);
271 }
272 }
273 else {
274 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
275 const t_pppSatObs* satObs = obsVector[jj];
276 if (satObs->prn() == par->prn()) {
277 par->setLastObsTime(epoTime);
278 if (par->firstObsTime().undef()) {
279 par->setFirstObsTime(epoTime);
280 }
281 break;
282 }
283 }
284 }
285 }
286
287 // Required Set of Parameters
288 // --------------------------
289 vector<t_pppParam*> required;
290
291 // Coordinates
292 // -----------
293 required.push_back(new t_pppParam(t_pppParam::crdX, t_prn(), t_lc::dummy));
294 required.push_back(new t_pppParam(t_pppParam::crdY, t_prn(), t_lc::dummy));
295 required.push_back(new t_pppParam(t_pppParam::crdZ, t_prn(), t_lc::dummy));
296
297 // Receiver Clock
298 // --------------
299 required.push_back(new t_pppParam(t_pppParam::clkR, t_prn(), t_lc::dummy));
300
301 // GPS-Glonass Clock Offset
302 // ------------------------
303 if (OPT->useGlonass()) {
304 required.push_back(new t_pppParam(t_pppParam::offGG, t_prn(), t_lc::dummy));
305 }
306
307 // Troposphere
308 // -----------
309 if (OPT->estTrp()) {
310 required.push_back(new t_pppParam(t_pppParam::trp, t_prn(), t_lc::dummy));
311 }
312
313 // Ambiguities
314 // -----------
315 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
316 const t_pppSatObs* satObs = obsVector[jj];
317 const vector<t_lc::type>& ambLCs = OPT->ambLCs(satObs->prn().system());
318 for (unsigned ii = 0; ii < ambLCs.size(); ii++) {
319 required.push_back(new t_pppParam(t_pppParam::amb, satObs->prn(), ambLCs[ii], &obsVector));
320 }
321 }
322
323 // Check if all required parameters are present
324 // --------------------------------------------
325 for (unsigned ii = 0; ii < required.size(); ii++) {
326 t_pppParam* parReq = required[ii];
327
328 bool found = false;
329 for (unsigned jj = 0; jj < _params.size(); jj++) {
330 t_pppParam* parOld = _params[jj];
331 if (parOld->isEqual(parReq)) {
332 found = true;
333 break;
334 }
335 }
336 if (found) {
337 delete parReq;
338 }
339 else {
340 _params.push_back(parReq);
341 }
342 }
343
344 // Set Parameter Indices
345 // ---------------------
346 sort(_params.begin(), _params.end(), t_pppParam::sortFunction);
347
348 for (unsigned ii = 0; ii < _params.size(); ii++) {
349 t_pppParam* par = _params[ii];
350 par->setIndex(ii);
351 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
352 const t_pppSatObs* satObs = obsVector[jj];
353 if (satObs->prn() == par->prn()) {
354 par->setAmbEleSat(satObs->eleSat());
355 par->stepAmbNumEpo();
356 }
357 }
358 }
359
360 return success;
361}
362
363//
364////////////////////////////////////////////////////////////////////////////
365void t_pppParlist::printResult(const bncTime& epoTime, const SymmetricMatrix& QQ,
366 const ColumnVector& xx) const {
367
368 string epoTimeStr = string(epoTime);
369
370 LOG << endl;
371
372 t_pppParam* parX = 0;
373 t_pppParam* parY = 0;
374 t_pppParam* parZ = 0;
375 for (unsigned ii = 0; ii < _params.size(); ii++) {
376 t_pppParam* par = _params[ii];
377 if (par->type() == t_pppParam::crdX) {
378 parX = par;
379 }
380 else if (par->type() == t_pppParam::crdY) {
381 parY = par;
382 }
383 else if (par->type() == t_pppParam::crdZ) {
384 parZ = par;
385 }
386 else {
387 int ind = par->indexNew();
388 LOG << epoTimeStr << ' ' << par->toString() << ' '
389 << setw(10) << setprecision(4) << par->x0() << ' '
390 << showpos << setw(10) << setprecision(4) << xx[ind] << noshowpos << " +- "
391 << setw(8) << setprecision(4) << sqrt(QQ[ind][ind]);
392 if (par->type() == t_pppParam::amb) {
393 LOG << " el = " << setw(6) << setprecision(2) << par->ambEleSat() * 180.0 / M_PI
394 << " epo = " << setw(4) << par->ambNumEpo();
395 }
396 LOG << endl;
397 }
398 }
399
400 if (parX && parY && parZ) {
401 const t_pppStation* sta = PPP_CLIENT->staRover();
402
403 ColumnVector xyz(3);
404 xyz[0] = xx[parX->indexNew()];
405 xyz[1] = xx[parY->indexNew()];
406 xyz[2] = xx[parZ->indexNew()];
407
408 ColumnVector neu(3);
409 xyz2neu(sta->ellApr().data(), xyz.data(), neu.data());
410
411 SymmetricMatrix QQxyz = QQ.SymSubMatrix(1,3);
412
413 SymmetricMatrix QQneu(3);
414 covariXYZ_NEU(QQxyz, sta->ellApr().data(), QQneu);
415
416 LOG << epoTimeStr << ' ' << sta->name()
417 << " X = " << setprecision(4) << sta->xyzApr()[0] + xyz[0] << " +- "
418 << setprecision(4) << sqrt(QQxyz[0][0])
419
420 << " Y = " << setprecision(4) << sta->xyzApr()[1] + xyz[1] << " +- "
421 << setprecision(4) << sqrt(QQxyz[1][1])
422
423 << " Z = " << setprecision(4) << sta->xyzApr()[2] + xyz[2] << " +- "
424 << setprecision(4) << sqrt(QQxyz[2][2])
425
426 << " dN = " << setprecision(4) << neu[0] << " +- "
427 << setprecision(4) << sqrt(QQneu[0][0])
428
429 << " dE = " << setprecision(4) << neu[1] << " +- "
430 << setprecision(4) << sqrt(QQneu[1][1])
431
432 << " dU = " << setprecision(4) << neu[2] << " +- "
433 << setprecision(4) << sqrt(QQneu[2][2])
434
435 << endl;
436 }
437}
438
Note: See TracBrowser for help on using the repository browser.