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