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

Last change on this file since 6800 was 6800, checked in by stuerze, 7 years ago

conversion from metric accuracy values into indices added for Galileo

File size: 21.7 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// Strip White Space
74////////////////////////////////////////////////////////////////////////////
75void stripWhiteSpace(string& str) {
76 if (!str.empty()) {
77 string::size_type beg = str.find_first_not_of(" \t\f\n\r\v");
78 string::size_type end = str.find_last_not_of(" \t\f\n\r\v");
79 if (beg > str.max_size())
80 str.erase();
81 else
82 str = str.substr(beg, end-beg+1);
83 }
84}
85
86//
87////////////////////////////////////////////////////////////////////////////
88QDateTime dateAndTimeFromGPSweek(int GPSWeek, double GPSWeeks) {
89
90 static const QDate zeroEpoch(1980, 1, 6);
91
92 QDate date(zeroEpoch);
93 QTime time(0,0,0,0);
94
95 int weekDays = int(GPSWeeks) / 86400;
96 date = date.addDays( GPSWeek * 7 + weekDays );
97 time = time.addMSecs( int( (GPSWeeks - 86400 * weekDays) * 1e3 ) );
98
99 return QDateTime(date,time);
100}
101
102//
103////////////////////////////////////////////////////////////////////////////
104void currentGPSWeeks(int& week, double& sec) {
105
106 QDateTime currDateTimeGPS;
107
108 if ( BNC_CORE->dateAndTimeGPSSet() ) {
109 currDateTimeGPS = BNC_CORE->dateAndTimeGPS();
110 }
111 else {
112 currDateTimeGPS = QDateTime::currentDateTime().toUTC();
113 QDate hlp = currDateTimeGPS.date();
114 currDateTimeGPS = currDateTimeGPS.addSecs(gnumleap(hlp.year(),
115 hlp.month(), hlp.day()));
116 }
117
118 QDate currDateGPS = currDateTimeGPS.date();
119 QTime currTimeGPS = currDateTimeGPS.time();
120
121 week = int( (double(currDateGPS.toJulianDay()) - 2444244.5) / 7 );
122
123 sec = (currDateGPS.dayOfWeek() % 7) * 24.0 * 3600.0 +
124 currTimeGPS.hour() * 3600.0 +
125 currTimeGPS.minute() * 60.0 +
126 currTimeGPS.second() +
127 currTimeGPS.msec() / 1000.0;
128}
129
130//
131////////////////////////////////////////////////////////////////////////////
132QDateTime currentDateAndTimeGPS() {
133 if ( BNC_CORE->dateAndTimeGPSSet() ) {
134 return BNC_CORE->dateAndTimeGPS();
135 }
136 else {
137 int GPSWeek;
138 double GPSWeeks;
139 currentGPSWeeks(GPSWeek, GPSWeeks);
140 return dateAndTimeFromGPSweek(GPSWeek, GPSWeeks);
141 }
142}
143
144//
145////////////////////////////////////////////////////////////////////////////
146QByteArray ggaString(const QByteArray& latitude,
147 const QByteArray& longitude,
148 const QByteArray& height,
149 const QString& ggaType) {
150
151 double lat = strtod(latitude,NULL);
152 double lon = strtod(longitude,NULL);
153 double hei = strtod(height,NULL);
154 QString sentences = "GPGGA,";
155 if (ggaType.contains("GNGGA")) {
156 sentences = "GNGGA,";
157 }
158
159 const char* flagN="N";
160 const char* flagE="E";
161 if (lon >180.) {lon=(lon-360.)*(-1.); flagE="W";}
162 if ((lon < 0.) && (lon >= -180.)) {lon=lon*(-1.); flagE="W";}
163 if (lon < -180.) {lon=(lon+360.); flagE="E";}
164 if (lat < 0.) {lat=lat*(-1.); flagN="S";}
165 QTime ttime(QDateTime::currentDateTime().toUTC().time());
166 int lat_deg = (int)lat;
167 double lat_min=(lat-lat_deg)*60.;
168 int lon_deg = (int)lon;
169 double lon_min=(lon-lon_deg)*60.;
170 int hh = 0 , mm = 0;
171 double ss = 0.0;
172 hh=ttime.hour();
173 mm=ttime.minute();
174 ss=(double)ttime.second()+0.001*ttime.msec();
175 QString gga;
176 gga += sentences;
177 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'));
178 gga += QString("%1%2,").arg((int)lat_deg,2, 10, QLatin1Char('0')).arg(lat_min, 7, 'f', 4, QLatin1Char('0'));
179 gga += flagN;
180 gga += QString(",%1%2,").arg((int)lon_deg,3, 10, QLatin1Char('0')).arg(lon_min, 7, 'f', 4, QLatin1Char('0'));
181 gga += flagE + QString(",1,05,1.00");
182 gga += QString(",%1,").arg(hei, 2, 'f', 1);
183 gga += QString("M,10.000,M,,");
184 int xori;
185 char XOR = 0;
186 char *Buff =gga.toAscii().data();
187 int iLen = strlen(Buff);
188 for (xori = 0; xori < iLen; xori++) {
189 XOR ^= (char)Buff[xori];
190 }
191 gga = "$" + gga + QString("*%1").arg(XOR, 2, 16, QLatin1Char('0'));
192
193 return gga.toAscii();
194}
195
196//
197////////////////////////////////////////////////////////////////////////////
198void RSW_to_XYZ(const ColumnVector& rr, const ColumnVector& vv,
199 const ColumnVector& rsw, ColumnVector& xyz) {
200
201 ColumnVector along = vv / vv.norm_Frobenius();
202 ColumnVector cross = crossproduct(rr, vv); cross /= cross.norm_Frobenius();
203 ColumnVector radial = crossproduct(along, cross);
204
205 Matrix RR(3,3);
206 RR.Column(1) = radial;
207 RR.Column(2) = along;
208 RR.Column(3) = cross;
209
210 xyz = RR * rsw;
211}
212
213// Transformation xyz --> radial, along track, out-of-plane
214////////////////////////////////////////////////////////////////////////////
215void XYZ_to_RSW(const ColumnVector& rr, const ColumnVector& vv,
216 const ColumnVector& xyz, ColumnVector& rsw) {
217
218 ColumnVector along = vv / vv.norm_Frobenius();
219 ColumnVector cross = crossproduct(rr, vv); cross /= cross.norm_Frobenius();
220 ColumnVector radial = crossproduct(along, cross);
221
222 rsw.ReSize(3);
223 rsw(1) = DotProduct(xyz, radial);
224 rsw(2) = DotProduct(xyz, along);
225 rsw(3) = DotProduct(xyz, cross);
226}
227
228// Rectangular Coordinates -> Ellipsoidal Coordinates
229////////////////////////////////////////////////////////////////////////////
230t_irc xyz2ell(const double* XYZ, double* Ell) {
231
232 const double bell = t_CST::aell*(1.0-1.0/t_CST::fInv) ;
233 const double e2 = (t_CST::aell*t_CST::aell-bell*bell)/(t_CST::aell*t_CST::aell) ;
234 const double e2c = (t_CST::aell*t_CST::aell-bell*bell)/(bell*bell) ;
235
236 double nn, ss, zps, hOld, phiOld, theta, sin3, cos3;
237
238 ss = sqrt(XYZ[0]*XYZ[0]+XYZ[1]*XYZ[1]) ;
239 zps = XYZ[2]/ss ;
240 theta = atan( (XYZ[2]*t_CST::aell) / (ss*bell) );
241 sin3 = sin(theta) * sin(theta) * sin(theta);
242 cos3 = cos(theta) * cos(theta) * cos(theta);
243
244 // Closed formula
245 Ell[0] = atan( (XYZ[2] + e2c * bell * sin3) / (ss - e2 * t_CST::aell * cos3) );
246 Ell[1] = atan2(XYZ[1],XYZ[0]) ;
247 nn = t_CST::aell/sqrt(1.0-e2*sin(Ell[0])*sin(Ell[0])) ;
248 Ell[2] = ss / cos(Ell[0]) - nn;
249
250 const int MAXITER = 100;
251 for (int ii = 1; ii <= MAXITER; ii++) {
252 nn = t_CST::aell/sqrt(1.0-e2*sin(Ell[0])*sin(Ell[0])) ;
253 hOld = Ell[2] ;
254 phiOld = Ell[0] ;
255 Ell[2] = ss/cos(Ell[0])-nn ;
256 Ell[0] = atan(zps/(1.0-e2*nn/(nn+Ell[2]))) ;
257 if ( fabs(phiOld-Ell[0]) <= 1.0e-11 && fabs(hOld-Ell[2]) <= 1.0e-5 ) {
258 return success;
259 }
260 }
261
262 return failure;
263}
264
265// Rectangular Coordinates -> North, East, Up Components
266////////////////////////////////////////////////////////////////////////////
267void xyz2neu(const double* Ell, const double* xyz, double* neu) {
268
269 double sinPhi = sin(Ell[0]);
270 double cosPhi = cos(Ell[0]);
271 double sinLam = sin(Ell[1]);
272 double cosLam = cos(Ell[1]);
273
274 neu[0] = - sinPhi*cosLam * xyz[0]
275 - sinPhi*sinLam * xyz[1]
276 + cosPhi * xyz[2];
277
278 neu[1] = - sinLam * xyz[0]
279 + cosLam * xyz[1];
280
281 neu[2] = + cosPhi*cosLam * xyz[0]
282 + cosPhi*sinLam * xyz[1]
283 + sinPhi * xyz[2];
284}
285
286// North, East, Up Components -> Rectangular Coordinates
287////////////////////////////////////////////////////////////////////////////
288void neu2xyz(const double* Ell, const double* neu, double* xyz) {
289
290 double sinPhi = sin(Ell[0]);
291 double cosPhi = cos(Ell[0]);
292 double sinLam = sin(Ell[1]);
293 double cosLam = cos(Ell[1]);
294
295 xyz[0] = - sinPhi*cosLam * neu[0]
296 - sinLam * neu[1]
297 + cosPhi*cosLam * neu[2];
298
299 xyz[1] = - sinPhi*sinLam * neu[0]
300 + cosLam * neu[1]
301 + cosPhi*sinLam * neu[2];
302
303 xyz[2] = + cosPhi * neu[0]
304 + sinPhi * neu[2];
305}
306
307//
308////////////////////////////////////////////////////////////////////////////
309double Frac (double x) {
310 return x-floor(x);
311}
312
313//
314////////////////////////////////////////////////////////////////////////////
315double Modulo (double x, double y) {
316 return y*Frac(x/y);
317}
318
319// Round to nearest integer
320////////////////////////////////////////////////////////////////////////////
321double nint(double val) {
322 return ((val < 0.0) ? -floor(fabs(val)+0.5) : floor(val+0.5));
323}
324
325// Jacobian XYZ --> NEU
326////////////////////////////////////////////////////////////////////////////
327void jacobiXYZ_NEU(const double* Ell, Matrix& jacobi) {
328
329 Tracer tracer("jacobiXYZ_NEU");
330
331 double sinPhi = sin(Ell[0]);
332 double cosPhi = cos(Ell[0]);
333 double sinLam = sin(Ell[1]);
334 double cosLam = cos(Ell[1]);
335
336 jacobi(1,1) = - sinPhi * cosLam;
337 jacobi(1,2) = - sinPhi * sinLam;
338 jacobi(1,3) = cosPhi;
339
340 jacobi(2,1) = - sinLam;
341 jacobi(2,2) = cosLam;
342 jacobi(2,3) = 0.0;
343
344 jacobi(3,1) = cosPhi * cosLam;
345 jacobi(3,2) = cosPhi * sinLam;
346 jacobi(3,3) = sinPhi;
347}
348
349// Jacobian Ell --> XYZ
350////////////////////////////////////////////////////////////////////////////
351void jacobiEll_XYZ(const double* Ell, Matrix& jacobi) {
352
353 Tracer tracer("jacobiEll_XYZ");
354
355 double sinPhi = sin(Ell[0]);
356 double cosPhi = cos(Ell[0]);
357 double sinLam = sin(Ell[1]);
358 double cosLam = cos(Ell[1]);
359 double hh = Ell[2];
360
361 double bell = t_CST::aell*(1.0-1.0/t_CST::fInv);
362 double e2 = (t_CST::aell*t_CST::aell-bell*bell)/(t_CST::aell*t_CST::aell) ;
363 double nn = t_CST::aell/sqrt(1.0-e2*sinPhi*sinPhi) ;
364
365 jacobi(1,1) = -(nn+hh) * sinPhi * cosLam;
366 jacobi(1,2) = -(nn+hh) * cosPhi * sinLam;
367 jacobi(1,3) = cosPhi * cosLam;
368
369 jacobi(2,1) = -(nn+hh) * sinPhi * sinLam;
370 jacobi(2,2) = (nn+hh) * cosPhi * cosLam;
371 jacobi(2,3) = cosPhi * sinLam;
372
373 jacobi(3,1) = (nn*(1.0-e2)+hh) * cosPhi;
374 jacobi(3,2) = 0.0;
375 jacobi(3,3) = sinPhi;
376}
377
378// Covariance Matrix in NEU
379////////////////////////////////////////////////////////////////////////////
380void covariXYZ_NEU(const SymmetricMatrix& QQxyz, const double* Ell,
381 SymmetricMatrix& Qneu) {
382
383 Tracer tracer("covariXYZ_NEU");
384
385 Matrix CC(3,3);
386 jacobiXYZ_NEU(Ell, CC);
387 Qneu << CC * QQxyz * CC.t();
388}
389
390// Covariance Matrix in XYZ
391////////////////////////////////////////////////////////////////////////////
392void covariNEU_XYZ(const SymmetricMatrix& QQneu, const double* Ell,
393 SymmetricMatrix& Qxyz) {
394
395 Tracer tracer("covariNEU_XYZ");
396
397 Matrix CC(3,3);
398 jacobiXYZ_NEU(Ell, CC);
399 Qxyz << CC.t() * QQneu * CC;
400}
401
402// Fourth order Runge-Kutta numerical integrator for ODEs
403////////////////////////////////////////////////////////////////////////////
404ColumnVector rungeKutta4(
405 double xi, // the initial x-value
406 const ColumnVector& yi, // vector of the initial y-values
407 double dx, // the step size for the integration
408 double* acc, // aditional acceleration
409 ColumnVector (*der)(double x, const ColumnVector& y, double* acc)
410 // A pointer to a function that computes the
411 // derivative of a function at a point (x,y)
412 ) {
413
414 ColumnVector k1 = der(xi , yi , acc) * dx;
415 ColumnVector k2 = der(xi+dx/2.0, yi+k1/2.0, acc) * dx;
416 ColumnVector k3 = der(xi+dx/2.0, yi+k2/2.0, acc) * dx;
417 ColumnVector k4 = der(xi+dx , yi+k3 , acc) * dx;
418
419 ColumnVector yf = yi + k1/6.0 + k2/3.0 + k3/3.0 + k4/6.0;
420
421 return yf;
422}
423//
424////////////////////////////////////////////////////////////////////////////
425double djul(long jj, long mm, double tt) {
426 long ii, kk;
427 double djul ;
428 if( mm <= 2 ) {
429 jj = jj - 1;
430 mm = mm + 12;
431 }
432 ii = jj/100;
433 kk = 2 - ii + ii/4;
434 djul = (365.25*jj - fmod( 365.25*jj, 1.0 )) - 679006.0;
435 djul = djul + floor( 30.6001*(mm + 1) ) + tt + kk;
436 return djul;
437}
438
439//
440////////////////////////////////////////////////////////////////////////////
441double gpjd(double second, int nweek) {
442 double deltat;
443 deltat = nweek*7.0 + second/86400.0 ;
444 return( 44244.0 + deltat) ;
445}
446
447//
448////////////////////////////////////////////////////////////////////////////
449void jdgp(double tjul, double & second, long & nweek) {
450 double deltat;
451 deltat = tjul - 44244.0 ;
452 nweek = (long) floor(deltat/7.0);
453 second = (deltat - (nweek)*7.0)*86400.0;
454}
455
456//
457////////////////////////////////////////////////////////////////////////////
458void jmt(double djul, long& jj, long& mm, double& dd) {
459 long ih, ih1, ih2 ;
460 double t1, t2, t3, t4;
461 t1 = 1.0 + djul - fmod( djul, 1.0 ) + 2400000.0;
462 t4 = fmod( djul, 1.0 );
463 ih = long( (t1 - 1867216.25)/36524.25 );
464 t2 = t1 + 1 + ih - ih/4;
465 t3 = t2 - 1720995.0;
466 ih1 = long( (t3 - 122.1)/365.25 );
467 t1 = 365.25*ih1 - fmod( 365.25*ih1, 1.0 );
468 ih2 = long( (t3 - t1)/30.6001 );
469 dd = t3 - t1 - (int)( 30.6001*ih2 ) + t4;
470 mm = ih2 - 1;
471 if ( ih2 > 13 ) mm = ih2 - 13;
472 jj = ih1;
473 if ( mm <= 2 ) jj = jj + 1;
474}
475
476//
477////////////////////////////////////////////////////////////////////////////
478void GPSweekFromDateAndTime(const QDateTime& dateTime,
479 int& GPSWeek, double& GPSWeeks) {
480
481 static const QDateTime zeroEpoch(QDate(1980, 1, 6),QTime(),Qt::UTC);
482
483 GPSWeek = zeroEpoch.daysTo(dateTime) / 7;
484
485 int weekDay = dateTime.date().dayOfWeek() + 1; // Qt: Monday = 1
486 if (weekDay > 7) weekDay = 1;
487
488 GPSWeeks = (weekDay - 1) * 86400.0
489 - dateTime.time().msecsTo(QTime()) / 1e3;
490}
491
492//
493////////////////////////////////////////////////////////////////////////////
494void GPSweekFromYMDhms(int year, int month, int day, int hour, int min,
495 double sec, int& GPSWeek, double& GPSWeeks) {
496
497 double mjd = djul(year, month, day);
498
499 long GPSWeek_long;
500 jdgp(mjd, GPSWeeks, GPSWeek_long);
501 GPSWeek = GPSWeek_long;
502 GPSWeeks += hour * 3600.0 + min * 60.0 + sec;
503}
504
505//
506////////////////////////////////////////////////////////////////////////////
507void mjdFromDateAndTime(const QDateTime& dateTime, int& mjd, double& dayfrac) {
508
509 static const QDate zeroDate(1858, 11, 17);
510
511 mjd = zeroDate.daysTo(dateTime.date());
512
513 dayfrac = (dateTime.time().hour() +
514 (dateTime.time().minute() +
515 (dateTime.time().second() +
516 dateTime.time().msec() / 1000.0) / 60.0) / 60.0) / 24.0;
517}
518
519//
520////////////////////////////////////////////////////////////////////////////
521bool findInVector(const vector<QString>& vv, const QString& str) {
522 std::vector<QString>::const_iterator it;
523 for (it = vv.begin(); it != vv.end(); ++it) {
524 if ( (*it) == str) {
525 return true;
526 }
527 }
528 return false;
529}
530
531//
532////////////////////////////////////////////////////////////////////////////
533int readInt(const QString& str, int pos, int len, int& value) {
534 bool ok;
535 value = str.mid(pos, len).toInt(&ok);
536 return ok ? 0 : 1;
537}
538
539//
540////////////////////////////////////////////////////////////////////////////
541int readDbl(const QString& str, int pos, int len, double& value) {
542 QString hlp = str.mid(pos, len);
543 for (int ii = 0; ii < hlp.length(); ii++) {
544 if (hlp[ii]=='D' || hlp[ii]=='d' || hlp[ii] == 'E') {
545 hlp[ii]='e';
546 }
547 }
548 bool ok;
549 value = hlp.toDouble(&ok);
550 return ok ? 0 : 1;
551}
552
553// Topocentrical Distance and Elevation
554////////////////////////////////////////////////////////////////////////////
555void topos(double xRec, double yRec, double zRec,
556 double xSat, double ySat, double zSat,
557 double& rho, double& eleSat, double& azSat) {
558
559 double dx[3];
560 dx[0] = xSat-xRec;
561 dx[1] = ySat-yRec;
562 dx[2] = zSat-zRec;
563
564 rho = sqrt( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
565
566 double xyzRec[3];
567 xyzRec[0] = xRec;
568 xyzRec[1] = yRec;
569 xyzRec[2] = zRec;
570
571 double Ell[3];
572 double neu[3];
573 xyz2ell(xyzRec, Ell);
574 xyz2neu(Ell, dx, neu);
575
576 eleSat = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / rho );
577 if (neu[2] < 0) {
578 eleSat *= -1.0;
579 }
580
581 azSat = atan2(neu[1], neu[0]);
582}
583
584// Degrees -> degrees, minutes, seconds
585////////////////////////////////////////////////////////////////////////////
586void deg2DMS(double decDeg, int& deg, int& min, double& sec) {
587 int sgn = (decDeg < 0.0 ? -1 : 1);
588 deg = sgn * static_cast<int>(decDeg);
589 min = static_cast<int>((decDeg - deg)*60);
590 sec = (decDeg - deg - min/60.0) * 3600.0;
591}
592
593//
594////////////////////////////////////////////////////////////////////////////
595QString fortranFormat(double value, int width, int prec) {
596 int expo = value == 0.0 ? 0 : int(log10(fabs(value)));
597 double mant = value == 0.0 ? 0 : value / pow(10, expo);
598 if (fabs(mant) >= 1.0) {
599 mant /= 10.0;
600 expo += 1;
601 }
602 if (expo >= 0) {
603 return QString("%1e+%2").arg(mant, width-4, 'f', prec).arg(expo, 2, 10, QChar('0'));
604 }
605 else {
606 return QString("%1e-%2").arg(mant, width-4, 'f', prec).arg(-expo, 2, 10, QChar('0'));
607 }
608}
609
610//
611//////////////////////////////////////////////////////////////////////////////
612void kalman(const Matrix& AA, const ColumnVector& ll, const DiagonalMatrix& PP,
613 SymmetricMatrix& QQ, ColumnVector& xx) {
614
615 Tracer tracer("kalman");
616
617 int nPar = AA.Ncols();
618 int nObs = AA.Nrows();
619 UpperTriangularMatrix SS = Cholesky(QQ).t();
620
621 Matrix SA = SS*AA.t();
622 Matrix SRF(nObs+nPar, nObs+nPar); SRF = 0;
623 for (int ii = 1; ii <= nObs; ++ii) {
624 SRF(ii,ii) = 1.0 / sqrt(PP(ii,ii));
625 }
626
627 SRF.SubMatrix (nObs+1, nObs+nPar, 1, nObs) = SA;
628 SRF.SymSubMatrix(nObs+1, nObs+nPar) = SS;
629
630 UpperTriangularMatrix UU;
631 QRZ(SRF, UU);
632
633 SS = UU.SymSubMatrix(nObs+1, nObs+nPar);
634 UpperTriangularMatrix SH_rt = UU.SymSubMatrix(1, nObs);
635 Matrix YY = UU.SubMatrix(1, nObs, nObs+1, nObs+nPar);
636
637 UpperTriangularMatrix SHi = SH_rt.i();
638
639 Matrix KT = SHi * YY;
640 SymmetricMatrix Hi; Hi << SHi * SHi.t();
641
642 xx += KT.t() * (ll - AA * xx);
643 QQ << (SS.t() * SS);
644}
645
646double accuracyFromIndex(int index, t_eph::e_type type) {
647
648 if (type == t_eph::GPS || type == t_eph::BDS || type == t_eph::SBAS
649 || type == t_eph::QZSS) {
650
651 if ((index >= 0) && (index <= 6)) {
652 if (index == 3) {
653 return ceil(10.0 * pow(2.0, (double(index) / 2.0) + 1.0)) / 10.0;
654 }
655 else {
656 return floor(10.0 * pow(2.0, (double(index) / 2.0) + 1.0)) / 10.0;
657 }
658 }
659 else if ((index > 6) && (index <= 15)) {
660 return (10.0 * pow(2.0, (double(index) - 2.0))) / 10.0;
661 }
662 else {
663 return 8192.0;
664 }
665 }
666
667 if (type == t_eph::Galileo) {
668
669 if ((index >= 0) && (index <= 49)) {
670 return (double(index) / 100.0);
671 }
672 else if ((index > 49) && (index <= 74)) {
673 return (50.0 + (double(index) - 50.0) * 2.0) / 100.0;
674 }
675 else if ((index > 74) && (index <= 99)) {
676 return 1.0 + (double(index) - 75.0) * 0.04;
677 }
678 else if ((index > 99) && (index <= 125)) {
679 return 2.0 + (double(index) - 100.0) * 0.16;
680 }
681 else {
682 return -1.0;
683 }
684 }
685
686 return double(index);
687}
688
689int indexFromAccuracy(double accuracy, t_eph::e_type type) {
690
691 if (type == t_eph::GPS || type == t_eph::BDS || type == t_eph::SBAS
692 || type == t_eph::QZSS) {
693
694 if (accuracy <= 2.40) {
695 return 0;
696 }
697 else if (accuracy <= 3.40) {
698 return 1;
699 }
700 else if (accuracy <= 4.85) {
701 return 2;
702 }
703 else if (accuracy <= 6.85) {
704 return 3;
705 }
706 else if (accuracy <= 9.65) {
707 return 4;
708 }
709 else if (accuracy <= 13.65) {
710 return 5;
711 }
712 else if (accuracy <= 24.00) {
713 return 6;
714 }
715 else if (accuracy <= 48.00) {
716 return 7;
717 }
718 else if (accuracy <= 96.00) {
719 return 8;
720 }
721 else if (accuracy <= 192.00) {
722 return 9;
723 }
724 else if (accuracy <= 384.00) {
725 return 10;
726 }
727 else if (accuracy <= 768.00) {
728 return 11;
729 }
730 else if (accuracy <= 1536.00) {
731 return 12;
732 }
733 else if (accuracy <= 3072.00) {
734 return 13;
735 }
736 else if (accuracy <= 6144.00) {
737 return 14;
738 }
739 else {
740 return 15;
741 }
742 }
743
744 if (type == t_eph::Galileo) {
745
746 if (accuracy <= 0.49) {
747 return int(ceil(accuracy * 100.0));
748 }
749 else if (accuracy <= 0.98) {
750 return int(50.0 + (((accuracy * 100.0) - 50) / 2.0));
751 }
752 else if (accuracy <= 2.0) {
753 return int(75.0 + ((accuracy - 1.0) / 0.04));
754 }
755 else if (accuracy <= 6.0) {
756 return int(100.0 + ((accuracy - 2.0) / 0.16));
757 }
758 else {
759 return 255;
760 }
761 }
762
763 return (type == t_eph::Galileo) ? 255 : 15;
764}
Note: See TracBrowser for help on using the repository browser.