source: ntrip/branches/BNC_2.12/src/bncutils.cpp@ 8855

Last change on this file since 8855 was 8855, checked in by stuerze, 5 years ago

changes regarding fit interval for GPS and QZSS completed

File size: 30.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: bncutils
30 *
31 * Purpose: Auxiliary Functions
32 *
33 * Author: L. Mervart
34 *
35 * Created: 30-Aug-2006
36 *
37 * Changes:
38 *
39 * -----------------------------------------------------------------------*/
40
41#include <iostream>
42#include <ctime>
43#include <math.h>
44
45#include <QRegExp>
46#include <QStringList>
47#include <QDateTime>
48
49#include <newmatap.h>
50
51#include "bncutils.h"
52#include "bnccore.h"
53
54using namespace std;
55
56struct leapseconds { /* specify the day of leap second */
57 int day; /* this is the day, where 23:59:59 exists 2 times */
58 int month; /* not the next day! */
59 int year;
60 int taicount;
61};
62static const int months[13] = {0,31,28,31,30,31,30,31,31,30,31,30,31};
63static const struct leapseconds leap[] = {
64/*{31, 12, 1971, 10},*/
65/*{30, 06, 1972, 11},*/
66/*{31, 12, 1972, 12},*/
67/*{31, 12, 1973, 13},*/
68/*{31, 12, 1974, 14},*/
69/*{31, 12, 1975, 15},*/
70/*{31, 12, 1976, 16},*/
71/*{31, 12, 1977, 17},*/
72/*{31, 12, 1978, 18},*/
73/*{31, 12, 1979, 19},*/
74{30, 06, 1981,20},
75{30, 06, 1982,21},
76{30, 06, 1983,22},
77{30, 06, 1985,23},
78{31, 12, 1987,24},
79{31, 12, 1989,25},
80{31, 12, 1990,26},
81{30, 06, 1992,27},
82{30, 06, 1993,28},
83{30, 06, 1994,29},
84{31, 12, 1995,30},
85{30, 06, 1997,31},
86{31, 12, 1998,32},
87{31, 12, 2005,33},
88{31, 12, 2008,34},
89{30, 06, 2012,35},
90{30, 06, 2015,36},
91{01, 01, 2017,37},
92{0,0,0,0} /* end marker */
93};
94
95#define GPSLEAPSTART 19 /* 19 leap seconds existed at 6.1.1980 */
96
97static int longyear(int year, int month)
98{
99 if(!(year % 4) && (!(year % 400) || (year % 100)))
100 {
101 if(!month || month == 2)
102 return 1;
103 }
104 return 0;
105}
106
107int gnumleap(int year, int month, int day)
108{
109 int ls = 0;
110 const struct leapseconds *l;
111
112 for(l = leap; l->taicount && year >= l->year; ++l)
113 {
114 if(year > l->year || month > l->month || (month == l->month && day > l->day))
115 ls = l->taicount - GPSLEAPSTART;
116 }
117 return ls;
118}
119
120/* Convert Moscow time into UTC (fixnumleap == 1) or GPS (fixnumleap == 0) */
121void updatetime(int *week, int *secOfWeek, int mSecOfWeek, bool fixnumleap)
122{
123 int y,m,d,k,l, nul;
124 unsigned int j = *week*(7*24*60*60) + *secOfWeek + 5*24*60*60+3*60*60;
125 int glo_daynumber = 0, glo_timeofday;
126 for(y = 1980; j >= (unsigned int)(k = (l = (365+longyear(y,0)))*24*60*60)
127 + gnumleap(y+1,1,1); ++y)
128 {
129 j -= k; glo_daynumber += l;
130 }
131 for(m = 1; j >= (unsigned int)(k = (l = months[m]+longyear(y, m))*24*60*60)
132 + gnumleap(y, m+1, 1); ++m)
133 {
134 j -= k; glo_daynumber += l;
135 }
136 for(d = 1; j >= 24UL*60UL*60UL + gnumleap(y, m, d+1); ++d)
137 j -= 24*60*60;
138 glo_daynumber -= 16*365+4-d;
139 nul = gnumleap(y, m, d);
140 glo_timeofday = j-nul;
141
142 // original version
143 // if(mSecOfWeek < 5*60*1000 && glo_timeofday > 23*60*60)
144 // *secOfWeek += 24*60*60;
145 // else if(glo_timeofday < 5*60 && mSecOfWeek > 23*60*60*1000)
146 // *secOfWeek -= 24*60*60;
147
148 // new version
149 if(mSecOfWeek < 4*60*60*1000 && glo_timeofday > 20*60*60)
150 *secOfWeek += 24*60*60;
151 else if(glo_timeofday < 4*60*60 && mSecOfWeek > 20*60*60*1000)
152 *secOfWeek -= 24*60*60;
153
154 *secOfWeek += mSecOfWeek/1000-glo_timeofday;
155 if(fixnumleap)
156 *secOfWeek -= nul;
157 if(*secOfWeek < 0) {*secOfWeek += 24*60*60*7; --*week; }
158 if(*secOfWeek >= 24*60*60*7) {*secOfWeek -= 24*60*60*7; ++*week; }
159}
160
161//
162////////////////////////////////////////////////////////////////////////////
163void expandEnvVar(QString& str) {
164
165 QRegExp rx("(\\$\\{.+\\})");
166
167 if (rx.indexIn(str) != -1) {
168 QStringListIterator it(rx.capturedTexts());
169 if (it.hasNext()) {
170 QString rxStr = it.next();
171 QString envVar = rxStr.mid(2,rxStr.length()-3);
172 str.replace(rxStr, qgetenv(envVar.toAscii()));
173 }
174 }
175}
176
177// Strip White Space
178////////////////////////////////////////////////////////////////////////////
179void stripWhiteSpace(string& str) {
180 if (!str.empty()) {
181 string::size_type beg = str.find_first_not_of(" \t\f\n\r\v");
182 string::size_type end = str.find_last_not_of(" \t\f\n\r\v");
183 if (beg > str.max_size())
184 str.erase();
185 else
186 str = str.substr(beg, end-beg+1);
187 }
188}
189
190//
191////////////////////////////////////////////////////////////////////////////
192QDateTime dateAndTimeFromGPSweek(int GPSWeek, double GPSWeeks) {
193
194 static const QDate zeroEpoch(1980, 1, 6);
195
196 QDate date(zeroEpoch);
197 QTime time(0,0,0,0);
198
199 int weekDays = int(GPSWeeks) / 86400;
200 date = date.addDays( GPSWeek * 7 + weekDays );
201 time = time.addMSecs( int( (GPSWeeks - 86400 * weekDays) * 1e3 ) );
202
203 return QDateTime(date,time);
204}
205
206//
207////////////////////////////////////////////////////////////////////////////
208void currentGPSWeeks(int& week, double& sec) {
209
210 QDateTime currDateTimeGPS;
211
212 if ( BNC_CORE->dateAndTimeGPSSet() ) {
213 currDateTimeGPS = BNC_CORE->dateAndTimeGPS();
214 }
215 else {
216 currDateTimeGPS = QDateTime::currentDateTime().toUTC();
217 QDate hlp = currDateTimeGPS.date();
218 currDateTimeGPS = currDateTimeGPS.addSecs(gnumleap(hlp.year(),
219 hlp.month(), hlp.day()));
220 }
221
222 QDate currDateGPS = currDateTimeGPS.date();
223 QTime currTimeGPS = currDateTimeGPS.time();
224
225 week = int( (double(currDateGPS.toJulianDay()) - 2444244.5) / 7 );
226
227 sec = (currDateGPS.dayOfWeek() % 7) * 24.0 * 3600.0 +
228 currTimeGPS.hour() * 3600.0 +
229 currTimeGPS.minute() * 60.0 +
230 currTimeGPS.second() +
231 currTimeGPS.msec() / 1000.0;
232}
233
234//
235////////////////////////////////////////////////////////////////////////////
236QDateTime currentDateAndTimeGPS() {
237 if ( BNC_CORE->dateAndTimeGPSSet() ) {
238 return BNC_CORE->dateAndTimeGPS();
239 }
240 else {
241 int GPSWeek;
242 double GPSWeeks;
243 currentGPSWeeks(GPSWeek, GPSWeeks);
244 return dateAndTimeFromGPSweek(GPSWeek, GPSWeeks);
245 }
246}
247
248//
249////////////////////////////////////////////////////////////////////////////
250bool checkForWrongObsEpoch(bncTime obsEpoch) {
251 const double maxDt = 600.0;
252 bncTime obsTime = obsEpoch;
253 int week;
254 double sec;
255 currentGPSWeeks(week, sec);
256 bncTime currTime(week, sec);
257 if (((currTime - obsTime) < 0.0) ||
258 (fabs(currTime - obsTime) > maxDt)) {
259 return true;
260 }
261 return false;
262}
263//
264////////////////////////////////////////////////////////////////////////////
265QByteArray ggaString(const QByteArray& latitude,
266 const QByteArray& longitude,
267 const QByteArray& height,
268 const QString& ggaType) {
269
270 double lat = strtod(latitude,NULL);
271 double lon = strtod(longitude,NULL);
272 double hei = strtod(height,NULL);
273 QString sentences = "GPGGA,";
274 if (ggaType.contains("GNGGA")) {
275 sentences = "GNGGA,";
276 }
277
278 const char* flagN="N";
279 const char* flagE="E";
280 if (lon >180.) {lon=(lon-360.)*(-1.); flagE="W";}
281 if ((lon < 0.) && (lon >= -180.)) {lon=lon*(-1.); flagE="W";}
282 if (lon < -180.) {lon=(lon+360.); flagE="E";}
283 if (lat < 0.) {lat=lat*(-1.); flagN="S";}
284 QTime ttime(QDateTime::currentDateTime().toUTC().time());
285 int lat_deg = (int)lat;
286 double lat_min=(lat-lat_deg)*60.;
287 int lon_deg = (int)lon;
288 double lon_min=(lon-lon_deg)*60.;
289 int hh = 0 , mm = 0;
290 double ss = 0.0;
291 hh=ttime.hour();
292 mm=ttime.minute();
293 ss=(double)ttime.second()+0.001*ttime.msec();
294 QString gga;
295 gga += sentences;
296 gga += QString("%1%2%3,").arg((int)hh, 2, 10, QLatin1Char('0')).arg((int)mm, 2, 10, QLatin1Char('0')).arg((int)ss, 2, 10, QLatin1Char('0'));
297 gga += QString("%1%2,").arg((int)lat_deg,2, 10, QLatin1Char('0')).arg(lat_min, 7, 'f', 4, QLatin1Char('0'));
298 gga += flagN;
299 gga += QString(",%1%2,").arg((int)lon_deg,3, 10, QLatin1Char('0')).arg(lon_min, 7, 'f', 4, QLatin1Char('0'));
300 gga += flagE + QString(",1,05,1.00");
301 gga += QString(",%1,").arg(hei, 2, 'f', 1);
302 gga += QString("M,10.000,M,,");
303
304 unsigned char XOR = 0;
305 for (int ii = 0; ii < gga.length(); ii++) {
306 XOR ^= (unsigned char) gga[ii].toAscii();
307 }
308 gga = "$" + gga + QString("*%1").arg(XOR, 2, 16, QLatin1Char('0')) + "\n";
309
310 return gga.toAscii();
311}
312
313//
314////////////////////////////////////////////////////////////////////////////
315void RSW_to_XYZ(const ColumnVector& rr, const ColumnVector& vv,
316 const ColumnVector& rsw, ColumnVector& xyz) {
317
318 ColumnVector along = vv / vv.norm_Frobenius();
319 ColumnVector cross = crossproduct(rr, vv); cross /= cross.norm_Frobenius();
320 ColumnVector radial = crossproduct(along, cross);
321
322 Matrix RR(3,3);
323 RR.Column(1) = radial;
324 RR.Column(2) = along;
325 RR.Column(3) = cross;
326
327 xyz = RR * rsw;
328}
329
330// Transformation xyz --> radial, along track, out-of-plane
331////////////////////////////////////////////////////////////////////////////
332void XYZ_to_RSW(const ColumnVector& rr, const ColumnVector& vv,
333 const ColumnVector& xyz, ColumnVector& rsw) {
334
335 ColumnVector along = vv / vv.norm_Frobenius();
336 ColumnVector cross = crossproduct(rr, vv); cross /= cross.norm_Frobenius();
337 ColumnVector radial = crossproduct(along, cross);
338
339 rsw.ReSize(3);
340 rsw(1) = DotProduct(xyz, radial);
341 rsw(2) = DotProduct(xyz, along);
342 rsw(3) = DotProduct(xyz, cross);
343}
344
345// Rectangular Coordinates -> Ellipsoidal Coordinates
346////////////////////////////////////////////////////////////////////////////
347t_irc xyz2ell(const double* XYZ, double* Ell) {
348
349 const double bell = t_CST::aell*(1.0-1.0/t_CST::fInv) ;
350 const double e2 = (t_CST::aell*t_CST::aell-bell*bell)/(t_CST::aell*t_CST::aell) ;
351 const double e2c = (t_CST::aell*t_CST::aell-bell*bell)/(bell*bell) ;
352
353 double nn, ss, zps, hOld, phiOld, theta, sin3, cos3;
354
355 ss = sqrt(XYZ[0]*XYZ[0]+XYZ[1]*XYZ[1]) ;
356 zps = XYZ[2]/ss ;
357 theta = atan( (XYZ[2]*t_CST::aell) / (ss*bell) );
358 sin3 = sin(theta) * sin(theta) * sin(theta);
359 cos3 = cos(theta) * cos(theta) * cos(theta);
360
361 // Closed formula
362 Ell[0] = atan( (XYZ[2] + e2c * bell * sin3) / (ss - e2 * t_CST::aell * cos3) );
363 Ell[1] = atan2(XYZ[1],XYZ[0]) ;
364 nn = t_CST::aell/sqrt(1.0-e2*sin(Ell[0])*sin(Ell[0])) ;
365 Ell[2] = ss / cos(Ell[0]) - nn;
366
367 const int MAXITER = 100;
368 for (int ii = 1; ii <= MAXITER; ii++) {
369 nn = t_CST::aell/sqrt(1.0-e2*sin(Ell[0])*sin(Ell[0])) ;
370 hOld = Ell[2] ;
371 phiOld = Ell[0] ;
372 Ell[2] = ss/cos(Ell[0])-nn ;
373 Ell[0] = atan(zps/(1.0-e2*nn/(nn+Ell[2]))) ;
374 if ( fabs(phiOld-Ell[0]) <= 1.0e-11 && fabs(hOld-Ell[2]) <= 1.0e-5 ) {
375 return success;
376 }
377 }
378
379 return failure;
380}
381
382// Rectangular Coordinates -> North, East, Up Components
383////////////////////////////////////////////////////////////////////////////
384void xyz2neu(const double* Ell, const double* xyz, double* neu) {
385
386 double sinPhi = sin(Ell[0]);
387 double cosPhi = cos(Ell[0]);
388 double sinLam = sin(Ell[1]);
389 double cosLam = cos(Ell[1]);
390
391 neu[0] = - sinPhi*cosLam * xyz[0]
392 - sinPhi*sinLam * xyz[1]
393 + cosPhi * xyz[2];
394
395 neu[1] = - sinLam * xyz[0]
396 + cosLam * xyz[1];
397
398 neu[2] = + cosPhi*cosLam * xyz[0]
399 + cosPhi*sinLam * xyz[1]
400 + sinPhi * xyz[2];
401}
402
403// North, East, Up Components -> Rectangular Coordinates
404////////////////////////////////////////////////////////////////////////////
405void neu2xyz(const double* Ell, const double* neu, double* xyz) {
406
407 double sinPhi = sin(Ell[0]);
408 double cosPhi = cos(Ell[0]);
409 double sinLam = sin(Ell[1]);
410 double cosLam = cos(Ell[1]);
411
412 xyz[0] = - sinPhi*cosLam * neu[0]
413 - sinLam * neu[1]
414 + cosPhi*cosLam * neu[2];
415
416 xyz[1] = - sinPhi*sinLam * neu[0]
417 + cosLam * neu[1]
418 + cosPhi*sinLam * neu[2];
419
420 xyz[2] = + cosPhi * neu[0]
421 + sinPhi * neu[2];
422}
423
424// Rectangular Coordinates -> Geocentric Coordinates
425////////////////////////////////////////////////////////////////////////////
426t_irc xyz2geoc(const double* XYZ, double* Geoc) {
427
428 const double bell = t_CST::aell*(1.0-1.0/t_CST::fInv) ;
429 const double e2 = (t_CST::aell*t_CST::aell-bell*bell)/(t_CST::aell*t_CST::aell) ;
430 double Ell[3];
431 if (xyz2ell(XYZ, Ell) != success) {
432 return failure;
433 }
434 double rho = sqrt(XYZ[0]*XYZ[0]+XYZ[1]*XYZ[1]+XYZ[2]*XYZ[2]);
435 double Rn = t_CST::aell/sqrt(1-e2*pow(sin(Ell[0]),2));
436
437 Geoc[0] = atan((1-e2 * Rn/(Rn + Ell[2])) * tan(Ell[0]));
438 Geoc[1] = Ell[1];
439 Geoc[2] = rho-t_CST::rgeoc;
440
441 return success;
442}
443
444//
445////////////////////////////////////////////////////////////////////////////
446double Frac (double x) {
447 return x-floor(x);
448}
449
450//
451////////////////////////////////////////////////////////////////////////////
452double Modulo (double x, double y) {
453 return y*Frac(x/y);
454}
455
456// Round to nearest integer
457////////////////////////////////////////////////////////////////////////////
458double nint(double val) {
459 return ((val < 0.0) ? -floor(fabs(val)+0.5) : floor(val+0.5));
460}
461
462//
463////////////////////////////////////////////////////////////////////////////
464double factorial(int n) {
465 if (n == 0) {
466 return 1;
467 }
468 else {
469 return (n * factorial(n - 1));
470 }
471}
472
473//
474////////////////////////////////////////////////////////////////////////////
475double associatedLegendreFunction(int n, int m, double t) {
476 double sum = 0.0;
477 int r = (int) floor((n - m) / 2);
478 for (int k = 0; k <= r; k++) {
479 sum += (pow(-1.0, (double)k) * factorial(2*n - 2*k)
480 / (factorial(k) * factorial(n-k) * factorial(n-m-2*k))
481 * pow(t, (double)n-m-2*k));
482 }
483 double fac = pow(2.0,(double) -n) * pow((1 - t*t), (double)m/2);
484 return sum *= fac;
485}
486
487
488// Jacobian XYZ --> NEU
489////////////////////////////////////////////////////////////////////////////
490void jacobiXYZ_NEU(const double* Ell, Matrix& jacobi) {
491
492 Tracer tracer("jacobiXYZ_NEU");
493
494 double sinPhi = sin(Ell[0]);
495 double cosPhi = cos(Ell[0]);
496 double sinLam = sin(Ell[1]);
497 double cosLam = cos(Ell[1]);
498
499 jacobi(1,1) = - sinPhi * cosLam;
500 jacobi(1,2) = - sinPhi * sinLam;
501 jacobi(1,3) = cosPhi;
502
503 jacobi(2,1) = - sinLam;
504 jacobi(2,2) = cosLam;
505 jacobi(2,3) = 0.0;
506
507 jacobi(3,1) = cosPhi * cosLam;
508 jacobi(3,2) = cosPhi * sinLam;
509 jacobi(3,3) = sinPhi;
510}
511
512// Jacobian Ell --> XYZ
513////////////////////////////////////////////////////////////////////////////
514void jacobiEll_XYZ(const double* Ell, Matrix& jacobi) {
515
516 Tracer tracer("jacobiEll_XYZ");
517
518 double sinPhi = sin(Ell[0]);
519 double cosPhi = cos(Ell[0]);
520 double sinLam = sin(Ell[1]);
521 double cosLam = cos(Ell[1]);
522 double hh = Ell[2];
523
524 double bell = t_CST::aell*(1.0-1.0/t_CST::fInv);
525 double e2 = (t_CST::aell*t_CST::aell-bell*bell)/(t_CST::aell*t_CST::aell) ;
526 double nn = t_CST::aell/sqrt(1.0-e2*sinPhi*sinPhi) ;
527
528 jacobi(1,1) = -(nn+hh) * sinPhi * cosLam;
529 jacobi(1,2) = -(nn+hh) * cosPhi * sinLam;
530 jacobi(1,3) = cosPhi * cosLam;
531
532 jacobi(2,1) = -(nn+hh) * sinPhi * sinLam;
533 jacobi(2,2) = (nn+hh) * cosPhi * cosLam;
534 jacobi(2,3) = cosPhi * sinLam;
535
536 jacobi(3,1) = (nn*(1.0-e2)+hh) * cosPhi;
537 jacobi(3,2) = 0.0;
538 jacobi(3,3) = sinPhi;
539}
540
541// Covariance Matrix in NEU
542////////////////////////////////////////////////////////////////////////////
543void covariXYZ_NEU(const SymmetricMatrix& QQxyz, const double* Ell,
544 SymmetricMatrix& Qneu) {
545
546 Tracer tracer("covariXYZ_NEU");
547
548 Matrix CC(3,3);
549 jacobiXYZ_NEU(Ell, CC);
550 Qneu << CC * QQxyz * CC.t();
551}
552
553// Covariance Matrix in XYZ
554////////////////////////////////////////////////////////////////////////////
555void covariNEU_XYZ(const SymmetricMatrix& QQneu, const double* Ell,
556 SymmetricMatrix& Qxyz) {
557
558 Tracer tracer("covariNEU_XYZ");
559
560 Matrix CC(3,3);
561 jacobiXYZ_NEU(Ell, CC);
562 Qxyz << CC.t() * QQneu * CC;
563}
564
565// Fourth order Runge-Kutta numerical integrator for ODEs
566////////////////////////////////////////////////////////////////////////////
567ColumnVector rungeKutta4(
568 double xi, // the initial x-value
569 const ColumnVector& yi, // vector of the initial y-values
570 double dx, // the step size for the integration
571 double* acc, // additional acceleration
572 ColumnVector (*der)(double x, const ColumnVector& y, double* acc)
573 // A pointer to a function that computes the
574 // derivative of a function at a point (x,y)
575 ) {
576
577 ColumnVector k1 = der(xi , yi , acc) * dx;
578 ColumnVector k2 = der(xi+dx/2.0, yi+k1/2.0, acc) * dx;
579 ColumnVector k3 = der(xi+dx/2.0, yi+k2/2.0, acc) * dx;
580 ColumnVector k4 = der(xi+dx , yi+k3 , acc) * dx;
581
582 ColumnVector yf = yi + k1/6.0 + k2/3.0 + k3/3.0 + k4/6.0;
583
584 return yf;
585}
586//
587////////////////////////////////////////////////////////////////////////////
588double djul(long jj, long mm, double tt) {
589 long ii, kk;
590 double djul ;
591 if( mm <= 2 ) {
592 jj = jj - 1;
593 mm = mm + 12;
594 }
595 ii = jj/100;
596 kk = 2 - ii + ii/4;
597 djul = (365.25*jj - fmod( 365.25*jj, 1.0 )) - 679006.0;
598 djul = djul + floor( 30.6001*(mm + 1) ) + tt + kk;
599 return djul;
600}
601
602//
603////////////////////////////////////////////////////////////////////////////
604double gpjd(double second, int nweek) {
605 double deltat;
606 deltat = nweek*7.0 + second/86400.0 ;
607 return( 44244.0 + deltat) ;
608}
609
610//
611////////////////////////////////////////////////////////////////////////////
612void jdgp(double tjul, double & second, long & nweek) {
613 double deltat;
614 deltat = tjul - 44244.0 ;
615 nweek = (long) floor(deltat/7.0);
616 second = (deltat - (nweek)*7.0)*86400.0;
617}
618
619//
620////////////////////////////////////////////////////////////////////////////
621void jmt(double djul, long& jj, long& mm, double& dd) {
622 long ih, ih1, ih2 ;
623 double t1, t2, t3, t4;
624 t1 = 1.0 + djul - fmod( djul, 1.0 ) + 2400000.0;
625 t4 = fmod( djul, 1.0 );
626 ih = long( (t1 - 1867216.25)/36524.25 );
627 t2 = t1 + 1 + ih - ih/4;
628 t3 = t2 - 1720995.0;
629 ih1 = long( (t3 - 122.1)/365.25 );
630 t1 = 365.25*ih1 - fmod( 365.25*ih1, 1.0 );
631 ih2 = long( (t3 - t1)/30.6001 );
632 dd = t3 - t1 - (int)( 30.6001*ih2 ) + t4;
633 mm = ih2 - 1;
634 if ( ih2 > 13 ) mm = ih2 - 13;
635 jj = ih1;
636 if ( mm <= 2 ) jj = jj + 1;
637}
638
639//
640////////////////////////////////////////////////////////////////////////////
641void GPSweekFromDateAndTime(const QDateTime& dateTime,
642 int& GPSWeek, double& GPSWeeks) {
643
644 static const QDateTime zeroEpoch(QDate(1980, 1, 6),QTime(),Qt::UTC);
645
646 GPSWeek = zeroEpoch.daysTo(dateTime) / 7;
647
648 int weekDay = dateTime.date().dayOfWeek() + 1; // Qt: Monday = 1
649 if (weekDay > 7) weekDay = 1;
650
651 GPSWeeks = (weekDay - 1) * 86400.0
652 - dateTime.time().msecsTo(QTime()) / 1e3;
653}
654
655//
656////////////////////////////////////////////////////////////////////////////
657void GPSweekFromYMDhms(int year, int month, int day, int hour, int min,
658 double sec, int& GPSWeek, double& GPSWeeks) {
659
660 double mjd = djul(year, month, day);
661
662 long GPSWeek_long;
663 jdgp(mjd, GPSWeeks, GPSWeek_long);
664 GPSWeek = GPSWeek_long;
665 GPSWeeks += hour * 3600.0 + min * 60.0 + sec;
666}
667
668//
669////////////////////////////////////////////////////////////////////////////
670void mjdFromDateAndTime(const QDateTime& dateTime, int& mjd, double& dayfrac) {
671
672 static const QDate zeroDate(1858, 11, 17);
673
674 mjd = zeroDate.daysTo(dateTime.date());
675
676 dayfrac = (dateTime.time().hour() +
677 (dateTime.time().minute() +
678 (dateTime.time().second() +
679 dateTime.time().msec() / 1000.0) / 60.0) / 60.0) / 24.0;
680}
681
682//
683////////////////////////////////////////////////////////////////////////////
684bool findInVector(const vector<QString>& vv, const QString& str) {
685 std::vector<QString>::const_iterator it;
686 for (it = vv.begin(); it != vv.end(); ++it) {
687 if ( (*it) == str) {
688 return true;
689 }
690 }
691 return false;
692}
693
694//
695////////////////////////////////////////////////////////////////////////////
696int readInt(const QString& str, int pos, int len, int& value) {
697 bool ok;
698 value = str.mid(pos, len).toInt(&ok);
699 return ok ? 0 : 1;
700}
701
702//
703////////////////////////////////////////////////////////////////////////////
704int readDbl(const QString& str, int pos, int len, double& value) {
705 QString hlp = str.mid(pos, len);
706 for (int ii = 0; ii < hlp.length(); ii++) {
707 if (hlp[ii]=='D' || hlp[ii]=='d' || hlp[ii] == 'E') {
708 hlp[ii]='e';
709 }
710 }
711 bool ok;
712 value = hlp.toDouble(&ok);
713 return ok ? 0 : 1;
714}
715
716// Topocentrical Distance and Elevation
717////////////////////////////////////////////////////////////////////////////
718void topos(double xRec, double yRec, double zRec,
719 double xSat, double ySat, double zSat,
720 double& rho, double& eleSat, double& azSat) {
721
722 double dx[3];
723 dx[0] = xSat-xRec;
724 dx[1] = ySat-yRec;
725 dx[2] = zSat-zRec;
726
727 rho = sqrt( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
728
729 double xyzRec[3];
730 xyzRec[0] = xRec;
731 xyzRec[1] = yRec;
732 xyzRec[2] = zRec;
733
734 double Ell[3];
735 double neu[3];
736 xyz2ell(xyzRec, Ell);
737 xyz2neu(Ell, dx, neu);
738
739 eleSat = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / rho );
740 if (neu[2] < 0) {
741 eleSat *= -1.0;
742 }
743
744 azSat = atan2(neu[1], neu[0]);
745}
746
747// Degrees -> degrees, minutes, seconds
748////////////////////////////////////////////////////////////////////////////
749void deg2DMS(double decDeg, int& deg, int& min, double& sec) {
750 int sgn = (decDeg < 0.0 ? -1 : 1);
751 deg = static_cast<int>(decDeg);
752 min = sgn * static_cast<int>((decDeg - deg)*60);
753 sec = (sgn* (decDeg - deg) - min/60.0) * 3600.0;
754}
755
756//
757////////////////////////////////////////////////////////////////////////////
758QString fortranFormat(double value, int width, int prec) {
759 int expo = value == 0.0 ? 0 : int(log10(fabs(value)));
760 double mant = value == 0.0 ? 0 : value / pow(10.0, double(expo));
761 if (fabs(mant) >= 1.0) {
762 mant /= 10.0;
763 expo += 1;
764 }
765 if (expo >= 0) {
766 return QString("%1e+%2").arg(mant, width-4, 'f', prec).arg(expo, 2, 10, QChar('0'));
767 }
768 else {
769 return QString("%1e-%2").arg(mant, width-4, 'f', prec).arg(-expo, 2, 10, QChar('0'));
770 }
771}
772
773//
774//////////////////////////////////////////////////////////////////////////////
775void kalman(const Matrix& AA, const ColumnVector& ll, const DiagonalMatrix& PP,
776 SymmetricMatrix& QQ, ColumnVector& xx) {
777
778 Tracer tracer("kalman");
779
780 int nPar = AA.Ncols();
781 int nObs = AA.Nrows();
782 UpperTriangularMatrix SS = Cholesky(QQ).t();
783
784 Matrix SA = SS*AA.t();
785 Matrix SRF(nObs+nPar, nObs+nPar); SRF = 0;
786 for (int ii = 1; ii <= nObs; ++ii) {
787 SRF(ii,ii) = 1.0 / sqrt(PP(ii,ii));
788 }
789
790 SRF.SubMatrix (nObs+1, nObs+nPar, 1, nObs) = SA;
791 SRF.SymSubMatrix(nObs+1, nObs+nPar) = SS;
792
793 UpperTriangularMatrix UU;
794 QRZ(SRF, UU);
795
796 SS = UU.SymSubMatrix(nObs+1, nObs+nPar);
797 UpperTriangularMatrix SH_rt = UU.SymSubMatrix(1, nObs);
798 Matrix YY = UU.SubMatrix(1, nObs, nObs+1, nObs+nPar);
799
800 UpperTriangularMatrix SHi = SH_rt.i();
801
802 Matrix KT = SHi * YY;
803 SymmetricMatrix Hi; Hi << SHi * SHi.t();
804
805 xx += KT.t() * (ll - AA * xx);
806 QQ << (SS.t() * SS);
807}
808
809//
810////////////////////////////////////////////////////////////////////////////
811double accuracyFromIndex(int index, t_eph::e_type type) {
812double accuracy = -1.0;
813
814 if (type == t_eph::GPS ||
815 type == t_eph::BDS ||
816 type == t_eph::SBAS||
817 type == t_eph::QZSS) {
818 if ((index >= 0) && (index <= 6)) {
819 if (index == 3) {
820 accuracy = ceil(10.0 * pow(2.0, (double(index) / 2.0) + 1.0)) / 10.0;
821 }
822 else {
823 accuracy = floor(10.0 * pow(2.0, (double(index) / 2.0) + 1.0)) / 10.0;
824 }
825 }
826 else if ((index > 6) && (index <= 15)) {
827 accuracy = (10.0 * pow(2.0, (double(index) - 2.0))) / 10.0;
828 }
829 else {
830 accuracy = 8192.0;
831 }
832 }
833 else if (type == t_eph::Galileo) {
834 if ((index >= 0) && (index <= 49)) {
835 accuracy = (double(index) / 100.0);
836 }
837 else if ((index > 49) && (index <= 74)) {
838 accuracy = (50.0 + (double(index) - 50.0) * 2.0) / 100.0;
839 }
840 else if ((index > 74) && (index <= 99)) {
841 accuracy = 1.0 + (double(index) - 75.0) * 0.04;
842 }
843 else if ((index > 99) && (index <= 125)) {
844 accuracy = 2.0 + (double(index) - 100.0) * 0.16;
845 }
846 else {
847 accuracy = -1.0;
848 }
849 }
850 else if (type == t_eph::IRNSS) {
851 if ((index >= 0) && (index <= 6)) {
852 if (index == 1) {
853 accuracy = 2.8;
854 }
855 else if (index == 3) {
856 accuracy = 5.7;
857 }
858 else if (index == 5) {
859 accuracy = 11.3;
860 }
861 else {
862 accuracy = pow(2, 1 + index / 2);
863 }
864 }
865 else if ((index > 6) && (index <= 15)) {
866 accuracy = pow(2, index - 2);
867 }
868 }
869 return accuracy;
870}
871
872//
873////////////////////////////////////////////////////////////////////////////
874int indexFromAccuracy(double accuracy, t_eph::e_type type) {
875
876 if (type == t_eph::GPS || type == t_eph::BDS || type == t_eph::SBAS
877 || type == t_eph::QZSS) {
878
879 if (accuracy <= 2.40) {
880 return 0;
881 }
882 else if (accuracy <= 3.40) {
883 return 1;
884 }
885 else if (accuracy <= 4.85) {
886 return 2;
887 }
888 else if (accuracy <= 6.85) {
889 return 3;
890 }
891 else if (accuracy <= 9.65) {
892 return 4;
893 }
894 else if (accuracy <= 13.65) {
895 return 5;
896 }
897 else if (accuracy <= 24.00) {
898 return 6;
899 }
900 else if (accuracy <= 48.00) {
901 return 7;
902 }
903 else if (accuracy <= 96.00) {
904 return 8;
905 }
906 else if (accuracy <= 192.00) {
907 return 9;
908 }
909 else if (accuracy <= 384.00) {
910 return 10;
911 }
912 else if (accuracy <= 768.00) {
913 return 11;
914 }
915 else if (accuracy <= 1536.00) {
916 return 12;
917 }
918 else if (accuracy <= 3072.00) {
919 return 13;
920 }
921 else if (accuracy <= 6144.00) {
922 return 14;
923 }
924 else {
925 return 15;
926 }
927 }
928
929 if (type == t_eph::Galileo) {
930
931 if (accuracy <= 0.49) {
932 return int(ceil(accuracy * 100.0));
933 }
934 else if (accuracy <= 0.98) {
935 return int(50.0 + (((accuracy * 100.0) - 50) / 2.0));
936 }
937 else if (accuracy <= 2.0) {
938 return int(75.0 + ((accuracy - 1.0) / 0.04));
939 }
940 else if (accuracy <= 6.0) {
941 return int(100.0 + ((accuracy - 2.0) / 0.16));
942 }
943 else {
944 return 255;
945 }
946 }
947
948 return (type == t_eph::Galileo) ? 255 : 15;
949}
950
951// Returns fit interval in hours from flag
952////////////////////////////////////////////////////////////////////////////
953double fitIntervalFromFlag(int flag, double iodc, t_eph::e_type type) {
954 double fitInterval = 0.0;
955
956 switch (flag) {
957 case 0:
958 if (type == t_eph::GPS) {
959 fitInterval = 4.0;
960 }
961 else if (type == t_eph::QZSS) {
962 fitInterval = 2.0;
963 }
964 break;
965 case 1:
966 if (type == t_eph::GPS) {
967 if (iodc >= 240 && iodc <= 247) {
968 fitInterval = 8.0;
969 }
970 else if ((iodc >= 248 && iodc <= 255) ||
971 (iodc == 496) ) {
972 fitInterval = 14.0;
973 }
974 else if ((iodc >= 497 && iodc <= 503) ||
975 (iodc >= 2021 && iodc <= 1023) ) {
976 fitInterval = 26.0;
977 }
978 else {
979 fitInterval = 6.0;
980 }
981 }
982 break;
983 }
984 return fitInterval;
985}
986
987// Returns CRC24
988////////////////////////////////////////////////////////////////////////////
989unsigned long CRC24(long size, const unsigned char *buf) {
990 unsigned long crc = 0;
991 int ii;
992 while (size--) {
993 crc ^= (*buf++) << (16);
994 for(ii = 0; ii < 8; ii++) {
995 crc <<= 1;
996 if (crc & 0x1000000) {
997 crc ^= 0x01864cfb;
998 }
999 }
1000 }
1001 return crc;
1002}
1003
1004// Convert RTCM3 lock-time indicator to lock time in seconds
1005////////////////////////////////////////////////////////////////////////////
1006double lti2sec(int type, int lti) {
1007
1008 if ( (type>=1001 && type<=1004) ||
1009 (type>=1009 && type<=1012) ) { // RTCM3 msg 100[1...4] and 10[09...12]
1010 if (lti< 0) return -1;
1011 else if (lti< 24) return 1*lti; // [ 0 1 23]
1012 else if (lti< 48) return 2*lti-24; // [ 24 2 70]
1013 else if (lti< 72) return 4*lti-120; // [ 72 4 164]
1014 else if (lti< 96) return 8*lti-408; // [168 8 352]
1015 else if (lti< 120) return 16*lti-1176; // [360 16 728]
1016 else if (lti< 127) return 32*lti-3096; // [744 32 905]
1017 else if (lti==127) return 937;
1018 else return -1;
1019 }
1020 else if (type%10==2 || type%10==3 ||
1021 type%10==4 || type%10==5) { // RTCM3 MSM-2/-3/-4/-5
1022 switch(lti) {
1023 case( 0) : return 0;
1024 case( 1) : return 32e-3;
1025 case( 2) : return 64e-3;
1026 case( 3) : return 128e-3;
1027 case( 4) : return 256e-3;
1028 case( 5) : return 512e-3;
1029 case( 6) : return 1024e-3;
1030 case( 7) : return 2048e-3;
1031 case( 8) : return 4096e-3;
1032 case( 9) : return 8192e-3;
1033 case(10) : return 16384e-3;
1034 case(11) : return 32768e-3;
1035 case(12) : return 65536e-3;
1036 case(13) : return 131072e-3;
1037 case(14) : return 262144e-3;
1038 case(15) : return 524288e-3;
1039 default : return -1;
1040 };
1041 }
1042 else if (type%10==6 || type%10==7) { // RTCM3 MSM-6 and MSM-7
1043 if (lti< 0) return ( -1 );
1044 else if (lti< 64) return ( 1*lti )*1e-3;
1045 else if (lti< 96) return ( 2*lti-64 )*1e-3;
1046 else if (lti< 128) return ( 4*lti-256 )*1e-3;
1047 else if (lti< 160) return ( 8*lti-768 )*1e-3;
1048 else if (lti< 192) return ( 16*lti-2048 )*1e-3;
1049 else if (lti< 224) return ( 32*lti-5120 )*1e-3;
1050 else if (lti< 256) return ( 64*lti-12288 )*1e-3;
1051 else if (lti< 288) return ( 128*lti-28672 )*1e-3;
1052 else if (lti< 320) return ( 256*lti-65536 )*1e-3;
1053 else if (lti< 352) return ( 512*lti-147456 )*1e-3;
1054 else if (lti< 384) return ( 1024*lti-327680 )*1e-3;
1055 else if (lti< 416) return ( 2048*lti-720896 )*1e-3;
1056 else if (lti< 448) return ( 4096*lti-1572864 )*1e-3;
1057 else if (lti< 480) return ( 8192*lti-3407872 )*1e-3;
1058 else if (lti< 512) return ( 16384*lti-7340032 )*1e-3;
1059 else if (lti< 544) return ( 32768*lti-15728640 )*1e-3;
1060 else if (lti< 576) return ( 65536*lti-33554432 )*1e-3;
1061 else if (lti< 608) return ( 131072*lti-71303168 )*1e-3;
1062 else if (lti< 640) return ( 262144*lti-150994944 )*1e-3;
1063 else if (lti< 672) return ( 524288*lti-318767104 )*1e-3;
1064 else if (lti< 704) return (1048576*lti-671088640 )*1e-3;
1065 else if (lti==704) return (2097152*lti-1409286144)*1e-3;
1066 else return ( -1 );
1067 }
1068 else {
1069 return -1;
1070 };
1071};
Note: See TracBrowser for help on using the repository browser.