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

Last change on this file since 9386 was 9386, checked in by stuerze, 3 years ago

update regarding PPP
  • Property svn:keywords set to Author Date Id Rev URL;svn:eol-style=native
  • Property svn:mime-type set to text/plain
File size: 19.8 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->_aprSigClk;
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 offGR = 0;
77 if (_prn.system() == 'R' && tLC != t_lc::MW) {
78 offGR = PPP_CLIENT->offGR();
79 }
80 double offGE = 0;
81 if (_prn.system() == 'E' && tLC != t_lc::MW) {
82 offGE = PPP_CLIENT->offGE();
83 }
84 double offGC = 0;
85 if (_prn.system() == 'C' && tLC != t_lc::MW) {
86 offGC = PPP_CLIENT->offGC();
87 }
88 _x0 = floor((obs->obsValue(tLC) - offGR -offGE - offGC - obs->cmpValue(tLC)) / obs->lambda(tLC) + 0.5);
89 break;
90 }
91 }
92 }
93 break;
94 case offGR:
95 _epoSpec = true;
96 _sigma0 = OPT->_aprSigOGR;
97 _x0 = PPP_CLIENT->offGR();
98 break;
99 case offGE:
100 _epoSpec = true;
101 _sigma0 = OPT->_aprSigOGE;
102 _x0 = PPP_CLIENT->offGE();
103 break;
104 case offGC:
105 _epoSpec = true;
106 _sigma0 = OPT->_aprSigOGC;
107 _x0 = PPP_CLIENT->offGC();
108 break;
109 case trp:
110 _epoSpec = false;
111 _sigma0 = OPT->_aprSigTrp;
112 _noise = OPT->_noiseTrp;
113 break;
114 case ion:
115 _epoSpec = true;
116 _sigma0 = OPT->_aprSigIon;
117 break;
118 case cBiasG1: case cBiasR1: case cBiasE1: case cBiasC1:
119 case cBiasG2: case cBiasR2: case cBiasE2: case cBiasC2:
120 _epoSpec = true;
121 _sigma0 = OPT->_aprSigCodeBias;
122 break;
123 case pBiasG1: case pBiasR1: case pBiasE1: case pBiasC1:
124 case pBiasG2: case pBiasR2: case pBiasE2: case pBiasC2:
125 _epoSpec = true;
126 _sigma0 = OPT->_aprSigPhaseBias;
127 break;
128 }
129}
130
131// Destructor
132////////////////////////////////////////////////////////////////////////////
133t_pppParam::~t_pppParam() {
134 delete _ambInfo;
135}
136
137//
138////////////////////////////////////////////////////////////////////////////
139double t_pppParam::partial(const bncTime& /* epoTime */, const t_pppSatObs* obs,
140 const t_lc::type& tLC, const t_prn refPrn) const {
141
142 // Special Case - Melbourne-Wuebbena
143 // ---------------------------------
144 if (tLC == t_lc::MW && _type != amb) {
145 return 0.0;
146 }
147
148 const t_pppStation* sta = PPP_CLIENT->staRover();
149 ColumnVector rhoV = sta->xyzApr() - obs->xc().Rows(1,3);
150
151 map<t_frequency::type, double> codeCoeff;
152 map<t_frequency::type, double> phaseCoeff;
153 map<t_frequency::type, double> ionoCoeff;
154 obs->lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
155
156 switch (_type) {
157 case crdX:
158 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
159 return (sta->xyzApr()[0] - obs->xc()[0]) / rhoV.NormFrobenius();
160 case crdY:
161 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
162 return (sta->xyzApr()[1] - obs->xc()[1]) / rhoV.NormFrobenius();
163 case crdZ:
164 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
165 return (sta->xyzApr()[2] - obs->xc()[2]) / rhoV.NormFrobenius();
166 case clkR:
167 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
168 return 1.0;
169 case offGR:
170 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
171 return (obs->prn().system() == 'R') ? 1.0 : 0.0;
172 case offGE:
173 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
174 return (obs->prn().system() == 'E') ? 1.0 : 0.0;
175 case offGC:
176 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
177 return (obs->prn().system() == 'C') ? 1.0 : 0.0;
178 case amb:
179 if (tLC == t_lc::GIM || tLC == t_lc::Tz0) {return 0.0;}
180 else if ((OPT->_obsModelType == OPT->IF) ||
181 (OPT->_obsModelType == OPT->PPPRTK) ||
182 (OPT->_obsModelType == OPT->UncombPPP) ||
183 (OPT->_obsModelType == OPT->DCMcodeBias && !obs->isReference()) ||
184 (OPT->_obsModelType == OPT->DCMphaseBias && !obs->isReference()) ) {
185
186 if (obs->prn() == _prn) {
187 if (tLC == _tLC) {
188 return (obs->lambda(tLC));
189 }
190 else if (tLC == t_lc::lIF && _tLC == t_lc::MW) {
191 return obs->lambda(t_lc::lIF) * obs->lambda(t_lc::MW) / obs->lambda(t_lc::l2);
192 }
193 else {
194 if (_tLC == t_lc::l1) {
195 return obs->lambda(t_lc::l1) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
196 }
197 else if (_tLC == t_lc::l2) {
198 return obs->lambda(t_lc::l2) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
199 }
200 }
201 }
202 }
203 break;
204 case trp:
205 if (tLC == t_lc::GIM) {
206 return 0.0;
207 }
208 else if (tLC == t_lc::Tz0) {
209 return 1.0;
210 }
211 else {
212 return 1.0 / sin(obs->eleSat());
213 }
214 case ion:
215 if (obs->prn() == _prn) {
216 if (tLC == t_lc::c1) {
217 return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::c1)];
218 }
219 else if (tLC == t_lc::c2) {
220 return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::c2)];
221 }
222 else if (tLC == t_lc::l1) {
223 return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
224 }
225 else if (tLC == t_lc::l2) {
226 return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
227 }
228 else if (tLC == t_lc::GIM) {
229 return -1.0;
230 }
231 }
232 if (tLC == t_lc::GIM && _prn == refPrn) {
233 return 1.0;
234 }
235 break;
236 case cBiasG1:
237 if ((obs->prn().system() == 'G') && (tLC == t_lc::c1)) {return 1.0;} else {return 0.0;}
238 break;
239 case cBiasR1:
240 if ((obs->prn().system() == 'R') && (tLC == t_lc::c1)) {return 1.0;} else {return 0.0;}
241 break;
242 case cBiasE1:
243 if ((obs->prn().system() == 'E') && (tLC == t_lc::c1)) {return 1.0;} else {return 0.0;}
244 break;
245 case cBiasC1:
246 if ((obs->prn().system() == 'C') && (tLC == t_lc::c1)) {return 1.0;} else {return 0.0;}
247 break;
248 case cBiasG2:
249 if ((obs->prn().system() == 'G') && (tLC == t_lc::c2)) {return 1.0;} else {return 0.0;}
250 break;
251 case cBiasR2:
252 if ((obs->prn().system() == 'R') && (tLC == t_lc::c2)) {return 1.0;} else {return 0.0;}
253 break;
254 case cBiasE2:
255 if ((obs->prn().system() == 'E') && (tLC == t_lc::c2)) {return 1.0;} else {return 0.0;}
256 break;
257 case cBiasC2:
258 if ((obs->prn().system() == 'C') && (tLC == t_lc::c2)) {return 1.0;} else {return 0.0;}
259 break;
260 case pBiasG1:
261 if ((obs->prn().system() == 'G') && (tLC == t_lc::l1)) {return 1.0;} else {return 0.0;}
262 break;
263 case pBiasR1:
264 if ((obs->prn().system() == 'R') && (tLC == t_lc::l1)) {return 1.0;} else {return 0.0;}
265 break;
266 case pBiasE1:
267 if ((obs->prn().system() == 'E') && (tLC == t_lc::l1)) {return 1.0;} else {return 0.0;}
268 break;
269 case pBiasC1:
270 if ((obs->prn().system() == 'C') && (tLC == t_lc::l1)) {return 1.0;} else {return 0.0;}
271 break;
272 case pBiasG2:
273 if ((obs->prn().system() == 'G') && (tLC == t_lc::l2)) {return 1.0;} else {return 0.0;}
274 break;
275 case pBiasR2:
276 if ((obs->prn().system() == 'R') && (tLC == t_lc::l2)) {return 1.0;} else {return 0.0;}
277 break;
278 case pBiasE2:
279 if ((obs->prn().system() == 'E') && (tLC == t_lc::l2)) {return 1.0;} else {return 0.0;}
280 break;
281 case pBiasC2:
282 if ((obs->prn().system() == 'C') && (tLC == t_lc::l2)) {return 1.0;} else {return 0.0;}
283 break;
284
285 }
286 return 0.0;
287}
288
289//
290////////////////////////////////////////////////////////////////////////////
291string t_pppParam::toString() const {
292 stringstream ss;
293 switch (_type) {
294 case crdX:
295 ss << "CRD_X";
296 break;
297 case crdY:
298 ss << "CRD_Y";
299 break;
300 case crdZ:
301 ss << "CRD_Z";
302 break;
303 case clkR:
304 ss << "REC_CLK ";
305 break;
306 case offGR:
307 ss << "OGR ";
308 break;
309 case offGE:
310 ss << "OGE ";
311 break;
312 case offGC:
313 ss << "OGC ";
314 break;
315 case trp:
316 ss << "TRP ";
317 break;
318 case amb:
319 ss << "AMB " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
320 break;
321 case ion:
322 ss << "ION " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
323 break;
324 case cBiasG1: case pBiasG1:
325 case cBiasG2: case pBiasG2:
326 ss << "BIA " << left << setw(3) << t_lc::toString(_tLC) << right << " G ";
327 break;
328 case cBiasR1: case pBiasR1:
329 case cBiasR2: case pBiasR2:
330 ss << "BIA " << left << setw(3) << t_lc::toString(_tLC) << right << " R ";
331 break;
332 case cBiasE1: case pBiasE1:
333 case cBiasE2: case pBiasE2:
334 ss << "BIA " << left << setw(3) << t_lc::toString(_tLC) << right << " E ";
335 break;
336 case cBiasC1: case pBiasC1:
337 case cBiasC2: case pBiasC2:
338 ss << "BIA " << left << setw(3) << t_lc::toString(_tLC) << right << " C ";
339 break;
340 }
341 return ss.str();
342}
343
344// Constructor
345////////////////////////////////////////////////////////////////////////////
346t_pppParlist::t_pppParlist() {
347}
348
349// Destructor
350////////////////////////////////////////////////////////////////////////////
351t_pppParlist::~t_pppParlist() {
352 for (unsigned ii = 0; ii < _params.size(); ii++) {
353 delete _params[ii];
354 }
355}
356
357//
358////////////////////////////////////////////////////////////////////////////
359t_irc t_pppParlist::set(const bncTime& epoTime,
360 const std::vector<t_pppSatObs*>& obsVector,
361 const QMap<char, t_pppRefSat*>& refSatMap) {
362
363 // Remove some Parameters
364 // ----------------------
365 vector<t_pppParam*>::iterator it = _params.begin();
366 while (it != _params.end()) {
367 t_pppParam* par = *it;
368
369 bool remove = false;
370
371 if (par->epoSpec()) {
372 remove = true;
373 }
374
375 else if (par->type() == t_pppParam::crdX ||
376 par->type() == t_pppParam::crdY ||
377 par->type() == t_pppParam::crdZ) {
378 if (par->lastObsTime().valid() && (epoTime - par->lastObsTime() > 60.0)) {
379 remove = true;
380 }
381 }
382 else if (par->type() == t_pppParam::amb) {
383 char system = par->prn().system();
384 t_prn refPrn = t_prn();
385 if (OPT->_refSatRequired) {
386 refPrn = (refSatMap[system])->prn();
387 }
388 if ((par->lastObsTime().valid() &&
389 (epoTime - par->lastObsTime() > 60.0)) ||
390 (refPrn == par->prn())) {
391 remove = true;
392 }
393 }
394
395 if (remove) {
396 delete par;
397 it = _params.erase(it);
398 }
399 else {
400 ++it;
401 }
402 }
403
404 // Check whether parameters have observations
405 // ------------------------------------------
406 for (unsigned ii = 0; ii < _params.size(); ii++) {
407 t_pppParam* par = _params[ii];
408 if (par->prn() == 0) {
409 par->setLastObsTime(epoTime);
410 if (par->firstObsTime().undef()) {
411 par->setFirstObsTime(epoTime);
412 }
413 }
414 else {
415 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
416 const t_pppSatObs* satObs = obsVector[jj];
417 if (satObs->prn() == par->prn()) {
418 par->setLastObsTime(epoTime);
419 if (par->firstObsTime().undef()) {
420 par->setFirstObsTime(epoTime);
421 }
422 break;
423 }
424 }
425 }
426 }
427
428 // Required Set of Parameters
429 // --------------------------
430 vector<t_pppParam*> required;
431
432 // Coordinates
433 // -----------
434 required.push_back(new t_pppParam(t_pppParam::crdX, t_prn(), t_lc::dummy));
435 required.push_back(new t_pppParam(t_pppParam::crdY, t_prn(), t_lc::dummy));
436 required.push_back(new t_pppParam(t_pppParam::crdZ, t_prn(), t_lc::dummy));
437
438 // Receiver Code Biases
439 // --------------------
440 if (OPT->_obsModelType == OPT->DCMcodeBias) {
441 if (OPT->useSystem('G')) {
442 required.push_back(new t_pppParam(t_pppParam::cBiasG1, t_prn(), t_lc::c1));
443 required.push_back(new t_pppParam(t_pppParam::cBiasG2, t_prn(), t_lc::c2));
444 }
445 if (OPT->useSystem('R')) {
446 required.push_back(new t_pppParam(t_pppParam::cBiasR1, t_prn(), t_lc::c1));
447 required.push_back(new t_pppParam(t_pppParam::cBiasR2, t_prn(), t_lc::c2));
448 }
449 if (OPT->useSystem('E')) {
450 required.push_back(new t_pppParam(t_pppParam::cBiasE1, t_prn(), t_lc::c1));
451 required.push_back(new t_pppParam(t_pppParam::cBiasE2, t_prn(), t_lc::c2));
452 }
453 if (OPT->useSystem('C')) {
454 required.push_back(new t_pppParam(t_pppParam::cBiasC1, t_prn(), t_lc::c1));
455 required.push_back(new t_pppParam(t_pppParam::cBiasC2, t_prn(), t_lc::c2));
456 }
457 }
458
459 // Receiver Phase Biases
460 // ---------------------
461 if ((OPT->_obsModelType == OPT->DCMphaseBias) ||
462 (OPT->_obsModelType == OPT->PPPRTK) ) {
463 if (OPT->useSystem('G')) {
464 required.push_back(new t_pppParam(t_pppParam::pBiasG1, t_prn(), t_lc::l1));
465 required.push_back(new t_pppParam(t_pppParam::pBiasG2, t_prn(), t_lc::l2));
466 }
467 if (OPT->useSystem('R')) {
468 required.push_back(new t_pppParam(t_pppParam::pBiasR1, t_prn(), t_lc::l1));
469 required.push_back(new t_pppParam(t_pppParam::pBiasR2, t_prn(), t_lc::l2));
470 }
471 if (OPT->useSystem('E')) {
472 required.push_back(new t_pppParam(t_pppParam::pBiasE1, t_prn(), t_lc::l1));
473 required.push_back(new t_pppParam(t_pppParam::pBiasE2, t_prn(), t_lc::l2));
474 }
475 if (OPT->useSystem('C')) {
476 required.push_back(new t_pppParam(t_pppParam::pBiasC1, t_prn(), t_lc::l1));
477 required.push_back(new t_pppParam(t_pppParam::pBiasC2, t_prn(), t_lc::l2));
478 }
479 }
480 // Receiver Clock
481 // --------------
482 required.push_back(new t_pppParam(t_pppParam::clkR, t_prn(), t_lc::dummy));
483
484 // GPS-GLONASS Clock Offset
485 // ------------------------
486 if (OPT->useSystem('R') && OPT->useSystem('G')) {
487 required.push_back(new t_pppParam(t_pppParam::offGR, t_prn(), t_lc::dummy));
488 }
489
490 // GPS-Galileo Clock Offset
491 // ------------------------
492 if (OPT->useSystem('E') && OPT->useSystem('G')) {
493 required.push_back(new t_pppParam(t_pppParam::offGE, t_prn(), t_lc::dummy));
494 }
495
496 // GPS-BDS Clock Offset
497 // ------------------------
498 if (OPT->useSystem('C') && OPT->useSystem('G')) {
499 required.push_back(new t_pppParam(t_pppParam::offGC, t_prn(), t_lc::dummy));
500 }
501
502 // Troposphere
503 // -----------
504 if (OPT->estTrp()) {
505 required.push_back(new t_pppParam(t_pppParam::trp, t_prn(), t_lc::dummy));
506 }
507
508 // Ionosphere
509 // ----------
510 if (OPT->_obsModelType == OPT->UncombPPP ||
511 OPT->_obsModelType == OPT->DCMcodeBias ||
512 OPT->_obsModelType == OPT->DCMphaseBias ) {
513 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
514 const t_pppSatObs* satObs = obsVector[jj];
515 required.push_back(new t_pppParam(t_pppParam::ion, satObs->prn(), t_lc::dummy));
516 }
517 }
518 // Ambiguities
519 // -----------
520 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
521 const t_pppSatObs* satObs = obsVector[jj];
522 if ((OPT->_obsModelType == OPT->IF) ||
523 (OPT->_obsModelType == OPT->PPPRTK) ||
524 (OPT->_obsModelType == OPT->UncombPPP) ||
525 (OPT->_obsModelType == OPT->DCMcodeBias && !satObs->isReference()) ||
526 (OPT->_obsModelType == OPT->DCMphaseBias && !satObs->isReference()) ) {
527 const vector<t_lc::type>& ambLCs = OPT->ambLCs(satObs->prn().system());
528 for (unsigned ii = 0; ii < ambLCs.size(); ii++) {
529 required.push_back(new t_pppParam(t_pppParam::amb, satObs->prn(), ambLCs[ii], &obsVector));
530 }
531 }
532 }
533
534 // Check if all required parameters are present
535 // --------------------------------------------
536 for (unsigned ii = 0; ii < required.size(); ii++) {
537 t_pppParam* parReq = required[ii];
538
539 bool found = false;
540 for (unsigned jj = 0; jj < _params.size(); jj++) {
541 t_pppParam* parOld = _params[jj];
542 if (parOld->isEqual(parReq)) {
543 found = true;
544 break;
545 }
546 }
547 if (found) {
548 delete parReq;
549 }
550 else {
551 _params.push_back(parReq);
552 }
553 }
554
555 // Set Parameter Indices
556 // ---------------------
557 sort(_params.begin(), _params.end(), t_pppParam::sortFunction);
558
559 for (unsigned ii = 0; ii < _params.size(); ii++) {
560 t_pppParam* par = _params[ii];
561 par->setIndex(ii);
562 for (unsigned jj = 0; jj < obsVector.size(); jj++) {
563 const t_pppSatObs* satObs = obsVector[jj];
564 if (satObs->prn() == par->prn()) {
565 par->setAmbEleSat(satObs->eleSat());
566 par->stepAmbNumEpo();
567 }
568 }
569 }
570
571 return success;
572}
573
574//
575////////////////////////////////////////////////////////////////////////////
576void t_pppParlist::printResult(const bncTime& epoTime, const SymmetricMatrix& QQ,
577 const ColumnVector& xx) const {
578
579 string epoTimeStr = string(epoTime);
580 const t_pppStation* sta = PPP_CLIENT->staRover();
581
582 LOG << endl;
583
584 t_pppParam* parX = 0;
585 t_pppParam* parY = 0;
586 t_pppParam* parZ = 0;
587 for (unsigned ii = 0; ii < _params.size(); ii++) {
588 t_pppParam* par = _params[ii];
589 if (par->type() == t_pppParam::crdX) {
590 parX = par;
591 }
592 else if (par->type() == t_pppParam::crdY) {
593 parY = par;
594 }
595 else if (par->type() == t_pppParam::crdZ) {
596 parZ = par;
597 }
598 else {
599 int ind = par->indexNew();
600 double apr = (par->type() == t_pppParam::trp) ?
601 t_tropo::delay_saast(sta->xyzApr(), M_PI/2.0) : par->x0();
602 LOG << epoTimeStr << ' ' << par->toString() << ' '
603 << setw(10) << setprecision(4) << apr << ' '
604 << showpos << setw(10) << setprecision(4) << xx[ind] << noshowpos << " +- "
605 << setw(8) << setprecision(4) << sqrt(QQ[ind][ind]);
606 if (par->type() == t_pppParam::amb) {
607 LOG << " el = " << setw(6) << setprecision(2) << par->ambEleSat() * 180.0 / M_PI
608 << " epo = " << setw(4) << par->ambNumEpo();
609 }
610 LOG << endl;
611 }
612 }
613
614 if (parX && parY && parZ) {
615
616 ColumnVector xyz(3);
617 xyz[0] = xx[parX->indexNew()];
618 xyz[1] = xx[parY->indexNew()];
619 xyz[2] = xx[parZ->indexNew()];
620
621 ColumnVector neu(3);
622 xyz2neu(sta->ellApr().data(), xyz.data(), neu.data());
623
624 SymmetricMatrix QQxyz = QQ.SymSubMatrix(1,3);
625
626 SymmetricMatrix QQneu(3);
627 covariXYZ_NEU(QQxyz, sta->ellApr().data(), QQneu);
628
629 LOG << epoTimeStr << ' ' << sta->name()
630 << " X = " << setprecision(4) << sta->xyzApr()[0] + xyz[0] << " +- "
631 << setprecision(4) << sqrt(QQxyz[0][0])
632
633 << " Y = " << setprecision(4) << sta->xyzApr()[1] + xyz[1] << " +- "
634 << setprecision(4) << sqrt(QQxyz[1][1])
635
636 << " Z = " << setprecision(4) << sta->xyzApr()[2] + xyz[2] << " +- "
637 << setprecision(4) << sqrt(QQxyz[2][2])
638
639 << " dN = " << setprecision(4) << neu[0] << " +- "
640 << setprecision(4) << sqrt(QQneu[0][0])
641
642 << " dE = " << setprecision(4) << neu[1] << " +- "
643 << setprecision(4) << sqrt(QQneu[1][1])
644
645 << " dU = " << setprecision(4) << neu[2] << " +- "
646 << setprecision(4) << sqrt(QQneu[2][2])
647
648 << endl;
649 }
650}
651
Note: See TracBrowser for help on using the repository browser.