[280] | 1 | // Part of BNC, a utility for retrieving decoding and
|
---|
[464] | 2 | // converting GNSS data streams from NTRIP broadcasters.
|
---|
[280] | 3 | //
|
---|
[464] | 4 | // Copyright (C) 2007
|
---|
[280] | 5 | // German Federal Agency for Cartography and Geodesy (BKG)
|
---|
| 6 | // http://www.bkg.bund.de
|
---|
[464] | 7 | // Czech Technical University Prague, Department of Geodesy
|
---|
[280] | 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.
|
---|
[83] | 24 |
|
---|
| 25 | /* -------------------------------------------------------------------------
|
---|
[93] | 26 | * BKG NTRIP Client
|
---|
[83] | 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 |
|
---|
[124] | 41 | #include <iostream>
|
---|
[218] | 42 | #include <ctime>
|
---|
[221] | 43 | #include <math.h>
|
---|
[124] | 44 |
|
---|
[83] | 45 | #include <QRegExp>
|
---|
| 46 | #include <QStringList>
|
---|
[271] | 47 | #include <QDateTime>
|
---|
[83] | 48 |
|
---|
[5807] | 49 | #include <newmatap.h>
|
---|
| 50 |
|
---|
[83] | 51 | #include "bncutils.h"
|
---|
[5070] | 52 | #include "bnccore.h"
|
---|
[83] | 53 |
|
---|
[124] | 54 | using namespace std;
|
---|
| 55 |
|
---|
[1381] | 56 | //
|
---|
| 57 | ////////////////////////////////////////////////////////////////////////////
|
---|
[83] | 58 | void 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 | }
|
---|
[124] | 72 |
|
---|
[1381] | 73 | //
|
---|
| 74 | ////////////////////////////////////////////////////////////////////////////
|
---|
[124] | 75 | QDateTime 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 | }
|
---|
[210] | 88 |
|
---|
[1381] | 89 | //
|
---|
| 90 | ////////////////////////////////////////////////////////////////////////////
|
---|
[218] | 91 | void currentGPSWeeks(int& week, double& sec) {
|
---|
[210] | 92 |
|
---|
[1942] | 93 | QDateTime currDateTimeGPS;
|
---|
[1155] | 94 |
|
---|
[5846] | 95 | if ( BNC_CORE->dateAndTimeGPSSet() ) {
|
---|
| 96 | currDateTimeGPS = BNC_CORE->dateAndTimeGPS();
|
---|
[1155] | 97 | }
|
---|
| 98 | else {
|
---|
[1942] | 99 | currDateTimeGPS = QDateTime::currentDateTime().toUTC();
|
---|
| 100 | QDate hlp = currDateTimeGPS.date();
|
---|
| 101 | currDateTimeGPS = currDateTimeGPS.addSecs(gnumleap(hlp.year(),
|
---|
| 102 | hlp.month(), hlp.day()));
|
---|
[1155] | 103 | }
|
---|
| 104 |
|
---|
[1942] | 105 | QDate currDateGPS = currDateTimeGPS.date();
|
---|
| 106 | QTime currTimeGPS = currDateTimeGPS.time();
|
---|
[210] | 107 |
|
---|
[1942] | 108 | week = int( (double(currDateGPS.toJulianDay()) - 2444244.5) / 7 );
|
---|
[1036] | 109 |
|
---|
[1942] | 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;
|
---|
[1036] | 115 | }
|
---|
[1154] | 116 |
|
---|
[1381] | 117 | //
|
---|
| 118 | ////////////////////////////////////////////////////////////////////////////
|
---|
[1154] | 119 | QDateTime currentDateAndTimeGPS() {
|
---|
[5846] | 120 | if ( BNC_CORE->dateAndTimeGPSSet() ) {
|
---|
| 121 | return BNC_CORE->dateAndTimeGPS();
|
---|
[2530] | 122 | }
|
---|
| 123 | else {
|
---|
| 124 | int GPSWeek;
|
---|
| 125 | double GPSWeeks;
|
---|
| 126 | currentGPSWeeks(GPSWeek, GPSWeeks);
|
---|
| 127 | return dateAndTimeFromGPSweek(GPSWeek, GPSWeeks);
|
---|
| 128 | }
|
---|
[1154] | 129 | }
|
---|
| 130 |
|
---|
[1381] | 131 | //
|
---|
| 132 | ////////////////////////////////////////////////////////////////////////////
|
---|
[1595] | 133 | QByteArray ggaString(const QByteArray& latitude,
|
---|
| 134 | const QByteArray& longitude,
|
---|
| 135 | const QByteArray& height) {
|
---|
[1381] | 136 |
|
---|
| 137 | double lat = strtod(latitude,NULL);
|
---|
| 138 | double lon = strtod(longitude,NULL);
|
---|
[1595] | 139 | double hei = strtod(height,NULL);
|
---|
[1381] | 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;
|
---|
[1506] | 158 | gga += "GPGGA,";
|
---|
[1381] | 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'));
|
---|
[1595] | 163 | gga += flagE + QString(",1,05,1.00");
|
---|
[1599] | 164 | gga += QString(",%1,").arg(hei, 2, 'f', 1);
|
---|
[1595] | 165 | gga += QString("M,10.000,M,,");
|
---|
[1381] | 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 | }
|
---|
[1506] | 173 | gga = "$" + gga + QString("*%1").arg(XOR, 2, 16, QLatin1Char('0'));
|
---|
[1381] | 174 |
|
---|
[1387] | 175 | return gga.toAscii();
|
---|
[1381] | 176 | }
|
---|
[2043] | 177 |
|
---|
| 178 | //
|
---|
| 179 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 180 | void 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 | }
|
---|
[2063] | 194 |
|
---|
[2988] | 195 | // Transformation xyz --> radial, along track, out-of-plane
|
---|
| 196 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 197 | void 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 |
|
---|
[2063] | 210 | // Rectangular Coordinates -> Ellipsoidal Coordinates
|
---|
| 211 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 212 | t_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 | }
|
---|
[2065] | 246 |
|
---|
| 247 | // Rectangular Coordinates -> North, East, Up Components
|
---|
| 248 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 249 | void 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 | }
|
---|
[2221] | 267 |
|
---|
[2582] | 268 | // North, East, Up Components -> Rectangular Coordinates
|
---|
| 269 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 270 | void 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 |
|
---|
[5807] | 289 | //
|
---|
| 290 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 291 | double Frac (double x) {
|
---|
| 292 | return x-floor(x);
|
---|
| 293 | }
|
---|
| 294 |
|
---|
| 295 | //
|
---|
| 296 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 297 | double Modulo (double x, double y) {
|
---|
| 298 | return y*Frac(x/y);
|
---|
| 299 | }
|
---|
| 300 |
|
---|
[5753] | 301 | // Round to nearest integer
|
---|
| 302 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 303 | double nint(double val) {
|
---|
| 304 | return ((val < 0.0) ? -floor(fabs(val)+0.5) : floor(val+0.5));
|
---|
| 305 | }
|
---|
| 306 |
|
---|
[5752] | 307 | // Jacobian XYZ --> NEU
|
---|
| 308 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 309 | void 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 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 333 | void 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 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 362 | void 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 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 374 | void 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 |
|
---|
[2221] | 384 | // Fourth order Runge-Kutta numerical integrator for ODEs
|
---|
| 385 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 386 | ColumnVector 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
|
---|
[2556] | 390 | double* acc, // aditional acceleration
|
---|
| 391 | ColumnVector (*der)(double x, const ColumnVector& y, double* acc)
|
---|
[2221] | 392 | // A pointer to a function that computes the
|
---|
| 393 | // derivative of a function at a point (x,y)
|
---|
| 394 | ) {
|
---|
| 395 |
|
---|
[2556] | 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;
|
---|
[2221] | 400 |
|
---|
| 401 | ColumnVector yf = yi + k1/6.0 + k2/3.0 + k3/3.0 + k4/6.0;
|
---|
| 402 |
|
---|
| 403 | return yf;
|
---|
| 404 | }
|
---|
[3044] | 405 | //
|
---|
| 406 | ////////////////////////////////////////////////////////////////////////////
|
---|
[5886] | 407 | double djul(long jj, long mm, double tt) {
|
---|
| 408 | long ii, kk;
|
---|
[3171] | 409 | double djul ;
|
---|
| 410 | if( mm <= 2 ) {
|
---|
| 411 | jj = jj - 1;
|
---|
| 412 | mm = mm + 12;
|
---|
| 413 | }
|
---|
| 414 | ii = jj/100;
|
---|
| 415 | kk = 2 - ii + ii/4;
|
---|
| 416 | djul = (365.25*jj - fmod( 365.25*jj, 1.0 )) - 679006.0;
|
---|
| 417 | djul = djul + floor( 30.6001*(mm + 1) ) + tt + kk;
|
---|
| 418 | return djul;
|
---|
| 419 | }
|
---|
| 420 |
|
---|
| 421 | //
|
---|
| 422 | ////////////////////////////////////////////////////////////////////////////
|
---|
[5886] | 423 | double gpjd(double second, int nweek) {
|
---|
| 424 | double deltat;
|
---|
| 425 | deltat = nweek*7.0 + second/86400.0 ;
|
---|
| 426 | return( 44244.0 + deltat) ;
|
---|
| 427 | }
|
---|
| 428 |
|
---|
| 429 | //
|
---|
| 430 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 431 | void jdgp(double tjul, double & second, long & nweek) {
|
---|
[3171] | 432 | double deltat;
|
---|
| 433 | deltat = tjul - 44244.0 ;
|
---|
[5886] | 434 | nweek = (long) floor(deltat/7.0);
|
---|
[3171] | 435 | second = (deltat - (nweek)*7.0)*86400.0;
|
---|
| 436 | }
|
---|
| 437 |
|
---|
| 438 | //
|
---|
| 439 | ////////////////////////////////////////////////////////////////////////////
|
---|
[5886] | 440 | void jmt(double djul, long& jj, long& mm, double& dd) {
|
---|
| 441 | long ih, ih1, ih2 ;
|
---|
| 442 | double t1, t2, t3, t4;
|
---|
| 443 | t1 = 1.0 + djul - fmod( djul, 1.0 ) + 2400000.0;
|
---|
| 444 | t4 = fmod( djul, 1.0 );
|
---|
| 445 | ih = long( (t1 - 1867216.25)/36524.25 );
|
---|
| 446 | t2 = t1 + 1 + ih - ih/4;
|
---|
| 447 | t3 = t2 - 1720995.0;
|
---|
| 448 | ih1 = long( (t3 - 122.1)/365.25 );
|
---|
| 449 | t1 = 365.25*ih1 - fmod( 365.25*ih1, 1.0 );
|
---|
| 450 | ih2 = long( (t3 - t1)/30.6001 );
|
---|
| 451 | dd = t3 - t1 - (int)( 30.6001*ih2 ) + t4;
|
---|
| 452 | mm = ih2 - 1;
|
---|
| 453 | if ( ih2 > 13 ) mm = ih2 - 13;
|
---|
| 454 | jj = ih1;
|
---|
| 455 | if ( mm <= 2 ) jj = jj + 1;
|
---|
| 456 | }
|
---|
| 457 |
|
---|
| 458 | //
|
---|
| 459 | ////////////////////////////////////////////////////////////////////////////
|
---|
[3044] | 460 | void GPSweekFromDateAndTime(const QDateTime& dateTime,
|
---|
| 461 | int& GPSWeek, double& GPSWeeks) {
|
---|
| 462 |
|
---|
| 463 | static const QDateTime zeroEpoch(QDate(1980, 1, 6),QTime(),Qt::UTC);
|
---|
| 464 |
|
---|
| 465 | GPSWeek = zeroEpoch.daysTo(dateTime) / 7;
|
---|
| 466 |
|
---|
| 467 | int weekDay = dateTime.date().dayOfWeek() + 1; // Qt: Monday = 1
|
---|
| 468 | if (weekDay > 7) weekDay = 1;
|
---|
| 469 |
|
---|
| 470 | GPSWeeks = (weekDay - 1) * 86400.0
|
---|
| 471 | - dateTime.time().msecsTo(QTime()) / 1e3;
|
---|
| 472 | }
|
---|
| 473 |
|
---|
| 474 | //
|
---|
| 475 | ////////////////////////////////////////////////////////////////////////////
|
---|
[3171] | 476 | void GPSweekFromYMDhms(int year, int month, int day, int hour, int min,
|
---|
| 477 | double sec, int& GPSWeek, double& GPSWeeks) {
|
---|
| 478 |
|
---|
| 479 | double mjd = djul(year, month, day);
|
---|
| 480 |
|
---|
[5888] | 481 | long GPSWeek_long;
|
---|
| 482 | jdgp(mjd, GPSWeeks, GPSWeek_long);
|
---|
| 483 | GPSWeek = GPSWeek_long;
|
---|
[3171] | 484 | GPSWeeks += hour * 3600.0 + min * 60.0 + sec;
|
---|
| 485 | }
|
---|
| 486 |
|
---|
| 487 | //
|
---|
| 488 | ////////////////////////////////////////////////////////////////////////////
|
---|
[3044] | 489 | void mjdFromDateAndTime(const QDateTime& dateTime, int& mjd, double& dayfrac) {
|
---|
| 490 |
|
---|
| 491 | static const QDate zeroDate(1858, 11, 17);
|
---|
| 492 |
|
---|
| 493 | mjd = zeroDate.daysTo(dateTime.date());
|
---|
| 494 |
|
---|
| 495 | dayfrac = (dateTime.time().hour() +
|
---|
| 496 | (dateTime.time().minute() +
|
---|
| 497 | (dateTime.time().second() +
|
---|
| 498 | dateTime.time().msec() / 1000.0) / 60.0) / 60.0) / 24.0;
|
---|
| 499 | }
|
---|
[3408] | 500 |
|
---|
| 501 | //
|
---|
| 502 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 503 | bool findInVector(const vector<QString>& vv, const QString& str) {
|
---|
| 504 | std::vector<QString>::const_iterator it;
|
---|
| 505 | for (it = vv.begin(); it != vv.end(); ++it) {
|
---|
| 506 | if ( (*it) == str) {
|
---|
| 507 | return true;
|
---|
| 508 | }
|
---|
| 509 | }
|
---|
| 510 | return false;
|
---|
| 511 | }
|
---|
| 512 |
|
---|
[3664] | 513 | //
|
---|
| 514 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 515 | int readInt(const QString& str, int pos, int len, int& value) {
|
---|
| 516 | bool ok;
|
---|
| 517 | value = str.mid(pos, len).toInt(&ok);
|
---|
| 518 | return ok ? 0 : 1;
|
---|
| 519 | }
|
---|
| 520 |
|
---|
| 521 | //
|
---|
| 522 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 523 | int readDbl(const QString& str, int pos, int len, double& value) {
|
---|
| 524 | QString hlp = str.mid(pos, len);
|
---|
| 525 | for (int ii = 0; ii < hlp.length(); ii++) {
|
---|
| 526 | if (hlp[ii]=='D' || hlp[ii]=='d' || hlp[ii] == 'E') {
|
---|
| 527 | hlp[ii]='e';
|
---|
| 528 | }
|
---|
| 529 | }
|
---|
| 530 | bool ok;
|
---|
| 531 | value = hlp.toDouble(&ok);
|
---|
| 532 | return ok ? 0 : 1;
|
---|
| 533 | }
|
---|
[4338] | 534 |
|
---|
| 535 | // Topocentrical Distance and Elevation
|
---|
| 536 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 537 | void topos(double xRec, double yRec, double zRec,
|
---|
| 538 | double xSat, double ySat, double zSat,
|
---|
| 539 | double& rho, double& eleSat, double& azSat) {
|
---|
| 540 |
|
---|
| 541 | double dx[3];
|
---|
| 542 | dx[0] = xSat-xRec;
|
---|
| 543 | dx[1] = ySat-yRec;
|
---|
| 544 | dx[2] = zSat-zRec;
|
---|
| 545 |
|
---|
| 546 | rho = sqrt( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
|
---|
| 547 |
|
---|
| 548 | double xyzRec[3];
|
---|
| 549 | xyzRec[0] = xRec;
|
---|
| 550 | xyzRec[1] = yRec;
|
---|
| 551 | xyzRec[2] = zRec;
|
---|
| 552 |
|
---|
| 553 | double Ell[3];
|
---|
| 554 | double neu[3];
|
---|
| 555 | xyz2ell(xyzRec, Ell);
|
---|
| 556 | xyz2neu(Ell, dx, neu);
|
---|
| 557 |
|
---|
| 558 | eleSat = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / rho );
|
---|
| 559 | if (neu[2] < 0) {
|
---|
| 560 | eleSat *= -1.0;
|
---|
| 561 | }
|
---|
| 562 |
|
---|
| 563 | azSat = atan2(neu[1], neu[0]);
|
---|
| 564 | }
|
---|
[5230] | 565 |
|
---|
| 566 | // Degrees -> degrees, minutes, seconds
|
---|
| 567 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 568 | void deg2DMS(double decDeg, int& deg, int& min, double& sec) {
|
---|
| 569 | int sgn = (decDeg < 0.0 ? -1 : 1);
|
---|
| 570 | deg = sgn * static_cast<int>(decDeg);
|
---|
| 571 | min = static_cast<int>((decDeg - deg)*60);
|
---|
| 572 | sec = (decDeg - deg - min/60.0) * 3600.0;
|
---|
| 573 | }
|
---|
[5310] | 574 |
|
---|
| 575 | //
|
---|
| 576 | ////////////////////////////////////////////////////////////////////////////
|
---|
| 577 | QString fortranFormat(double value, int width, int prec) {
|
---|
| 578 | int expo = value == 0.0 ? 0 : log10(fabs(value));
|
---|
| 579 | double mant = value == 0.0 ? 0 : value / pow(10, expo);
|
---|
| 580 | if (fabs(mant) >= 1.0) {
|
---|
| 581 | mant /= 10.0;
|
---|
| 582 | expo += 1;
|
---|
| 583 | }
|
---|
| 584 | if (expo >= 0) {
|
---|
| 585 | return QString("%1e+%2").arg(mant, width-4, 'f', prec).arg(expo, 2, 10, QChar('0'));
|
---|
| 586 | }
|
---|
| 587 | else {
|
---|
| 588 | return QString("%1e-%2").arg(mant, width-4, 'f', prec).arg(-expo, 2, 10, QChar('0'));
|
---|
| 589 | }
|
---|
| 590 | }
|
---|
[5807] | 591 |
|
---|
| 592 | //
|
---|
| 593 | //////////////////////////////////////////////////////////////////////////////
|
---|
| 594 | void kalman(const Matrix& AA, const ColumnVector& ll, const DiagonalMatrix& PP,
|
---|
[5866] | 595 | SymmetricMatrix& QQ, ColumnVector& xx) {
|
---|
[5807] | 596 |
|
---|
| 597 | Tracer tracer("kalman");
|
---|
| 598 |
|
---|
| 599 | int nPar = AA.Ncols();
|
---|
| 600 | int nObs = AA.Nrows();
|
---|
| 601 | UpperTriangularMatrix SS = Cholesky(QQ).t();
|
---|
| 602 |
|
---|
| 603 | Matrix SA = SS*AA.t();
|
---|
| 604 | Matrix SRF(nObs+nPar, nObs+nPar); SRF = 0;
|
---|
| 605 | for (int ii = 1; ii <= nObs; ++ii) {
|
---|
| 606 | SRF(ii,ii) = 1.0 / sqrt(PP(ii,ii));
|
---|
| 607 | }
|
---|
| 608 |
|
---|
| 609 | SRF.SubMatrix (nObs+1, nObs+nPar, 1, nObs) = SA;
|
---|
| 610 | SRF.SymSubMatrix(nObs+1, nObs+nPar) = SS;
|
---|
| 611 |
|
---|
| 612 | UpperTriangularMatrix UU;
|
---|
| 613 | QRZ(SRF, UU);
|
---|
| 614 |
|
---|
| 615 | SS = UU.SymSubMatrix(nObs+1, nObs+nPar);
|
---|
| 616 | UpperTriangularMatrix SH_rt = UU.SymSubMatrix(1, nObs);
|
---|
| 617 | Matrix YY = UU.SubMatrix(1, nObs, nObs+1, nObs+nPar);
|
---|
| 618 |
|
---|
| 619 | UpperTriangularMatrix SHi = SH_rt.i();
|
---|
| 620 |
|
---|
| 621 | Matrix KT = SHi * YY;
|
---|
| 622 | SymmetricMatrix Hi; Hi << SHi * SHi.t();
|
---|
| 623 |
|
---|
[5866] | 624 | xx += KT.t() * (ll - AA * xx);
|
---|
[5807] | 625 | QQ << (SS.t() * SS);
|
---|
| 626 | }
|
---|
| 627 |
|
---|