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 <iomanip>
|
---|
20 | #include <cmath>
|
---|
21 | #include <newmatio.h>
|
---|
22 |
|
---|
23 | #include "pppSatObs.h"
|
---|
24 | #include "bncconst.h"
|
---|
25 | #include "pppEphPool.h"
|
---|
26 | #include "pppStation.h"
|
---|
27 | #include "bncutils.h"
|
---|
28 | #include "bncantex.h"
|
---|
29 | #include "pppObsPool.h"
|
---|
30 | #include "pppClient.h"
|
---|
31 |
|
---|
32 | using namespace BNC_PPP;
|
---|
33 | using namespace std;
|
---|
34 |
|
---|
35 | const double GLO_WEIGHT_FACTOR = 5.0;
|
---|
36 | const double BDS_WEIGHT_FACTOR = 2.0;
|
---|
37 |
|
---|
38 | // Constructor
|
---|
39 | ////////////////////////////////////////////////////////////////////////////
|
---|
40 | t_pppSatObs::t_pppSatObs(const t_satObs& pppSatObs) {
|
---|
41 | _prn = pppSatObs._prn;
|
---|
42 | _time = pppSatObs._time;
|
---|
43 | _outlier = false;
|
---|
44 | _valid = true;
|
---|
45 | _reference = false;
|
---|
46 | _stecRefSat = 0.0;
|
---|
47 | _stecSat = 0.0;
|
---|
48 | _signalPriorities = QString::fromStdString(OPT->_signalPriorities);
|
---|
49 | if (!_signalPriorities.size()) {
|
---|
50 | if (OPT->_obsModelType == OPT->DCMcodeBias ||
|
---|
51 | OPT->_obsModelType == OPT->DCMphaseBias) {
|
---|
52 | // at the moment only one code or phase bias per system (G,R,E,C)/modulation considered
|
---|
53 | _signalPriorities = "G:12&CW R:12&CP E:1&CX E:5&QX C:26&I";
|
---|
54 | }
|
---|
55 | else {
|
---|
56 | _signalPriorities = "G:12&CWPSLX R:12&CP E:1&CBX E:5&QIX C:26&IQX";
|
---|
57 | }
|
---|
58 | }
|
---|
59 |
|
---|
60 | for (unsigned ii = 0; ii < t_frequency::max; ii++) {
|
---|
61 | _obs[ii] = 0;
|
---|
62 | }
|
---|
63 | prepareObs(pppSatObs);
|
---|
64 | }
|
---|
65 |
|
---|
66 | // Destructor
|
---|
67 | ////////////////////////////////////////////////////////////////////////////
|
---|
68 | t_pppSatObs::~t_pppSatObs() {
|
---|
69 | for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
|
---|
70 | delete _obs[iFreq];
|
---|
71 | }
|
---|
72 | }
|
---|
73 |
|
---|
74 | //
|
---|
75 | ////////////////////////////////////////////////////////////////////////////
|
---|
76 | void t_pppSatObs::prepareObs(const t_satObs& pppSatObs) {
|
---|
77 |
|
---|
78 | _model.reset();
|
---|
79 |
|
---|
80 | // Select pseudo-ranges and phase observations
|
---|
81 | // -------------------------------------------
|
---|
82 | QStringList priorList = _signalPriorities.split(" ", QString::SkipEmptyParts);
|
---|
83 | string preferredAttrib;
|
---|
84 | char obsSys = pppSatObs._prn.system(); //cout << "SATELLITE: " << pppSatObs._prn.toString() << endl;
|
---|
85 | for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
|
---|
86 | t_frequency::type frqType = static_cast<t_frequency::type>(iFreq);
|
---|
87 | char frqSys = t_frequency::toString(frqType)[0]; //cout << "frqSys: " << frqSys << endl;
|
---|
88 | char frqNum = t_frequency::toString(frqType)[1]; //cout << "frqNum: " << frqNum << endl;
|
---|
89 | if (obsSys != frqSys) {
|
---|
90 | continue;
|
---|
91 | }
|
---|
92 | QStringList hlp;
|
---|
93 | for (int ii = 0; ii < priorList.size(); ii++) {
|
---|
94 | if (priorList[ii].indexOf(":") != -1) {
|
---|
95 | hlp = priorList[ii].split(":", QString::SkipEmptyParts);
|
---|
96 | if (hlp.size() == 2 && hlp[0].length() == 1 && hlp[0][0] == frqSys) {
|
---|
97 | hlp = hlp[1].split("&", QString::SkipEmptyParts);
|
---|
98 | }
|
---|
99 | if (hlp.size() == 2 && hlp[0].indexOf(frqNum) != -1) {
|
---|
100 | preferredAttrib = hlp[1].toStdString(); //cout << "preferredAttrib: " << preferredAttrib << endl;
|
---|
101 | }
|
---|
102 | }
|
---|
103 | for (unsigned iPref = 0; iPref < preferredAttrib.length(); iPref++) {
|
---|
104 | QString obsType = QString("%1").arg(frqNum) + preferredAttrib[iPref]; //cout << "obstype: " << obsType.toStdString().c_str() << endl;
|
---|
105 | if (_obs[iFreq] == 0) {
|
---|
106 | for (unsigned ii = 0; ii < pppSatObs._obs.size(); ii++) {
|
---|
107 | const t_frqObs* obs = pppSatObs._obs[ii];
|
---|
108 | //cout << "observation2char: " << obs->_rnxType2ch << " vs. " << obsType.toStdString().c_str()<< endl;
|
---|
109 | if (obs->_rnxType2ch == obsType.toStdString() &&
|
---|
110 | obs->_codeValid && obs->_code &&
|
---|
111 | obs->_phaseValid && obs->_phase) {
|
---|
112 | _obs[iFreq] = new t_frqObs(*obs); //cout << "================> newObs: " << obs->_rnxType2ch <<endl;
|
---|
113 | }
|
---|
114 | }
|
---|
115 | }
|
---|
116 | }
|
---|
117 | }
|
---|
118 | }
|
---|
119 |
|
---|
120 | // Used frequency types
|
---|
121 | // --------------------
|
---|
122 | _fType1 = t_lc::toFreq(_prn.system(),t_lc::l1);
|
---|
123 | _fType2 = t_lc::toFreq(_prn.system(),t_lc::l2);
|
---|
124 |
|
---|
125 | // Check whether all required frequencies available
|
---|
126 | // ------------------------------------------------
|
---|
127 | for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) {
|
---|
128 | t_lc::type tLC = OPT->LCs(_prn.system())[ii];
|
---|
129 | if (tLC == t_lc::GIM) {continue;}
|
---|
130 | if (!isValid(tLC)) {
|
---|
131 | _valid = false;
|
---|
132 | return;
|
---|
133 | }
|
---|
134 | }
|
---|
135 |
|
---|
136 | // Find GLONASS Channel Number
|
---|
137 | // ---------------------------
|
---|
138 | if (_prn.system() == 'R') {
|
---|
139 | _channel = PPP_CLIENT->ephPool()->getChannel(_prn);
|
---|
140 | }
|
---|
141 | else {
|
---|
142 | _channel = 0;
|
---|
143 | }
|
---|
144 |
|
---|
145 | // Compute Satellite Coordinates at Time of Transmission
|
---|
146 | // -----------------------------------------------------
|
---|
147 | _xcSat.ReSize(6); _xcSat = 0.0;
|
---|
148 | _vvSat.ReSize(3); _vvSat = 0.0;
|
---|
149 | bool totOK = false;
|
---|
150 | ColumnVector satPosOld(6); satPosOld = 0.0;
|
---|
151 | t_lc::type tLC = t_lc::dummy;
|
---|
152 | if (isValid(t_lc::cIF)) {
|
---|
153 | tLC = t_lc::cIF;
|
---|
154 | }
|
---|
155 | if (tLC == t_lc::dummy && isValid(t_lc::c1)) {
|
---|
156 | tLC = t_lc::c1;
|
---|
157 | }
|
---|
158 | if (tLC == t_lc::dummy && isValid(t_lc::c2)) {
|
---|
159 | tLC = t_lc::c2;
|
---|
160 | }
|
---|
161 | if (tLC == t_lc::dummy) {
|
---|
162 | _valid = false;
|
---|
163 | return;
|
---|
164 | }
|
---|
165 | double prange = obsValue(tLC);
|
---|
166 | for (int ii = 1; ii <= 10; ii++) {
|
---|
167 | bncTime ToT = _time - prange / t_CST::c - _xcSat[3];
|
---|
168 | if (PPP_CLIENT->ephPool()->getCrd(_prn, ToT, _xcSat, _vvSat) != success) {
|
---|
169 | _valid = false;
|
---|
170 | return;
|
---|
171 | }
|
---|
172 | ColumnVector dx = _xcSat - satPosOld;
|
---|
173 | dx[3] *= t_CST::c;
|
---|
174 | if (dx.NormFrobenius() < 1.e-4) {
|
---|
175 | totOK = true;
|
---|
176 | break;
|
---|
177 | }
|
---|
178 | satPosOld = _xcSat;
|
---|
179 | }
|
---|
180 | if (totOK) {
|
---|
181 | _signalPropagationTime = prange / t_CST::c - _xcSat[3];
|
---|
182 | _model._satClkM = _xcSat[3] * t_CST::c;
|
---|
183 | }
|
---|
184 | else {
|
---|
185 | _valid = false;
|
---|
186 | }
|
---|
187 | }
|
---|
188 |
|
---|
189 | //
|
---|
190 | ////////////////////////////////////////////////////////////////////////////
|
---|
191 | void t_pppSatObs::lcCoeff(t_lc::type tLC,
|
---|
192 | map<t_frequency::type, double>& codeCoeff,
|
---|
193 | map<t_frequency::type, double>& phaseCoeff,
|
---|
194 | map<t_frequency::type, double>& ionoCoeff) const {
|
---|
195 |
|
---|
196 | codeCoeff.clear();
|
---|
197 | phaseCoeff.clear();
|
---|
198 | ionoCoeff.clear();
|
---|
199 |
|
---|
200 | double f1 = t_CST::freq(_fType1, _channel);
|
---|
201 | double f2 = t_CST::freq(_fType2, _channel);
|
---|
202 | double f1GPS = t_CST::freq(t_frequency::G1, 0);
|
---|
203 |
|
---|
204 | switch (tLC) {
|
---|
205 | case t_lc::l1:
|
---|
206 | phaseCoeff[_fType1] = 1.0;
|
---|
207 | ionoCoeff [_fType1] = -1.0 * pow(f1GPS, 2) / pow(f1, 2);
|
---|
208 | return;
|
---|
209 | case t_lc::l2:
|
---|
210 | phaseCoeff[_fType2] = 1.0;
|
---|
211 | ionoCoeff [_fType2] = -1.0 * pow(f1GPS, 2) / pow(f2, 2);
|
---|
212 | return;
|
---|
213 | case t_lc::lIF:
|
---|
214 | phaseCoeff[_fType1] = f1 * f1 / (f1 * f1 - f2 * f2);
|
---|
215 | phaseCoeff[_fType2] = -f2 * f2 / (f1 * f1 - f2 * f2);
|
---|
216 | return;
|
---|
217 | case t_lc::MW:
|
---|
218 | phaseCoeff[_fType1] = f1 / (f1 - f2);
|
---|
219 | phaseCoeff[_fType2] = -f2 / (f1 - f2);
|
---|
220 | codeCoeff[_fType1] = -f1 / (f1 + f2);
|
---|
221 | codeCoeff[_fType2] = -f2 / (f1 + f2);
|
---|
222 | return;
|
---|
223 | case t_lc::CL:
|
---|
224 | phaseCoeff[_fType1] = 0.5;
|
---|
225 | codeCoeff [_fType1] = 0.5;
|
---|
226 | return;
|
---|
227 | case t_lc::c1:
|
---|
228 | codeCoeff[_fType1] = 1.0;
|
---|
229 | ionoCoeff[_fType1] = pow(f1GPS, 2) / pow(f1, 2);
|
---|
230 | return;
|
---|
231 | case t_lc::c2:
|
---|
232 | codeCoeff[_fType2] = 1.0;
|
---|
233 | ionoCoeff[_fType2] = pow(f1GPS, 2) / pow(f2, 2);
|
---|
234 | return;
|
---|
235 | case t_lc::cIF:
|
---|
236 | codeCoeff[_fType1] = f1 * f1 / (f1 * f1 - f2 * f2);
|
---|
237 | codeCoeff[_fType2] = -f2 * f2 / (f1 * f1 - f2 * f2);
|
---|
238 | return;
|
---|
239 | case t_lc::GIM:
|
---|
240 | case t_lc::dummy:
|
---|
241 | case t_lc::maxLc:
|
---|
242 | return;
|
---|
243 | }
|
---|
244 | }
|
---|
245 |
|
---|
246 | //
|
---|
247 | ////////////////////////////////////////////////////////////////////////////
|
---|
248 | bool t_pppSatObs::isValid(t_lc::type tLC) const {
|
---|
249 | bool valid = true;
|
---|
250 | obsValue(tLC, &valid);
|
---|
251 |
|
---|
252 | return valid;
|
---|
253 | }
|
---|
254 | //
|
---|
255 | ////////////////////////////////////////////////////////////////////////////
|
---|
256 | double t_pppSatObs::obsValue(t_lc::type tLC, bool* valid) const {
|
---|
257 |
|
---|
258 | double retVal = 0.0;
|
---|
259 | if (valid) *valid = true;
|
---|
260 |
|
---|
261 | // Pseudo observations
|
---|
262 | if (tLC == t_lc::GIM) {
|
---|
263 | if (_stecRefSat == 0.0 || _stecSat == 0.0) {
|
---|
264 | if (valid) *valid = false;
|
---|
265 | return 0.0;
|
---|
266 | }
|
---|
267 | else {
|
---|
268 | return _stecRefSat;
|
---|
269 | }
|
---|
270 | }
|
---|
271 |
|
---|
272 | map<t_frequency::type, double> codeCoeff;
|
---|
273 | map<t_frequency::type, double> phaseCoeff;
|
---|
274 | map<t_frequency::type, double> ionoCoeff;
|
---|
275 | lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
|
---|
276 |
|
---|
277 | map<t_frequency::type, double>::const_iterator it;
|
---|
278 |
|
---|
279 | // Code observations
|
---|
280 | for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
|
---|
281 | t_frequency::type tFreq = it->first;
|
---|
282 | if (_obs[tFreq] == 0) {
|
---|
283 | if (valid) *valid = false;
|
---|
284 | return 0.0;
|
---|
285 | }
|
---|
286 | else {
|
---|
287 | retVal += it->second * _obs[tFreq]->_code;
|
---|
288 | }
|
---|
289 | }
|
---|
290 | // Phase observations
|
---|
291 | for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
|
---|
292 | t_frequency::type tFreq = it->first;
|
---|
293 | if (_obs[tFreq] == 0) {
|
---|
294 | if (valid) *valid = false;
|
---|
295 | return 0.0;
|
---|
296 | }
|
---|
297 | else {
|
---|
298 | retVal += it->second * _obs[tFreq]->_phase * t_CST::lambda(tFreq, _channel);
|
---|
299 | }
|
---|
300 | }
|
---|
301 | return retVal;
|
---|
302 | }
|
---|
303 |
|
---|
304 | //
|
---|
305 | ////////////////////////////////////////////////////////////////////////////
|
---|
306 | double t_pppSatObs::lambda(t_lc::type tLC) const {
|
---|
307 |
|
---|
308 | double f1 = t_CST::freq(_fType1, _channel);
|
---|
309 | double f2 = t_CST::freq(_fType2, _channel);
|
---|
310 |
|
---|
311 | if (tLC == t_lc::l1) {
|
---|
312 | return t_CST::c / f1;
|
---|
313 | }
|
---|
314 | else if (tLC == t_lc::l2) {
|
---|
315 | return t_CST::c / f2;
|
---|
316 | }
|
---|
317 | else if (tLC == t_lc::lIF) {
|
---|
318 | return t_CST::c / (f1 + f2);
|
---|
319 | }
|
---|
320 | else if (tLC == t_lc::MW) {
|
---|
321 | return t_CST::c / (f1 - f2);
|
---|
322 | }
|
---|
323 | else if (tLC == t_lc::CL) {
|
---|
324 | return t_CST::c / f1 / 2.0;
|
---|
325 | }
|
---|
326 |
|
---|
327 | return 0.0;
|
---|
328 | }
|
---|
329 |
|
---|
330 | //
|
---|
331 | ////////////////////////////////////////////////////////////////////////////
|
---|
332 | double t_pppSatObs::sigma(t_lc::type tLC) const {
|
---|
333 |
|
---|
334 | double retVal = 0.0;
|
---|
335 | map<t_frequency::type, double> codeCoeff;
|
---|
336 | map<t_frequency::type, double> phaseCoeff;
|
---|
337 | map<t_frequency::type, double> ionoCoeff;
|
---|
338 | lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
|
---|
339 |
|
---|
340 | if (tLC == t_lc::GIM) {
|
---|
341 | retVal = OPT->_sigmaGIM * OPT->_sigmaGIM + OPT->_sigmaGIM * OPT->_sigmaGIM;
|
---|
342 | }
|
---|
343 |
|
---|
344 | map<t_frequency::type, double>::const_iterator it;
|
---|
345 | for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
|
---|
346 | retVal += it->second * it->second * OPT->_sigmaC1 * OPT->_sigmaC1;
|
---|
347 | }
|
---|
348 |
|
---|
349 | for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
|
---|
350 | retVal += it->second * it->second * OPT->_sigmaL1 * OPT->_sigmaL1;
|
---|
351 | }
|
---|
352 |
|
---|
353 | retVal = sqrt(retVal);
|
---|
354 |
|
---|
355 | // De-Weight GLO+BDS
|
---|
356 | // -----------------
|
---|
357 | if (_prn.system() == 'R') {
|
---|
358 | retVal *= GLO_WEIGHT_FACTOR;
|
---|
359 | }
|
---|
360 | if (_prn.system() == 'C') {
|
---|
361 | retVal *= BDS_WEIGHT_FACTOR;
|
---|
362 | }
|
---|
363 |
|
---|
364 |
|
---|
365 | // Elevation-Dependent Weighting
|
---|
366 | // -----------------------------
|
---|
367 | double cEle = 1.0;
|
---|
368 | if ( (OPT->_eleWgtCode && t_lc::includesCode(tLC)) ||
|
---|
369 | (OPT->_eleWgtPhase && t_lc::includesPhase(tLC)) ) {
|
---|
370 | double eleD = eleSat()*180.0/M_PI;
|
---|
371 | double hlp = fabs(90.0 - eleD);
|
---|
372 | cEle = (1.0 + hlp*hlp*hlp*0.000004);
|
---|
373 | }
|
---|
374 |
|
---|
375 | return cEle * retVal;
|
---|
376 | }
|
---|
377 |
|
---|
378 | //
|
---|
379 | ////////////////////////////////////////////////////////////////////////////
|
---|
380 | double t_pppSatObs::maxRes(t_lc::type tLC) const {
|
---|
381 | double retVal = 0.0;
|
---|
382 |
|
---|
383 | map<t_frequency::type, double> codeCoeff;
|
---|
384 | map<t_frequency::type, double> phaseCoeff;
|
---|
385 | map<t_frequency::type, double> ionoCoeff;
|
---|
386 | lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
|
---|
387 |
|
---|
388 | map<t_frequency::type, double>::const_iterator it;
|
---|
389 | for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
|
---|
390 | retVal += it->second * it->second * OPT->_maxResC1 * OPT->_maxResC1;
|
---|
391 | }
|
---|
392 | for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
|
---|
393 | retVal += it->second * it->second * OPT->_maxResL1 * OPT->_maxResL1;
|
---|
394 | }
|
---|
395 | if (tLC == t_lc::GIM) {
|
---|
396 | retVal = OPT->_maxResGIM * OPT->_maxResGIM + OPT->_maxResGIM * OPT->_maxResGIM;
|
---|
397 | }
|
---|
398 |
|
---|
399 | retVal = sqrt(retVal);
|
---|
400 |
|
---|
401 | return retVal;
|
---|
402 | }
|
---|
403 |
|
---|
404 |
|
---|
405 | //
|
---|
406 | ////////////////////////////////////////////////////////////////////////////
|
---|
407 | t_irc t_pppSatObs::cmpModel(const t_pppStation* station) {
|
---|
408 |
|
---|
409 | // Reset all model values
|
---|
410 | // ----------------------
|
---|
411 | _model.reset();
|
---|
412 |
|
---|
413 | // Topocentric Satellite Position
|
---|
414 | // ------------------------------
|
---|
415 | ColumnVector rSat = _xcSat.Rows(1,3);
|
---|
416 | ColumnVector rRec = station->xyzApr();
|
---|
417 | ColumnVector rhoV = rSat - rRec;
|
---|
418 | _model._rho = rhoV.NormFrobenius();
|
---|
419 |
|
---|
420 | ColumnVector vSat = _vvSat;
|
---|
421 |
|
---|
422 | ColumnVector neu(3);
|
---|
423 | xyz2neu(station->ellApr().data(), rhoV.data(), neu.data());
|
---|
424 |
|
---|
425 | _model._eleSat = acos(sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / _model._rho);
|
---|
426 | if (neu[2] < 0) {
|
---|
427 | _model._eleSat *= -1.0;
|
---|
428 | }
|
---|
429 | _model._azSat = atan2(neu[1], neu[0]);
|
---|
430 |
|
---|
431 | // Sun unit vector
|
---|
432 | ColumnVector xSun = t_astro::Sun(_time.mjddec());
|
---|
433 | xSun /= xSun.norm_Frobenius();
|
---|
434 |
|
---|
435 | // Satellite unit vectors sz, sy, sx
|
---|
436 | ColumnVector sz = -rSat / rSat.norm_Frobenius();
|
---|
437 | ColumnVector sy = crossproduct(sz, xSun);
|
---|
438 | ColumnVector sx = crossproduct(sy, sz);
|
---|
439 |
|
---|
440 | sx /= sx.norm_Frobenius();
|
---|
441 | sy /= sy.norm_Frobenius();
|
---|
442 |
|
---|
443 | // LOS unit vector satellite --> receiver
|
---|
444 | ColumnVector rho = rRec - rSat;
|
---|
445 | rho /= rho.norm_Frobenius();
|
---|
446 |
|
---|
447 | // LOS vector in satellite frame
|
---|
448 | ColumnVector u(3);
|
---|
449 | u(1) = dotproduct(sx, rho);
|
---|
450 | u(2) = dotproduct(sy, rho);
|
---|
451 | u(3) = dotproduct(sz, rho);
|
---|
452 |
|
---|
453 | // Azimuth and elevation in satellite antenna frame
|
---|
454 | _model._elTx = atan2(u(3),sqrt(pow(u(2),2)+pow(u(1),2)));
|
---|
455 | _model._azTx = atan2(u(2),u(1));
|
---|
456 |
|
---|
457 |
|
---|
458 | // Satellite Clocks
|
---|
459 | // ----------------
|
---|
460 | _model._satClkM = _xcSat[3] * t_CST::c;
|
---|
461 |
|
---|
462 | // Receiver Clocks
|
---|
463 | // ---------------
|
---|
464 | _model._recClkM = station->dClk() * t_CST::c;
|
---|
465 |
|
---|
466 | // Sagnac Effect (correction due to Earth rotation)
|
---|
467 | // ------------------------------------------------
|
---|
468 | ColumnVector Omega(3);
|
---|
469 | Omega[0] = 0.0;
|
---|
470 | Omega[1] = 0.0;
|
---|
471 | Omega[2] = t_CST::omega / t_CST::c;
|
---|
472 | _model._sagnac = DotProduct(Omega, crossproduct(rSat, rRec));
|
---|
473 |
|
---|
474 | // Antenna Eccentricity
|
---|
475 | // --------------------
|
---|
476 | _model._antEcc = -DotProduct(station->xyzEcc(), rhoV) / _model._rho;
|
---|
477 |
|
---|
478 | // Antenna Phase Center Offsets and Variations
|
---|
479 | // -------------------------------------------
|
---|
480 | if (PPP_CLIENT->antex()) {
|
---|
481 | for (unsigned ii = 0; ii < t_frequency::max; ii++) {
|
---|
482 | t_frequency::type frqType = static_cast<t_frequency::type>(ii);
|
---|
483 | string frqStr = t_frequency::toString(frqType);
|
---|
484 | if (frqStr[0] != _prn.system()) {continue;}
|
---|
485 | bool found;
|
---|
486 | QString prn(_prn.toString().c_str());
|
---|
487 | _model._antPCO[ii] = PPP_CLIENT->antex()->rcvCorr(station->antName(), frqType, _model._eleSat, _model._azSat, found);
|
---|
488 | _model._antPCO[ii] += PPP_CLIENT->antex()->satCorr(prn, frqType, _model._elTx, _model._azTx, found);
|
---|
489 | if (OPT->_isAPC && found) {
|
---|
490 | // the PCOs as given in the satellite antenna correction for all frequencies
|
---|
491 | // have to be reduced by the PCO of the respective reference frequency
|
---|
492 | if (_prn.system() == 'G') {
|
---|
493 | _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::G1, _model._elTx, _model._azTx, found);
|
---|
494 | }
|
---|
495 | else if (_prn.system() == 'R') {
|
---|
496 | _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::R1, _model._elTx, _model._azTx, found);
|
---|
497 | }
|
---|
498 | else if (_prn.system() == 'E') {
|
---|
499 | _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::E1, _model._elTx, _model._azTx, found);
|
---|
500 | }
|
---|
501 | else if (_prn.system() == 'C') {
|
---|
502 | _model._antPCO[ii] -= PPP_CLIENT->antex()->satCorr(prn, t_frequency::C2, _model._elTx, _model._azTx, found);
|
---|
503 | }
|
---|
504 | }
|
---|
505 | }
|
---|
506 | }
|
---|
507 |
|
---|
508 | // Tropospheric Delay
|
---|
509 | // ------------------
|
---|
510 | _model._tropo = t_tropo::delay_saast(rRec, _model._eleSat);
|
---|
511 |
|
---|
512 | // Code Biases
|
---|
513 | // -----------
|
---|
514 | const t_satCodeBias* satCodeBias = PPP_CLIENT->obsPool()->satCodeBias(_prn);
|
---|
515 | if (satCodeBias) {
|
---|
516 | for (unsigned ii = 0; ii < satCodeBias->_bias.size(); ii++) {
|
---|
517 | const t_frqCodeBias& bias = satCodeBias->_bias[ii];
|
---|
518 | for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
|
---|
519 | string frqStr = t_frequency::toString(t_frequency::type(iFreq));
|
---|
520 | if (frqStr[0] != _prn.system()) {
|
---|
521 | continue;
|
---|
522 | }
|
---|
523 | const t_frqObs* obs = _obs[iFreq];
|
---|
524 | if (obs && obs->_rnxType2ch == bias._rnxType2ch) {
|
---|
525 | _model._codeBias[iFreq] = bias._value;
|
---|
526 | }
|
---|
527 | }
|
---|
528 | }
|
---|
529 | }
|
---|
530 |
|
---|
531 | // Phase Biases
|
---|
532 | // -----------
|
---|
533 | const t_satPhaseBias* satPhaseBias = PPP_CLIENT->obsPool()->satPhaseBias(_prn);
|
---|
534 | double yaw = 0.0;
|
---|
535 | bool ssr = false;
|
---|
536 | if (satPhaseBias) {
|
---|
537 | double dt = station->epochTime() - satPhaseBias->_time;
|
---|
538 | if (satPhaseBias->_updateInt) {
|
---|
539 | dt -= (0.5 * ssrUpdateInt[satPhaseBias->_updateInt]);
|
---|
540 | }
|
---|
541 | yaw = satPhaseBias->_yaw + satPhaseBias->_yawRate * dt;
|
---|
542 | ssr = true;
|
---|
543 | for (unsigned ii = 0; ii < satPhaseBias->_bias.size(); ii++) {
|
---|
544 | const t_frqPhaseBias& bias = satPhaseBias->_bias[ii];
|
---|
545 | for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
|
---|
546 | string frqStr = t_frequency::toString(t_frequency::type(iFreq));
|
---|
547 | if (frqStr[0] != _prn.system()) {
|
---|
548 | continue;
|
---|
549 | }
|
---|
550 | const t_frqObs* obs = _obs[iFreq];
|
---|
551 | if (obs && obs->_rnxType2ch == bias._rnxType2ch) {
|
---|
552 | _model._phaseBias[iFreq] = bias._value;
|
---|
553 | }
|
---|
554 | }
|
---|
555 | }
|
---|
556 | }
|
---|
557 |
|
---|
558 | // Phase Wind-Up
|
---|
559 | // -------------
|
---|
560 | _model._windUp = station->windUp(_time, _prn, rSat, ssr, yaw, vSat) ;
|
---|
561 |
|
---|
562 | // Relativistic effect due to earth gravity
|
---|
563 | // ----------------------------------------
|
---|
564 | double a = rSat.NormFrobenius() + rRec.NormFrobenius();
|
---|
565 | double b = (rSat - rRec).NormFrobenius();
|
---|
566 | double gm = 3.986004418e14; // m3/s2
|
---|
567 | _model._rel = 2 * gm / t_CST::c / t_CST::c * log((a + b) / (a - b));
|
---|
568 |
|
---|
569 | // Tidal Correction
|
---|
570 | // ----------------
|
---|
571 | _model._tideEarth = -DotProduct(station->tideDsplEarth(), rhoV) / _model._rho;
|
---|
572 | _model._tideOcean = -DotProduct(station->tideDsplOcean(), rhoV) / _model._rho;
|
---|
573 |
|
---|
574 | // Ionospheric Delay
|
---|
575 | // -----------------
|
---|
576 | const t_vTec* vTec = PPP_CLIENT->obsPool()->vTec();
|
---|
577 | bool vTecUsage = true;
|
---|
578 | for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) {
|
---|
579 | t_lc::type tLC = OPT->LCs(_prn.system())[ii];
|
---|
580 | if (tLC == t_lc::cIF || tLC == t_lc::lIF) {
|
---|
581 | vTecUsage = false;
|
---|
582 | }
|
---|
583 | }
|
---|
584 |
|
---|
585 | if (vTecUsage && vTec) {
|
---|
586 | double stec = station->stec(vTec, _signalPropagationTime, rSat);
|
---|
587 | double f1GPS = t_CST::freq(t_frequency::G1, 0);
|
---|
588 | for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
|
---|
589 | if (OPT->_pseudoObsIono) { // DCMcodeBias, DCMphaseBias
|
---|
590 | // For scaling the slant ionospheric delays the trick is to be consistent with units!
|
---|
591 | // The conversion of TECU into meters requires the frequency of the signal.
|
---|
592 | // Hence, GPS L1 frequency is used for all systems. The same is true for mu_i in lcCoeff().
|
---|
593 | _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(f1GPS, 2) * stec;
|
---|
594 | }
|
---|
595 | else { // PPP-RTK
|
---|
596 | t_frequency::type frqType = static_cast<t_frequency::type>(iFreq);
|
---|
597 | _model._ionoCodeDelay[iFreq] = 40.3E16 / pow(t_CST::freq(frqType, _channel), 2) * stec;
|
---|
598 | }
|
---|
599 | }
|
---|
600 | }
|
---|
601 |
|
---|
602 | // Set Model Set Flag
|
---|
603 | // ------------------
|
---|
604 | _model._set = true;
|
---|
605 |
|
---|
606 | //printModel();
|
---|
607 |
|
---|
608 | return success;
|
---|
609 | }
|
---|
610 |
|
---|
611 | //
|
---|
612 | ////////////////////////////////////////////////////////////////////////////
|
---|
613 | void t_pppSatObs::printModel() const {
|
---|
614 |
|
---|
615 | LOG.setf(ios::fixed);
|
---|
616 | LOG << "\nMODEL for Satellite " << _prn.toString() << (isReference() ? " (Reference Satellite)" : "")
|
---|
617 |
|
---|
618 | << "\n======================= " << endl
|
---|
619 | << "PPP STRATEGY : " << OPT->_obsmodelTypeStr.at((int)OPT->_obsModelType).toLocal8Bit().constData()
|
---|
620 | << ((OPT->_pseudoObsIono) ? " with pseudo-observations for STEC" : "") << endl
|
---|
621 | << "RHO : " << setw(12) << setprecision(3) << _model._rho << endl
|
---|
622 | << "ELE : " << setw(12) << setprecision(3) << _model._eleSat * RHO_DEG << endl
|
---|
623 | << "AZI : " << setw(12) << setprecision(3) << _model._azSat * RHO_DEG << endl
|
---|
624 | << "SATCLK : " << setw(12) << setprecision(3) << _model._satClkM << endl
|
---|
625 | << "RECCLK : " << setw(12) << setprecision(3) << _model._recClkM << endl
|
---|
626 | << "SAGNAC : " << setw(12) << setprecision(3) << _model._sagnac << endl
|
---|
627 | << "ANTECC : " << setw(12) << setprecision(3) << _model._antEcc << endl
|
---|
628 | << "TROPO : " << setw(12) << setprecision(3) << _model._tropo << endl
|
---|
629 | << "WINDUP : " << setw(12) << setprecision(3) << _model._windUp << endl
|
---|
630 | << "REL : " << setw(12) << setprecision(3) << _model._rel << endl
|
---|
631 | << "EARTH TIDES : " << setw(12) << setprecision(3) << _model._tideEarth << endl
|
---|
632 | << "OCEAN TIDES : " << setw(12) << setprecision(3) << _model._tideOcean << endl
|
---|
633 | << endl
|
---|
634 | << "FREQUENCY DEPENDENT CORRECTIONS:" << endl
|
---|
635 | << "-------------------------------" << endl;
|
---|
636 | for (unsigned iFreq = 1; iFreq < t_frequency::max; iFreq++) {
|
---|
637 | if (_obs[iFreq]) {
|
---|
638 | string frqStr = t_frequency::toString(t_frequency::type(iFreq));
|
---|
639 | if (_prn.system() == frqStr[0]) {
|
---|
640 | LOG << "PCO : " << frqStr << setw(12) << setprecision(3) << _model._antPCO[iFreq] << endl
|
---|
641 | << "BIAS CODE : " << frqStr << setw(12) << setprecision(3) << _model._codeBias[iFreq] << "\t(" << _obs[iFreq]->_rnxType2ch[1] << ") " << endl
|
---|
642 | << "BIAS PHASE : " << frqStr << setw(12) << setprecision(3) << _model._phaseBias[iFreq] << "\t(" << _obs[iFreq]->_rnxType2ch[1] << ") " << endl
|
---|
643 | << "IONO CODEDELAY: " << frqStr << setw(12) << setprecision(3) << _model._ionoCodeDelay[iFreq]<< endl;
|
---|
644 | }
|
---|
645 | }
|
---|
646 | }
|
---|
647 | }
|
---|
648 |
|
---|
649 | //
|
---|
650 | ////////////////////////////////////////////////////////////////////////////
|
---|
651 | void t_pppSatObs::printObsMinusComputed() const {
|
---|
652 | // TODO: cout should be LOG
|
---|
653 | LOG.setf(ios::fixed);
|
---|
654 | LOG << "\nOBS-COMP for Satellite " << _prn.toString() << (isReference() ? " (Reference Satellite)" : "") << endl
|
---|
655 | << "========================== " << endl;
|
---|
656 | for (unsigned ii = 0; ii < OPT->LCs(_prn.system()).size(); ii++) {
|
---|
657 | t_lc::type tLC = OPT->LCs(_prn.system())[ii];
|
---|
658 | LOG << "OBS-CMP " << setw(4) << t_lc::toString(tLC) << ": " << _prn.toString() << " "
|
---|
659 | << setw(12) << setprecision(3) << obsValue(tLC) << " "
|
---|
660 | << setw(12) << setprecision(3) << cmpValue(tLC) << " "
|
---|
661 | << setw(12) << setprecision(3) << obsValue(tLC) - cmpValue(tLC) << endl;
|
---|
662 | }
|
---|
663 | }
|
---|
664 |
|
---|
665 | //
|
---|
666 | ////////////////////////////////////////////////////////////////////////////
|
---|
667 | double t_pppSatObs::cmpValueForBanc(t_lc::type tLC) const {
|
---|
668 | return cmpValue(tLC) - _model._rho - _model._sagnac - _model._recClkM;
|
---|
669 | }
|
---|
670 |
|
---|
671 | //
|
---|
672 | ////////////////////////////////////////////////////////////////////////////
|
---|
673 | double t_pppSatObs::cmpValue(t_lc::type tLC) const {
|
---|
674 | double cmpValue;
|
---|
675 |
|
---|
676 | if (!isValid(tLC)) {
|
---|
677 | cmpValue = 0.0;
|
---|
678 | }
|
---|
679 | else if (tLC == t_lc::GIM) {
|
---|
680 | cmpValue = _stecSat;
|
---|
681 | }
|
---|
682 | else {
|
---|
683 | // Non-Dispersive Part
|
---|
684 | // -------------------
|
---|
685 | double nonDisp = _model._rho
|
---|
686 | + _model._recClkM - _model._satClkM
|
---|
687 | + _model._sagnac + _model._antEcc + _model._tropo
|
---|
688 | + _model._tideEarth + _model._tideOcean + _model._rel;
|
---|
689 |
|
---|
690 | // Add Dispersive Part
|
---|
691 | // -------------------
|
---|
692 | double dispPart = 0.0;
|
---|
693 | map<t_frequency::type, double> codeCoeff;
|
---|
694 | map<t_frequency::type, double> phaseCoeff;
|
---|
695 | map<t_frequency::type, double> ionoCoeff;
|
---|
696 | lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
|
---|
697 | map<t_frequency::type, double>::const_iterator it;
|
---|
698 | for (it = codeCoeff.begin(); it != codeCoeff.end(); it++) {
|
---|
699 | t_frequency::type tFreq = it->first;
|
---|
700 | dispPart += it->second * (_model._antPCO[tFreq] - _model._codeBias[tFreq]);
|
---|
701 | if (OPT->PPPRTK) {
|
---|
702 | dispPart += it->second * (_model._ionoCodeDelay[tFreq]);
|
---|
703 | }
|
---|
704 | }
|
---|
705 | for (it = phaseCoeff.begin(); it != phaseCoeff.end(); it++) {
|
---|
706 | t_frequency::type tFreq = it->first;
|
---|
707 | dispPart += it->second * (_model._antPCO[tFreq] - _model._phaseBias[tFreq] +
|
---|
708 | _model._windUp * t_CST::lambda(tFreq, _channel));
|
---|
709 | if (OPT->PPPRTK) {
|
---|
710 | dispPart += it->second * (- _model._ionoCodeDelay[tFreq]);
|
---|
711 | }
|
---|
712 | }
|
---|
713 | cmpValue = nonDisp + dispPart;
|
---|
714 | }
|
---|
715 |
|
---|
716 | return cmpValue;
|
---|
717 | }
|
---|
718 |
|
---|
719 | //
|
---|
720 | ////////////////////////////////////////////////////////////////////////////
|
---|
721 | void t_pppSatObs::setRes(t_lc::type tLC, double res) {
|
---|
722 | _res[tLC] = res;
|
---|
723 | }
|
---|
724 |
|
---|
725 | //
|
---|
726 | ////////////////////////////////////////////////////////////////////////////
|
---|
727 | double t_pppSatObs::getRes(t_lc::type tLC) const {
|
---|
728 | map<t_lc::type, double>::const_iterator it = _res.find(tLC);
|
---|
729 | if (it != _res.end()) {
|
---|
730 | return it->second;
|
---|
731 | }
|
---|
732 | else {
|
---|
733 | return 0.0;
|
---|
734 | }
|
---|
735 | }
|
---|
736 |
|
---|
737 | //
|
---|
738 | ////////////////////////////////////////////////////////////////////////////
|
---|
739 | void t_pppSatObs::setPseudoObsIono(t_frequency::type freq, double stecRefSat) {
|
---|
740 | _stecSat = _model._ionoCodeDelay[freq];
|
---|
741 | _stecRefSat = stecRefSat;
|
---|
742 | }
|
---|