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

Last change on this file since 7971 was 7626, checked in by stuerze, 9 years ago

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