source: ntrip/trunk/BNC/src/bncutils.cpp@ 7288

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

minor changes

File size: 26.1 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{0,0,0,0} /* end marker */
92};
93
94#define GPSLEAPSTART 19 /* 19 leap seconds existed at 6.1.1980 */
95
96static int longyear(int year, int month)
97{
98 if(!(year % 4) && (!(year % 400) || (year % 100)))
99 {
100 if(!month || month == 2)
101 return 1;
102 }
103 return 0;
104}
105
106int gnumleap(int year, int month, int day)
107{
108 int ls = 0;
109 const struct leapseconds *l;
110
111 for(l = leap; l->taicount && year >= l->year; ++l)
112 {
113 if(year > l->year || month > l->month || (month == l->month && day > l->day))
114 ls = l->taicount - GPSLEAPSTART;
115 }
116 return ls;
117}
118
119/* Convert Moscow time into UTC (fixnumleap == 1) or GPS (fixnumleap == 0) */
120void updatetime(int *week, int *secOfWeek, int mSecOfWeek, bool fixnumleap)
121{
122 int y,m,d,k,l, nul;
123 unsigned int j = *week*(7*24*60*60) + *secOfWeek + 5*24*60*60+3*60*60;
124 int glo_daynumber = 0, glo_timeofday;
125 for(y = 1980; j >= (unsigned int)(k = (l = (365+longyear(y,0)))*24*60*60)
126 + gnumleap(y+1,1,1); ++y)
127 {
128 j -= k; glo_daynumber += l;
129 }
130 for(m = 1; j >= (unsigned int)(k = (l = months[m]+longyear(y, m))*24*60*60)
131 + gnumleap(y, m+1, 1); ++m)
132 {
133 j -= k; glo_daynumber += l;
134 }
135 for(d = 1; j >= 24UL*60UL*60UL + gnumleap(y, m, d+1); ++d)
136 j -= 24*60*60;
137 glo_daynumber -= 16*365+4-d;
138 nul = gnumleap(y, m, d);
139 glo_timeofday = j-nul;
140
141 // original version
142 // if(mSecOfWeek < 5*60*1000 && glo_timeofday > 23*60*60)
143 // *secOfWeek += 24*60*60;
144 // else if(glo_timeofday < 5*60 && mSecOfWeek > 23*60*60*1000)
145 // *secOfWeek -= 24*60*60;
146
147 // new version
148 if(mSecOfWeek < 4*60*60*1000 && glo_timeofday > 20*60*60)
149 *secOfWeek += 24*60*60;
150 else if(glo_timeofday < 4*60*60 && mSecOfWeek > 20*60*60*1000)
151 *secOfWeek -= 24*60*60;
152
153 *secOfWeek += mSecOfWeek/1000-glo_timeofday;
154 if(fixnumleap)
155 *secOfWeek -= nul;
156 if(*secOfWeek < 0) {*secOfWeek += 24*60*60*7; --*week; }
157 if(*secOfWeek >= 24*60*60*7) {*secOfWeek -= 24*60*60*7; ++*week; }
158}
159
160//
161////////////////////////////////////////////////////////////////////////////
162void expandEnvVar(QString& str) {
163
164 QRegExp rx("(\\$\\{.+\\})");
165
166 if (rx.indexIn(str) != -1) {
167 QStringListIterator it(rx.capturedTexts());
168 if (it.hasNext()) {
169 QString rxStr = it.next();
170 QString envVar = rxStr.mid(2,rxStr.length()-3);
171 str.replace(rxStr, qgetenv(envVar.toAscii()));
172 }
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////////////////////////////////////////////////////////////////////////////
250QByteArray ggaString(const QByteArray& latitude,
251 const QByteArray& longitude,
252 const QByteArray& height,
253 const QString& ggaType) {
254
255 double lat = strtod(latitude,NULL);
256 double lon = strtod(longitude,NULL);
257 double hei = strtod(height,NULL);
258 QString sentences = "GPGGA,";
259 if (ggaType.contains("GNGGA")) {
260 sentences = "GNGGA,";
261 }
262
263 const char* flagN="N";
264 const char* flagE="E";
265 if (lon >180.) {lon=(lon-360.)*(-1.); flagE="W";}
266 if ((lon < 0.) && (lon >= -180.)) {lon=lon*(-1.); flagE="W";}
267 if (lon < -180.) {lon=(lon+360.); flagE="E";}
268 if (lat < 0.) {lat=lat*(-1.); flagN="S";}
269 QTime ttime(QDateTime::currentDateTime().toUTC().time());
270 int lat_deg = (int)lat;
271 double lat_min=(lat-lat_deg)*60.;
272 int lon_deg = (int)lon;
273 double lon_min=(lon-lon_deg)*60.;
274 int hh = 0 , mm = 0;
275 double ss = 0.0;
276 hh=ttime.hour();
277 mm=ttime.minute();
278 ss=(double)ttime.second()+0.001*ttime.msec();
279 QString gga;
280 gga += sentences;
281 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'));
282 gga += QString("%1%2,").arg((int)lat_deg,2, 10, QLatin1Char('0')).arg(lat_min, 7, 'f', 4, QLatin1Char('0'));
283 gga += flagN;
284 gga += QString(",%1%2,").arg((int)lon_deg,3, 10, QLatin1Char('0')).arg(lon_min, 7, 'f', 4, QLatin1Char('0'));
285 gga += flagE + QString(",1,05,1.00");
286 gga += QString(",%1,").arg(hei, 2, 'f', 1);
287 gga += QString("M,10.000,M,,");
288 int xori;
289 char XOR = 0;
290 char *Buff =gga.toAscii().data();
291 int iLen = strlen(Buff);
292 for (xori = 0; xori < iLen; xori++) {
293 XOR ^= (char)Buff[xori];
294 }
295 gga = "$" + gga + QString("*%1").arg(XOR, 2, 16, QLatin1Char('0'));
296
297 return gga.toAscii();
298}
299
300//
301////////////////////////////////////////////////////////////////////////////
302void RSW_to_XYZ(const ColumnVector& rr, const ColumnVector& vv,
303 const ColumnVector& rsw, ColumnVector& xyz) {
304
305 ColumnVector along = vv / vv.norm_Frobenius();
306 ColumnVector cross = crossproduct(rr, vv); cross /= cross.norm_Frobenius();
307 ColumnVector radial = crossproduct(along, cross);
308
309 Matrix RR(3,3);
310 RR.Column(1) = radial;
311 RR.Column(2) = along;
312 RR.Column(3) = cross;
313
314 xyz = RR * rsw;
315}
316
317// Transformation xyz --> radial, along track, out-of-plane
318////////////////////////////////////////////////////////////////////////////
319void XYZ_to_RSW(const ColumnVector& rr, const ColumnVector& vv,
320 const ColumnVector& xyz, ColumnVector& rsw) {
321
322 ColumnVector along = vv / vv.norm_Frobenius();
323 ColumnVector cross = crossproduct(rr, vv); cross /= cross.norm_Frobenius();
324 ColumnVector radial = crossproduct(along, cross);
325
326 rsw.ReSize(3);
327 rsw(1) = DotProduct(xyz, radial);
328 rsw(2) = DotProduct(xyz, along);
329 rsw(3) = DotProduct(xyz, cross);
330}
331
332// Rectangular Coordinates -> Ellipsoidal Coordinates
333////////////////////////////////////////////////////////////////////////////
334t_irc xyz2ell(const double* XYZ, double* Ell) {
335
336 const double bell = t_CST::aell*(1.0-1.0/t_CST::fInv) ;
337 const double e2 = (t_CST::aell*t_CST::aell-bell*bell)/(t_CST::aell*t_CST::aell) ;
338 const double e2c = (t_CST::aell*t_CST::aell-bell*bell)/(bell*bell) ;
339
340 double nn, ss, zps, hOld, phiOld, theta, sin3, cos3;
341
342 ss = sqrt(XYZ[0]*XYZ[0]+XYZ[1]*XYZ[1]) ;
343 zps = XYZ[2]/ss ;
344 theta = atan( (XYZ[2]*t_CST::aell) / (ss*bell) );
345 sin3 = sin(theta) * sin(theta) * sin(theta);
346 cos3 = cos(theta) * cos(theta) * cos(theta);
347
348 // Closed formula
349 Ell[0] = atan( (XYZ[2] + e2c * bell * sin3) / (ss - e2 * t_CST::aell * cos3) );
350 Ell[1] = atan2(XYZ[1],XYZ[0]) ;
351 nn = t_CST::aell/sqrt(1.0-e2*sin(Ell[0])*sin(Ell[0])) ;
352 Ell[2] = ss / cos(Ell[0]) - nn;
353
354 const int MAXITER = 100;
355 for (int ii = 1; ii <= MAXITER; ii++) {
356 nn = t_CST::aell/sqrt(1.0-e2*sin(Ell[0])*sin(Ell[0])) ;
357 hOld = Ell[2] ;
358 phiOld = Ell[0] ;
359 Ell[2] = ss/cos(Ell[0])-nn ;
360 Ell[0] = atan(zps/(1.0-e2*nn/(nn+Ell[2]))) ;
361 if ( fabs(phiOld-Ell[0]) <= 1.0e-11 && fabs(hOld-Ell[2]) <= 1.0e-5 ) {
362 return success;
363 }
364 }
365
366 return failure;
367}
368
369// Rectangular Coordinates -> North, East, Up Components
370////////////////////////////////////////////////////////////////////////////
371void xyz2neu(const double* Ell, const double* xyz, double* neu) {
372
373 double sinPhi = sin(Ell[0]);
374 double cosPhi = cos(Ell[0]);
375 double sinLam = sin(Ell[1]);
376 double cosLam = cos(Ell[1]);
377
378 neu[0] = - sinPhi*cosLam * xyz[0]
379 - sinPhi*sinLam * xyz[1]
380 + cosPhi * xyz[2];
381
382 neu[1] = - sinLam * xyz[0]
383 + cosLam * xyz[1];
384
385 neu[2] = + cosPhi*cosLam * xyz[0]
386 + cosPhi*sinLam * xyz[1]
387 + sinPhi * xyz[2];
388}
389
390// North, East, Up Components -> Rectangular Coordinates
391////////////////////////////////////////////////////////////////////////////
392void neu2xyz(const double* Ell, const double* neu, double* xyz) {
393
394 double sinPhi = sin(Ell[0]);
395 double cosPhi = cos(Ell[0]);
396 double sinLam = sin(Ell[1]);
397 double cosLam = cos(Ell[1]);
398
399 xyz[0] = - sinPhi*cosLam * neu[0]
400 - sinLam * neu[1]
401 + cosPhi*cosLam * neu[2];
402
403 xyz[1] = - sinPhi*sinLam * neu[0]
404 + cosLam * neu[1]
405 + cosPhi*sinLam * neu[2];
406
407 xyz[2] = + cosPhi * neu[0]
408 + sinPhi * neu[2];
409}
410
411// Rectangular Coordinates -> Geocentric Coordinates
412////////////////////////////////////////////////////////////////////////////
413t_irc xyz2geoc(const double* XYZ, double* Geoc) {
414
415 const double bell = t_CST::aell*(1.0-1.0/t_CST::fInv) ;
416 const double e2 = (t_CST::aell*t_CST::aell-bell*bell)/(t_CST::aell*t_CST::aell) ;
417 double Ell[3];
418 if (xyz2ell(XYZ, Ell) != success) {
419 return failure;
420 }
421 double rho = sqrt(XYZ[0]*XYZ[0]+XYZ[1]*XYZ[1]+XYZ[2]*XYZ[2]);
422 double Rn = t_CST::aell/sqrt(1-e2*pow(sin(Ell[0]),2));
423
424 Geoc[0] = atan((1-e2 * Rn/(Rn + Ell[2])) * tan(Ell[0]));
425 Geoc[1] = Ell[1];
426 Geoc[2] = rho-t_CST::rgeoc;
427
428 return success;
429}
430
431//
432////////////////////////////////////////////////////////////////////////////
433double Frac (double x) {
434 return x-floor(x);
435}
436
437//
438////////////////////////////////////////////////////////////////////////////
439double Modulo (double x, double y) {
440 return y*Frac(x/y);
441}
442
443// Round to nearest integer
444////////////////////////////////////////////////////////////////////////////
445double nint(double val) {
446 return ((val < 0.0) ? -floor(fabs(val)+0.5) : floor(val+0.5));
447}
448
449//
450////////////////////////////////////////////////////////////////////////////
451int factorial(int n) {
452 if (n == 0) {
453 return 1;
454 }
455 else {
456 return (n * factorial(n - 1));
457 }
458}
459
460//
461////////////////////////////////////////////////////////////////////////////
462double associatedLegendreFunction(int n, int m, double t) {
463 double sum = 0.0;
464 int r = (int) floor((n - m) / 2);
465 for (int k = 0; k <= r; k++) {
466 sum += (pow(-1.0, (double)k) * factorial(2*n - 2*k)
467 / (factorial(k) * factorial(n-k) * factorial(n-m-2*k))
468 * pow(t, (double)n-m-2*k));
469 }
470 double fac = pow(2.0,(double) -n) * pow((1 - t*t), (double)m/2);
471 return sum *= fac;
472}
473
474
475// Jacobian XYZ --> NEU
476////////////////////////////////////////////////////////////////////////////
477void jacobiXYZ_NEU(const double* Ell, Matrix& jacobi) {
478
479 Tracer tracer("jacobiXYZ_NEU");
480
481 double sinPhi = sin(Ell[0]);
482 double cosPhi = cos(Ell[0]);
483 double sinLam = sin(Ell[1]);
484 double cosLam = cos(Ell[1]);
485
486 jacobi(1,1) = - sinPhi * cosLam;
487 jacobi(1,2) = - sinPhi * sinLam;
488 jacobi(1,3) = cosPhi;
489
490 jacobi(2,1) = - sinLam;
491 jacobi(2,2) = cosLam;
492 jacobi(2,3) = 0.0;
493
494 jacobi(3,1) = cosPhi * cosLam;
495 jacobi(3,2) = cosPhi * sinLam;
496 jacobi(3,3) = sinPhi;
497}
498
499// Jacobian Ell --> XYZ
500////////////////////////////////////////////////////////////////////////////
501void jacobiEll_XYZ(const double* Ell, Matrix& jacobi) {
502
503 Tracer tracer("jacobiEll_XYZ");
504
505 double sinPhi = sin(Ell[0]);
506 double cosPhi = cos(Ell[0]);
507 double sinLam = sin(Ell[1]);
508 double cosLam = cos(Ell[1]);
509 double hh = Ell[2];
510
511 double bell = t_CST::aell*(1.0-1.0/t_CST::fInv);
512 double e2 = (t_CST::aell*t_CST::aell-bell*bell)/(t_CST::aell*t_CST::aell) ;
513 double nn = t_CST::aell/sqrt(1.0-e2*sinPhi*sinPhi) ;
514
515 jacobi(1,1) = -(nn+hh) * sinPhi * cosLam;
516 jacobi(1,2) = -(nn+hh) * cosPhi * sinLam;
517 jacobi(1,3) = cosPhi * cosLam;
518
519 jacobi(2,1) = -(nn+hh) * sinPhi * sinLam;
520 jacobi(2,2) = (nn+hh) * cosPhi * cosLam;
521 jacobi(2,3) = cosPhi * sinLam;
522
523 jacobi(3,1) = (nn*(1.0-e2)+hh) * cosPhi;
524 jacobi(3,2) = 0.0;
525 jacobi(3,3) = sinPhi;
526}
527
528// Covariance Matrix in NEU
529////////////////////////////////////////////////////////////////////////////
530void covariXYZ_NEU(const SymmetricMatrix& QQxyz, const double* Ell,
531 SymmetricMatrix& Qneu) {
532
533 Tracer tracer("covariXYZ_NEU");
534
535 Matrix CC(3,3);
536 jacobiXYZ_NEU(Ell, CC);
537 Qneu << CC * QQxyz * CC.t();
538}
539
540// Covariance Matrix in XYZ
541////////////////////////////////////////////////////////////////////////////
542void covariNEU_XYZ(const SymmetricMatrix& QQneu, const double* Ell,
543 SymmetricMatrix& Qxyz) {
544
545 Tracer tracer("covariNEU_XYZ");
546
547 Matrix CC(3,3);
548 jacobiXYZ_NEU(Ell, CC);
549 Qxyz << CC.t() * QQneu * CC;
550}
551
552// Fourth order Runge-Kutta numerical integrator for ODEs
553////////////////////////////////////////////////////////////////////////////
554ColumnVector rungeKutta4(
555 double xi, // the initial x-value
556 const ColumnVector& yi, // vector of the initial y-values
557 double dx, // the step size for the integration
558 double* acc, // aditional acceleration
559 ColumnVector (*der)(double x, const ColumnVector& y, double* acc)
560 // A pointer to a function that computes the
561 // derivative of a function at a point (x,y)
562 ) {
563
564 ColumnVector k1 = der(xi , yi , acc) * dx;
565 ColumnVector k2 = der(xi+dx/2.0, yi+k1/2.0, acc) * dx;
566 ColumnVector k3 = der(xi+dx/2.0, yi+k2/2.0, acc) * dx;
567 ColumnVector k4 = der(xi+dx , yi+k3 , acc) * dx;
568
569 ColumnVector yf = yi + k1/6.0 + k2/3.0 + k3/3.0 + k4/6.0;
570
571 return yf;
572}
573//
574////////////////////////////////////////////////////////////////////////////
575double djul(long jj, long mm, double tt) {
576 long ii, kk;
577 double djul ;
578 if( mm <= 2 ) {
579 jj = jj - 1;
580 mm = mm + 12;
581 }
582 ii = jj/100;
583 kk = 2 - ii + ii/4;
584 djul = (365.25*jj - fmod( 365.25*jj, 1.0 )) - 679006.0;
585 djul = djul + floor( 30.6001*(mm + 1) ) + tt + kk;
586 return djul;
587}
588
589//
590////////////////////////////////////////////////////////////////////////////
591double gpjd(double second, int nweek) {
592 double deltat;
593 deltat = nweek*7.0 + second/86400.0 ;
594 return( 44244.0 + deltat) ;
595}
596
597//
598////////////////////////////////////////////////////////////////////////////
599void jdgp(double tjul, double & second, long & nweek) {
600 double deltat;
601 deltat = tjul - 44244.0 ;
602 nweek = (long) floor(deltat/7.0);
603 second = (deltat - (nweek)*7.0)*86400.0;
604}
605
606//
607////////////////////////////////////////////////////////////////////////////
608void jmt(double djul, long& jj, long& mm, double& dd) {
609 long ih, ih1, ih2 ;
610 double t1, t2, t3, t4;
611 t1 = 1.0 + djul - fmod( djul, 1.0 ) + 2400000.0;
612 t4 = fmod( djul, 1.0 );
613 ih = long( (t1 - 1867216.25)/36524.25 );
614 t2 = t1 + 1 + ih - ih/4;
615 t3 = t2 - 1720995.0;
616 ih1 = long( (t3 - 122.1)/365.25 );
617 t1 = 365.25*ih1 - fmod( 365.25*ih1, 1.0 );
618 ih2 = long( (t3 - t1)/30.6001 );
619 dd = t3 - t1 - (int)( 30.6001*ih2 ) + t4;
620 mm = ih2 - 1;
621 if ( ih2 > 13 ) mm = ih2 - 13;
622 jj = ih1;
623 if ( mm <= 2 ) jj = jj + 1;
624}
625
626//
627////////////////////////////////////////////////////////////////////////////
628void GPSweekFromDateAndTime(const QDateTime& dateTime,
629 int& GPSWeek, double& GPSWeeks) {
630
631 static const QDateTime zeroEpoch(QDate(1980, 1, 6),QTime(),Qt::UTC);
632
633 GPSWeek = zeroEpoch.daysTo(dateTime) / 7;
634
635 int weekDay = dateTime.date().dayOfWeek() + 1; // Qt: Monday = 1
636 if (weekDay > 7) weekDay = 1;
637
638 GPSWeeks = (weekDay - 1) * 86400.0
639 - dateTime.time().msecsTo(QTime()) / 1e3;
640}
641
642//
643////////////////////////////////////////////////////////////////////////////
644void GPSweekFromYMDhms(int year, int month, int day, int hour, int min,
645 double sec, int& GPSWeek, double& GPSWeeks) {
646
647 double mjd = djul(year, month, day);
648
649 long GPSWeek_long;
650 jdgp(mjd, GPSWeeks, GPSWeek_long);
651 GPSWeek = GPSWeek_long;
652 GPSWeeks += hour * 3600.0 + min * 60.0 + sec;
653}
654
655//
656////////////////////////////////////////////////////////////////////////////
657void mjdFromDateAndTime(const QDateTime& dateTime, int& mjd, double& dayfrac) {
658
659 static const QDate zeroDate(1858, 11, 17);
660
661 mjd = zeroDate.daysTo(dateTime.date());
662
663 dayfrac = (dateTime.time().hour() +
664 (dateTime.time().minute() +
665 (dateTime.time().second() +
666 dateTime.time().msec() / 1000.0) / 60.0) / 60.0) / 24.0;
667}
668
669//
670////////////////////////////////////////////////////////////////////////////
671bool findInVector(const vector<QString>& vv, const QString& str) {
672 std::vector<QString>::const_iterator it;
673 for (it = vv.begin(); it != vv.end(); ++it) {
674 if ( (*it) == str) {
675 return true;
676 }
677 }
678 return false;
679}
680
681//
682////////////////////////////////////////////////////////////////////////////
683int readInt(const QString& str, int pos, int len, int& value) {
684 bool ok;
685 value = str.mid(pos, len).toInt(&ok);
686 return ok ? 0 : 1;
687}
688
689//
690////////////////////////////////////////////////////////////////////////////
691int readDbl(const QString& str, int pos, int len, double& value) {
692 QString hlp = str.mid(pos, len);
693 for (int ii = 0; ii < hlp.length(); ii++) {
694 if (hlp[ii]=='D' || hlp[ii]=='d' || hlp[ii] == 'E') {
695 hlp[ii]='e';
696 }
697 }
698 bool ok;
699 value = hlp.toDouble(&ok);
700 return ok ? 0 : 1;
701}
702
703// Topocentrical Distance and Elevation
704////////////////////////////////////////////////////////////////////////////
705void topos(double xRec, double yRec, double zRec,
706 double xSat, double ySat, double zSat,
707 double& rho, double& eleSat, double& azSat) {
708
709 double dx[3];
710 dx[0] = xSat-xRec;
711 dx[1] = ySat-yRec;
712 dx[2] = zSat-zRec;
713
714 rho = sqrt( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
715
716 double xyzRec[3];
717 xyzRec[0] = xRec;
718 xyzRec[1] = yRec;
719 xyzRec[2] = zRec;
720
721 double Ell[3];
722 double neu[3];
723 xyz2ell(xyzRec, Ell);
724 xyz2neu(Ell, dx, neu);
725
726 eleSat = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / rho );
727 if (neu[2] < 0) {
728 eleSat *= -1.0;
729 }
730
731 azSat = atan2(neu[1], neu[0]);
732}
733
734// Degrees -> degrees, minutes, seconds
735////////////////////////////////////////////////////////////////////////////
736void deg2DMS(double decDeg, int& deg, int& min, double& sec) {
737 int sgn = (decDeg < 0.0 ? -1 : 1);
738 deg = sgn * static_cast<int>(decDeg);
739 min = static_cast<int>((decDeg - deg)*60);
740 sec = (decDeg - deg - min/60.0) * 3600.0;
741}
742
743//
744////////////////////////////////////////////////////////////////////////////
745QString fortranFormat(double value, int width, int prec) {
746 int expo = value == 0.0 ? 0 : int(log10(fabs(value)));
747 double mant = value == 0.0 ? 0 : value / pow(10, expo);
748 if (fabs(mant) >= 1.0) {
749 mant /= 10.0;
750 expo += 1;
751 }
752 if (expo >= 0) {
753 return QString("%1e+%2").arg(mant, width-4, 'f', prec).arg(expo, 2, 10, QChar('0'));
754 }
755 else {
756 return QString("%1e-%2").arg(mant, width-4, 'f', prec).arg(-expo, 2, 10, QChar('0'));
757 }
758}
759
760//
761//////////////////////////////////////////////////////////////////////////////
762void kalman(const Matrix& AA, const ColumnVector& ll, const DiagonalMatrix& PP,
763 SymmetricMatrix& QQ, ColumnVector& xx) {
764
765 Tracer tracer("kalman");
766
767 int nPar = AA.Ncols();
768 int nObs = AA.Nrows();
769 UpperTriangularMatrix SS = Cholesky(QQ).t();
770
771 Matrix SA = SS*AA.t();
772 Matrix SRF(nObs+nPar, nObs+nPar); SRF = 0;
773 for (int ii = 1; ii <= nObs; ++ii) {
774 SRF(ii,ii) = 1.0 / sqrt(PP(ii,ii));
775 }
776
777 SRF.SubMatrix (nObs+1, nObs+nPar, 1, nObs) = SA;
778 SRF.SymSubMatrix(nObs+1, nObs+nPar) = SS;
779
780 UpperTriangularMatrix UU;
781 QRZ(SRF, UU);
782
783 SS = UU.SymSubMatrix(nObs+1, nObs+nPar);
784 UpperTriangularMatrix SH_rt = UU.SymSubMatrix(1, nObs);
785 Matrix YY = UU.SubMatrix(1, nObs, nObs+1, nObs+nPar);
786
787 UpperTriangularMatrix SHi = SH_rt.i();
788
789 Matrix KT = SHi * YY;
790 SymmetricMatrix Hi; Hi << SHi * SHi.t();
791
792 xx += KT.t() * (ll - AA * xx);
793 QQ << (SS.t() * SS);
794}
795
796double accuracyFromIndex(int index, t_eph::e_type type) {
797
798 if (type == t_eph::GPS || type == t_eph::BDS || type == t_eph::SBAS
799 || type == t_eph::QZSS) {
800
801 if ((index >= 0) && (index <= 6)) {
802 if (index == 3) {
803 return ceil(10.0 * pow(2.0, (double(index) / 2.0) + 1.0)) / 10.0;
804 }
805 else {
806 return floor(10.0 * pow(2.0, (double(index) / 2.0) + 1.0)) / 10.0;
807 }
808 }
809 else if ((index > 6) && (index <= 15)) {
810 return (10.0 * pow(2.0, (double(index) - 2.0))) / 10.0;
811 }
812 else {
813 return 8192.0;
814 }
815 }
816
817 if (type == t_eph::Galileo) {
818
819 if ((index >= 0) && (index <= 49)) {
820 return (double(index) / 100.0);
821 }
822 else if ((index > 49) && (index <= 74)) {
823 return (50.0 + (double(index) - 50.0) * 2.0) / 100.0;
824 }
825 else if ((index > 74) && (index <= 99)) {
826 return 1.0 + (double(index) - 75.0) * 0.04;
827 }
828 else if ((index > 99) && (index <= 125)) {
829 return 2.0 + (double(index) - 100.0) * 0.16;
830 }
831 else {
832 return -1.0;
833 }
834 }
835
836 return double(index);
837}
838
839int indexFromAccuracy(double accuracy, t_eph::e_type type) {
840
841 if (type == t_eph::GPS || type == t_eph::BDS || type == t_eph::SBAS
842 || type == t_eph::QZSS) {
843
844 if (accuracy <= 2.40) {
845 return 0;
846 }
847 else if (accuracy <= 3.40) {
848 return 1;
849 }
850 else if (accuracy <= 4.85) {
851 return 2;
852 }
853 else if (accuracy <= 6.85) {
854 return 3;
855 }
856 else if (accuracy <= 9.65) {
857 return 4;
858 }
859 else if (accuracy <= 13.65) {
860 return 5;
861 }
862 else if (accuracy <= 24.00) {
863 return 6;
864 }
865 else if (accuracy <= 48.00) {
866 return 7;
867 }
868 else if (accuracy <= 96.00) {
869 return 8;
870 }
871 else if (accuracy <= 192.00) {
872 return 9;
873 }
874 else if (accuracy <= 384.00) {
875 return 10;
876 }
877 else if (accuracy <= 768.00) {
878 return 11;
879 }
880 else if (accuracy <= 1536.00) {
881 return 12;
882 }
883 else if (accuracy <= 3072.00) {
884 return 13;
885 }
886 else if (accuracy <= 6144.00) {
887 return 14;
888 }
889 else {
890 return 15;
891 }
892 }
893
894 if (type == t_eph::Galileo) {
895
896 if (accuracy <= 0.49) {
897 return int(ceil(accuracy * 100.0));
898 }
899 else if (accuracy <= 0.98) {
900 return int(50.0 + (((accuracy * 100.0) - 50) / 2.0));
901 }
902 else if (accuracy <= 2.0) {
903 return int(75.0 + ((accuracy - 1.0) / 0.04));
904 }
905 else if (accuracy <= 6.0) {
906 return int(100.0 + ((accuracy - 2.0) / 0.16));
907 }
908 else {
909 return 255;
910 }
911 }
912
913 return (type == t_eph::Galileo) ? 255 : 15;
914}
915
916// Returns CRC24
917////////////////////////////////////////////////////////////////////////////
918unsigned long CRC24(long size, const unsigned char *buf) {
919 unsigned long crc = 0;
920 int ii;
921 while (size--) {
922 crc ^= (*buf++) << (16);
923 for(ii = 0; ii < 8; ii++) {
924 crc <<= 1;
925 if (crc & 0x1000000) {
926 crc ^= 0x01864cfb;
927 }
928 }
929 }
930 return crc;
931}
932
Note: See TracBrowser for help on using the repository browser.