source: ntrip/trunk/BNS/bnsutils.cpp@ 798

Last change on this file since 798 was 798, checked in by mervart, 16 years ago

* empty log message *

File size: 4.0 KB
Line 
1/* -------------------------------------------------------------------------
2 * BKG NTRIP Server
3 * -------------------------------------------------------------------------
4 *
5 * Class: bncutils
6 *
7 * Purpose: Auxiliary Functions
8 *
9 * Author: L. Mervart
10 *
11 * Created: 08-Apr-2008
12 *
13 * Changes:
14 *
15 * -----------------------------------------------------------------------*/
16
17#include <iostream>
18#include <ctime>
19#include <math.h>
20
21#include <QRegExp>
22#include <QStringList>
23#include <QDateTime>
24
25#include "bnsutils.h"
26#include "bnseph.h"
27
28using namespace std;
29
30//
31////////////////////////////////////////////////////////////////////////////
32void expandEnvVar(QString& str) {
33
34 QRegExp rx("(\\$\\{.+\\})");
35
36 if (rx.indexIn(str) != -1) {
37 QStringListIterator it(rx.capturedTexts());
38 if (it.hasNext()) {
39 QString rxStr = it.next();
40 QString envVar = rxStr.mid(2,rxStr.length()-3);
41 str.replace(rxStr, qgetenv(envVar.toAscii()));
42 }
43 }
44
45}
46
47//
48////////////////////////////////////////////////////////////////////////////
49QDateTime dateAndTimeFromGPSweek(int GPSWeek, double GPSWeeks) {
50
51 static const QDate zeroEpoch(1980, 1, 6);
52
53 QDate date(zeroEpoch);
54 QTime time(0,0,0,0);
55
56 int weekDays = int(GPSWeeks) / 86400;
57 date = date.addDays( GPSWeek * 7 + weekDays );
58 time = time.addMSecs( int( (GPSWeeks - 86400 * weekDays) * 1e3 ) );
59
60 return QDateTime(date,time);
61}
62
63//
64////////////////////////////////////////////////////////////////////////////
65void GPSweekFromDateAndTime(const QDateTime& dateTime,
66 int& GPSWeek, double& GPSWeeks) {
67
68 static const QDateTime zeroEpoch(QDate(1980, 1, 6));
69
70 GPSWeek = zeroEpoch.daysTo(dateTime) / 7;
71
72 int weekDay = dateTime.date().dayOfWeek() + 1; // Qt: Monday = 1
73 if (weekDay > 7) weekDay = 1;
74
75 GPSWeeks = (weekDay - 1) * 86400.0
76 - dateTime.time().msecsTo(QTime()) / 1e3;
77}
78
79//
80////////////////////////////////////////////////////////////////////////////
81void currentGPSWeeks(int& week, double& sec) {
82
83 QDateTime currDateTime = QDateTime::currentDateTime().toUTC();
84 QDate currDate = currDateTime.date();
85 QTime currTime = currDateTime.time();
86
87 week = int( (double(currDate.toJulianDay()) - 2444244.5) / 7 );
88
89 sec = (currDate.dayOfWeek() % 7) * 24.0 * 3600.0 +
90 currTime.hour() * 3600.0 +
91 currTime.minute() * 60.0 +
92 currTime.second() +
93 currTime.msec() / 1000.0;
94}
95
96//
97////////////////////////////////////////////////////////////////////////////
98void satellitePosition(double GPSweeks, const gpsEph* ep,
99 double& X, double& Y, double& Z, double& dt) {
100
101 const static double omegaEarth = 7292115.1467e-11;
102 const static double gmWGS = 398.6005e12;
103
104 X = Y = Z = dt = 0.0;
105
106 double a0 = ep->sqrt_A * ep->sqrt_A;
107 if (a0 == 0) {
108 return;
109 }
110
111 double n0 = sqrt(gmWGS/(a0*a0*a0));
112 double tk = GPSweeks - ep->TOE;
113 double n = n0 + ep->Delta_n;
114 double M = ep->M0 + n*tk;
115 double E = M;
116 double E_last;
117 do {
118 E_last = E;
119 E = M + ep->e*sin(E);
120 } while ( fabs(E-E_last)*a0 > 0.001 );
121 double v = 2.0*atan( sqrt( (1.0 + ep->e)/(1.0 - ep->e) )*tan( E/2 ) );
122 double u0 = v + ep->omega;
123 double sin2u0 = sin(2*u0);
124 double cos2u0 = cos(2*u0);
125 double r = a0*(1 - ep->e*cos(E)) + ep->Crc*cos2u0 + ep->Crs*sin2u0;
126 double i = ep->i0 + ep->IDOT*tk + ep->Cic*cos2u0 + ep->Cis*sin2u0;
127 double u = u0 + ep->Cuc*cos2u0 + ep->Cus*sin2u0;
128 double xp = r*cos(u);
129 double yp = r*sin(u);
130 double OM = ep->OMEGA0 + (ep->OMEGADOT - omegaEarth)*tk -
131 omegaEarth*ep->TOE;
132
133 double sinom = sin(OM);
134 double cosom = cos(OM);
135 double sini = sin(i);
136 double cosi = cos(i);
137 X = xp*cosom - yp*cosi*sinom;
138 Y = xp*sinom + yp*cosi*cosom;
139 Z = yp*sini;
140
141 double tc = GPSweeks - ep->TOC;
142 dt = ep->clock_bias + ep->clock_drift*tc + ep->clock_driftrate*tc*tc
143 - 4.442807633e-10 * ep->e * sqrt(a0) *sin(E);
144}
Note: See TracBrowser for help on using the repository browser.