1 /********************************************************************
2 ** @source JEEPS arithmetic/conversion functions
3 **
4 ** @author Copyright (C) 1999 Alan Bleasby
5 ** @version 1.0
6 ** @modified Dec 28 1999 Alan Bleasby. First version
7 ** @@
8 **
9 ** This library is free software; you can redistribute it and/or
10 ** modify it under the terms of the GNU Library General Public
11 ** License as published by the Free Software Foundation; either
12 ** version 2 of the License, or (at your option) any later version.
13 **
14 ** This library is distributed in the hope that it will be useful,
15 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
16 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 ** Library General Public License for more details.
18 **
19 ** You should have received a copy of the GNU Library General Public
20 ** License along with this library; if not, write to the
21 ** Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
22 ** Boston, MA 02110-1301, USA.
23 ********************************************************************/
24 #include "garmin_gps.h"
25 #include <stdlib.h>
26 #include <math.h>
27 #include <string.h>
28 #include "gpsdatum.h"
29
30
31 static int32 GPS_Math_LatLon_To_UTM_Param(double lat, double lon, int32 *zone,
32 char *zc, double *Mc, double *E0,
33 double *N0, double *F0);
34 static int32 GPS_Math_UTM_Param_To_Mc(int32 zone, char zc, double *Mc,
35 double *E0, double *N0, double *F0);
36
37
38
39 /* @func GPS_Math_Deg_To_Rad *******************************************
40 **
41 ** Convert degrees to radians
42 **
43 ** @param [r] v [double] degrees
44 **
45 ** @return [double] radians
46 ************************************************************************/
47
GPS_Math_Deg_To_Rad(double v)48 double GPS_Math_Deg_To_Rad(double v)
49 {
50 return v*(double)((double)GPS_PI/(double)180.);
51 }
52
53
54
55 /* @func GPS_Math_Rad_To_Deg *******************************************
56 **
57 ** Convert radians to degrees
58 **
59 ** @param [r] v [double] radians
60 **
61 ** @return [double] degrees
62 ************************************************************************/
63
GPS_Math_Rad_To_Deg(double v)64 double GPS_Math_Rad_To_Deg(double v)
65 {
66 return v*(double)((double)180./(double)GPS_PI);
67 }
68
69
70
71 /* @func GPS_Math_Deg_To_DegMin *****************************************
72 **
73 ** Convert degrees to degrees and minutes
74 **
75 ** @param [r] v [double] degrees
76 ** @param [w] d [int32 *] whole degrees
77 ** @param [w] m [double *] minutes
78 **
79 ** @return [void]
80 ************************************************************************/
81
GPS_Math_Deg_To_DegMin(double v,int32 * d,double * m)82 void GPS_Math_Deg_To_DegMin(double v, int32 *d, double *m)
83 {
84 int32 sign;
85
86 if(v<(double)0.)
87 {
88 v *= (double)-1.;
89 sign = 1;
90 }
91 else
92 sign = 0;
93
94 *d = (int32)v;
95 *m = (v-(double)*d) * (double)60.0;
96 if(*m>(double)59.999)
97 {
98 ++*d;
99 *m = (double)0.0;
100 }
101
102 if(sign)
103 *d = -*d;
104
105 return;
106 }
107
108
109
110 /* @func GPS_Math_DegMin_To_Deg *****************************************
111 **
112 ** Convert degrees and minutes to degrees
113 **
114 ** @param [r] d [int32] whole degrees
115 ** @param [r] m [double] minutes
116 ** @param [w] deg [double *] degrees
117 **
118 ** @return [void]
119 ************************************************************************/
120
GPS_Math_DegMin_To_Deg(int32 d,double m,double * deg)121 void GPS_Math_DegMin_To_Deg(int32 d, double m, double *deg)
122 {
123
124 *deg = ((double)abs(d)) + m / (double)60.0;
125 if(d<0)
126 *deg = -*deg;
127
128 return;
129 }
130
131
132
133 /* @func GPS_Math_Deg_To_DegMinSec *************************************
134 **
135 ** Convert degrees to degrees, minutes and seconds
136 **
137 ** @param [r] v [double] degrees
138 ** @param [w] d [int32 *] whole degrees
139 ** @param [w] m [int32 *] whole minutes
140 ** @param [w] s [double *] seconds
141 **
142 ** @return [void]
143 ************************************************************************/
144
GPS_Math_Deg_To_DegMinSec(double v,int32 * d,int32 * m,double * s)145 void GPS_Math_Deg_To_DegMinSec(double v, int32 *d, int32 *m, double *s)
146 {
147 int32 sign;
148 double t;
149
150 if(v<(double)0.)
151 {
152 v *= (double)-1.;
153 sign = 1;
154 }
155 else
156 sign = 0;
157
158 *d = (int32)v;
159 t = (v -(double)*d) * (double)60.0;
160 *m = (int32)((v-(double)*d) * (double)60.0);
161 *s = (t - (int32)t) * (double)60.0;
162
163 if(*s>(double)59.999)
164 {
165 ++t;
166 *s = (double)0.0;
167 }
168
169
170 if(t>(double)59.999)
171 {
172 ++*d;
173 t = 0;
174 }
175
176 *m = (int32)t;
177
178 if(sign)
179 *d = -*d;
180
181 return;
182 }
183
184
185
186 /* @func GPS_Math_DegMinSec_To_Deg *****************************************
187 **
188 ** Convert degrees, minutes and seconds to degrees
189 **
190 ** @param [r] d [int32] whole degrees
191 ** @param [r] m [int32] whole minutes
192 ** @param [r] s [double] seconds
193 ** @param [w] deg [double *] degrees
194 **
195 ** @return [void]
196 ************************************************************************/
197
GPS_Math_DegMinSec_To_Deg(int32 d,int32 m,double s,double * deg)198 void GPS_Math_DegMinSec_To_Deg(int32 d, int32 m, double s, double *deg)
199 {
200
201 *deg = ((double)abs(d)) + ((double)m + s / (double)60.0) / (double)60.0;
202 if(d<0)
203 *deg = -*deg;
204
205 return;
206 }
207
208
209
210 /* @func GPS_Math_Metres_To_Feet *******************************************
211 **
212 ** Convert metres to feet
213 **
214 ** @param [r] v [double] metres
215 **
216 ** @return [double] feet
217 ************************************************************************/
218
GPS_Math_Metres_To_Feet(double v)219 double GPS_Math_Metres_To_Feet(double v)
220 {
221 return v/0.3048;
222 }
223
224
225
226 /* @func GPS_Math_Feet_To_Metres *******************************************
227 **
228 ** Convert feet to metres
229 **
230 ** @param [r] v [double] feet
231 **
232 ** @return [double] metres
233 ************************************************************************/
234
GPS_Math_Feet_To_Metres(double v)235 double GPS_Math_Feet_To_Metres(double v)
236 {
237 return v * 0.3048;
238 }
239
240
241
242 /* @func GPS_Math_Deg_To_Semi *******************************************
243 **
244 ** Convert degrees to semi-circles
245 **
246 ** @param [r] v [double] degrees
247 **
248 ** @return [int32] semicircles
249 ************************************************************************/
250
GPS_Math_Deg_To_Semi(double v)251 int32 GPS_Math_Deg_To_Semi(double v)
252 {
253 return (int32)(( (double)(1U<<31) / (double)180) * v);
254 }
255
256
257
258 /* @func GPS_Math_Semi_To_Deg *******************************************
259 **
260 ** Convert semi-circles to degrees
261 **
262 ** @param [r] v [int32] semi-circles
263 **
264 ** @return [double] degrees
265 ************************************************************************/
266
GPS_Math_Semi_To_Deg(int32 v)267 double GPS_Math_Semi_To_Deg(int32 v)
268 {
269 return (double) (((double)v / (double)(1U<<31)) * (double)180);
270 }
271
272
273
274 /* @func GPS_Math_Utime_To_Gtime *******************************************
275 **
276 ** Convert Unix time (1970) to GPS time (1990)
277 **
278 ** @param [r] v [time_t] Unix time
279 **
280 ** @return [time_t] GPS time
281 ************************************************************************/
282
GPS_Math_Utime_To_Gtime(time_t v)283 time_t GPS_Math_Utime_To_Gtime(time_t v)
284 {
285 return v-631065600;
286 }
287
288
289
290 /* @func GPS_Math_Gtime_To_Utime *******************************************
291 **
292 ** Convert GPS time (1990) to Unix time (1970)
293 **
294 ** @param [r] v [time_t] GPS time
295 **
296 ** @return [time_t] Unix time
297 ************************************************************************/
298
GPS_Math_Gtime_To_Utime(time_t v)299 time_t GPS_Math_Gtime_To_Utime(time_t v)
300 {
301 return v+631065600;
302 }
303
304
305
306
307 /* @func GPS_Math_LatLonH_To_XYZ **********************************
308 **
309 ** Convert latitude and longitude and ellipsoidal height to
310 ** X, Y & Z coordinates
311 **
312 ** @param [r] phi [double] latitude (deg)
313 ** @param [r] lambda [double] longitude (deg)
314 ** @param [r] H [double] ellipsoidal height (metres)
315 ** @param [w] x [double *] X
316 ** @param [w] y [double *] Y
317 ** @param [w] z [double *] Z
318 ** @param [r] a [double] semi-major axis (metres)
319 ** @param [r] b [double] semi-minor axis (metres)
320 **
321 ** @return [void]
322 ************************************************************************/
GPS_Math_LatLonH_To_XYZ(double phi,double lambda,double H,double * x,double * y,double * z,double a,double b)323 void GPS_Math_LatLonH_To_XYZ(double phi, double lambda, double H,
324 double *x, double *y, double *z,
325 double a, double b)
326 {
327 double esq;
328 double nu;
329
330 phi = GPS_Math_Deg_To_Rad(phi);
331 lambda = GPS_Math_Deg_To_Rad(lambda);
332
333
334 esq = ((a*a)-(b*b)) / (a*a);
335
336 nu = a / pow(((double)1.0-esq*sin(phi)*sin(phi)),(double)0.5);
337 *x = (nu+H) * cos(phi) * cos(lambda);
338 *y = (nu+H) * cos(phi) * sin(lambda);
339 *z = (((double)1.0-esq)*nu+H) * sin(phi);
340
341 return;
342 }
343
344
345
346
347 /* @func GPS_Math_XYX_To_LatLonH ***************************************
348 **
349 ** Convert XYZ coordinates to latitude and longitude and ellipsoidal height
350 **
351 ** @param [w] phi [double] latitude (deg)
352 ** @param [w] lambda [double] longitude (deg)
353 ** @param [w] H [double] ellipsoidal height (metres)
354 ** @param [r] x [double *] X
355 ** @param [r] y [double *] Y
356 ** @param [r] z [double *] Z
357 ** @param [r] a [double] semi-major axis (metres)
358 ** @param [r] b [double] semi-minor axis (metres)
359 **
360 ** @return [void]
361 ************************************************************************/
GPS_Math_XYZ_To_LatLonH(double * phi,double * lambda,double * H,double x,double y,double z,double a,double b)362 void GPS_Math_XYZ_To_LatLonH(double *phi, double *lambda, double *H,
363 double x, double y, double z,
364 double a, double b)
365 {
366 double esq;
367 double nu=0.0;
368 double phix;
369 double nphi;
370 double rho;
371 int32 t1=0;
372 int32 t2=0;
373
374 if(x<(double)0 && y>=(double)0) t1=1;
375 if(x<(double)0 && y<(double)0) t2=1;
376
377 rho = pow(((x*x)+(y*y)),(double)0.5);
378 esq = ((a*a)-(b*b)) / (a*a);
379 phix = atan(z/(((double)1.0 - esq) * rho));
380 nphi = (double)1e20;
381
382 while(fabs(phix-nphi)>(double)0.00000000001)
383 {
384 nphi = phix;
385 nu = a / pow(((double)1.0-esq*sin(nphi)*sin(nphi)),(double)0.5);
386 phix = atan((z+esq*nu*sin(nphi))/rho);
387 }
388
389 *phi = GPS_Math_Rad_To_Deg(phix);
390 *H = (rho / cos(phix)) - nu;
391 *lambda = GPS_Math_Rad_To_Deg(atan(y/x));
392
393 if(t1) *lambda += (double)180.0;
394 if(t2) *lambda -= (double)180.0;
395
396 return;
397 }
398
399
400
401 /* @func GPS_Math_Airy1830LatLonH_To_XYZ **********************************
402 **
403 ** Convert Airy 1830 latitude and longitude and ellipsoidal height to
404 ** X, Y & Z coordinates
405 **
406 ** @param [r] phi [double] latitude (deg)
407 ** @param [r] lambda [double] longitude (deg)
408 ** @param [r] H [double] ellipsoidal height (metres)
409 ** @param [w] x [double *] X
410 ** @param [w] y [double *] Y
411 ** @param [w] z [double *] Z
412 **
413 ** @return [void]
414 ************************************************************************/
GPS_Math_Airy1830LatLonH_To_XYZ(double phi,double lambda,double H,double * x,double * y,double * z)415 void GPS_Math_Airy1830LatLonH_To_XYZ(double phi, double lambda, double H,
416 double *x, double *y, double *z)
417 {
418 double a = 6377563.396;
419 double b = 6356256.910;
420
421 GPS_Math_LatLonH_To_XYZ(phi,lambda,H,x,y,z,a,b);
422
423 return;
424 }
425
426
427
428 /* @func GPS_Math_WGS84LatLonH_To_XYZ **********************************
429 **
430 ** Convert WGS84 latitude and longitude and ellipsoidal height to
431 ** X, Y & Z coordinates
432 **
433 ** @param [r] phi [double] latitude (deg)
434 ** @param [r] lambda [double] longitude (deg)
435 ** @param [r] H [double] ellipsoidal height (metres)
436 ** @param [w] x [double *] X
437 ** @param [w] y [double *] Y
438 ** @param [w] z [double *] Z
439 **
440 ** @return [void]
441 ************************************************************************/
GPS_Math_WGS84LatLonH_To_XYZ(double phi,double lambda,double H,double * x,double * y,double * z)442 void GPS_Math_WGS84LatLonH_To_XYZ(double phi, double lambda, double H,
443 double *x, double *y, double *z)
444 {
445 double a = 6378137.000;
446 double b = 6356752.3141;
447
448 GPS_Math_LatLonH_To_XYZ(phi,lambda,H,x,y,z,a,b);
449
450 return;
451 }
452
453
454
455
456 /* @func GPS_Math_XYZ_To_Airy1830LatLonH **********************************
457 **
458 ** Convert XYZ to Airy 1830 latitude and longitude and ellipsoidal height
459 **
460 ** @param [r] phi [double] latitude (deg)
461 ** @param [r] lambda [double] longitude (deg)
462 ** @param [r] H [double] ellipsoidal height (metres)
463 ** @param [w] x [double *] X
464 ** @param [w] y [double *] Y
465 ** @param [w] z [double *] Z
466 **
467 ** @return [void]
468 ************************************************************************/
GPS_Math_XYZ_To_Airy1830LatLonH(double * phi,double * lambda,double * H,double x,double y,double z)469 void GPS_Math_XYZ_To_Airy1830LatLonH(double *phi, double *lambda, double *H,
470 double x, double y, double z)
471 {
472 double a = 6377563.396;
473 double b = 6356256.910;
474
475 GPS_Math_XYZ_To_LatLonH(phi,lambda,H,x,y,z,a,b);
476
477 return;
478 }
479
480
481
482 /* @func GPS_Math_XYZ_To_WGS84LatLonH **********************************
483 **
484 ** Convert XYZ to WGS84 latitude and longitude and ellipsoidal height
485 **
486 ** @param [r] phi [double] latitude (deg)
487 ** @param [r] lambda [double] longitude (deg)
488 ** @param [r] H [double] ellipsoidal height (metres)
489 ** @param [w] x [double *] X
490 ** @param [w] y [double *] Y
491 ** @param [w] z [double *] Z
492 **
493 ** @return [void]
494 ************************************************************************/
GPS_Math_XYZ_To_WGS84LatLonH(double * phi,double * lambda,double * H,double x,double y,double z)495 void GPS_Math_XYZ_To_WGS84LatLonH(double *phi, double *lambda, double *H,
496 double x, double y, double z)
497 {
498 double a = 6378137.000;
499 double b = 66356752.3141;
500
501 GPS_Math_XYZ_To_LatLonH(phi,lambda,H,x,y,z,a,b);
502
503 return;
504 }
505
506
507 /* @func GPS_Math_LatLon_To_EN **********************************
508 **
509 ** Convert latitude and longitude to eastings and northings
510 ** Standard Gauss-Kruger Transverse Mercator
511 **
512 ** @param [w] E [double *] easting (metres)
513 ** @param [w] N [double *] northing (metres)
514 ** @param [r] phi [double] latitude (deg)
515 ** @param [r] lambda [double] longitude (deg)
516 ** @param [r] N0 [double] true northing origin (metres)
517 ** @param [r] E0 [double] true easting origin (metres)
518 ** @param [r] phi0 [double] true latitude origin (deg)
519 ** @param [r] lambda0 [double] true longitude origin (deg)
520 ** @param [r] F0 [double] scale factor on central meridian
521 ** @param [r] a [double] semi-major axis (metres)
522 ** @param [r] b [double] semi-minor axis (metres)
523 **
524 ** @return [void]
525 ************************************************************************/
GPS_Math_LatLon_To_EN(double * E,double * N,double phi,double lambda,double N0,double E0,double phi0,double lambda0,double F0,double a,double b)526 void GPS_Math_LatLon_To_EN(double *E, double *N, double phi,
527 double lambda, double N0, double E0,
528 double phi0, double lambda0,
529 double F0, double a, double b)
530 {
531 double esq;
532 double n;
533 double etasq;
534 double nu;
535 double rho;
536 double M;
537 double I;
538 double II;
539 double III;
540 double IIIA;
541 double IV;
542 double V;
543 double VI;
544
545 double tmp;
546 double tmp2;
547 double fdf;
548 double fde;
549
550 phi0 = GPS_Math_Deg_To_Rad(phi0);
551 lambda0 = GPS_Math_Deg_To_Rad(lambda0);
552 phi = GPS_Math_Deg_To_Rad(phi);
553 lambda = GPS_Math_Deg_To_Rad(lambda);
554
555 esq = ((a*a)-(b*b)) / (a*a);
556 n = (a-b) / (a+b);
557
558 tmp = (double)1.0 - (esq * sin(phi) * sin(phi));
559 nu = a * F0 * pow(tmp,(double)-0.5);
560 rho = a * F0 * ((double)1.0 - esq) * pow(tmp,(double)-1.5);
561 etasq = (nu / rho) - (double)1.0;
562
563 fdf = (double)5.0 / (double)4.0;
564 tmp = (double)1.0 + n + (fdf * n * n) + (fdf * n * n * n);
565 tmp *= (phi - phi0);
566 tmp2 = (double)3.0*n + (double)3.0*n*n + ((double)21./(double)8.)*n*n*n;
567 tmp2 *= (sin(phi-phi0) * cos(phi+phi0));
568 tmp -= tmp2;
569
570 fde = ((double)15.0 / (double)8.0);
571 tmp2 = ((fde*n*n) + (fde*n*n*n)) * sin((double)2.0 * (phi-phi0));
572 tmp2 *= cos((double)2.0 * (phi+phi0));
573 tmp += tmp2;
574
575 tmp2 = ((double)35.0/(double)24.0) * n * n * n;
576 tmp2 *= sin((double)3.0 * (phi-phi0));
577 tmp2 *= cos((double)3.0 * (phi+phi0));
578 tmp -= tmp2;
579
580 M = b * F0 * tmp;
581 I = M + N0;
582 II = (nu / (double)2.0) * sin(phi) * cos(phi);
583 III = (nu / (double)24.0) * sin(phi) * cos(phi) * cos(phi) * cos(phi);
584 III *= ((double)5.0 - (tan(phi) * tan(phi)) + ((double)9.0 * etasq));
585 IIIA = (nu / (double)720.0) * sin(phi) * pow(cos(phi),(double)5.0);
586 IIIA *= ((double)61.0 - ((double)58.0*tan(phi)*tan(phi)) +
587 pow(tan(phi),(double)4.0));
588 IV = nu * cos(phi);
589
590 tmp = pow(cos(phi),(double)3.0);
591 tmp *= ((nu/rho) - tan(phi) * tan(phi));
592 V = (nu/(double)6.0) * tmp;
593
594 tmp = (double)5.0 - ((double)18.0 * tan(phi) * tan(phi));
595 tmp += tan(phi)*tan(phi)*tan(phi)*tan(phi) + ((double)14.0 * etasq);
596 tmp -= ((double)58.0 * tan(phi) * tan(phi) * etasq);
597 tmp2 = cos(phi)*cos(phi)*cos(phi)*cos(phi)*cos(phi) * tmp;
598 VI = (nu / (double)120.0) * tmp2;
599
600 *N = I + II*(lambda-lambda0)*(lambda-lambda0) +
601 III*pow((lambda-lambda0),(double)4.0) +
602 IIIA*pow((lambda-lambda0),(double)6.0);
603
604 *E = E0 + IV*(lambda-lambda0) + V*pow((lambda-lambda0),(double)3.0) +
605 VI * pow((lambda-lambda0),(double)5.0);
606
607 return;
608 }
609
610
611
612 /* @func GPS_Math_Airy1830MLatLonToINGEN ************************************
613 **
614 ** Convert Modified Airy 1830 datum latitude and longitude to Irish
615 ** National Grid Eastings and Northings
616 **
617 ** @param [r] phi [double] modified Airy latitude (deg)
618 ** @param [r] lambda [double] modified Airy longitude (deg)
619 ** @param [w] E [double *] NG easting (metres)
620 ** @param [w] N [double *] NG northing (metres)
621 **
622 ** @return [void]
623 ************************************************************************/
GPS_Math_Airy1830M_LatLonToINGEN(double phi,double lambda,double * E,double * N)624 void GPS_Math_Airy1830M_LatLonToINGEN(double phi, double lambda, double *E,
625 double *N)
626 {
627 double N0 = 250000;
628 double E0 = 200000;
629 double F0 = 1.000035;
630 double phi0 = 53.5;
631 double lambda0 = -8.;
632 double a = 6377340.189;
633 double b = 6356034.447;
634
635 GPS_Math_LatLon_To_EN(E,N,phi,lambda,N0,E0,phi0,lambda0,F0,a,b);
636
637 return;
638 }
639
640
641
642
643 /* @func GPS_Math_Airy1830LatLonToNGEN **************************************
644 **
645 ** Convert Airy 1830 datum latitude and longitude to UK Ordnance Survey
646 ** National Grid Eastings and Northings
647 **
648 ** @param [r] phi [double] WGS84 latitude (deg)
649 ** @param [r] lambda [double] WGS84 longitude (deg)
650 ** @param [w] E [double *] NG easting (metres)
651 ** @param [w] N [double *] NG northing (metres)
652 **
653 ** @return [void]
654 ************************************************************************/
GPS_Math_Airy1830LatLonToNGEN(double phi,double lambda,double * E,double * N)655 void GPS_Math_Airy1830LatLonToNGEN(double phi, double lambda, double *E,
656 double *N)
657 {
658 double N0 = -100000;
659 double E0 = 400000;
660 double F0 = 0.9996012717;
661 double phi0 = 49.;
662 double lambda0 = -2.;
663 double a = 6377563.396;
664 double b = 6356256.910;
665
666 GPS_Math_LatLon_To_EN(E,N,phi,lambda,N0,E0,phi0,lambda0,F0,a,b);
667
668 return;
669 }
670
671
672 /* @func int32 GPS_Math_WGS84_To_Swiss_EN ******************************
673 **
674 ** Convert WGS84 latitude and longitude to
675 ** Swiss CH-1903 National Grid Eastings and Northings
676 ** ( Oblique Mercator Projection )
677 **
678 ** @param [r] phi [double] WGS84 latitude (deg)
679 ** @param [r] lambda [double] WGS84 longitude (deg)
680 ** @param [w] E [double *] Swiss-NG easting (metres)
681 ** @param [w] N [double *] Swiss-NG northing (metres)
682 **
683 ** @return [void]
684 ************************************************************************/
685
GPS_Math_WGS84_To_Swiss_EN(double lat,double lon,double * E,double * N)686 int32 GPS_Math_WGS84_To_Swiss_EN(double lat, double lon, double *E,
687 double *N)
688 {
689 const double phi0 = 46.95240556;
690 const double lambda0 = 7.43958333;
691 const double E0 = 600000.0;
692 const double N0 = 200000.0;
693 double phi, lambda, alt, a, b;
694
695 if (lat < 44.89022757) return 0;
696 if (lon < -0.16386312) return 0;
697
698 a = GPS_Ellipse[4].a;
699 b = a - (a / GPS_Ellipse[4].invf);
700
701 GPS_Math_WGS84_To_Known_Datum_M(lat, lon, 0, &phi, &lambda, &alt, 123);
702 GPS_Math_Swiss_LatLon_To_EN(phi, lambda, E, N, phi0, lambda0, E0, N0, a, b);
703
704 return 1;
705 }
706
707
708 /* @GPS_Math_Swiss_EN_To_WGS84 *****************************************
709 **
710 ** Convert WGS84 latitude and longitude to
711 ** Swiss CH-1903 National Grid Eastings and Northings
712 **
713 ** @param [r] E [double] Swiss-NG easting (metres)
714 ** @param [r] N [double] Swiss-NG northing (metres)
715 ** @param [w] lat [double *] WGS84 latitude (deg)
716 ** @param [w] lon [double *] WGS84 longitude (deg)
717 **
718 ** @return [void]
719 ************************************************************************/
GPS_Math_Swiss_EN_To_WGS84(double E,double N,double * lat,double * lon)720 void GPS_Math_Swiss_EN_To_WGS84(double E, double N, double *lat, double *lon)
721 {
722 const double phi0 = 46.95240556;
723 const double lambda0 = 7.43958333;
724 const double E0 = 600000.0;
725 const double N0 = 200000.0;
726 double phi, lambda, alt, a, b;
727
728 a = GPS_Ellipse[4].a;
729 b = a - (a / GPS_Ellipse[4].invf);
730
731 GPS_Math_Swiss_EN_To_LatLon(E, N, &phi, &lambda, phi0, lambda0, E0, N0, a, b);
732 GPS_Math_Known_Datum_To_WGS84_M(phi, lambda, 0, lat, lon, &alt, 123);
733 }
734
735
736 /* @func GPS_Math_Cassini_LatLon_To_EN **********************************
737 **
738 ** Convert latitude and longitude to Cassini transverse cylindrical projection
739 ** easting and northing
740 **
741 ** @param [r] phi [double] latitude (deg)
742 ** @param [r] lambda [double] longitude (deg)
743 ** @param [w] E [double *] easting (metre)
744 ** @param [w] N [double *] northing (metre)
745 ** @param [r] phi0 [double] latitude of origin (deg)
746 ** @param [r] M0 [double] central meridian (deg)
747 ** @param [r] E0 [double] false easting
748 ** @param [r] N0 [double] false northing
749 ** @param [r] a [double] semi-major axis
750 ** @param [r] b [double] semi-minor axis
751 **
752 ** @return [void]
753 ************************************************************************/
GPS_Math_Cassini_LatLon_To_EN(double phi,double lambda,double * E,double * N,double phi0,double M0,double E0,double N0,double a,double b)754 void GPS_Math_Cassini_LatLon_To_EN(double phi, double lambda, double *E,
755 double *N, double phi0, double M0,
756 double E0, double N0, double a, double b)
757 {
758 double p2;
759 double po2;
760 double a2;
761 double b2;
762 double e2;
763 double e4;
764 double e6;
765 double AM0;
766 double c0;
767 double c1;
768 double c2;
769 double c3;
770 double om0;
771 double A0;
772 double A1;
773 double A2;
774 double A3;
775 double j;
776 double te4;
777 double phi0s2;
778 double phi0s4;
779 double phi0s6;
780 double lat;
781 double x;
782 double E1;
783 double E2;
784 double E3;
785 double E4;
786
787 double phis;
788 double phic;
789 double phit;
790 double phis2;
791 double phis4;
792 double phis6;
793 double RD;
794 double dlam;
795 double NN;
796 double TT;
797 double WW;
798 double WW2;
799 double WW3;
800 double WW4;
801 double WW5;
802 double CC;
803 double MM;
804
805
806 lambda = GPS_Math_Deg_To_Rad(lambda);
807 phi0 = GPS_Math_Deg_To_Rad(phi0);
808 phi = GPS_Math_Deg_To_Rad(phi);
809 M0 = GPS_Math_Deg_To_Rad(M0);
810
811
812 p2 = (double)GPS_PI * (double)2.;
813 po2 = (double)GPS_PI / (double)2.;
814
815 a2 = a*a;
816 b2 = b*b;
817 e2 = (a2-b2)/a2;
818 e4 = e2*e2;
819 e6 = e2*e4;
820
821 te4 = (double)3.0 * e4;
822 j = (double)45. * e6 / (double)1024.;
823 c0 = (double)1.0-e2/(double)4.-te4/(double)64.-(double)5.*e6/(double)256.;
824 c1 = (double)3.*e2/(double)8.+te4/(double)32.+j;
825 c2 = (double)15.*e4/(double)256.+j;
826 c3 = (double)35.*e6/(double)3072.;
827
828 lat = c0*phi0;
829 phi0s2 = c1 * sin((double)2.*phi0);
830 phi0s4 = c2 * sin((double)4.*phi0);
831 phi0s6 = c3 * sin((double)6.*phi0);
832 AM0 = a * (lat-phi0s2+phi0s4-phi0s6);
833
834 om0 = (double)1.0 - e2;
835 x = pow(om0,(double)0.5);
836 E1 = ((double)1.0 - x) / ((double)1.0 + x);
837 E2 = E1*E1;
838 E3 = E1*E2;
839 E4 = E1*E3;
840 A0 = (double)3.*E1/(double)2.-(double)27.*E3/(double)32.;
841 A1 = (double)21.*E2/(double)16.-(double)55.*E4/(double)32.;
842 A2 = (double)151.*E3/(double)96.;
843 A3 = (double)1097.*E4/(double)512.;
844
845
846 dlam = lambda - M0;
847 if(dlam>GPS_PI)
848 dlam -= p2;
849 if(dlam<-GPS_PI)
850 dlam += p2;
851
852 phis = sin(phi);
853 phic = cos(phi);
854 phit = tan(phi);
855 RD = pow((double)1.-e2*phis*phis,(double).5);
856 NN = a/RD;
857 TT = phit*phit;
858 WW = dlam*phic;
859 WW2 = WW*WW;
860 WW3 = WW*WW2;
861 WW4 = WW*WW3;
862 WW5 = WW*WW4;
863 CC = e2*phic*phic/om0;
864 lat = c0*phi;
865 phis2 = c1 * sin((double)2.*phi);
866 phis4 = c2 * sin((double)4.*phi);
867 phis6 = c3 * sin((double)6.*phi);
868 MM = a * (lat-phis2+phis4-phis6);
869
870 *E = NN*(WW-(TT*WW3/(double)6.)-((double)8.-TT+(double)8.*CC)*
871 (TT*WW5/(double)120.)) + E0;
872 *N = MM-AM0+NN*phit*((WW2/(double)2.)+((double)5.-TT+(double)6.*CC) *
873 WW4/(double)24.) + N0;
874 return;
875 }
876
877
878
879
880 /* @func GPS_Math_Cassini_EN_To_LatLon **********************************
881 **
882 ** Convert Cassini transverse cylindrical easting and northing projection
883 ** to latitude and longitude
884 **
885 ** @param [r] E [double] easting (metre)
886 ** @param [r] N [double] northing (metre)
887 ** @param [w] phi [double *] latitude (deg)
888 ** @param [w] lambda [double *] longitude (deg)
889 ** @param [r] phi0 [double] latitude of origin (deg)
890 ** @param [r] M0 [double] central meridian (deg)
891 ** @param [r] E0 [double] false easting
892 ** @param [r] N0 [double] false northing
893 ** @param [r] a [double] semi-major axis
894 ** @param [r] b [double] semi-minor axis
895 **
896 ** @return [void]
897 ************************************************************************/
GPS_Math_Cassini_EN_To_LatLon(double E,double N,double * phi,double * lambda,double phi0,double M0,double E0,double N0,double a,double b)898 void GPS_Math_Cassini_EN_To_LatLon(double E, double N, double *phi,
899 double *lambda, double phi0, double M0,
900 double E0, double N0, double a, double b)
901
902 {
903 double p2;
904 double po2;
905 double a2;
906 double b2;
907 double e2;
908 double e4;
909 double e6;
910 double AM0;
911 double c0;
912 double c1;
913 double c2;
914 double c3;
915 double om0;
916 double A0;
917 double A1;
918 double A2;
919 double A3;
920 double j;
921 double te4;
922 double phi0s2;
923 double phi0s4;
924 double phi0s6;
925 double lat;
926 double x;
927 double E1;
928 double E2;
929 double E3;
930 double E4;
931
932 double dx;
933 double dy;
934 double mu;
935 double mus2;
936 double mus4;
937 double mus6;
938 double mus8;
939 double M1;
940 double phi1;
941 double phi1s;
942 double phi1c;
943 double phi1t;
944 double T;
945 double T1;
946 double N1;
947 double R1;
948 double RD;
949 double DD;
950 double D2;
951 double D3;
952 double D4;
953 double D5;
954 double tol;
955
956 M0 = GPS_Math_Deg_To_Rad(M0);
957 phi0 = GPS_Math_Deg_To_Rad(phi0);
958
959 p2 = (double)GPS_PI * (double)2.;
960 po2 = (double)GPS_PI / (double)2.;
961
962 a2 = a*a;
963 b2 = b*b;
964 e2 = (a2-b2)/a2;
965 e4 = e2*e2;
966 e6 = e2*e4;
967
968 te4 = (double)3.0 * e4;
969 j = (double)45. * e6 / (double)1024.;
970 c0 = (double)1.0-e2/(double)4.-te4/(double)64.-(double)5.*e6/(double)256.;
971 c1 = (double)3.*e2/(double)8.+te4/(double)32.+j;
972 c2 = (double)15.*e4/(double)256.+j;
973 c3 = (double)35.*e6/(double)3072.;
974
975 lat = c0*phi0;
976 phi0s2 = c1 * sin((double)2.*phi0);
977 phi0s4 = c2 * sin((double)4.*phi0);
978 phi0s6 = c3 * sin((double)6.*phi0);
979 AM0 = a * (lat-phi0s2+phi0s4-phi0s6);
980
981 om0 = (double)1.0 - e2;
982 x = pow(om0,(double)0.5);
983 E1 = ((double)1.0 - x) / ((double)1.0 + x);
984 E2 = E1*E1;
985 E3 = E1*E2;
986 E4 = E1*E3;
987 A0 = (double)3.*E1/(double)2.-(double)27.*E3/(double)32.;
988 A1 = (double)21.*E2/(double)16.-(double)55.*E4/(double)32.;
989 A2 = (double)151.*E3/(double)96.;
990 A3 = (double)1097.*E4/(double)512.;
991
992
993
994 tol = (double)1.e-5;
995
996 dx = E - E0;
997 dy = N - N0;
998 M1 = AM0 + dy;
999 mu = M1 / (a*c0);
1000 mus2 = A0 * sin((double)2.*mu);
1001 mus4 = A1 * sin((double)4.*mu);
1002 mus6 = A2 * sin((double)6.*mu);
1003 mus8 = A3 * sin((double)8.*mu);
1004 phi1 = mu + mus2 + mus4 + mus6 + mus8;
1005
1006 if((((po2-tol)<phi1)&&(phi1<(po2+tol))))
1007 {
1008 *phi = po2;
1009 *lambda = M0;
1010 }
1011 else if((((-po2-tol)<phi1)&&(phi1<(-po2+tol))))
1012 {
1013 *phi = -po2;
1014 *lambda = M0;
1015 }
1016 else
1017 {
1018 phi1s = sin(phi1);
1019 phi1c = cos(phi1);
1020 phi1t = tan(phi1);
1021 T1 = phi1t*phi1t;
1022 RD = pow((double)1.-e2*phi1s*phi1s,(double).5);
1023 N1 = a/RD;
1024 R1 = N1 * om0 / (RD*RD);
1025 DD = dx/N1;
1026 D2 = DD*DD;
1027 D3 = DD*D2;
1028 D4 = DD*D3;
1029 D5 = DD*D4;
1030 T = (double)1. + (double)3.*T1;
1031 *phi = phi1-(N1*phi1t/R1)*(D2/(double)2.-T*D4/(double)24.);
1032 *lambda = M0+(DD-T1*D3/(double)3.+T*T1*D5/(double)15.)/phi1c;
1033
1034 if(*phi>po2)
1035 *phi=po2;
1036 else if(*phi<-po2)
1037 *phi=-po2;
1038
1039 if(*lambda>GPS_PI)
1040 *lambda -= p2;
1041 if(*lambda<-GPS_PI)
1042 *lambda += p2;
1043
1044 if(*lambda>GPS_PI)
1045 *lambda=GPS_PI;
1046 else if(*lambda<-GPS_PI)
1047 *lambda=-GPS_PI;
1048 }
1049
1050 *lambda = GPS_Math_Rad_To_Deg(*lambda);
1051 *phi = GPS_Math_Rad_To_Deg(*phi);
1052
1053 return;
1054 }
1055
1056
1057
1058
1059 /* @func int32 GPS_Math_WGS84_To_ICS_EN ******************************
1060 **
1061 ** Convert WGS84 latitude and longitude to
1062 ** Israeli old Grid Eastings and Northings
1063 ** ( Israeli Cassini Soldner )
1064 **
1065 ** @param [r] phi [double] WGS84 latitude (deg)
1066 ** @param [r] lambda [double] WGS84 longitude (deg)
1067 ** @param [w] E [double *] ICS easting (metres)
1068 ** @param [w] N [double *] ICS northing (metres)
1069 **
1070 ** @return [void]
1071 ************************************************************************/
1072
GPS_Math_WGS84_To_ICS_EN(double lat,double lon,double * E,double * N)1073 int32 GPS_Math_WGS84_To_ICS_EN(double lat, double lon, double *E,
1074 double *N)
1075 {
1076 double const phi0 = 31.73409694444; // 31 44 2.749
1077 double const lambda0 = 35.21208055556; // 35 12 43.49
1078 double const E0 = 170251.555;
1079 double const N0 = 1126867.909;
1080 double phi, lambda, alt, a, b;
1081
1082 int32 datum = GPS_Lookup_Datum_Index("Palestine 1923");
1083 int32 ellipse = GPS_Datum[datum].ellipse;
1084
1085 a = GPS_Ellipse[ellipse].a;
1086 b = a - (a / GPS_Ellipse[ellipse].invf);
1087
1088 GPS_Math_WGS84_To_Known_Datum_M(lat, lon, 0, &phi, &lambda, &alt, datum);
1089 GPS_Math_Cassini_LatLon_To_EN(phi, lambda, E, N,
1090 phi0, lambda0, E0, N0, a, b);
1091
1092 return 1;
1093 }
1094
1095
1096 /* @GPS_Math_ICS_EN_To_WGS84 *****************************************
1097 **
1098 ** Convert WGS84 latitude and longitude to
1099 ** Israeli Oldl Grid Eastings and Northings
1100 **
1101 ** @param [r] E [double] ICS easting (metres)
1102 ** @param [r] N [double] ICS northing (metres)
1103 ** @param [w] lat [double *] WGS84 latitude (deg)
1104 ** @param [w] lon [double *] WGS84 longitude (deg)
1105 **
1106 ** @return [void]
1107 ************************************************************************/
GPS_Math_ICS_EN_To_WGS84(double E,double N,double * lat,double * lon)1108 void GPS_Math_ICS_EN_To_WGS84(double E, double N, double *lat, double *lon)
1109 {
1110 double const phi0 = 31.73409694444; // 31 44 2.749
1111 double const lambda0 = 35.21208055556; // 35 12 43.49
1112 double const E0 = 170251.555;
1113 double const N0 = 1126867.909;
1114 double phi, lambda, alt, a, b;
1115 int32 datum = GPS_Lookup_Datum_Index("Palestine 1923");
1116 int32 ellipse = GPS_Datum[datum].ellipse;
1117
1118 a = GPS_Ellipse[ellipse].a;
1119 b = a - (a / GPS_Ellipse[ellipse].invf);
1120
1121 GPS_Math_Cassini_EN_To_LatLon(E, N, &phi, &lambda, phi0, lambda0,
1122 E0, N0, a, b);
1123 GPS_Math_Known_Datum_To_WGS84_M(phi, lambda, 0, lat, lon, &alt, datum);
1124 }
1125
1126
1127
1128 /* @func GPS_Math_EN_To_LatLon **************************************
1129 **
1130 ** Convert Eastings and Northings to latitude and longitude
1131 **
1132 ** @param [w] E [double] NG easting (metres)
1133 ** @param [w] N [double] NG northing (metres)
1134 ** @param [r] phi [double *] Airy latitude (deg)
1135 ** @param [r] lambda [double *] Airy longitude (deg)
1136 ** @param [r] N0 [double] true northing origin (metres)
1137 ** @param [r] E0 [double] true easting origin (metres)
1138 ** @param [r] phi0 [double] true latitude origin (deg)
1139 ** @param [r] lambda0 [double] true longitude origin (deg)
1140 ** @param [r] F0 [double] scale factor on central meridian
1141 ** @param [r] a [double] semi-major axis (metres)
1142 ** @param [r] b [double] semi-minor axis (metres)
1143 **
1144 ** @return [void]
1145 ************************************************************************/
GPS_Math_EN_To_LatLon(double E,double N,double * phi,double * lambda,double N0,double E0,double phi0,double lambda0,double F0,double a,double b)1146 void GPS_Math_EN_To_LatLon(double E, double N, double *phi,
1147 double *lambda, double N0, double E0,
1148 double phi0, double lambda0,
1149 double F0, double a, double b)
1150 {
1151 double esq;
1152 double n;
1153 double etasq;
1154 double nu;
1155 double rho;
1156 double M;
1157 double VII;
1158 double VIII;
1159 double IX;
1160 double X;
1161 double XI;
1162 double XII;
1163 double XIIA;
1164 double phix;
1165 double nphi=0.0;
1166
1167 double tmp;
1168 double tmp2;
1169 double fdf;
1170 double fde;
1171
1172 phi0 = GPS_Math_Deg_To_Rad(phi0);
1173 lambda0 = GPS_Math_Deg_To_Rad(lambda0);
1174
1175 n = (a-b) / (a+b);
1176 fdf = (double)5.0 / (double)4.0;
1177 fde = ((double)15.0 / (double)8.0);
1178
1179 esq = ((a*a)-(b*b)) / (a*a);
1180
1181
1182 phix = ((N-N0)/(a*F0)) + phi0;
1183
1184 tmp = (double)1.0 - (esq * sin(phix) * sin(phix));
1185 nu = a * F0 * pow(tmp,(double)-0.5);
1186 rho = a * F0 * ((double)1.0 - esq) * pow(tmp,(double)-1.5);
1187 etasq = (nu / rho) - (double)1.0;
1188
1189 M = (double)-1e20;
1190
1191 while(N-N0-M > (double)0.000001)
1192 {
1193 nphi = phix;
1194
1195 tmp = (double)1.0 + n + (fdf * n * n) + (fdf * n * n * n);
1196 tmp *= (nphi - phi0);
1197 tmp2 = (double)3.0*n + (double)3.0*n*n +
1198 ((double)21./(double)8.)*n*n*n;
1199 tmp2 *= (sin(nphi-phi0) * cos(nphi+phi0));
1200 tmp -= tmp2;
1201
1202
1203 tmp2 = ((fde*n*n) + (fde*n*n*n)) * sin((double)2.0 * (nphi-phi0));
1204 tmp2 *= cos((double)2.0 * (nphi+phi0));
1205 tmp += tmp2;
1206
1207 tmp2 = ((double)35.0/(double)24.0) * n * n * n;
1208 tmp2 *= sin((double)3.0 * (nphi-phi0));
1209 tmp2 *= cos((double)3.0 * (nphi+phi0));
1210 tmp -= tmp2;
1211
1212 M = b * F0 * tmp;
1213
1214 if(N-N0-M > (double)0.000001)
1215 phix = ((N-N0-M)/(a*F0)) + nphi;
1216 }
1217
1218
1219 VII = tan(nphi) / ((double)2.0 * rho * nu);
1220
1221 tmp = (double)5.0 + (double)3.0 * tan(nphi) * tan(nphi) + etasq;
1222 tmp -= (double)9.0 * tan(nphi) * tan(nphi) * etasq;
1223 VIII = (tan(nphi)*tmp) / ((double)24.0 * rho * nu*nu*nu);
1224
1225 tmp = (double)61.0 + (double)90.0 * tan(nphi) * tan(nphi);
1226 tmp += (double)45.0 * pow(tan(nphi),(double)4.0);
1227 IX = tan(nphi) / ((double)720.0 * rho * pow(nu,(double)5.0)) * tmp;
1228
1229 X = (double)1.0 / (cos(nphi) * nu);
1230
1231 tmp = (nu / rho) + (double)2.0 * tan(nphi) * tan(nphi);
1232 XI = ((double)1.0 / (cos(nphi) * (double)6.0 * nu*nu*nu)) * tmp;
1233
1234 tmp = (double)5.0 + (double)28.0 * tan(nphi)*tan(nphi);
1235 tmp += (double)24.0 * pow(tan(nphi),(double)4.0);
1236 XII = ((double)1.0 / ((double)120.0 * pow(nu,(double)5.0) * cos(nphi)))
1237 * tmp;
1238
1239 tmp = (double)61.0 + (double)662.0 * tan(nphi) * tan(nphi);
1240 tmp += (double)1320.0 * pow(tan(nphi),(double)4.0);
1241 tmp += (double)720.0 * pow(tan(nphi),(double)6.0);
1242 XIIA = ((double)1.0 / (cos(nphi) * (double)5040.0 * pow(nu,(double)7.0)))
1243 * tmp;
1244
1245 *phi = nphi - VII*pow((E-E0),(double)2.0) + VIII*pow((E-E0),(double)4.0) -
1246 IX*pow((E-E0),(double)6.0);
1247
1248 *lambda = lambda0 + X*(E-E0) - XI*pow((E-E0),(double)3.0) +
1249 XII*pow((E-E0),(double)5.0) - XIIA*pow((E-E0),(double)7.0);
1250
1251 *phi = GPS_Math_Rad_To_Deg(*phi);
1252 *lambda = GPS_Math_Rad_To_Deg(*lambda);
1253
1254 return;
1255 }
1256
1257
1258
1259
1260 /* @func GPS_Math_NGENToAiry1830LatLon **************************************
1261 **
1262 ** Convert to UK Ordnance Survey National Grid Eastings and Northings to
1263 ** Airy 1830 datum latitude and longitude
1264 **
1265 ** @param [r] E [double] NG easting (metres)
1266 ** @param [r] N [double] NG northing (metres)
1267 ** @param [w] phi [double *] Airy latitude (deg)
1268 ** @param [w] lambda [double *] Airy longitude (deg)
1269 **
1270 ** @return [void]
1271 ************************************************************************/
GPS_Math_NGENToAiry1830LatLon(double E,double N,double * phi,double * lambda)1272 void GPS_Math_NGENToAiry1830LatLon(double E, double N, double *phi,
1273 double *lambda)
1274 {
1275 double N0 = -100000;
1276 double E0 = 400000;
1277 double F0 = 0.9996012717;
1278 double phi0 = 49.;
1279 double lambda0 = -2.;
1280 double a = 6377563.396;
1281 double b = 6356256.910;
1282
1283 GPS_Math_EN_To_LatLon(E,N,phi,lambda,N0,E0,phi0,lambda0,F0,a,b);
1284
1285 return;
1286 }
1287
1288
1289
1290 /* @func GPS_Math_INGENToAiry1830MLatLon **************************************
1291 **
1292 ** Convert Irish National Grid Eastings and Northings to modified
1293 ** Airy 1830 datum latitude and longitude
1294 **
1295 ** @param [r] E [double] ING easting (metres)
1296 ** @param [r] N [double] ING northing (metres)
1297 ** @param [w] phi [double *] modified Airy latitude (deg)
1298 ** @param [w] lambda [double *] modified Airy longitude (deg)
1299 **
1300 ** @return [void]
1301 ************************************************************************/
GPS_Math_INGENToAiry1830MLatLon(double E,double N,double * phi,double * lambda)1302 void GPS_Math_INGENToAiry1830MLatLon(double E, double N, double *phi,
1303 double *lambda)
1304 {
1305 double N0 = 250000;
1306 double E0 = 200000;
1307 double F0 = 1.000035;
1308 double phi0 = 53.5;
1309 double lambda0 = -8.;
1310 double a = 6377340.189;
1311 double b = 6356034.447;
1312
1313 GPS_Math_EN_To_LatLon(E,N,phi,lambda,N0,E0,phi0,lambda0,F0,a,b);
1314
1315 return;
1316 }
1317
1318
1319
1320 /* @func GPS_Math_EN_To_UKOSNG_Map *************************************
1321 **
1322 ** Convert Airy 1830 eastings and northings to Ordnance Survey map
1323 ** two letter code plus modified eastings and northings
1324 **
1325 ** @param [r] E [double] NG easting (metres)
1326 ** @param [r] N [double] NG northing (metres)
1327 ** @param [w] mE [double *] modified easting (metres)
1328 ** @param [w] mN [double *] modified northing (metres)
1329 ** @param [w] map [char *] map code
1330 **
1331 ** @return [int32] success
1332 ************************************************************************/
GPS_Math_EN_To_UKOSNG_Map(double E,double N,double * mE,double * mN,char * map)1333 int32 GPS_Math_EN_To_UKOSNG_Map(double E, double N, double *mE,
1334 double *mN, char *map)
1335 {
1336 int32 t;
1337 int32 idx;
1338
1339 if(E>=(double)700000. || E<(double)0.0 || N<(double)0.0 ||
1340 N>=(double)1300000.0)
1341 return 0;
1342
1343 idx = ((int32)N/100000)*7 + (int32)E/100000;
1344 (void) strcpy(map,UKNG[idx]);
1345
1346 t = ((int32)E / 100000) * 100000;
1347 *mE = E - (double)t;
1348
1349 t = ((int32)N / 100000) * 100000;
1350 *mN = N - (double)t;
1351
1352 return 1;
1353 }
1354
1355
1356
1357 /* @func GPS_Math_UKOSNG_Map_To_EN *************************************
1358 **
1359 ** Convert Ordnance Survey map eastings and northings plus
1360 ** two letter code to Airy 1830 eastings and northings
1361 **
1362 ** @param [w] map [char *] map code
1363 ** @param [r] mapE [double] easting (metres)
1364 ** @param [r] mapN [double] northing (metres)
1365 ** @param [w] E [double *] full Airy easting (metres)
1366 ** @param [w] N [double *] full Airy northing (metres)
1367
1368 **
1369 ** @return [int32] success
1370 ************************************************************************/
GPS_Math_UKOSNG_Map_To_EN(char * map,double mapE,double mapN,double * E,double * N)1371 int32 GPS_Math_UKOSNG_Map_To_EN(char *map, double mapE, double mapN, double *E,
1372 double *N)
1373 {
1374 int32 t;
1375 int32 idx;
1376
1377 if(mapE>=(double)100000.0 || mapE<(double)0.0 || mapN<(double)0.0 ||
1378 mapN>(double)100000.0)
1379 return 0;
1380
1381 idx=0;
1382 while(*UKNG[idx])
1383 {
1384 if(!strcmp(UKNG[idx],map)) break;
1385 ++idx;
1386 }
1387 if(!*UKNG[idx])
1388 return 0;
1389
1390
1391 t = (idx / 7) * 100000;
1392 *N = mapN + (double)t;
1393
1394 t = (idx % 7) * 100000;
1395 *E = mapE + (double)t;
1396
1397 return 1;
1398 }
1399
1400
1401
1402 /* @func GPS_Math_Molodensky *******************************************
1403 **
1404 ** Transform one datum to another
1405 **
1406 ** @param [r] Sphi [double] source latitude (deg)
1407 ** @param [r] Slam [double] source longitude (deg)
1408 ** @param [r] SH [double] source height (metres)
1409 ** @param [r] Sa [double] source semi-major axis (metres)
1410 ** @param [r] Sif [double] source inverse flattening
1411 ** @param [w] Dphi [double *] dest latitude (deg)
1412 ** @param [w] Dlam [double *] dest longitude (deg)
1413 ** @param [w] DH [double *] dest height (metres)
1414 ** @param [r] Da [double] dest semi-major axis (metres)
1415 ** @param [r] Dif [double] dest inverse flattening
1416 ** @param [r] dx [double] dx
1417 ** @param [r] dy [double] dy
1418 ** @param [r] dz [double] dz
1419 **
1420 ** @return [void]
1421 ************************************************************************/
GPS_Math_Molodensky(double Sphi,double Slam,double SH,double Sa,double Sif,double * Dphi,double * Dlam,double * DH,double Da,double Dif,double dx,double dy,double dz)1422 void GPS_Math_Molodensky(double Sphi, double Slam, double SH, double Sa,
1423 double Sif, double *Dphi, double *Dlam,
1424 double *DH, double Da, double Dif, double dx,
1425 double dy, double dz)
1426 {
1427 double Sf;
1428 double Df;
1429 double esq;
1430 double bda;
1431 double da;
1432 double df;
1433 double N;
1434 double M;
1435 double tmp;
1436 double tmp2;
1437 double dphi;
1438 double dlambda;
1439 double dheight;
1440 double phis;
1441 double phic;
1442 double lams;
1443 double lamc;
1444
1445 Sf = (double)1.0 / Sif;
1446 Df = (double)1.0 / Dif;
1447
1448 esq = (double)2.0*Sf - pow(Sf,(double)2.0);
1449 bda = (double)1.0 - Sf;
1450 Sphi = GPS_Math_Deg_To_Rad(Sphi);
1451 Slam = GPS_Math_Deg_To_Rad(Slam);
1452
1453 da = Da - Sa;
1454 df = Df - Sf;
1455
1456 phis = sin(Sphi);
1457 phic = cos(Sphi);
1458 lams = sin(Slam);
1459 lamc = cos(Slam);
1460
1461 N = Sa / sqrt((double)1.0 - esq*pow(phis,(double)2.0));
1462
1463 tmp = ((double)1.0-esq) /pow(((double)1.0-esq*pow(phis,(double)2.0)),1.5);
1464 M = Sa * tmp;
1465
1466 tmp = df * ((M/bda)+N*bda) * phis * phic;
1467 tmp2 = da * N * esq * phis * phic / Sa;
1468 tmp2 += ((-dx*phis*lamc-dy*phis*lams) + dz*phic);
1469 dphi = (tmp2 + tmp) / (M + SH);
1470
1471 dlambda = (-dx*lams+dy*lamc) / ((N+SH)*phic);
1472
1473 dheight = dx*phic*lamc + dy*phic*lams + dz*phis - da*(Sa/N) +
1474 df*bda*N*phis*phis;
1475
1476 *Dphi = Sphi + dphi;
1477 *Dlam = Slam + dlambda;
1478 *DH = SH + dheight;
1479
1480 *Dphi = GPS_Math_Rad_To_Deg(*Dphi);
1481 *Dlam = GPS_Math_Rad_To_Deg(*Dlam);
1482
1483 return;
1484 }
1485
1486
1487
1488 /* @func GPS_Math_Known_Datum_To_WGS84_M **********************************
1489 **
1490 ** Transform datum to WGS84 using Molodensky
1491 **
1492 ** @param [r] Sphi [double] source latitude (deg)
1493 ** @param [r] Slam [double] source longitude (deg)
1494 ** @param [r] SH [double] source height (metres)
1495 ** @param [w] Dphi [double *] dest latitude (deg)
1496 ** @param [w] Dlam [double *] dest longitude (deg)
1497 ** @param [w] DH [double *] dest height (metres)
1498 ** @param [r] n [int32] datum number from GPS_Datum structure
1499 **
1500 ** @return [void]
1501 ************************************************************************/
GPS_Math_Known_Datum_To_WGS84_M(double Sphi,double Slam,double SH,double * Dphi,double * Dlam,double * DH,int32 n)1502 void GPS_Math_Known_Datum_To_WGS84_M(double Sphi, double Slam, double SH,
1503 double *Dphi, double *Dlam, double *DH,
1504 int32 n)
1505 {
1506 double Sa;
1507 double Sif;
1508 double Da;
1509 double Dif;
1510 double x;
1511 double y;
1512 double z;
1513 int32 idx;
1514
1515 Da = (double) 6378137.0;
1516 Dif = (double) 298.257223563;
1517
1518 idx = GPS_Datum[n].ellipse;
1519 Sa = GPS_Ellipse[idx].a;
1520 Sif = GPS_Ellipse[idx].invf;
1521 x = GPS_Datum[n].dx;
1522 y = GPS_Datum[n].dy;
1523 z = GPS_Datum[n].dz;
1524
1525 GPS_Math_Molodensky(Sphi,Slam,SH,Sa,Sif,Dphi,Dlam,DH,Da,Dif,x,y,z);
1526
1527 return;
1528 }
1529
1530
1531
1532 /* @func GPS_Math_WGS84_To_Known_Datum_M ********************************
1533 **
1534 ** Transform WGS84 to other datum using Molodensky
1535 **
1536 ** @param [r] Sphi [double] source latitude (deg)
1537 ** @param [r] Slam [double] source longitude (deg)
1538 ** @param [r] SH [double] source height (metres)
1539 ** @param [w] Dphi [double *] dest latitude (deg)
1540 ** @param [w] Dlam [double *] dest longitude (deg)
1541 ** @param [w] DH [double *] dest height (metres)
1542 ** @param [r] n [int32] datum number from GPS_Datum structure
1543 **
1544 ** @return [void]
1545 ************************************************************************/
GPS_Math_WGS84_To_Known_Datum_M(double Sphi,double Slam,double SH,double * Dphi,double * Dlam,double * DH,int32 n)1546 void GPS_Math_WGS84_To_Known_Datum_M(double Sphi, double Slam, double SH,
1547 double *Dphi, double *Dlam, double *DH,
1548 int32 n)
1549 {
1550 double Sa;
1551 double Sif;
1552 double Da;
1553 double Dif;
1554 double x;
1555 double y;
1556 double z;
1557 int32 idx;
1558
1559 Sa = (double) 6378137.0;
1560 Sif = (double) 298.257223563;
1561
1562 idx = GPS_Datum[n].ellipse;
1563 Da = GPS_Ellipse[idx].a;
1564 Dif = GPS_Ellipse[idx].invf;
1565 x = -GPS_Datum[n].dx;
1566 y = -GPS_Datum[n].dy;
1567 z = -GPS_Datum[n].dz;
1568
1569 GPS_Math_Molodensky(Sphi,Slam,SH,Sa,Sif,Dphi,Dlam,DH,Da,Dif,x,y,z);
1570
1571 return;
1572 }
1573
1574
1575
1576 /* @func GPS_Math_Known_Datum_To_WGS84_C **********************************
1577 **
1578 ** Transform datum to WGS84 using Cartesian coordinates
1579 **
1580 ** @param [r] Sphi [double] source latitude (deg)
1581 ** @param [r] Slam [double] source longitude (deg)
1582 ** @param [r] SH [double] source height (metres)
1583 ** @param [w] Dphi [double *] dest latitude (deg)
1584 ** @param [w] Dlam [double *] dest longitude (deg)
1585 ** @param [w] DH [double *] dest height (metres)
1586 ** @param [r] n [int32] datum number from GPS_Datum structure
1587 **
1588 ** @return [void]
1589 ************************************************************************/
GPS_Math_Known_Datum_To_WGS84_C(double Sphi,double Slam,double SH,double * Dphi,double * Dlam,double * DH,int32 n)1590 void GPS_Math_Known_Datum_To_WGS84_C(double Sphi, double Slam, double SH,
1591 double *Dphi, double *Dlam, double *DH,
1592 int32 n)
1593 {
1594 double Sa;
1595 double Sif;
1596 double Sb;
1597 double Da;
1598 double Dif;
1599 double Db;
1600 double x;
1601 double y;
1602 double z;
1603 int32 idx;
1604 double sx;
1605 double sy;
1606 double sz;
1607
1608 Da = (double) 6378137.0;
1609 Dif = (double) 298.257223563;
1610 Db = Da - (Da / Dif);
1611
1612 idx = GPS_Datum[n].ellipse;
1613 Sa = GPS_Ellipse[idx].a;
1614 Sif = GPS_Ellipse[idx].invf;
1615 Sb = Sa - (Sa / Sif);
1616
1617 x = GPS_Datum[n].dx;
1618 y = GPS_Datum[n].dy;
1619 z = GPS_Datum[n].dz;
1620
1621 GPS_Math_LatLonH_To_XYZ(Sphi,Slam,SH,&sx,&sy,&sz,Sa,Sb);
1622 sx += x;
1623 sy += y;
1624 sz += z;
1625
1626 GPS_Math_XYZ_To_LatLonH(Dphi,Dlam,DH,sx,sy,sz,Da,Db);
1627
1628 return;
1629 }
1630
1631
1632
1633 /* @func GPS_Math_WGS84_To_Known_Datum_C ********************************
1634 **
1635 ** Transform WGS84 to other datum using Cartesian coordinates
1636 **
1637 ** @param [r] Sphi [double] source latitude (deg)
1638 ** @param [r] Slam [double] source longitude (deg)
1639 ** @param [r] SH [double] source height (metres)
1640 ** @param [w] Dphi [double *] dest latitude (deg)
1641 ** @param [w] Dlam [double *] dest longitude (deg)
1642 ** @param [w] DH [double *] dest height (metres)
1643 ** @param [r] n [int32] datum number from GPS_Datum structure
1644 **
1645 ** @return [void]
1646 ************************************************************************/
GPS_Math_WGS84_To_Known_Datum_C(double Sphi,double Slam,double SH,double * Dphi,double * Dlam,double * DH,int32 n)1647 void GPS_Math_WGS84_To_Known_Datum_C(double Sphi, double Slam, double SH,
1648 double *Dphi, double *Dlam, double *DH,
1649 int32 n)
1650 {
1651 double Sa;
1652 double Sif;
1653 double Da;
1654 double Dif;
1655 double x;
1656 double y;
1657 double z;
1658 int32 idx;
1659 double Sb;
1660 double Db;
1661 double dx;
1662 double dy;
1663 double dz;
1664
1665 Sa = (double) 6378137.0;
1666 Sif = (double) 298.257223563;
1667 Sb = Sa - (Sa / Sif);
1668
1669 idx = GPS_Datum[n].ellipse;
1670 Da = GPS_Ellipse[idx].a;
1671 Dif = GPS_Ellipse[idx].invf;
1672 Db = Da - (Da / Dif);
1673
1674 x = -GPS_Datum[n].dx;
1675 y = -GPS_Datum[n].dy;
1676 z = -GPS_Datum[n].dz;
1677
1678 GPS_Math_LatLonH_To_XYZ(Sphi,Slam,SH,&dx,&dy,&dz,Sa,Sb);
1679 dx += x;
1680 dy += y;
1681 dz += z;
1682
1683 GPS_Math_XYZ_To_LatLonH(Dphi,Dlam,DH,dx,dy,dz,Da,Db);
1684
1685 return;
1686 }
1687
1688
1689
1690 /* @func GPS_Math_Known_Datum_To_Known_Datum_M *************************
1691 **
1692 ** Transform WGS84 to other datum using Molodensky
1693 **
1694 ** @param [r] Sphi [double] source latitude (deg)
1695 ** @param [r] Slam [double] source longitude (deg)
1696 ** @param [r] SH [double] source height (metres)
1697 ** @param [w] Dphi [double *] dest latitude (deg)
1698 ** @param [w] Dlam [double *] dest longitude (deg)
1699 ** @param [w] DH [double *] dest height (metres)
1700 ** @param [r] n1 [int32] source datum number from GPS_Datum structure
1701 ** @param [r] n2 [int32] dest datum number from GPS_Datum structure
1702 **
1703 ** @return [void]
1704 ************************************************************************/
GPS_Math_Known_Datum_To_Known_Datum_M(double Sphi,double Slam,double SH,double * Dphi,double * Dlam,double * DH,int32 n1,int32 n2)1705 void GPS_Math_Known_Datum_To_Known_Datum_M(double Sphi, double Slam, double SH,
1706 double *Dphi, double *Dlam,
1707 double *DH, int32 n1, int32 n2)
1708 {
1709 double Sa;
1710 double Sif;
1711 double Da;
1712 double Dif;
1713 double x1;
1714 double y1;
1715 double z1;
1716 double x2;
1717 double y2;
1718 double z2;
1719 double x;
1720 double y;
1721 double z;
1722
1723 int32 idx1;
1724 int32 idx2;
1725
1726
1727 idx1 = GPS_Datum[n1].ellipse;
1728 Sa = GPS_Ellipse[idx1].a;
1729 Sif = GPS_Ellipse[idx1].invf;
1730 x1 = GPS_Datum[n1].dx;
1731 y1 = GPS_Datum[n1].dy;
1732 z1 = GPS_Datum[n1].dz;
1733
1734 idx2 = GPS_Datum[n2].ellipse;
1735 Da = GPS_Ellipse[idx2].a;
1736 Dif = GPS_Ellipse[idx2].invf;
1737 x2 = GPS_Datum[n2].dx;
1738 y2 = GPS_Datum[n2].dy;
1739 z2 = GPS_Datum[n2].dz;
1740
1741 x = -(x2-x1);
1742 y = -(y2-y1);
1743 z = -(z2-z1);
1744
1745 GPS_Math_Molodensky(Sphi,Slam,SH,Sa,Sif,Dphi,Dlam,DH,Da,Dif,x,y,z);
1746
1747 return;
1748 }
1749
1750
1751
1752 /* @func GPS_Math_Known_Datum_To_Known_Datum_C *************************
1753 **
1754 ** Transform known datum to other datum using Cartesian coordinates
1755 **
1756 ** @param [r] Sphi [double] source latitude (deg)
1757 ** @param [r] Slam [double] source longitude (deg)
1758 ** @param [r] SH [double] source height (metres)
1759 ** @param [w] Dphi [double *] dest latitude (deg)
1760 ** @param [w] Dlam [double *] dest longitude (deg)
1761 ** @param [w] DH [double *] dest height (metres)
1762 ** @param [r] n1 [int32] source datum number from GPS_Datum structure
1763 ** @param [r] n2 [int32] dest datum number from GPS_Datum structure
1764 **
1765 ** @return [void]
1766 ************************************************************************/
GPS_Math_Known_Datum_To_Known_Datum_C(double Sphi,double Slam,double SH,double * Dphi,double * Dlam,double * DH,int32 n1,int32 n2)1767 void GPS_Math_Known_Datum_To_Known_Datum_C(double Sphi, double Slam, double SH,
1768 double *Dphi, double *Dlam,
1769 double *DH, int32 n1, int32 n2)
1770 {
1771 double Sa;
1772 double Sif;
1773 double Da;
1774 double Dif;
1775 double x1;
1776 double y1;
1777 double z1;
1778 double x2;
1779 double y2;
1780 double z2;
1781
1782 int32 idx1;
1783 int32 idx2;
1784
1785 double Sb;
1786 double Db;
1787 double dx;
1788 double dy;
1789 double dz;
1790
1791 idx1 = GPS_Datum[n1].ellipse;
1792 Sa = GPS_Ellipse[idx1].a;
1793 Sif = GPS_Ellipse[idx1].invf;
1794 Sb = Sa - (Sa / Sif);
1795
1796 x1 = GPS_Datum[n1].dx;
1797 y1 = GPS_Datum[n1].dy;
1798 z1 = GPS_Datum[n1].dz;
1799
1800 idx2 = GPS_Datum[n2].ellipse;
1801 Da = GPS_Ellipse[idx2].a;
1802 Dif = GPS_Ellipse[idx2].invf;
1803 Db = Da - (Da / Dif);
1804
1805 x2 = GPS_Datum[n2].dx;
1806 y2 = GPS_Datum[n2].dy;
1807 z2 = GPS_Datum[n2].dz;
1808
1809 GPS_Math_LatLonH_To_XYZ(Sphi,Slam,SH,&dx,&dy,&dz,Sa,Sb);
1810 dx += -(x2-x1);
1811 dy += -(y2-y1);
1812 dz += -(z2-z1);
1813
1814 GPS_Math_XYZ_To_LatLonH(Dphi,Dlam,DH,dx,dy,dz,Da,Db);
1815
1816 return;
1817 }
1818
1819
1820
1821 /* @func GPS_Math_WGS84_To_UKOSMap_M ***********************************
1822 **
1823 ** Convert WGS84 lat/lon to Ordnance survey map code and easting and
1824 ** northing. Uses Molodensky
1825 **
1826 ** @param [r] lat [double] WGS84 latitude (deg)
1827 ** @param [r] lon [double] WGS84 longitude (deg)
1828 ** @param [w] mE [double *] map easting (metres)
1829 ** @param [w] mN [double *] map northing (metres)
1830 ** @param [w] map [char *] map two letter code
1831 **
1832 ** @return [int32] success
1833 ************************************************************************/
GPS_Math_WGS84_To_UKOSMap_M(double lat,double lon,double * mE,double * mN,char * map)1834 int32 GPS_Math_WGS84_To_UKOSMap_M(double lat, double lon, double *mE,
1835 double *mN, char *map)
1836 {
1837 double alat;
1838 double alon;
1839 double aht;
1840 double aE;
1841 double aN;
1842
1843
1844 GPS_Math_WGS84_To_Known_Datum_M(lat,lon,30,&alat,&alon,&aht,86);
1845
1846 GPS_Math_Airy1830LatLonToNGEN(alat,alon,&aE,&aN);
1847
1848 if(!GPS_Math_EN_To_UKOSNG_Map(aE,aN,mE,mN,map))
1849 return 0;
1850
1851 return 1;
1852 }
1853
1854
1855
1856 /* @func GPS_Math_UKOSMap_To_WGS84_M ***********************************
1857 **
1858 ** Transform UK Ordnance survey map position to WGS84 lat/lon
1859 ** Uses Molodensky transformation
1860 **
1861 ** @param [r] map [char *] map two letter code
1862 ** @param [r] mE [double] map easting (metres)
1863 ** @param [r] mN [double] map northing (metres)
1864 ** @param [w] lat [double *] WGS84 latitude (deg)
1865 ** @param [w] lon [double *] WGS84 longitude (deg)
1866 **
1867 ** @return [int32] success
1868 ************************************************************************/
GPS_Math_UKOSMap_To_WGS84_M(char * map,double mE,double mN,double * lat,double * lon)1869 int32 GPS_Math_UKOSMap_To_WGS84_M(char *map, double mE, double mN,
1870 double *lat, double *lon)
1871 {
1872 double E;
1873 double N;
1874 double alat;
1875 double alon;
1876 double ht;
1877
1878 if(!GPS_Math_UKOSNG_Map_To_EN(map,mE,mN,&E,&N))
1879 return 0;
1880
1881 GPS_Math_NGENToAiry1830LatLon(E,N,&alat,&alon);
1882
1883 GPS_Math_Known_Datum_To_WGS84_M(alat,alon,0,lat,lon,&ht,86);
1884
1885 return 1;
1886 }
1887
1888
1889
1890 /* @func GPS_Math_WGS84_To_UKOSMap_C ***********************************
1891 **
1892 ** Convert WGS84 lat/lon to Ordnance survey map code and easting and
1893 ** northing. Uses cartesian transformation
1894 **
1895 ** @param [r] lat [double] WGS84 latitude (deg)
1896 ** @param [r] lon [double] WGS84 longitude (deg)
1897 ** @param [w] mE [double *] map easting (metres)
1898 ** @param [w] mN [double *] map northing (metres)
1899 ** @param [w] map [char *] map two letter code
1900 **
1901 ** @return [int32] success
1902 ************************************************************************/
GPS_Math_WGS84_To_UKOSMap_C(double lat,double lon,double * mE,double * mN,char * map)1903 int32 GPS_Math_WGS84_To_UKOSMap_C(double lat, double lon, double *mE,
1904 double *mN, char *map)
1905 {
1906 double alat;
1907 double alon;
1908 double aht;
1909 double aE;
1910 double aN;
1911
1912
1913 GPS_Math_WGS84_To_Known_Datum_C(lat,lon,30,&alat,&alon,&aht,86);
1914
1915 GPS_Math_Airy1830LatLonToNGEN(alat,alon,&aE,&aN);
1916
1917 if(!GPS_Math_EN_To_UKOSNG_Map(aE,aN,mE,mN,map))
1918 return 0;
1919
1920 return 1;
1921 }
1922
1923
1924
1925 /* @func GPS_Math_UKOSMap_To_WGS84_C ***********************************
1926 **
1927 ** Transform UK Ordnance survey map position to WGS84 lat/lon
1928 ** Uses cartesian transformation
1929 **
1930 ** @param [r] map [char *] map two letter code
1931 ** @param [r] mE [double] map easting (metres)
1932 ** @param [r] mN [double] map northing (metres)
1933 ** @param [w] lat [double *] WGS84 latitude (deg)
1934 ** @param [w] lon [double *] WGS84 longitude (deg)
1935 **
1936 ** @return [int32] success
1937 ************************************************************************/
GPS_Math_UKOSMap_To_WGS84_C(char * map,double mE,double mN,double * lat,double * lon)1938 int32 GPS_Math_UKOSMap_To_WGS84_C(char *map, double mE, double mN,
1939 double *lat, double *lon)
1940 {
1941 double E;
1942 double N;
1943 double alat;
1944 double alon;
1945 double ht;
1946
1947 if(!GPS_Math_UKOSNG_Map_To_EN(map,mE,mN,&E,&N))
1948 return 0;
1949
1950 GPS_Math_NGENToAiry1830LatLon(E,N,&alat,&alon);
1951
1952 GPS_Math_Known_Datum_To_WGS84_C(alat,alon,0,lat,lon,&ht,86);
1953
1954 return 1;
1955 }
1956
1957
1958 /* @funcstatic GPS_Math_LatLon_To_UTM_Param *****************************
1959 **
1960 ** Transform NAD33
1961 **
1962 ** @param [r] lat [double] NAD latitude (deg)
1963 ** @param [r] lon [double] NAD longitude (deg)
1964 ** @param [w] zone [int32 *] zone number
1965 ** @param [w] zc [char *] zone character
1966 ** @param [w] Mc [double *] central meridian
1967 ** @param [w] E0 [double *] false easting
1968 ** @param [w] N0 [double *] false northing
1969 ** @param [w] F0 [double *] scale factor
1970 **
1971 ** @return [int32] success
1972 ************************************************************************/
GPS_Math_LatLon_To_UTM_Param(double lat,double lon,int32 * zone,char * zc,double * Mc,double * E0,double * N0,double * F0)1973 static int32 GPS_Math_LatLon_To_UTM_Param(double lat, double lon, int32 *zone,
1974 char *zc, double *Mc, double *E0,
1975 double *N0, double *F0)
1976 {
1977 int32 ilon;
1978 int32 ilat;
1979 int32 psign;
1980 int32 lsign;
1981
1982 if(lat >= (double)84.0 || lat < (double)-80.0)
1983 return 0;
1984
1985 psign = lsign = 0;
1986 if(lon < (double)0.0)
1987 lsign=1;
1988 if(lat < (double)0.0)
1989 psign=1;
1990
1991 ilon = abs((int32)lon);
1992 ilat = abs((int32)lat);
1993
1994 if(!lsign)
1995 {
1996 *zone = 31 + (ilon / 6);
1997 *Mc = (double)((ilon / 6) * 6 + 3);
1998 }
1999 else
2000 {
2001 *zone = 30 - (ilon / 6);
2002 *Mc = -(double)((ilon / 6) * 6 + 3);
2003 }
2004
2005 if(!psign)
2006 {
2007 *zc = 'N' + (UC)ilat / 8;
2008 if(*zc > 'N') ++*zc;
2009 }
2010 else
2011 {
2012 *zc = 'M' - (UC)(ilat / 8);
2013 if(*zc <= 'I') --*zc;
2014 }
2015
2016
2017 if(lat>=(double)56.0 && lat<(double)64.0 && lon>=(double)3.0 &&
2018 lon<(double)12.0)
2019 {
2020 *zone = 32;
2021 *zc = 'V';
2022 *Mc = (double)9.0;
2023 }
2024
2025 if(*zc=='X' && lon>=(double)0.0 && lon<(double)42.0)
2026 {
2027 if(lon<(double)9.0)
2028 {
2029 *zone = 31;
2030 *Mc = (double)3.0;
2031 }
2032 else if(lon<(double)21.0)
2033 {
2034 *zone = 33;
2035 *Mc = (double)15.0;
2036 }
2037 else if(lon<(double)33.0)
2038 {
2039 *zone = 35;
2040 *Mc = (double)27.0;
2041 }
2042 else
2043 {
2044 *zone = 37;
2045 *Mc = (double)39.0;
2046 }
2047 }
2048
2049 if(!psign)
2050 *N0 = (double)0.0;
2051 else
2052 *N0 = (double)10000000;
2053
2054 *E0 = (double)500000;
2055 *F0 = (double)0.9996;
2056
2057 return 1;
2058 }
2059
2060
2061
2062 /* @func GPS_Math_NAD83_To_UTM_EN **************************************
2063 **
2064 ** Transform NAD33 lat/lon to UTM zone, easting and northing
2065 **
2066 ** @param [r] lat [double] NAD latitude (deg)
2067 ** @param [r] lon [double] NAD longitude (deg)
2068 ** @param [w] E [double *] easting (metres)
2069 ** @param [w] N [double *] northing (metres)
2070 ** @param [w] zone [int32 *] zone number
2071 ** @param [w] zc [char *] zone character
2072 **
2073 ** @return [int32] success
2074 ************************************************************************/
GPS_Math_NAD83_To_UTM_EN(double lat,double lon,double * E,double * N,int32 * zone,char * zc)2075 int32 GPS_Math_NAD83_To_UTM_EN(double lat, double lon, double *E,
2076 double *N, int32 *zone, char *zc)
2077 {
2078 double phi0;
2079 double lambda0;
2080 double N0;
2081 double E0;
2082 double F0;
2083 double a;
2084 double b;
2085
2086 if(!GPS_Math_LatLon_To_UTM_Param(lat,lon,zone,zc,&lambda0,&E0,
2087 &N0,&F0))
2088 return 0;
2089
2090 phi0 = (double)0.0;
2091
2092 a = (double) GPS_Ellipse[21].a;
2093 b = a - (a/GPS_Ellipse[21].invf);
2094
2095 GPS_Math_LatLon_To_EN(E,N,lat,lon,N0,E0,phi0,lambda0,F0,a,b);
2096
2097 return 1;
2098 }
2099
2100
2101
2102 /* @func GPS_Math_WGS84_To_UTM_EN **************************************
2103 **
2104 ** Transform WGS84 lat/lon to UTM zone, easting and northing
2105 **
2106 ** @param [r] lat [double] WGS84 latitude (deg)
2107 ** @param [r] lon [double] WGS84 longitude (deg)
2108 ** @param [w] E [double *] easting (metres)
2109 ** @param [w] N [double *] northing (metres)
2110 ** @param [w] zone [int32 *] zone number
2111 ** @param [w] zc [char *] zone character
2112 **
2113 ** @return [int32] success
2114 ************************************************************************/
GPS_Math_WGS84_To_UTM_EN(double lat,double lon,double * E,double * N,int32 * zone,char * zc)2115 int32 GPS_Math_WGS84_To_UTM_EN(double lat, double lon, double *E,
2116 double *N, int32 *zone, char *zc)
2117 {
2118 double phi;
2119 double lambda;
2120 double H;
2121
2122 GPS_Math_WGS84_To_Known_Datum_M(lat,lon,0,&phi,&lambda,&H,77);
2123 if(!GPS_Math_NAD83_To_UTM_EN(phi,lambda,E,N,zone,zc))
2124 return 0;
2125
2126 return 1;
2127 }
2128
2129
2130
2131 /* @funcstatic GPS_Math_UTM_Param_To_Mc ********************************
2132 **
2133 ** Convert UTM zone and zone character to central meridian value.
2134 ** Also return false eastings, northings and scale factor
2135 **
2136 ** @param [w] zone [int32] zone number
2137 ** @param [w] zc [char] zone character
2138 ** @param [w] Mc [double *] central meridian
2139 ** @param [w] E0 [double *] false easting
2140 ** @param [w] N0 [double *] false northing
2141 ** @param [w] F0 [double *] scale factor
2142 **
2143 ** @return [int32] success
2144 ************************************************************************/
GPS_Math_UTM_Param_To_Mc(int32 zone,char zc,double * Mc,double * E0,double * N0,double * F0)2145 static int32 GPS_Math_UTM_Param_To_Mc(int32 zone, char zc, double *Mc,
2146 double *E0, double *N0, double *F0)
2147 {
2148
2149 if(zone>60 || zone<0 || zc<'C' || zc>'X')
2150 return 0;
2151
2152 if(zone > 30)
2153 *Mc = (double)((zone-31)*6) + (double)3.0;
2154 else
2155 *Mc = (double) -(((30-zone)*6)+3);
2156
2157 if(zone==32 && zc=='V')
2158 *Mc = (double)9.0;
2159
2160 if(zone==31 && zc=='X')
2161 *Mc = (double)3.0;
2162 if(zone==33 && zc=='X')
2163 *Mc = (double)15.0;
2164 if(zone==35 && zc=='X')
2165 *Mc = (double)27.0;
2166 if(zone==37 && zc=='X')
2167 *Mc = (double)39.0;
2168
2169 if(zc>'M')
2170 *N0 = (double)0.0;
2171 else
2172 *N0 = (double)10000000;
2173
2174 *E0 = (double)500000;
2175 *F0 = (double)0.9996;
2176
2177 return 1;
2178 }
2179
2180
2181
2182 /* @func GPS_Math_UTM_EN_To_NAD83 **************************************
2183 **
2184 ** Transform UTM zone, easting and northing to NAD83 lat/lon
2185 **
2186 ** @param [r] lat [double *] NAD latitude (deg)
2187 ** @param [r] lon [double *] NAD longitude (deg)
2188 ** @param [w] E [double] easting (metres)
2189 ** @param [w] N [double] northing (metres)
2190 ** @param [w] zone [int32] zone number
2191 ** @param [w] zc [char] zone character
2192 **
2193 ** @return [int32] success
2194 ************************************************************************/
GPS_Math_UTM_EN_To_NAD83(double * lat,double * lon,double E,double N,int32 zone,char zc)2195 int32 GPS_Math_UTM_EN_To_NAD83(double *lat, double *lon, double E,
2196 double N, int32 zone, char zc)
2197 {
2198 return GPS_Math_UTM_EN_To_Known_Datum(lat, lon, E, N, zone, zc, 77);
2199 }
2200
2201
2202
2203 /* @func GPS_Math_UTM_EN_To_WGS84 **************************************
2204 **
2205 ** Transform UTM zone, easting and northing to WGS84 lat/lon
2206 **
2207 ** @param [w] lat [double *] WGS84 latitude (deg)
2208 ** @param [r] lon [double *] WGS84 longitude (deg)
2209 ** @param [w] E [double] easting (metres)
2210 ** @param [w] N [double] northing (metres)
2211 ** @param [w] zone [int32] zone number
2212 ** @param [w] zc [char] zone character
2213 **
2214 ** @return [int32] success
2215 ************************************************************************/
GPS_Math_UTM_EN_To_WGS84(double * lat,double * lon,double E,double N,int32 zone,char zc)2216 int32 GPS_Math_UTM_EN_To_WGS84(double *lat, double *lon, double E,
2217 double N, int32 zone, char zc)
2218 {
2219 double lambda0;
2220 double N0;
2221 double E0;
2222 double F0;
2223
2224 if (!GPS_Math_UTM_Param_To_Mc(zone,zc,&lambda0,&E0,&N0,&F0)) return 0;
2225
2226 GPS_Math_UTM_EN_to_LatLon(GPS_Datum[118].ellipse, N, E, lat, lon, lambda0, E0, N0);
2227 return 1;
2228 }
2229
2230
2231 /* @func GPS_Math_Known_Datum_To_UTM_EN *********************************
2232 **
2233 ** Transform known datum lat/lon to UTM zone, easting and northing
2234 **
2235 ** @param [r] lat [double] WGS84 latitude (deg)
2236 ** @param [r] lon [double] WGS84 longitude (deg)
2237 ** @param [w] E [double *] easting (metres)
2238 ** @param [w] N [double *] northing (metres)
2239 ** @param [w] zone [int32 *] zone number
2240 ** @param [w] zc [char *] zone character
2241 ** @param [r] n [int32] datum number from GPS_Datum structure
2242 **
2243 ** @return [int32] success
2244 ************************************************************************/
GPS_Math_Known_Datum_To_UTM_EN(double lat,double lon,double * E,double * N,int32 * zone,char * zc,const int n)2245 int32 GPS_Math_Known_Datum_To_UTM_EN(double lat, double lon, double *E,
2246 double *N, int32 *zone, char *zc, const int n)
2247 {
2248 double phi0;
2249 double lambda0;
2250 double N0;
2251 double E0;
2252 double F0;
2253 double a;
2254 double b;
2255 int32 idx;
2256
2257 if(!GPS_Math_LatLon_To_UTM_Param(lat,lon,zone,zc,&lambda0,&E0,
2258 &N0,&F0))
2259 return 0;
2260
2261 phi0 = (double)0.0;
2262
2263 idx = GPS_Datum[n].ellipse;
2264 a = (double) GPS_Ellipse[idx].a;
2265 b = a - (a/GPS_Ellipse[idx].invf);
2266
2267 GPS_Math_LatLon_To_EN(E,N,lat,lon,N0,E0,phi0,lambda0,F0,a,b);
2268
2269 return 1;
2270 }
2271
2272 /* @func GPS_Math_UTM_EN_To_Known_Datum *********************************
2273 **
2274 ** Transform UTM zone, easting and northing to known datum lat/lon
2275 **
2276 ** @param [w] lat [double *] WGS84 latitude (deg)
2277 ** @param [r] lon [double *] WGS84 longitude (deg)
2278 ** @param [w] E [double] easting (metres)
2279 ** @param [w] N [double] northing (metres)
2280 ** @param [w] zone [int32] zone number
2281 ** @param [w] zc [char] zone character
2282 ** @param [r] n [int32] datum number from GPS_Datum structure
2283 **
2284 ** @return [int32] success
2285 ************************************************************************/
GPS_Math_UTM_EN_To_Known_Datum(double * lat,double * lon,double E,double N,int32 zone,char zc,const int n)2286 int32 GPS_Math_UTM_EN_To_Known_Datum(double *lat, double *lon, double E,
2287 double N, int32 zone, char zc, const int n)
2288 {
2289 double lambda0;
2290 double N0;
2291 double E0;
2292 double F0;
2293
2294 if (!GPS_Math_UTM_Param_To_Mc(zone,zc,&lambda0,&E0,&N0,&F0)) return 0;
2295
2296 GPS_Math_UTM_EN_to_LatLon(GPS_Datum[n].ellipse, N, E, lat, lon, lambda0, E0, N0);
2297 return 1;
2298 }
2299
2300 /* !!! copied from unused gpsproj.c !!! */
2301
2302 /* @func GPS_Math_Swiss_LatLon_To_EN ***********************************
2303 **
2304 ** Convert latitude and longitude to Swiss grid easting and northing
2305 **
2306 ** @param [r] phi [double] latitude (deg)
2307 ** @param [r] lambda [double] longitude (deg)
2308 ** @param [w] E [double *] easting (metre)
2309 ** @param [w] N [double *] northing (metre)
2310 ** @param [r] phi0 [double] latitude origin (deg) [normally 46.95240556]
2311 ** @param [r] lambda0 [double] longitude origin (deg) [normally 7.43958333]
2312 ** @param [r] E0 [double] false easting (metre) [normally 600000.0]
2313 ** @param [r] N0 [double] false northing (metre) [normally 200000.0]
2314 ** @param [r] a [double] semi-major axis [normally 6377397.000]
2315 ** @param [r] b [double] semi-minor axis [normally 6356078.823]
2316 **
2317 ** @return [void]
2318 ***************************************************************************/
GPS_Math_Swiss_LatLon_To_EN(double phi,double lambda,double * E,double * N,double phi0,double lambda0,double E0,double N0,double a,double b)2319 void GPS_Math_Swiss_LatLon_To_EN(double phi, double lambda, double *E,
2320 double *N,double phi0,double lambda0,
2321 double E0, double N0, double a, double b)
2322
2323 {
2324 double a2;
2325 double b2;
2326 double esq;
2327 double e;
2328 double c;
2329 double ephi0p;
2330 double phip;
2331 double sphip;
2332 double phid;
2333 double slambda2;
2334 double lambda1;
2335 double lambda2;
2336 double K;
2337 double po4;
2338 double w;
2339 double R;
2340
2341 lambda0 = GPS_Math_Deg_To_Rad(lambda0);
2342 phi0 = GPS_Math_Deg_To_Rad(phi0);
2343 lambda = GPS_Math_Deg_To_Rad(lambda);
2344 phi = GPS_Math_Deg_To_Rad(phi);
2345
2346 po4=GPS_PI/(double)4.0;
2347
2348 a2 = a*a;
2349 b2 = b*b;
2350 esq = (a2-b2)/a2;
2351 e = pow(esq,(double)0.5);
2352
2353 c = sqrt(1+((esq*pow(cos(phi0),(double)4.))/((double)1.-esq)));
2354
2355 ephi0p = asin(sin(phi0)/c);
2356
2357 K = log(tan(po4+ephi0p/(double)2.)) - c*(log(tan(po4+phi0/(double)2.)) -
2358 e/(double)2. * log(((double)1.+e*sin(phi0)) /
2359 ((double)1.-e*sin(phi0))));
2360 lambda1 = c*(lambda-lambda0);
2361 w = c*(log(tan(po4+phi/(double)2.)) - e/(double)2. *
2362 log(((double)1.+e*sin(phi)) / ((double)1.-e*sin(phi)))) + K;
2363
2364
2365 phip = (double)2. * (atan(exp(w)) - po4);
2366
2367 sphip = cos(ephi0p) * sin(phip) - sin(ephi0p) * cos(phip) * cos(lambda1);
2368 phid = asin(sphip);
2369
2370 slambda2 = cos(phip)*sin(lambda1) / cos(phid);
2371 lambda2 = asin(slambda2);
2372
2373 R = a*sqrt((double)1.-esq) / ((double)1.-esq*sin(phi0) * sin(phi0));
2374
2375 *N = R*log(tan(po4 + phid/(double)2.)) + N0;
2376 *E = R*lambda2 + E0;
2377 return;
2378 }
2379
2380 /* !!! copied from unused gpsproj.c !!! */
2381
2382 /* @func GPS_Math_Swiss_EN_To_LatLon ************************************
2383 **
2384 ** Convert Swiss Grid easting and northing to latitude and longitude
2385 **
2386 ** @param [r] E [double] easting (metre)
2387 ** @param [r] N [double] northing (metre)
2388 ** @param [w] phi [double *] latitude (deg)
2389 ** @param [w] lambda [double *] longitude (deg)
2390 ** @param [r] phi0 [double] latitude origin (deg) [normally 46.95240556]
2391 ** @param [r] lambda0 [double] longitude origin (deg) [normally 7.43958333]
2392 ** @param [r] E0 [double] false easting (metre) [normally 600000.0]
2393 ** @param [r] N0 [double] false northing (metre) [normally 200000.0]
2394 ** @param [r] a [double] semi-major axis [normally 6377397.000]
2395 ** @param [r] b [double] semi-minor axis [normally 6356078.823]
2396 **
2397 ** @return [void]
2398 *************************************************************************/
2399
GPS_Math_Swiss_EN_To_LatLon(double E,double N,double * phi,double * lambda,double phi0,double lambda0,double E0,double N0,double a,double b)2400 void GPS_Math_Swiss_EN_To_LatLon(double E, double N, double *phi,
2401 double *lambda, double phi0, double lambda0,
2402 double E0, double N0, double a, double b)
2403 {
2404 double a2;
2405 double b2;
2406 double esq;
2407 double e;
2408 double R;
2409 double c;
2410 double po4;
2411 double phid;
2412 double phi1;
2413 double lambdad;
2414 double lambda1;
2415 double slambda1;
2416 double ephi0p;
2417 double sphip;
2418 double tol;
2419 double cr;
2420 double C;
2421 double K;
2422
2423 lambda0 = GPS_Math_Deg_To_Rad(lambda0);
2424 phi0 = GPS_Math_Deg_To_Rad(phi0);
2425
2426 po4=GPS_PI/(double)4.0;
2427 tol=(double)0.00001;
2428
2429 a2 = a*a;
2430 b2 = b*b;
2431 esq = (a2-b2)/a2;
2432 e = pow(esq,(double)0.5);
2433
2434 R = a*sqrt((double)1.-esq) / ((double)1.-esq*sin(phi0) * sin(phi0));
2435
2436 phid = (double)2.*(atan(exp((N - N0)/R)) - po4);
2437 lambdad = (E - E0)/R;
2438
2439 c = sqrt((double)1.+((esq * pow(cos(phi0), (double)4.)) /
2440 ((double)1.-esq)));
2441 ephi0p = asin(sin(phi0) / c);
2442
2443 sphip = cos(ephi0p)*sin(phid) + sin(ephi0p)*cos(phid)*cos(lambdad);
2444 phi1 = asin(sphip);
2445
2446 slambda1 = cos(phid)*sin(lambdad)/cos(phi1);
2447 lambda1 = asin(slambda1);
2448
2449 *lambda = GPS_Math_Rad_To_Deg((lambda1/c + lambda0));
2450
2451 K = log(tan(po4 + ephi0p/(double)2.)) -c*(log(tan(po4 + phi0/(double)2.))
2452 - e/(double)2. * log(((double)1.+e*sin(phi0)) /
2453 ((double)1.-e*sin(phi0))));
2454 C = (K - log(tan(po4 + phi1/(double)2.)))/c;
2455
2456 do
2457 {
2458 cr = (C + log(tan(po4 + phi1/(double)2.)) - e/(double)2. *
2459 log(((double)1.+e*sin(phi1)) / ((double)1.-e*sin(phi1)))) *
2460 ((((double)1.-esq*sin(phi1)*sin(phi1)) * cos(phi1)) /
2461 ((double)1.-esq));
2462 phi1 -= cr;
2463 }
2464 while (fabs(cr) > tol);
2465
2466 *phi = GPS_Math_Rad_To_Deg(phi1);
2467
2468 return;
2469 }
2470
2471 /********************************************************************/
2472
GPS_Math_UTM_EN_to_LatLon(int ReferenceEllipsoid,const double UTMNorthing,const double UTMEasting,double * Lat,double * Lon,const double lambda0,const double E0,const double N0)2473 void GPS_Math_UTM_EN_to_LatLon(int ReferenceEllipsoid,
2474 const double UTMNorthing, const double UTMEasting,
2475 double *Lat, double *Lon,
2476 const double lambda0,
2477 const double E0,
2478 const double N0)
2479 {
2480 //converts UTM coords to lat/long. Equations from USGS Bulletin 1532
2481 //East Longitudes are positive, West longitudes are negative.
2482 //North latitudes are positive, South latitudes are negative
2483 //Lat and Long are in decimal degrees.
2484 //based on code witten by Chuck Gantz- chuck.gantz@globalstar.com
2485 //found at http://www.gpsy.com/gpsinfo/geotoutm/index.html
2486
2487 double k0 = 0.9996;
2488 double a, b;
2489 double eccSquared;
2490 double eccPrimeSquared;
2491 double e1;
2492 double N1, T1, C1, R1, D, M;
2493 double mu, phi1, phi1Rad;
2494 double x, y;
2495
2496 a = GPS_Ellipse[ReferenceEllipsoid].a;
2497 b = 1 / GPS_Ellipse[ReferenceEllipsoid].invf;
2498 eccSquared = b * (2.0 - b);
2499 e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared));
2500
2501 x = UTMEasting - E0; //remove false easting
2502 y = UTMNorthing - N0; //remove false northing
2503
2504 eccPrimeSquared = (eccSquared)/(1-eccSquared);
2505
2506 M = y / k0;
2507 mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256));
2508
2509 phi1Rad = mu+ (3*e1/2-27*e1*e1*e1/32)*sin(2*mu) +
2510 (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu) +
2511 (151*e1*e1*e1/96)*sin(6*mu);
2512 phi1 = GPS_Math_Rad_To_Deg(phi1Rad);
2513
2514 N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad));
2515 T1 = tan(phi1Rad)*tan(phi1Rad);
2516 C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad);
2517 R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5);
2518 D = x/(N1*k0);
2519
2520 *Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24
2521 +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720);
2522 *Lat = GPS_Math_Rad_To_Deg(*Lat);
2523
2524 *Lon = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)*D*D*D*D*D/120)/cos(phi1Rad);
2525 *Lon = lambda0 + GPS_Math_Rad_To_Deg(*Lon);
2526 }
2527
2528 /********************************************************************/
2529
GPS_Lookup_Datum_Index(const char * n)2530 int32 GPS_Lookup_Datum_Index(const char *n)
2531 {
2532 GPS_PDatum dp;
2533 GPS_PDatum_Alias al;
2534
2535 for (al = GPS_DatumAlias; al->alias; al++) {
2536 if (case_ignore_strcmp(al->alias, n) == 0) {
2537 return al->datum;
2538 }
2539 }
2540
2541 for (dp = GPS_Datum; dp->name; dp++) {
2542 if (0 == case_ignore_strcmp(dp->name, n)) {
2543 return dp - GPS_Datum;
2544 }
2545 }
2546
2547 return -1;
2548 }
2549
2550 const char *
GPS_Math_Get_Datum_Name(const int datum_index)2551 GPS_Math_Get_Datum_Name(const int datum_index)
2552 {
2553 return GPS_Datum[datum_index].name;
2554 }
2555
2556
2557