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

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