source: ntrip/trunk/BNC/src/ephemeris.cpp@ 6109

Last change on this file since 6109 was 6109, checked in by mervart, 10 years ago
File size: 29.0 KB
Line 
1#include <math.h>
2#include <sstream>
3#include <iostream>
4#include <iomanip>
5#include <cstring>
6
7#include <newmatio.h>
8
9#include "ephemeris.h"
10#include "bncutils.h"
11#include "bnctime.h"
12#include "bnccore.h"
13#include "bncutils.h"
14#include "pppInclude.h"
15
16using namespace std;
17
18// Constructor
19////////////////////////////////////////////////////////////////////////////
20t_eph::t_eph() {
21 _ok = false;
22 _orbCorr = 0;
23 _clkCorr = 0;
24}
25
26//
27////////////////////////////////////////////////////////////////////////////
28void t_eph::setOrbCorr(const BNC_PPP::t_orbCorr* orbCorr) {
29 delete _orbCorr;
30 _orbCorr = new BNC_PPP::t_orbCorr(*orbCorr);
31}
32
33//
34////////////////////////////////////////////////////////////////////////////
35void t_eph::setClkCorr(const BNC_PPP::t_clkCorr* clkCorr) {
36 delete _clkCorr;
37 _clkCorr = new BNC_PPP::t_clkCorr(*clkCorr);
38}
39
40//
41////////////////////////////////////////////////////////////////////////////
42t_irc t_eph::getCrd(const bncTime& tt, ColumnVector& xc, ColumnVector& vv, bool useCorr) const {
43 xc.ReSize(4);
44 vv.ReSize(3);
45 position(tt.gpsw(), tt.gpssec(), xc.data(), vv.data());
46 if (useCorr) {
47 if (_orbCorr && _clkCorr) {
48
49 double dtO = tt - _orbCorr->_time;
50 ColumnVector dx(3);
51 dx[0] = _orbCorr->_xr[0] + _orbCorr->_dotXr[0] * dtO;
52 dx[1] = _orbCorr->_xr[1] + _orbCorr->_dotXr[1] * dtO;
53 dx[2] = _orbCorr->_xr[2] + _orbCorr->_dotXr[2] * dtO;
54
55 if (_orbCorr->_system == 'R') {
56 RSW_to_XYZ(xc.Rows(1,3), vv.Rows(1,3), dx, dx);
57 }
58
59 xc[0] -= dx[0];
60 xc[1] -= dx[1];
61 xc[2] -= dx[2];
62
63 double dtC = tt - _clkCorr->_time;
64 xc[3] += _clkCorr->_dClk + _clkCorr->_dotDClk * dtC + _clkCorr->_dotDotDClk * dtC * dtC;
65 }
66 else {
67 return failure;
68 }
69 }
70 return success;
71}
72
73// Set GPS Satellite Position
74////////////////////////////////////////////////////////////////////////////
75void t_ephGPS::set(const gpsephemeris* ee) {
76
77 _receptDateTime = currentDateAndTimeGPS();
78
79 _prn.set('G', ee->satellite);
80
81 _TOC.set(ee->GPSweek, ee->TOC);
82 _clock_bias = ee->clock_bias;
83 _clock_drift = ee->clock_drift;
84 _clock_driftrate = ee->clock_driftrate;
85
86 _IODE = ee->IODE;
87 _Crs = ee->Crs;
88 _Delta_n = ee->Delta_n;
89 _M0 = ee->M0;
90
91 _Cuc = ee->Cuc;
92 _e = ee->e;
93 _Cus = ee->Cus;
94 _sqrt_A = ee->sqrt_A;
95
96 _TOEsec = ee->TOE;
97 _Cic = ee->Cic;
98 _OMEGA0 = ee->OMEGA0;
99 _Cis = ee->Cis;
100
101 _i0 = ee->i0;
102 _Crc = ee->Crc;
103 _omega = ee->omega;
104 _OMEGADOT = ee->OMEGADOT;
105
106 _IDOT = ee->IDOT;
107 _L2Codes = 0.0;
108 _TOEweek = ee->GPSweek;
109 _L2PFlag = 0.0;
110
111 if (ee->URAindex <= 6) {
112 _ura = ceil(10.0*pow(2.0, 1.0+((double)ee->URAindex)/2.0))/10.0;
113 }
114 else {
115 _ura = ceil(10.0*pow(2.0, ((double)ee->URAindex)/2.0))/10.0;
116 }
117 _health = ee->SVhealth;
118 _TGD = ee->TGD;
119 _IODC = ee->IODC;
120
121 _TOT = 0.9999e9;
122 _fitInterval = 0.0;
123
124 _ok = true;
125}
126
127// Compute GPS Satellite Position (virtual)
128////////////////////////////////////////////////////////////////////////////
129void t_ephGPS::position(int GPSweek, double GPSweeks,
130 double* xc,
131 double* vv) const {
132
133
134 static const double omegaEarth = 7292115.1467e-11;
135 static const double gmGRS = 398.6005e12;
136
137 memset(xc, 0, 4*sizeof(double));
138 memset(vv, 0, 3*sizeof(double));
139
140 double a0 = _sqrt_A * _sqrt_A;
141 if (a0 == 0) {
142 return;
143 }
144
145 double n0 = sqrt(gmGRS/(a0*a0*a0));
146
147 bncTime tt(GPSweek, GPSweeks);
148 double tk = tt - bncTime(int(_TOEweek), _TOEsec);
149
150 double n = n0 + _Delta_n;
151 double M = _M0 + n*tk;
152 double E = M;
153 double E_last;
154 do {
155 E_last = E;
156 E = M + _e*sin(E);
157 } while ( fabs(E-E_last)*a0 > 0.001 );
158 double v = 2.0*atan( sqrt( (1.0 + _e)/(1.0 - _e) )*tan( E/2 ) );
159 double u0 = v + _omega;
160 double sin2u0 = sin(2*u0);
161 double cos2u0 = cos(2*u0);
162 double r = a0*(1 - _e*cos(E)) + _Crc*cos2u0 + _Crs*sin2u0;
163 double i = _i0 + _IDOT*tk + _Cic*cos2u0 + _Cis*sin2u0;
164 double u = u0 + _Cuc*cos2u0 + _Cus*sin2u0;
165 double xp = r*cos(u);
166 double yp = r*sin(u);
167 double OM = _OMEGA0 + (_OMEGADOT - omegaEarth)*tk -
168 omegaEarth*_TOEsec;
169
170 double sinom = sin(OM);
171 double cosom = cos(OM);
172 double sini = sin(i);
173 double cosi = cos(i);
174 xc[0] = xp*cosom - yp*cosi*sinom;
175 xc[1] = xp*sinom + yp*cosi*cosom;
176 xc[2] = yp*sini;
177
178 double tc = tt - _TOC;
179 xc[3] = _clock_bias + _clock_drift*tc + _clock_driftrate*tc*tc;
180
181 // Velocity
182 // --------
183 double tanv2 = tan(v/2);
184 double dEdM = 1 / (1 - _e*cos(E));
185 double dotv = sqrt((1.0 + _e)/(1.0 - _e)) / cos(E/2)/cos(E/2) / (1 + tanv2*tanv2)
186 * dEdM * n;
187 double dotu = dotv + (-_Cuc*sin2u0 + _Cus*cos2u0)*2*dotv;
188 double dotom = _OMEGADOT - omegaEarth;
189 double doti = _IDOT + (-_Cic*sin2u0 + _Cis*cos2u0)*2*dotv;
190 double dotr = a0 * _e*sin(E) * dEdM * n
191 + (-_Crc*sin2u0 + _Crs*cos2u0)*2*dotv;
192 double dotx = dotr*cos(u) - r*sin(u)*dotu;
193 double doty = dotr*sin(u) + r*cos(u)*dotu;
194
195 vv[0] = cosom *dotx - cosi*sinom *doty // dX / dr
196 - xp*sinom*dotom - yp*cosi*cosom*dotom // dX / dOMEGA
197 + yp*sini*sinom*doti; // dX / di
198
199 vv[1] = sinom *dotx + cosi*cosom *doty
200 + xp*cosom*dotom - yp*cosi*sinom*dotom
201 - yp*sini*cosom*doti;
202
203 vv[2] = sini *doty + yp*cosi *doti;
204
205 // Relativistic Correction
206 // -----------------------
207 xc[3] -= 2.0 * (xc[0]*vv[0] + xc[1]*vv[1] + xc[2]*vv[2]) / t_CST::c / t_CST::c;
208}
209
210// Derivative of the state vector using a simple force model (static)
211////////////////////////////////////////////////////////////////////////////
212ColumnVector t_ephGlo::glo_deriv(double /* tt */, const ColumnVector& xv,
213 double* acc) {
214
215 // State vector components
216 // -----------------------
217 ColumnVector rr = xv.rows(1,3);
218 ColumnVector vv = xv.rows(4,6);
219
220 // Acceleration
221 // ------------
222 static const double gmWGS = 398.60044e12;
223 static const double AE = 6378136.0;
224 static const double OMEGA = 7292115.e-11;
225 static const double C20 = -1082.6257e-6;
226
227 double rho = rr.norm_Frobenius();
228 double t1 = -gmWGS/(rho*rho*rho);
229 double t2 = 3.0/2.0 * C20 * (gmWGS*AE*AE) / (rho*rho*rho*rho*rho);
230 double t3 = OMEGA * OMEGA;
231 double t4 = 2.0 * OMEGA;
232 double z2 = rr(3) * rr(3);
233
234 // Vector of derivatives
235 // ---------------------
236 ColumnVector va(6);
237 va(1) = vv(1);
238 va(2) = vv(2);
239 va(3) = vv(3);
240 va(4) = (t1 + t2*(1.0-5.0*z2/(rho*rho)) + t3) * rr(1) + t4*vv(2) + acc[0];
241 va(5) = (t1 + t2*(1.0-5.0*z2/(rho*rho)) + t3) * rr(2) - t4*vv(1) + acc[1];
242 va(6) = (t1 + t2*(3.0-5.0*z2/(rho*rho)) ) * rr(3) + acc[2];
243
244 return va;
245}
246
247// Compute Glonass Satellite Position (virtual)
248////////////////////////////////////////////////////////////////////////////
249void t_ephGlo::position(int GPSweek, double GPSweeks,
250 double* xc, double* vv) const {
251
252 static const double nominalStep = 10.0;
253
254 memset(xc, 0, 4*sizeof(double));
255 memset(vv, 0, 3*sizeof(double));
256
257 double dtPos = bncTime(GPSweek, GPSweeks) - _tt;
258
259 int nSteps = int(fabs(dtPos) / nominalStep) + 1;
260 double step = dtPos / nSteps;
261
262 double acc[3];
263 acc[0] = _x_acceleration * 1.e3;
264 acc[1] = _y_acceleration * 1.e3;
265 acc[2] = _z_acceleration * 1.e3;
266 for (int ii = 1; ii <= nSteps; ii++) {
267 _xv = rungeKutta4(_tt.gpssec(), _xv, step, acc, glo_deriv);
268 _tt = _tt + step;
269 }
270
271 // Position and Velocity
272 // ---------------------
273 xc[0] = _xv(1);
274 xc[1] = _xv(2);
275 xc[2] = _xv(3);
276
277 vv[0] = _xv(4);
278 vv[1] = _xv(5);
279 vv[2] = _xv(6);
280
281 // Clock Correction
282 // ----------------
283 double dtClk = bncTime(GPSweek, GPSweeks) - _TOC;
284 xc[3] = -_tau + _gamma * dtClk;
285}
286
287// IOD of Glonass Ephemeris (virtual)
288////////////////////////////////////////////////////////////////////////////
289int t_ephGlo::IOD() const {
290 bncTime tMoscow = _TOC - _gps_utc + 3 * 3600.0;
291 return int(tMoscow.daysec() / 900);
292}
293
294// Set Glonass Ephemeris
295////////////////////////////////////////////////////////////////////////////
296void t_ephGlo::set(const glonassephemeris* ee) {
297
298 _receptDateTime = currentDateAndTimeGPS();
299
300 _prn.set('R', ee->almanac_number);
301
302 int ww = ee->GPSWeek;
303 int tow = ee->GPSTOW;
304 updatetime(&ww, &tow, ee->tb*1000, 0); // Moscow -> GPS
305
306 // Check the day once more
307 // -----------------------
308 bool timeChanged = false;
309 {
310 const double secPerDay = 24 * 3600.0;
311 const double secPerWeek = 7 * secPerDay;
312 int ww_old = ww;
313 int tow_old = tow;
314 int currentWeek;
315 double currentSec;
316 currentGPSWeeks(currentWeek, currentSec);
317 bncTime currentTime(currentWeek, currentSec);
318 bncTime hTime(ww, (double) tow);
319
320 if (hTime - currentTime > secPerDay/2.0) {
321 timeChanged = true;
322 tow -= int(secPerDay);
323 if (tow < 0) {
324 tow += int(secPerWeek);
325 ww -= 1;
326 }
327 }
328 else if (hTime - currentTime < -secPerDay/2.0) {
329 timeChanged = true;
330 tow += int(secPerDay);
331 if (tow > secPerWeek) {
332 tow -= int(secPerWeek);
333 ww += 1;
334 }
335 }
336
337 if (false && timeChanged && BNC_CORE->mode() == t_bncCore::batchPostProcessing) {
338 bncTime newHTime(ww, (double) tow);
339 cout << "GLONASS " << ee->almanac_number << " Time Changed at "
340 << currentTime.datestr() << " " << currentTime.timestr()
341 << endl
342 << "old: " << hTime.datestr() << " " << hTime.timestr()
343 << endl
344 << "new: " << newHTime.datestr() << " " << newHTime.timestr()
345 << endl
346 << "eph: " << ee->GPSWeek << " " << ee->GPSTOW << " " << ee->tb
347 << endl
348 << "ww, tow (old): " << ww_old << " " << tow_old
349 << endl
350 << "ww, tow (new): " << ww << " " << tow
351 << endl << endl;
352 }
353 }
354
355 bncTime hlpTime(ww, (double) tow);
356 unsigned year, month, day;
357 hlpTime.civil_date(year, month, day);
358 _gps_utc = gnumleap(year, month, day);
359
360 _TOC.set(ww, tow);
361 _E = ee->E;
362 _tau = ee->tau;
363 _gamma = ee->gamma;
364 _x_pos = ee->x_pos;
365 _x_velocity = ee->x_velocity;
366 _x_acceleration = ee->x_acceleration;
367 _y_pos = ee->y_pos;
368 _y_velocity = ee->y_velocity;
369 _y_acceleration = ee->y_acceleration;
370 _z_pos = ee->z_pos;
371 _z_velocity = ee->z_velocity;
372 _z_acceleration = ee->z_acceleration;
373 _health = 0;
374 _frequency_number = ee->frequency_number;
375 _tki = ee->tk-3*60*60; if (_tki < 0) _tki += 86400;
376
377 // Initialize status vector
378 // ------------------------
379 _tt = _TOC;
380
381 _xv(1) = _x_pos * 1.e3;
382 _xv(2) = _y_pos * 1.e3;
383 _xv(3) = _z_pos * 1.e3;
384 _xv(4) = _x_velocity * 1.e3;
385 _xv(5) = _y_velocity * 1.e3;
386 _xv(6) = _z_velocity * 1.e3;
387
388 _ok = true;
389}
390
391// Set Galileo Satellite Position
392////////////////////////////////////////////////////////////////////////////
393void t_ephGal::set(const galileoephemeris* ee) {
394
395 _receptDateTime = currentDateAndTimeGPS();
396
397 _prn.set('E', ee->satellite);
398
399 _TOC.set(ee->Week, ee->TOC);
400 _clock_bias = ee->clock_bias;
401 _clock_drift = ee->clock_drift;
402 _clock_driftrate = ee->clock_driftrate;
403
404 _IODnav = ee->IODnav;
405 _Crs = ee->Crs;
406 _Delta_n = ee->Delta_n;
407 _M0 = ee->M0;
408
409 _Cuc = ee->Cuc;
410 _e = ee->e;
411 _Cus = ee->Cus;
412 _sqrt_A = ee->sqrt_A;
413
414 _TOEsec = _TOC.gpssec();
415 //// _TOEsec = ee->TOE; //// TODO:
416 _Cic = ee->Cic;
417 _OMEGA0 = ee->OMEGA0;
418 _Cis = ee->Cis;
419
420 _i0 = ee->i0;
421 _Crc = ee->Crc;
422 _omega = ee->omega;
423 _OMEGADOT = ee->OMEGADOT;
424
425 _IDOT = ee->IDOT;
426 _TOEweek = ee->Week;
427
428 _SISA = ee->SISA;
429 _E5aHS = ee->E5aHS;
430 _E5bHS = ee->E5bHS;
431 _BGD_1_5A = ee->BGD_1_5A;
432 _BGD_1_5B = ee->BGD_1_5B;
433
434 _TOT = 0.9999e9;
435
436 _flags = ee->flags;
437
438 _ok = true;
439}
440
441// Compute Galileo Satellite Position (virtual)
442////////////////////////////////////////////////////////////////////////////
443void t_ephGal::position(int GPSweek, double GPSweeks,
444 double* xc,
445 double* vv) const {
446
447 static const double omegaEarth = 7292115.1467e-11;
448 static const double gmWGS = 398.60044e12;
449
450 memset(xc, 0, 4*sizeof(double));
451 memset(vv, 0, 3*sizeof(double));
452
453 double a0 = _sqrt_A * _sqrt_A;
454 if (a0 == 0) {
455 return;
456 }
457
458 double n0 = sqrt(gmWGS/(a0*a0*a0));
459
460 bncTime tt(GPSweek, GPSweeks);
461 double tk = tt - bncTime(_TOC.gpsw(), _TOEsec);
462
463 double n = n0 + _Delta_n;
464 double M = _M0 + n*tk;
465 double E = M;
466 double E_last;
467 do {
468 E_last = E;
469 E = M + _e*sin(E);
470 } while ( fabs(E-E_last)*a0 > 0.001 );
471 double v = 2.0*atan( sqrt( (1.0 + _e)/(1.0 - _e) )*tan( E/2 ) );
472 double u0 = v + _omega;
473 double sin2u0 = sin(2*u0);
474 double cos2u0 = cos(2*u0);
475 double r = a0*(1 - _e*cos(E)) + _Crc*cos2u0 + _Crs*sin2u0;
476 double i = _i0 + _IDOT*tk + _Cic*cos2u0 + _Cis*sin2u0;
477 double u = u0 + _Cuc*cos2u0 + _Cus*sin2u0;
478 double xp = r*cos(u);
479 double yp = r*sin(u);
480 double OM = _OMEGA0 + (_OMEGADOT - omegaEarth)*tk -
481 omegaEarth*_TOEsec;
482
483 double sinom = sin(OM);
484 double cosom = cos(OM);
485 double sini = sin(i);
486 double cosi = cos(i);
487 xc[0] = xp*cosom - yp*cosi*sinom;
488 xc[1] = xp*sinom + yp*cosi*cosom;
489 xc[2] = yp*sini;
490
491 double tc = tt - _TOC;
492 xc[3] = _clock_bias + _clock_drift*tc + _clock_driftrate*tc*tc;
493
494 // Velocity
495 // --------
496 double tanv2 = tan(v/2);
497 double dEdM = 1 / (1 - _e*cos(E));
498 double dotv = sqrt((1.0 + _e)/(1.0 - _e)) / cos(E/2)/cos(E/2) / (1 + tanv2*tanv2)
499 * dEdM * n;
500 double dotu = dotv + (-_Cuc*sin2u0 + _Cus*cos2u0)*2*dotv;
501 double dotom = _OMEGADOT - omegaEarth;
502 double doti = _IDOT + (-_Cic*sin2u0 + _Cis*cos2u0)*2*dotv;
503 double dotr = a0 * _e*sin(E) * dEdM * n
504 + (-_Crc*sin2u0 + _Crs*cos2u0)*2*dotv;
505 double dotx = dotr*cos(u) - r*sin(u)*dotu;
506 double doty = dotr*sin(u) + r*cos(u)*dotu;
507
508 vv[0] = cosom *dotx - cosi*sinom *doty // dX / dr
509 - xp*sinom*dotom - yp*cosi*cosom*dotom // dX / dOMEGA
510 + yp*sini*sinom*doti; // dX / di
511
512 vv[1] = sinom *dotx + cosi*cosom *doty
513 + xp*cosom*dotom - yp*cosi*sinom*dotom
514 - yp*sini*cosom*doti;
515
516 vv[2] = sini *doty + yp*cosi *doti;
517
518 // Relativistic Correction
519 // -----------------------
520 // xc(4) -= 4.442807633e-10 * _e * sqrt(a0) *sin(E);
521 xc[3] -= 2.0 * (xc[0]*vv[0] + xc[1]*vv[1] + xc[2]*vv[2]) / t_CST::c / t_CST::c;
522}
523
524// Constructor
525//////////////////////////////////////////////////////////////////////////////
526t_ephGPS::t_ephGPS(float rnxVersion, const QStringList& lines) {
527
528 const int nLines = 8;
529
530 _ok = false;
531
532 if (lines.size() != nLines) {
533 return;
534 }
535
536 // RINEX Format
537 // ------------
538 int fieldLen = 19;
539
540 int pos[4];
541 pos[0] = (rnxVersion <= 2.12) ? 3 : 4;
542 pos[1] = pos[0] + fieldLen;
543 pos[2] = pos[1] + fieldLen;
544 pos[3] = pos[2] + fieldLen;
545
546 // Read eight lines
547 // ----------------
548 for (int iLine = 0; iLine < nLines; iLine++) {
549 QString line = lines[iLine];
550
551 if ( iLine == 0 ) {
552 QTextStream in(line.left(pos[1]).toAscii());
553
554 int year, month, day, hour, min;
555 double sec;
556
557 QString prnStr;
558 in >> prnStr >> year >> month >> day >> hour >> min >> sec;
559 if (prnStr.at(0) == 'G') {
560 _prn.set('G', prnStr.mid(1).toInt());
561 }
562 else {
563 _prn.set('G', prnStr.toInt());
564 }
565
566 if (year < 80) {
567 year += 2000;
568 }
569 else if (year < 100) {
570 year += 1900;
571 }
572
573 _TOC.set(year, month, day, hour, min, sec);
574
575 if ( readDbl(line, pos[1], fieldLen, _clock_bias ) ||
576 readDbl(line, pos[2], fieldLen, _clock_drift ) ||
577 readDbl(line, pos[3], fieldLen, _clock_driftrate) ) {
578 return;
579 }
580 }
581
582 else if ( iLine == 1 ) {
583 if ( readDbl(line, pos[0], fieldLen, _IODE ) ||
584 readDbl(line, pos[1], fieldLen, _Crs ) ||
585 readDbl(line, pos[2], fieldLen, _Delta_n) ||
586 readDbl(line, pos[3], fieldLen, _M0 ) ) {
587 return;
588 }
589 }
590
591 else if ( iLine == 2 ) {
592 if ( readDbl(line, pos[0], fieldLen, _Cuc ) ||
593 readDbl(line, pos[1], fieldLen, _e ) ||
594 readDbl(line, pos[2], fieldLen, _Cus ) ||
595 readDbl(line, pos[3], fieldLen, _sqrt_A) ) {
596 return;
597 }
598 }
599
600 else if ( iLine == 3 ) {
601 if ( readDbl(line, pos[0], fieldLen, _TOEsec) ||
602 readDbl(line, pos[1], fieldLen, _Cic ) ||
603 readDbl(line, pos[2], fieldLen, _OMEGA0) ||
604 readDbl(line, pos[3], fieldLen, _Cis ) ) {
605 return;
606 }
607 }
608
609 else if ( iLine == 4 ) {
610 if ( readDbl(line, pos[0], fieldLen, _i0 ) ||
611 readDbl(line, pos[1], fieldLen, _Crc ) ||
612 readDbl(line, pos[2], fieldLen, _omega ) ||
613 readDbl(line, pos[3], fieldLen, _OMEGADOT) ) {
614 return;
615 }
616 }
617
618 else if ( iLine == 5 ) {
619 if ( readDbl(line, pos[0], fieldLen, _IDOT ) ||
620 readDbl(line, pos[1], fieldLen, _L2Codes) ||
621 readDbl(line, pos[2], fieldLen, _TOEweek ) ||
622 readDbl(line, pos[3], fieldLen, _L2PFlag) ) {
623 return;
624 }
625 }
626
627 else if ( iLine == 6 ) {
628 if ( readDbl(line, pos[0], fieldLen, _ura ) ||
629 readDbl(line, pos[1], fieldLen, _health) ||
630 readDbl(line, pos[2], fieldLen, _TGD ) ||
631 readDbl(line, pos[3], fieldLen, _IODC ) ) {
632 return;
633 }
634 }
635
636 else if ( iLine == 7 ) {
637 if ( readDbl(line, pos[0], fieldLen, _TOT) ) {
638 return;
639 }
640 readDbl(line, pos[1], fieldLen, _fitInterval); // _fitInterval optional
641 }
642 }
643
644 _ok = true;
645}
646
647// Constructor
648//////////////////////////////////////////////////////////////////////////////
649t_ephGlo::t_ephGlo(float rnxVersion, const QStringList& lines) {
650
651 const int nLines = 4;
652
653 _ok = false;
654
655 if (lines.size() != nLines) {
656 return;
657 }
658
659 // RINEX Format
660 // ------------
661 int fieldLen = 19;
662
663 int pos[4];
664 pos[0] = (rnxVersion <= 2.12) ? 3 : 4;
665 pos[1] = pos[0] + fieldLen;
666 pos[2] = pos[1] + fieldLen;
667 pos[3] = pos[2] + fieldLen;
668
669 // Read four lines
670 // ---------------
671 for (int iLine = 0; iLine < nLines; iLine++) {
672 QString line = lines[iLine];
673
674 if ( iLine == 0 ) {
675 QTextStream in(line.left(pos[1]).toAscii());
676
677 int year, month, day, hour, min;
678 double sec;
679
680 QString prnStr;
681 in >> prnStr >> year >> month >> day >> hour >> min >> sec;
682 if (prnStr.at(0) == 'R') {
683 _prn.set('R', prnStr.mid(1).toInt());
684 }
685 else {
686 _prn.set('R', prnStr.toInt());
687 }
688
689 if (year < 80) {
690 year += 2000;
691 }
692 else if (year < 100) {
693 year += 1900;
694 }
695
696 _gps_utc = gnumleap(year, month, day);
697
698 _TOC.set(year, month, day, hour, min, sec);
699 _TOC = _TOC + _gps_utc;
700
701 if ( readDbl(line, pos[1], fieldLen, _tau ) ||
702 readDbl(line, pos[2], fieldLen, _gamma) ||
703 readDbl(line, pos[3], fieldLen, _tki ) ) {
704 return;
705 }
706
707 _tau = -_tau;
708 }
709
710 else if ( iLine == 1 ) {
711 if ( readDbl(line, pos[0], fieldLen, _x_pos ) ||
712 readDbl(line, pos[1], fieldLen, _x_velocity ) ||
713 readDbl(line, pos[2], fieldLen, _x_acceleration) ||
714 readDbl(line, pos[3], fieldLen, _health ) ) {
715 return;
716 }
717 }
718
719 else if ( iLine == 2 ) {
720 if ( readDbl(line, pos[0], fieldLen, _y_pos ) ||
721 readDbl(line, pos[1], fieldLen, _y_velocity ) ||
722 readDbl(line, pos[2], fieldLen, _y_acceleration ) ||
723 readDbl(line, pos[3], fieldLen, _frequency_number) ) {
724 return;
725 }
726 }
727
728 else if ( iLine == 3 ) {
729 if ( readDbl(line, pos[0], fieldLen, _z_pos ) ||
730 readDbl(line, pos[1], fieldLen, _z_velocity ) ||
731 readDbl(line, pos[2], fieldLen, _z_acceleration) ||
732 readDbl(line, pos[3], fieldLen, _E ) ) {
733 return;
734 }
735 }
736 }
737
738 // Initialize status vector
739 // ------------------------
740 _tt = _TOC;
741 _xv.ReSize(6);
742 _xv(1) = _x_pos * 1.e3;
743 _xv(2) = _y_pos * 1.e3;
744 _xv(3) = _z_pos * 1.e3;
745 _xv(4) = _x_velocity * 1.e3;
746 _xv(5) = _y_velocity * 1.e3;
747 _xv(6) = _z_velocity * 1.e3;
748
749 _ok = true;
750}
751
752// Constructor
753//////////////////////////////////////////////////////////////////////////////
754t_ephGal::t_ephGal(float rnxVersion, const QStringList& lines) {
755
756 const int nLines = 8;
757
758 _ok = false;
759
760 if (lines.size() != nLines) {
761 return;
762 }
763
764 // RINEX Format
765 // ------------
766 int fieldLen = 19;
767
768 int pos[4];
769 pos[0] = (rnxVersion <= 2.12) ? 3 : 4;
770 pos[1] = pos[0] + fieldLen;
771 pos[2] = pos[1] + fieldLen;
772 pos[3] = pos[2] + fieldLen;
773
774 // Read eight lines
775 // ----------------
776 for (int iLine = 0; iLine < nLines; iLine++) {
777 QString line = lines[iLine];
778
779 if ( iLine == 0 ) {
780 QTextStream in(line.left(pos[1]).toAscii());
781
782 int year, month, day, hour, min;
783 double sec;
784
785 QString prnStr;
786 in >> prnStr >> year >> month >> day >> hour >> min >> sec;
787 if (prnStr.at(0) == 'E') {
788 _prn.set('E', prnStr.mid(1).toInt());
789 }
790 else {
791 _prn.set('E', prnStr.toInt());
792 }
793
794 if (year < 80) {
795 year += 2000;
796 }
797 else if (year < 100) {
798 year += 1900;
799 }
800
801 _TOC.set(year, month, day, hour, min, sec);
802
803 if ( readDbl(line, pos[1], fieldLen, _clock_bias ) ||
804 readDbl(line, pos[2], fieldLen, _clock_drift ) ||
805 readDbl(line, pos[3], fieldLen, _clock_driftrate) ) {
806 return;
807 }
808 }
809
810 else if ( iLine == 1 ) {
811 if ( readDbl(line, pos[0], fieldLen, _IODnav ) ||
812 readDbl(line, pos[1], fieldLen, _Crs ) ||
813 readDbl(line, pos[2], fieldLen, _Delta_n) ||
814 readDbl(line, pos[3], fieldLen, _M0 ) ) {
815 return;
816 }
817 }
818
819 else if ( iLine == 2 ) {
820 if ( readDbl(line, pos[0], fieldLen, _Cuc ) ||
821 readDbl(line, pos[1], fieldLen, _e ) ||
822 readDbl(line, pos[2], fieldLen, _Cus ) ||
823 readDbl(line, pos[3], fieldLen, _sqrt_A) ) {
824 return;
825 }
826 }
827
828 else if ( iLine == 3 ) {
829 if ( readDbl(line, pos[0], fieldLen, _TOEsec) ||
830 readDbl(line, pos[1], fieldLen, _Cic ) ||
831 readDbl(line, pos[2], fieldLen, _OMEGA0) ||
832 readDbl(line, pos[3], fieldLen, _Cis ) ) {
833 return;
834 }
835 }
836
837 else if ( iLine == 4 ) {
838 if ( readDbl(line, pos[0], fieldLen, _i0 ) ||
839 readDbl(line, pos[1], fieldLen, _Crc ) ||
840 readDbl(line, pos[2], fieldLen, _omega ) ||
841 readDbl(line, pos[3], fieldLen, _OMEGADOT) ) {
842 return;
843 }
844 }
845
846 else if ( iLine == 5 ) {
847 if ( readDbl(line, pos[0], fieldLen, _IDOT ) ||
848 readDbl(line, pos[2], fieldLen, _TOEweek) ) {
849 return;
850 }
851 }
852
853 else if ( iLine == 6 ) {
854 if ( readDbl(line, pos[0], fieldLen, _SISA ) ||
855 readDbl(line, pos[1], fieldLen, _E5aHS ) ||
856 readDbl(line, pos[2], fieldLen, _BGD_1_5A) ||
857 readDbl(line, pos[3], fieldLen, _BGD_1_5B) ) {
858 return;
859 }
860 }
861
862 else if ( iLine == 7 ) {
863 if ( readDbl(line, pos[0], fieldLen, _TOT) ) {
864 return;
865 }
866 }
867 }
868
869 _ok = true;
870}
871
872//
873//////////////////////////////////////////////////////////////////////////////
874QString t_eph::rinexDateStr(const bncTime& tt, const t_prn& prn, double version) {
875 QString prnStr(prn.toString().c_str());
876 return rinexDateStr(tt, prnStr, version);
877}
878
879//
880//////////////////////////////////////////////////////////////////////////////
881QString t_eph::rinexDateStr(const bncTime& tt, const QString& prnStr, double version) {
882
883 QString datStr;
884
885 unsigned year, month, day, hour, min;
886 double sec;
887 tt.civil_date(year, month, day);
888 tt.civil_time(hour, min, sec);
889
890 QTextStream out(&datStr);
891
892 if (version < 3.0) {
893 QString prnHlp = prnStr.mid(1,2); if (prnHlp[0] == '0') prnHlp[0] = ' ';
894 out << prnHlp << QString(" %1 %2 %3 %4 %5%6")
895 .arg(year % 100, 2, 10, QChar('0'))
896 .arg(month, 2)
897 .arg(day, 2)
898 .arg(hour, 2)
899 .arg(min, 2)
900 .arg(sec, 5, 'f',1);
901 }
902 else {
903 out << prnStr << QString(" %1 %2 %3 %4 %5 %6")
904 .arg(year, 4)
905 .arg(month, 2, 10, QChar('0'))
906 .arg(day, 2, 10, QChar('0'))
907 .arg(hour, 2, 10, QChar('0'))
908 .arg(min, 2, 10, QChar('0'))
909 .arg(int(sec), 2, 10, QChar('0'));
910 }
911
912 return datStr;
913}
914
915// RINEX Format String
916//////////////////////////////////////////////////////////////////////////////
917QString t_ephGPS::toString(double version) const {
918
919 QString rnxStr = rinexDateStr(_TOC, _prn, version);
920
921 QTextStream out(&rnxStr);
922
923 out << QString("%1%2%3\n")
924 .arg(_clock_bias, 19, 'e', 12)
925 .arg(_clock_drift, 19, 'e', 12)
926 .arg(_clock_driftrate, 19, 'e', 12);
927
928 QString fmt = version < 3.0 ? " %1%2%3%4\n" : " %1%2%3%4\n";
929
930 out << QString(fmt)
931 .arg(_IODE, 19, 'e', 12)
932 .arg(_Crs, 19, 'e', 12)
933 .arg(_Delta_n, 19, 'e', 12)
934 .arg(_M0, 19, 'e', 12);
935
936 out << QString(fmt)
937 .arg(_Cuc, 19, 'e', 12)
938 .arg(_e, 19, 'e', 12)
939 .arg(_Cus, 19, 'e', 12)
940 .arg(_sqrt_A, 19, 'e', 12);
941
942 out << QString(fmt)
943 .arg(_TOEsec, 19, 'e', 12)
944 .arg(_Cic, 19, 'e', 12)
945 .arg(_OMEGA0, 19, 'e', 12)
946 .arg(_Cis, 19, 'e', 12);
947
948 out << QString(fmt)
949 .arg(_i0, 19, 'e', 12)
950 .arg(_Crc, 19, 'e', 12)
951 .arg(_omega, 19, 'e', 12)
952 .arg(_OMEGADOT, 19, 'e', 12);
953
954 out << QString(fmt)
955 .arg(_IDOT, 19, 'e', 12)
956 .arg(_L2Codes, 19, 'e', 12)
957 .arg(_TOEweek, 19, 'e', 12)
958 .arg(_L2PFlag, 19, 'e', 12);
959
960 out << QString(fmt)
961 .arg(_ura, 19, 'e', 12)
962 .arg(_health, 19, 'e', 12)
963 .arg(_TGD, 19, 'e', 12)
964 .arg(_IODC, 19, 'e', 12);
965
966 out << QString(fmt)
967 .arg(_TOT, 19, 'e', 12)
968 .arg(_fitInterval, 19, 'e', 12)
969 .arg("", 19, QChar(' '))
970 .arg("", 19, QChar(' '));
971
972 return rnxStr;
973}
974
975// RINEX Format String
976//////////////////////////////////////////////////////////////////////////////
977QString t_ephGlo::toString(double version) const {
978
979 QString rnxStr = rinexDateStr(_TOC-_gps_utc, _prn, version);
980
981 QTextStream out(&rnxStr);
982
983 out << QString("%1%2%3\n")
984 .arg(-_tau, 19, 'e', 12)
985 .arg(_gamma, 19, 'e', 12)
986 .arg(_tki, 19, 'e', 12);
987
988 QString fmt = version < 3.0 ? " %1%2%3%4\n" : " %1%2%3%4\n";
989
990 out << QString(fmt)
991 .arg(_x_pos, 19, 'e', 12)
992 .arg(_x_velocity, 19, 'e', 12)
993 .arg(_x_acceleration, 19, 'e', 12)
994 .arg(_health, 19, 'e', 12);
995
996 out << QString(fmt)
997 .arg(_y_pos, 19, 'e', 12)
998 .arg(_y_velocity, 19, 'e', 12)
999 .arg(_y_acceleration, 19, 'e', 12)
1000 .arg(_frequency_number, 19, 'e', 12);
1001
1002 out << QString(fmt)
1003 .arg(_z_pos, 19, 'e', 12)
1004 .arg(_z_velocity, 19, 'e', 12)
1005 .arg(_z_acceleration, 19, 'e', 12)
1006 .arg(_E, 19, 'e', 12);
1007
1008 return rnxStr;
1009}
1010
1011// RINEX Format String
1012//////////////////////////////////////////////////////////////////////////////
1013QString t_ephGal::toString(double version) const {
1014
1015 QString rnxStr = rinexDateStr(_TOC, _prn, version);
1016
1017 QTextStream out(&rnxStr);
1018
1019 out << QString("%1%2%3\n")
1020 .arg(_clock_bias, 19, 'e', 12)
1021 .arg(_clock_drift, 19, 'e', 12)
1022 .arg(_clock_driftrate, 19, 'e', 12);
1023
1024 QString fmt = version < 3.0 ? " %1%2%3%4\n" : " %1%2%3%4\n";
1025
1026 out << QString(fmt)
1027 .arg(_IODnav, 19, 'e', 12)
1028 .arg(_Crs, 19, 'e', 12)
1029 .arg(_Delta_n, 19, 'e', 12)
1030 .arg(_M0, 19, 'e', 12);
1031
1032 out << QString(fmt)
1033 .arg(_Cuc, 19, 'e', 12)
1034 .arg(_e, 19, 'e', 12)
1035 .arg(_Cus, 19, 'e', 12)
1036 .arg(_sqrt_A, 19, 'e', 12);
1037
1038 out << QString(fmt)
1039 .arg(_TOEsec, 19, 'e', 12)
1040 .arg(_Cic, 19, 'e', 12)
1041 .arg(_OMEGA0, 19, 'e', 12)
1042 .arg(_Cis, 19, 'e', 12);
1043
1044 out << QString(fmt)
1045 .arg(_i0, 19, 'e', 12)
1046 .arg(_Crc, 19, 'e', 12)
1047 .arg(_omega, 19, 'e', 12)
1048 .arg(_OMEGADOT, 19, 'e', 12);
1049
1050 int dataSource = 0;
1051 double HS = 0.0;
1052 if ( (_flags & GALEPHF_INAV) == GALEPHF_INAV ) {
1053 dataSource |= (1<<0);
1054 dataSource |= (1<<9);
1055 HS = _E5bHS;
1056 }
1057 else if ( (_flags & GALEPHF_FNAV) == GALEPHF_FNAV ) {
1058 dataSource |= (1<<1);
1059 dataSource |= (1<<8);
1060 HS = _E5aHS;
1061 }
1062 out << QString(fmt)
1063 .arg(_IDOT, 19, 'e', 12)
1064 .arg(double(dataSource), 19, 'e', 12)
1065 .arg(_TOEweek, 19, 'e', 12)
1066 .arg(0.0, 19, 'e', 12);
1067
1068 out << QString(fmt)
1069 .arg(_SISA, 19, 'e', 12)
1070 .arg(HS, 19, 'e', 12)
1071 .arg(_BGD_1_5A, 19, 'e', 12)
1072 .arg(_BGD_1_5B, 19, 'e', 12);
1073
1074 out << QString(fmt)
1075 .arg(_TOT, 19, 'e', 12)
1076 .arg("", 19, QChar(' '))
1077 .arg("", 19, QChar(' '))
1078 .arg("", 19, QChar(' '));
1079
1080 return rnxStr;
1081}
1082
Note: See TracBrowser for help on using the repository browser.