source: ntrip/trunk/BNC/src/PPP/pppSatObs.cpp@ 7302

Last change on this file since 7302 was 7288, checked in by stuerze, 9 years ago

phase biases 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: 16.4 KB
Line 
1/* -------------------------------------------------------------------------
2 * BKG NTRIP Client
3 * -------------------------------------------------------------------------
4 *
5 * Class: t_pppSatObs
6 *
7 * Purpose: Satellite observations
8 *
9 * Author: L. Mervart
10 *
11 * Created: 29-Jul-2014
12 *
13 * Changes:
14 *
15 * -----------------------------------------------------------------------*/
16
17
18#include <iostream>
19#include <cmath>
20#include <newmatio.h>
21
22#include "pppSatObs.h"
23#include "bncconst.h"
24#include "pppEphPool.h"
25#include "pppStation.h"
26#include "bncutils.h"
27#include "bncantex.h"
28#include "pppObsPool.h"
29#include "pppClient.h"
30
31using namespace BNC_PPP;
32using namespace std;
33
34// Constructor
35////////////////////////////////////////////////////////////////////////////
36t_pppSatObs::t_pppSatObs(const t_satObs& pppSatObs) {
37 _prn = pppSatObs._prn;
38 _time = pppSatObs._time;
39 _outlier = false;
40 _valid = true;
41 for (unsigned ii = 0; ii < t_frequency::max; ii++) {
42 _obs[ii] = 0;
43 }
44 prepareObs(pppSatObs);
45}
46
47// Destructor
48////////////////////////////////////////////////////////////////////////////
49t_pppSatObs::~t_pppSatObs() {
50 for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
51 delete _obs[iFreq];
52 }
53}
54
55//
56////////////////////////////////////////////////////////////////////////////
57void t_pppSatObs::prepareObs(const t_satObs& pppSatObs) {
58
59 _model.reset();
60
61 // Select pseudoranges and phase observations
62 // ------------------------------------------
63 const string preferredAttrib = "CWPXI_";
64
65 for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
66 string frqNum = t_frequency::toString(t_frequency::type(iFreq)).substr(1);
67 for (unsigned iPref = 0; iPref < preferredAttrib.length(); iPref++) {
68 string obsType = (preferredAttrib[iPref] == '_') ? frqNum : frqNum + preferredAttrib[iPref];
69 if (_obs[iFreq] == 0) {
70 for (unsigned ii = 0; ii < pppSatObs._obs.size(); ii++) {
71 const t_frqObs* obs = pppSatObs._obs[ii];
72 if (obs->_rnxType2ch == obsType && obs->_codeValid && obs->_phaseValid) {
73 _obs[iFreq] = new t_frqObs(*obs);
74 }
75 }
76 }
77 }
78 }
79
80 // Used frequency types
81 // --------------------
82 _fType1 = t_lc::toFreq(_prn.system(),t_lc::l1);
83 _fType2 = t_lc::toFreq(_prn.system(),t_lc::l2);
84
85 // Check whether all required frequencies available
86 // ------------------------------------------------
87 for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) {
88 t_lc::type tLC = OPT->LCs(_prn.system())[ii];
89 if (!isValid(tLC)) {
90 _valid = false;
91 return;
92 }
93 }
94
95 // Find Glonass Channel Number
96 // ---------------------------
97 if (_prn.system() == 'R') {
98 _channel = PPP_CLIENT->ephPool()->getChannel(_prn);
99 }
100 else {
101 _channel = 0;
102 }
103
104 // Compute Satellite Coordinates at Time of Transmission
105 // -----------------------------------------------------
106 _xcSat.ReSize(4); _xcSat = 0.0;
107 _vvSat.ReSize(4); _vvSat = 0.0;
108 bool totOK = false;
109 ColumnVector satPosOld(4); satPosOld = 0.0;
110 t_lc::type tLC = isValid(t_lc::cIF) ? t_lc::cIF : t_lc::c1;
111 double prange = obsValue(tLC);
112 for (int ii = 1; ii <= 10; ii++) {
113 bncTime ToT = _time - prange / t_CST::c - _xcSat[3];
114 if (PPP_CLIENT->ephPool()->getCrd(_prn, ToT, _xcSat, _vvSat) != success) {
115 _valid = false;
116 return;
117 }
118 ColumnVector dx = _xcSat - satPosOld;
119 dx[3] *= t_CST::c;
120 if (dx.norm_Frobenius() < 1.e-4) {
121 totOK = true;
122 break;
123 }
124 satPosOld = _xcSat;
125 }
126 if (totOK) {
127 _signalPropagationTime = prange / t_CST::c - _xcSat[3];
128 _model._satClkM = _xcSat[3] * t_CST::c;
129 }
130 else {
131 _valid = false;
132 }
133}
134
135//
136////////////////////////////////////////////////////////////////////////////
137void t_pppSatObs::lcCoeff(t_lc::type tLC,
138 map<t_frequency::type, double>& codeCoeff,
139 map<t_frequency::type, double>& phaseCoeff) const {
140
141 codeCoeff.clear();
142 phaseCoeff.clear();
143
144 double f1 = t_CST::freq(_fType1, _channel);
145 double f2 = t_CST::freq(_fType2, _channel);
146
147 switch (tLC) {
148 case t_lc::l1:
149 phaseCoeff[_fType1] = 1.0;
150 return;
151 case t_lc::l2:
152 phaseCoeff[_fType2] = 1.0;
153 return;
154 case t_lc::lIF:
155 phaseCoeff[_fType1] = f1 * f1 / (f1 * f1 - f2 * f2);
156 phaseCoeff[_fType2] = -f2 * f2 / (f1 * f1 - f2 * f2);
157 return;
158 case t_lc::MW:
159 phaseCoeff[_fType1] = f1 / (f1 - f2);
160 phaseCoeff[_fType2] = -f2 / (f1 - f2);
161 codeCoeff[_fType1] = -f1 / (f1 + f2);
162 codeCoeff[_fType2] = -f2 / (f1 + f2);
163 return;
164 case t_lc::CL:
165 phaseCoeff[_fType1] = 0.5;
166 codeCoeff[_fType1] = 0.5;
167 return;
168 case t_lc::c1:
169 codeCoeff[_fType1] = 1.0;
170 return;
171 case t_lc::c2:
172 codeCoeff[_fType2] = 1.0;
173 return;
174 case t_lc::cIF:
175 codeCoeff[_fType1] = f1 * f1 / (f1 * f1 - f2 * f2);
176 codeCoeff[_fType2] = -f2 * f2 / (f1 * f1 - f2 * f2);
177 return;
178 case t_lc::dummy:
179 case t_lc::maxLc:
180 return;
181 }
182}
183
184//
185////////////////////////////////////////////////////////////////////////////
186bool t_pppSatObs::isValid(t_lc::type tLC) const {
187 bool valid = true;
188 obsValue(tLC, &valid);
189 return valid;
190}
191//
192////////////////////////////////////////////////////////////////////////////
193double t_pppSatObs::obsValue(t_lc::type tLC, bool* valid) const {
194
195 map<t_frequency::type, double> codeCoeff;
196 map<t_frequency::type, double> phaseCoeff;
197 lcCoeff(tLC, codeCoeff, phaseCoeff);
198
199 double retVal = 0.0;
200 if (valid) *valid = true;
201
202 map<t_frequency::type, double>::const_iterator it;
203 for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
204 t_frequency::type tFreq = it->first;
205 if (_obs[tFreq] == 0) {
206 if (valid) *valid = false;
207 return 0.0;
208 }
209 else {
210 retVal += it->second * _obs[tFreq]->_code;
211 }
212 }
213 for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
214 t_frequency::type tFreq = it->first;
215 if (_obs[tFreq] == 0) {
216 if (valid) *valid = false;
217 return 0.0;
218 }
219 else {
220 retVal += it->second * _obs[tFreq]->_phase * t_CST::lambda(tFreq, _channel);
221 }
222 }
223
224 return retVal;
225}
226
227//
228////////////////////////////////////////////////////////////////////////////
229double t_pppSatObs::lambda(t_lc::type tLC) const {
230
231 double f1 = t_CST::freq(_fType1, _channel);
232 double f2 = t_CST::freq(_fType2, _channel);
233
234 if (tLC == t_lc::l1) {
235 return t_CST::c / f1;
236 }
237 else if (tLC == t_lc::l2) {
238 return t_CST::c / f2;
239 }
240 else if (tLC == t_lc::lIF) {
241 return t_CST::c / (f1 + f2);
242 }
243 else if (tLC == t_lc::MW) {
244 return t_CST::c / (f1 - f2);
245 }
246 else if (tLC == t_lc::CL) {
247 return t_CST::c / f1 / 2.0;
248 }
249
250 return 0.0;
251}
252
253//
254////////////////////////////////////////////////////////////////////////////
255double t_pppSatObs::sigma(t_lc::type tLC) const {
256
257 map<t_frequency::type, double> codeCoeff;
258 map<t_frequency::type, double> phaseCoeff;
259 lcCoeff(tLC, codeCoeff, phaseCoeff);
260
261 double retVal = 0.0;
262
263 map<t_frequency::type, double>::const_iterator it;
264 for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
265 retVal += it->second * it->second * OPT->_sigmaC1 * OPT->_sigmaC1;
266 }
267 for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
268 retVal += it->second * it->second * OPT->_sigmaL1 * OPT->_sigmaL1;
269 }
270
271 retVal = sqrt(retVal);
272
273 // De-Weight GLONASS
274 // -----------------
275 if (_prn.system() == 'R') {
276 retVal *= 5.0;
277 }
278
279 // Elevation-Dependent Weighting
280 // -----------------------------
281 double cEle = 1.0;
282 if ( (OPT->_eleWgtCode && t_lc::includesCode(tLC)) ||
283 (OPT->_eleWgtPhase && t_lc::includesPhase(tLC)) ) {
284 double eleD = eleSat()*180.0/M_PI;
285 double hlp = fabs(90.0 - eleD);
286 cEle = (1.0 + hlp*hlp*hlp*0.000004);
287 }
288
289 return cEle * retVal;
290}
291
292//
293////////////////////////////////////////////////////////////////////////////
294double t_pppSatObs::maxRes(t_lc::type tLC) const {
295
296 map<t_frequency::type, double> codeCoeff;
297 map<t_frequency::type, double> phaseCoeff;
298 lcCoeff(tLC, codeCoeff, phaseCoeff);
299
300 double retVal = 0.0;
301
302 map<t_frequency::type, double>::const_iterator it;
303 for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
304 retVal += it->second * it->second * OPT->_maxResC1 * OPT->_maxResC1;
305 }
306 for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
307 retVal += it->second * it->second * OPT->_maxResL1 * OPT->_maxResL1;
308 }
309
310 return sqrt(retVal);
311}
312
313
314//
315////////////////////////////////////////////////////////////////////////////
316t_irc t_pppSatObs::cmpModel(const t_pppStation* station) {
317
318 // Reset all model values
319 // ----------------------
320 _model.reset();
321
322 // Topocentric Satellite Position
323 // ------------------------------
324 ColumnVector rSat = _xcSat.Rows(1,3);
325 ColumnVector rhoV = rSat - station->xyzApr();
326 _model._rho = rhoV.norm_Frobenius();
327
328 ColumnVector neu(3);
329 xyz2neu(station->ellApr().data(), rhoV.data(), neu.data());
330
331 _model._eleSat = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / _model._rho );
332 if (neu[2] < 0) {
333 _model._eleSat *= -1.0;
334 }
335 _model._azSat = atan2(neu[1], neu[0]);
336
337 // Satellite Clocks
338 // ----------------
339 _model._satClkM = _xcSat[3] * t_CST::c;
340
341 // Receiver Clocks
342 // ---------------
343 _model._recClkM = station->dClk() * t_CST::c;
344
345 // Sagnac Effect (correction due to Earth rotation)
346 // ------------------------------------------------
347 ColumnVector Omega(3);
348 Omega[0] = 0.0;
349 Omega[1] = 0.0;
350 Omega[2] = t_CST::omega / t_CST::c;
351 _model._sagnac = DotProduct(Omega, crossproduct(rSat, station->xyzApr()));
352
353 // Antenna Eccentricity
354 // --------------------
355 _model._antEcc = -DotProduct(station->xyzEcc(), rhoV) / _model._rho;
356
357 // Antenna Phase Center Offsets and Variations
358 // -------------------------------------------
359 if (PPP_CLIENT->antex()) {
360 for (unsigned ii = 0; ii < t_frequency::max; ii++) {
361 t_frequency::type frqType = static_cast<t_frequency::type>(ii);
362 bool found;
363 _model._antPCO[ii] = PPP_CLIENT->antex()->rcvCorr(station->antName(), frqType,
364 _model._eleSat, _model._azSat, found);
365 }
366 }
367
368 // Tropospheric Delay
369 // ------------------
370 _model._tropo = t_tropo::delay_saast(station->xyzApr(), _model._eleSat);
371
372 // Phase Wind-Up
373 // -------------
374 _model._windUp = station->windUp(_time, _prn, rSat);
375
376 // Code Biases
377 // -----------
378 const t_satCodeBias* satCodeBias = PPP_CLIENT->obsPool()->satCodeBias(_prn);
379 if (satCodeBias) {
380 for (unsigned ii = 0; ii < satCodeBias->_bias.size(); ii++) {
381 const t_frqCodeBias& bias = satCodeBias->_bias[ii];
382 for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
383 const t_frqObs* obs = _obs[iFreq];
384 if (obs && obs->_rnxType2ch == bias._rnxType2ch) {
385 _model._codeBias[iFreq] = bias._value;
386 }
387 }
388 }
389 }
390
391 // Phase Biases
392 // -----------
393 // TODO: consideration of fix indicators, yaw angle and jump counter
394 const t_satPhaseBias* satPhaseBias = PPP_CLIENT->obsPool()->satPhaseBias(_prn);
395 if (satPhaseBias) {
396 for (unsigned ii = 0; ii < satPhaseBias->_bias.size(); ii++) {
397 const t_frqPhaseBias& bias = satPhaseBias->_bias[ii];
398 for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
399 const t_frqObs* obs = _obs[iFreq];
400 if (obs && obs->_rnxType2ch == bias._rnxType2ch) {
401 _model._phaseBias[iFreq] = bias._value;
402 }
403 }
404 }
405 }
406
407 // Tidal Correction
408 // ----------------
409 _model._tide = -DotProduct(station->tideDspl(), rhoV) / _model._rho;
410
411 // Ionospheric Delay
412 // -----------------
413 const t_vTec* vTec = PPP_CLIENT->obsPool()->vTec();
414 bool vTecUsage = true;
415 for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) {
416 t_lc::type tLC = OPT->LCs(_prn.system())[ii];
417 if (tLC == t_lc::cIF || tLC == t_lc::lIF) {
418 vTecUsage = false;
419 }
420 }
421 if (vTecUsage && vTec) {
422 double stec = station->stec(vTec, _signalPropagationTime, rSat);
423 for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
424 t_frequency::type frqType = static_cast<t_frequency::type>(iFreq);
425 _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(t_CST::freq(frqType, _channel), 2) * stec;
426 }
427 }
428
429 // Ocean Loading
430 // -------------
431 // TODO
432
433 // Set Model Set Flag
434 // ------------------
435 _model._set = true;
436
437 //printModel();
438
439 return success;
440}
441
442//
443////////////////////////////////////////////////////////////////////////////
444void t_pppSatObs::printModel() const {
445
446 LOG.setf(ios::fixed);
447 LOG << "MODEL for Satellite " << _prn.toString() << endl
448 << "RHO: " << setw(12) << setprecision(3) << _model._rho << endl
449 << "ELE: " << setw(12) << setprecision(3) << _model._eleSat * 180.0 / M_PI << endl
450 << "AZI: " << setw(12) << setprecision(3) << _model._azSat * 180.0 / M_PI << endl
451 << "SATCLK: " << setw(12) << setprecision(3) << _model._satClkM << endl
452 << "RECCLK: " << setw(12) << setprecision(3) << _model._recClkM << endl
453 << "SAGNAC: " << setw(12) << setprecision(3) << _model._sagnac << endl
454 << "ANTECC: " << setw(12) << setprecision(3) << _model._antEcc << endl
455 << "TROPO: " << setw(12) << setprecision(3) << _model._tropo << endl
456 << "WINDUP: " << setw(12) << setprecision(3) << _model._windUp << endl
457 << "TIDES: " << setw(12) << setprecision(3) << _model._tide << endl;
458 for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
459 if (_obs[iFreq]) {
460 string frqStr = t_frequency::toString(t_frequency::type(iFreq));
461 if (_prn.system() == frqStr[0]) {
462 LOG << "PCO : " << frqStr << setw(12) << setprecision(3) << _model._antPCO[iFreq] << endl
463 << "BIAS CODE : " << frqStr << setw(12) << setprecision(3) << _model._codeBias[iFreq] << endl
464 << "BIAS PHASE : " << frqStr << setw(12) << setprecision(3) << _model._phaseBias[iFreq] << endl
465 << "IONO CODEDELAY: " << frqStr << setw(12) << setprecision(3) << _model._ionoCodeDelay[iFreq] << endl;
466 }
467 }
468 }
469 for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) {
470 t_lc::type tLC = OPT->LCs(_prn.system())[ii];
471 LOG << "OBS-CMP " << t_lc::toString(tLC) << ": " << _prn.toString() << " "
472 << setw(12) << setprecision(3) << obsValue(tLC) << " "
473 << setw(12) << setprecision(3) << cmpValue(tLC) << " "
474 << setw(12) << setprecision(3) << obsValue(tLC) - cmpValue(tLC) << endl;
475
476 }
477 LOG << "OBS-CMP MW: " << _prn.toString() << " "
478 << setw(12) << setprecision(3) << obsValue(t_lc::MW) << " "
479 << setw(12) << setprecision(3) << cmpValue(t_lc::MW) << " "
480 << setw(12) << setprecision(3) << obsValue(t_lc::MW) - cmpValue(t_lc::MW) << endl;
481}
482
483//
484////////////////////////////////////////////////////////////////////////////
485double t_pppSatObs::cmpValueForBanc(t_lc::type tLC) const {
486 return cmpValue(tLC) - _model._rho - _model._sagnac - _model._recClkM;
487}
488
489//
490////////////////////////////////////////////////////////////////////////////
491double t_pppSatObs::cmpValue(t_lc::type tLC) const {
492
493 if (!isValid(tLC)) {
494 return 0.0;
495 }
496
497 // Non-Dispersive Part
498 // -------------------
499 double nonDisp = _model._rho + _model._recClkM - _model._satClkM
500 + _model._sagnac + _model._antEcc + _model._tropo
501 + _model._tide;
502
503 // Add Dispersive Part
504 // -------------------
505 map<t_frequency::type, double> codeCoeff;
506 map<t_frequency::type, double> phaseCoeff;
507 lcCoeff(tLC, codeCoeff, phaseCoeff);
508
509 double dispPart = 0.0;
510
511 map<t_frequency::type, double>::const_iterator it;
512 for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
513 t_frequency::type tFreq = it->first;
514 dispPart += it->second * (_model._antPCO[tFreq] - _model._codeBias[tFreq] +
515 _model._ionoCodeDelay[tFreq]);
516 }
517 for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
518 t_frequency::type tFreq = it->first;
519 dispPart += it->second * (_model._antPCO[tFreq] - _model._phaseBias[tFreq] +
520 _model._windUp * t_CST::lambda(tFreq, _channel) -
521 _model._ionoCodeDelay[tFreq]);
522 }
523
524 return nonDisp + dispPart;
525}
526
527//
528////////////////////////////////////////////////////////////////////////////
529void t_pppSatObs::setRes(t_lc::type tLC, double res) {
530 _res[tLC] = res;
531}
532
533//
534////////////////////////////////////////////////////////////////////////////
535double t_pppSatObs::getRes(t_lc::type tLC) const {
536 map<t_lc::type, double>::const_iterator it = _res.find(tLC);
537 if (it != _res.end()) {
538 return it->second;
539 }
540 else {
541 return 0.0;
542 }
543}
Note: See TracBrowser for help on using the repository browser.