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

Last change on this file since 5846 was 5846, checked in by mervart, 8 years ago
File size: 17.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
56//
57////////////////////////////////////////////////////////////////////////////
58void expandEnvVar(QString& str) {
59
60 QRegExp rx("(\\$\\{.+\\})");
61
62 if (rx.indexIn(str) != -1) {
63 QStringListIterator it(rx.capturedTexts());
64 if (it.hasNext()) {
65 QString rxStr = it.next();
66 QString envVar = rxStr.mid(2,rxStr.length()-3);
67 str.replace(rxStr, qgetenv(envVar.toAscii()));
68 }
69 }
70
71}
72
73//
74////////////////////////////////////////////////////////////////////////////
75QDateTime dateAndTimeFromGPSweek(int GPSWeek, double GPSWeeks) {
76
77 static const QDate zeroEpoch(1980, 1, 6);
78
79 QDate date(zeroEpoch);
80 QTime time(0,0,0,0);
81
82 int weekDays = int(GPSWeeks) / 86400;
83 date = date.addDays( GPSWeek * 7 + weekDays );
84 time = time.addMSecs( int( (GPSWeeks - 86400 * weekDays) * 1e3 ) );
85
86 return QDateTime(date,time);
87}
88
89//
90////////////////////////////////////////////////////////////////////////////
91void currentGPSWeeks(int& week, double& sec) {
92
93 QDateTime currDateTimeGPS;
94
95 if ( BNC_CORE->dateAndTimeGPSSet() ) {
96 currDateTimeGPS = BNC_CORE->dateAndTimeGPS();
97 }
98 else {
99 currDateTimeGPS = QDateTime::currentDateTime().toUTC();
100 QDate hlp = currDateTimeGPS.date();
101 currDateTimeGPS = currDateTimeGPS.addSecs(gnumleap(hlp.year(),
102 hlp.month(), hlp.day()));
103 }
104
105 QDate currDateGPS = currDateTimeGPS.date();
106 QTime currTimeGPS = currDateTimeGPS.time();
107
108 week = int( (double(currDateGPS.toJulianDay()) - 2444244.5) / 7 );
109
110 sec = (currDateGPS.dayOfWeek() % 7) * 24.0 * 3600.0 +
111 currTimeGPS.hour() * 3600.0 +
112 currTimeGPS.minute() * 60.0 +
113 currTimeGPS.second() +
114 currTimeGPS.msec() / 1000.0;
115}
116
117//
118////////////////////////////////////////////////////////////////////////////
119QDateTime currentDateAndTimeGPS() {
120 if ( BNC_CORE->dateAndTimeGPSSet() ) {
121 return BNC_CORE->dateAndTimeGPS();
122 }
123 else {
124 int GPSWeek;
125 double GPSWeeks;
126 currentGPSWeeks(GPSWeek, GPSWeeks);
127 return dateAndTimeFromGPSweek(GPSWeek, GPSWeeks);
128 }
129}
130
131//
132////////////////////////////////////////////////////////////////////////////
133QByteArray ggaString(const QByteArray& latitude,
134 const QByteArray& longitude,
135 const QByteArray& height) {
136
137 double lat = strtod(latitude,NULL);
138 double lon = strtod(longitude,NULL);
139 double hei = strtod(height,NULL);
140
141 const char* flagN="N";
142 const char* flagE="E";
143 if (lon >180.) {lon=(lon-360.)*(-1.); flagE="W";}
144 if ((lon < 0.) && (lon >= -180.)) {lon=lon*(-1.); flagE="W";}
145 if (lon < -180.) {lon=(lon+360.); flagE="E";}
146 if (lat < 0.) {lat=lat*(-1.); flagN="S";}
147 QTime ttime(QDateTime::currentDateTime().toUTC().time());
148 int lat_deg = (int)lat;
149 double lat_min=(lat-lat_deg)*60.;
150 int lon_deg = (int)lon;
151 double lon_min=(lon-lon_deg)*60.;
152 int hh = 0 , mm = 0;
153 double ss = 0.0;
154 hh=ttime.hour();
155 mm=ttime.minute();
156 ss=(double)ttime.second()+0.001*ttime.msec();
157 QString gga;
158 gga += "GPGGA,";
159 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'));
160 gga += QString("%1%2,").arg((int)lat_deg,2, 10, QLatin1Char('0')).arg(lat_min, 7, 'f', 4, QLatin1Char('0'));
161 gga += flagN;
162 gga += QString(",%1%2,").arg((int)lon_deg,3, 10, QLatin1Char('0')).arg(lon_min, 7, 'f', 4, QLatin1Char('0'));
163 gga += flagE + QString(",1,05,1.00");
164 gga += QString(",%1,").arg(hei, 2, 'f', 1);
165 gga += QString("M,10.000,M,,");
166 int xori;
167 char XOR = 0;
168 char *Buff =gga.toAscii().data();
169 int iLen = strlen(Buff);
170 for (xori = 0; xori < iLen; xori++) {
171 XOR ^= (char)Buff[xori];
172 }
173 gga = "$" + gga + QString("*%1").arg(XOR, 2, 16, QLatin1Char('0'));
174
175 return gga.toAscii();
176}
177
178//
179////////////////////////////////////////////////////////////////////////////
180void RSW_to_XYZ(const ColumnVector& rr, const ColumnVector& vv,
181 const ColumnVector& rsw, ColumnVector& xyz) {
182
183 ColumnVector along = vv / vv.norm_Frobenius();
184 ColumnVector cross = crossproduct(rr, vv); cross /= cross.norm_Frobenius();
185 ColumnVector radial = crossproduct(along, cross);
186
187 Matrix RR(3,3);
188 RR.Column(1) = radial;
189 RR.Column(2) = along;
190 RR.Column(3) = cross;
191
192 xyz = RR * rsw;
193}
194
195// Transformation xyz --> radial, along track, out-of-plane
196////////////////////////////////////////////////////////////////////////////
197void XYZ_to_RSW(const ColumnVector& rr, const ColumnVector& vv,
198 const ColumnVector& xyz, ColumnVector& rsw) {
199
200 ColumnVector along = vv / vv.norm_Frobenius();
201 ColumnVector cross = crossproduct(rr, vv); cross /= cross.norm_Frobenius();
202 ColumnVector radial = crossproduct(along, cross);
203
204 rsw.ReSize(3);
205 rsw(1) = DotProduct(xyz, radial);
206 rsw(2) = DotProduct(xyz, along);
207 rsw(3) = DotProduct(xyz, cross);
208}
209
210// Rectangular Coordinates -> Ellipsoidal Coordinates
211////////////////////////////////////////////////////////////////////////////
212t_irc xyz2ell(const double* XYZ, double* Ell) {
213
214 const double bell = t_CST::aell*(1.0-1.0/t_CST::fInv) ;
215 const double e2 = (t_CST::aell*t_CST::aell-bell*bell)/(t_CST::aell*t_CST::aell) ;
216 const double e2c = (t_CST::aell*t_CST::aell-bell*bell)/(bell*bell) ;
217
218 double nn, ss, zps, hOld, phiOld, theta, sin3, cos3;
219
220 ss = sqrt(XYZ[0]*XYZ[0]+XYZ[1]*XYZ[1]) ;
221 zps = XYZ[2]/ss ;
222 theta = atan( (XYZ[2]*t_CST::aell) / (ss*bell) );
223 sin3 = sin(theta) * sin(theta) * sin(theta);
224 cos3 = cos(theta) * cos(theta) * cos(theta);
225
226 // Closed formula
227 Ell[0] = atan( (XYZ[2] + e2c * bell * sin3) / (ss - e2 * t_CST::aell * cos3) );
228 Ell[1] = atan2(XYZ[1],XYZ[0]) ;
229 nn = t_CST::aell/sqrt(1.0-e2*sin(Ell[0])*sin(Ell[0])) ;
230 Ell[2] = ss / cos(Ell[0]) - nn;
231
232 const int MAXITER = 100;
233 for (int ii = 1; ii <= MAXITER; ii++) {
234 nn = t_CST::aell/sqrt(1.0-e2*sin(Ell[0])*sin(Ell[0])) ;
235 hOld = Ell[2] ;
236 phiOld = Ell[0] ;
237 Ell[2] = ss/cos(Ell[0])-nn ;
238 Ell[0] = atan(zps/(1.0-e2*nn/(nn+Ell[2]))) ;
239 if ( fabs(phiOld-Ell[0]) <= 1.0e-11 && fabs(hOld-Ell[2]) <= 1.0e-5 ) {
240 return success;
241 }
242 }
243
244 return failure;
245}
246
247// Rectangular Coordinates -> North, East, Up Components
248////////////////////////////////////////////////////////////////////////////
249void xyz2neu(const double* Ell, const double* xyz, double* neu) {
250
251 double sinPhi = sin(Ell[0]);
252 double cosPhi = cos(Ell[0]);
253 double sinLam = sin(Ell[1]);
254 double cosLam = cos(Ell[1]);
255
256 neu[0] = - sinPhi*cosLam * xyz[0]
257 - sinPhi*sinLam * xyz[1]
258 + cosPhi * xyz[2];
259
260 neu[1] = - sinLam * xyz[0]
261 + cosLam * xyz[1];
262
263 neu[2] = + cosPhi*cosLam * xyz[0]
264 + cosPhi*sinLam * xyz[1]
265 + sinPhi * xyz[2];
266}
267
268// North, East, Up Components -> Rectangular Coordinates
269////////////////////////////////////////////////////////////////////////////
270void neu2xyz(const double* Ell, const double* neu, double* xyz) {
271
272 double sinPhi = sin(Ell[0]);
273 double cosPhi = cos(Ell[0]);
274 double sinLam = sin(Ell[1]);
275 double cosLam = cos(Ell[1]);
276
277 xyz[0] = - sinPhi*cosLam * neu[0]
278 - sinLam * neu[1]
279 + cosPhi*cosLam * neu[2];
280
281 xyz[1] = - sinPhi*sinLam * neu[0]
282 + cosLam * neu[1]
283 + cosPhi*sinLam * neu[2];
284
285 xyz[2] = + cosPhi * neu[0]
286 + sinPhi * neu[2];
287}
288
289//
290////////////////////////////////////////////////////////////////////////////
291double Frac (double x) {
292 return x-floor(x);
293}
294
295//
296////////////////////////////////////////////////////////////////////////////
297double Modulo (double x, double y) {
298 return y*Frac(x/y);
299}
300
301// Round to nearest integer
302////////////////////////////////////////////////////////////////////////////
303double nint(double val) {
304 return ((val < 0.0) ? -floor(fabs(val)+0.5) : floor(val+0.5));
305}
306
307// Jacobian XYZ --> NEU
308////////////////////////////////////////////////////////////////////////////
309void jacobiXYZ_NEU(const double* Ell, Matrix& jacobi) {
310
311 Tracer tracer("jacobiXYZ_NEU");
312
313 double sinPhi = sin(Ell[0]);
314 double cosPhi = cos(Ell[0]);
315 double sinLam = sin(Ell[1]);
316 double cosLam = cos(Ell[1]);
317
318 jacobi(1,1) = - sinPhi * cosLam;
319 jacobi(1,2) = - sinPhi * sinLam;
320 jacobi(1,3) = cosPhi;
321
322 jacobi(2,1) = - sinLam;
323 jacobi(2,2) = cosLam;
324 jacobi(2,3) = 0.0;
325
326 jacobi(3,1) = cosPhi * cosLam;
327 jacobi(3,2) = cosPhi * sinLam;
328 jacobi(3,3) = sinPhi;
329}
330
331// Jacobian Ell --> XYZ
332////////////////////////////////////////////////////////////////////////////
333void jacobiEll_XYZ(const double* Ell, Matrix& jacobi) {
334
335 Tracer tracer("jacobiEll_XYZ");
336
337 double sinPhi = sin(Ell[0]);
338 double cosPhi = cos(Ell[0]);
339 double sinLam = sin(Ell[1]);
340 double cosLam = cos(Ell[1]);
341 double hh = Ell[2];
342
343 double bell = t_CST::aell*(1.0-1.0/t_CST::fInv);
344 double e2 = (t_CST::aell*t_CST::aell-bell*bell)/(t_CST::aell*t_CST::aell) ;
345 double nn = t_CST::aell/sqrt(1.0-e2*sinPhi*sinPhi) ;
346
347 jacobi(1,1) = -(nn+hh) * sinPhi * cosLam;
348 jacobi(1,2) = -(nn+hh) * cosPhi * sinLam;
349 jacobi(1,3) = cosPhi * cosLam;
350
351 jacobi(2,1) = -(nn+hh) * sinPhi * sinLam;
352 jacobi(2,2) = (nn+hh) * cosPhi * cosLam;
353 jacobi(2,3) = cosPhi * sinLam;
354
355 jacobi(3,1) = (nn*(1.0-e2)+hh) * cosPhi;
356 jacobi(3,2) = 0.0;
357 jacobi(3,3) = sinPhi;
358}
359
360// Covariance Matrix in NEU
361////////////////////////////////////////////////////////////////////////////
362void covariXYZ_NEU(const SymmetricMatrix& QQxyz, const double* Ell,
363 SymmetricMatrix& Qneu) {
364
365 Tracer tracer("covariXYZ_NEU");
366
367 Matrix CC(3,3);
368 jacobiXYZ_NEU(Ell, CC);
369 Qneu << CC * QQxyz * CC.t();
370}
371
372// Covariance Matrix in XYZ
373////////////////////////////////////////////////////////////////////////////
374void covariNEU_XYZ(const SymmetricMatrix& QQneu, const double* Ell,
375 SymmetricMatrix& Qxyz) {
376
377 Tracer tracer("covariNEU_XYZ");
378
379 Matrix CC(3,3);
380 jacobiXYZ_NEU(Ell, CC);
381 Qxyz << CC.t() * QQneu * CC;
382}
383
384// Fourth order Runge-Kutta numerical integrator for ODEs
385////////////////////////////////////////////////////////////////////////////
386ColumnVector rungeKutta4(
387 double xi, // the initial x-value
388 const ColumnVector& yi, // vector of the initial y-values
389 double dx, // the step size for the integration
390 double* acc, // aditional acceleration
391 ColumnVector (*der)(double x, const ColumnVector& y, double* acc)
392 // A pointer to a function that computes the
393 // derivative of a function at a point (x,y)
394 ) {
395
396 ColumnVector k1 = der(xi , yi , acc) * dx;
397 ColumnVector k2 = der(xi+dx/2.0, yi+k1/2.0, acc) * dx;
398 ColumnVector k3 = der(xi+dx/2.0, yi+k2/2.0, acc) * dx;
399 ColumnVector k4 = der(xi+dx , yi+k3 , acc) * dx;
400
401 ColumnVector yf = yi + k1/6.0 + k2/3.0 + k3/3.0 + k4/6.0;
402
403 return yf;
404}
405
406//
407////////////////////////////////////////////////////////////////////////////
408double djul(int jj, int mm, double tt) {
409 int ii, kk;
410 double djul ;
411 if( mm <= 2 ) {
412 jj = jj - 1;
413 mm = mm + 12;
414 }
415 ii = jj/100;
416 kk = 2 - ii + ii/4;
417 djul = (365.25*jj - fmod( 365.25*jj, 1.0 )) - 679006.0;
418 djul = djul + floor( 30.6001*(mm + 1) ) + tt + kk;
419 return djul;
420}
421
422//
423////////////////////////////////////////////////////////////////////////////
424void jdgp(double tjul, double & second, int & nweek) {
425 double deltat;
426 deltat = tjul - 44244.0 ;
427 // current gps week
428 nweek = (int) floor(deltat/7.0);
429 // seconds past midnight of last weekend
430 second = (deltat - (nweek)*7.0)*86400.0;
431}
432
433//
434////////////////////////////////////////////////////////////////////////////
435void GPSweekFromDateAndTime(const QDateTime& dateTime,
436 int& GPSWeek, double& GPSWeeks) {
437
438 static const QDateTime zeroEpoch(QDate(1980, 1, 6),QTime(),Qt::UTC);
439
440 GPSWeek = zeroEpoch.daysTo(dateTime) / 7;
441
442 int weekDay = dateTime.date().dayOfWeek() + 1; // Qt: Monday = 1
443 if (weekDay > 7) weekDay = 1;
444
445 GPSWeeks = (weekDay - 1) * 86400.0
446 - dateTime.time().msecsTo(QTime()) / 1e3;
447}
448
449//
450////////////////////////////////////////////////////////////////////////////
451void GPSweekFromYMDhms(int year, int month, int day, int hour, int min,
452 double sec, int& GPSWeek, double& GPSWeeks) {
453
454 double mjd = djul(year, month, day);
455
456 jdgp(mjd, GPSWeeks, GPSWeek);
457 GPSWeeks += hour * 3600.0 + min * 60.0 + sec;
458}
459
460//
461////////////////////////////////////////////////////////////////////////////
462void mjdFromDateAndTime(const QDateTime& dateTime, int& mjd, double& dayfrac) {
463
464 static const QDate zeroDate(1858, 11, 17);
465
466 mjd = zeroDate.daysTo(dateTime.date());
467
468 dayfrac = (dateTime.time().hour() +
469 (dateTime.time().minute() +
470 (dateTime.time().second() +
471 dateTime.time().msec() / 1000.0) / 60.0) / 60.0) / 24.0;
472}
473
474//
475////////////////////////////////////////////////////////////////////////////
476bool findInVector(const vector<QString>& vv, const QString& str) {
477 std::vector<QString>::const_iterator it;
478 for (it = vv.begin(); it != vv.end(); ++it) {
479 if ( (*it) == str) {
480 return true;
481 }
482 }
483 return false;
484}
485
486//
487////////////////////////////////////////////////////////////////////////////
488int readInt(const QString& str, int pos, int len, int& value) {
489 bool ok;
490 value = str.mid(pos, len).toInt(&ok);
491 return ok ? 0 : 1;
492}
493
494//
495////////////////////////////////////////////////////////////////////////////
496int readDbl(const QString& str, int pos, int len, double& value) {
497 QString hlp = str.mid(pos, len);
498 for (int ii = 0; ii < hlp.length(); ii++) {
499 if (hlp[ii]=='D' || hlp[ii]=='d' || hlp[ii] == 'E') {
500 hlp[ii]='e';
501 }
502 }
503 bool ok;
504 value = hlp.toDouble(&ok);
505 return ok ? 0 : 1;
506}
507
508// Topocentrical Distance and Elevation
509////////////////////////////////////////////////////////////////////////////
510void topos(double xRec, double yRec, double zRec,
511 double xSat, double ySat, double zSat,
512 double& rho, double& eleSat, double& azSat) {
513
514 double dx[3];
515 dx[0] = xSat-xRec;
516 dx[1] = ySat-yRec;
517 dx[2] = zSat-zRec;
518
519 rho = sqrt( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
520
521 double xyzRec[3];
522 xyzRec[0] = xRec;
523 xyzRec[1] = yRec;
524 xyzRec[2] = zRec;
525
526 double Ell[3];
527 double neu[3];
528 xyz2ell(xyzRec, Ell);
529 xyz2neu(Ell, dx, neu);
530
531 eleSat = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / rho );
532 if (neu[2] < 0) {
533 eleSat *= -1.0;
534 }
535
536 azSat = atan2(neu[1], neu[0]);
537}
538
539// Degrees -> degrees, minutes, seconds
540////////////////////////////////////////////////////////////////////////////
541void deg2DMS(double decDeg, int& deg, int& min, double& sec) {
542 int sgn = (decDeg < 0.0 ? -1 : 1);
543 deg = sgn * static_cast<int>(decDeg);
544 min = static_cast<int>((decDeg - deg)*60);
545 sec = (decDeg - deg - min/60.0) * 3600.0;
546}
547
548//
549////////////////////////////////////////////////////////////////////////////
550QString fortranFormat(double value, int width, int prec) {
551 int expo = value == 0.0 ? 0 : log10(fabs(value));
552 double mant = value == 0.0 ? 0 : value / pow(10, expo);
553 if (fabs(mant) >= 1.0) {
554 mant /= 10.0;
555 expo += 1;
556 }
557 if (expo >= 0) {
558 return QString("%1e+%2").arg(mant, width-4, 'f', prec).arg(expo, 2, 10, QChar('0'));
559 }
560 else {
561 return QString("%1e-%2").arg(mant, width-4, 'f', prec).arg(-expo, 2, 10, QChar('0'));
562 }
563}
564
565//
566//////////////////////////////////////////////////////////////////////////////
567void kalman(const Matrix& AA, const ColumnVector& ll, const DiagonalMatrix& PP,
568 SymmetricMatrix& QQ, ColumnVector& dx) {
569
570 Tracer tracer("kalman");
571
572 int nPar = AA.Ncols();
573 int nObs = AA.Nrows();
574 UpperTriangularMatrix SS = Cholesky(QQ).t();
575
576 Matrix SA = SS*AA.t();
577 Matrix SRF(nObs+nPar, nObs+nPar); SRF = 0;
578 for (int ii = 1; ii <= nObs; ++ii) {
579 SRF(ii,ii) = 1.0 / sqrt(PP(ii,ii));
580 }
581
582 SRF.SubMatrix (nObs+1, nObs+nPar, 1, nObs) = SA;
583 SRF.SymSubMatrix(nObs+1, nObs+nPar) = SS;
584
585 UpperTriangularMatrix UU;
586 QRZ(SRF, UU);
587
588 SS = UU.SymSubMatrix(nObs+1, nObs+nPar);
589 UpperTriangularMatrix SH_rt = UU.SymSubMatrix(1, nObs);
590 Matrix YY = UU.SubMatrix(1, nObs, nObs+1, nObs+nPar);
591
592 UpperTriangularMatrix SHi = SH_rt.i();
593
594 Matrix KT = SHi * YY;
595 SymmetricMatrix Hi; Hi << SHi * SHi.t();
596
597 dx = KT.t() * ll;
598 QQ << (SS.t() * SS);
599}
600
Note: See TracBrowser for help on using the repository browser.