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

Last change on this file since 7787 was 7787, checked in by stuerze, 4 years ago

small bug fixed in deg2DMS

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