Geo Stars Library  0.9.3
Geodetic and Astrometry library
geoEfg2Llh.c
Go to the documentation of this file.
1 /*! \file geoEfg2Llh.c
2  \brief This file contains the geo library routines for
3  conversion between Geocentric and Geodetic coordinates
4 
5  Methods avalailable:
6  \li geoStarLib method
7  \li Hirvonen & Moritz method
8  \li Torge method
9  \li Astronomical Almanac 2002 method
10  \li Borkowski method
11  \li Bowring method
12 
13 */
14 #include <stdio.h>
15 #include <math.h>
16 #include <string.h>
17 #include <stdlib.h>
18 #include <ctype.h>
19 #include "geoStars.h"
20 
21 double geoAccuracy = GEO_EFG2LLH_ACCURACY;
22 
23 static int its; // keep track of the number of iterations
24 
25 int DLL_CALLCONV geoEfg2LlhGetIts (void)
26 {
27  return(its);
28 }
29 
30 //-----------------------------------------------------------------
31 /*! \brief This routine will set the accuracy of the iterative
32 Efg2Llh routines (geoEfg2Llh_hm, geoEfg2Llh_torge, geoEfg2Llh_bowring,
33  geoEfg2Llh_aa)
34 
35  \param double acc : Accuracy in degrees
36 \retval GEO_OK on success
37 \retval GEO_ERROR on error
38 
39  Several definitions may be used with this function:
40  GEO_EFG2LLH_ACCURACY_METER = 0.00001 degrees
41  GEO_EFG2LLH_ACCURACY_CM = 0.0000001 degrees
42  GEO_EFG2LLH_ACCURACY_MM = 0.00000001 degrees
43 
44 */
45 int DLL_CALLCONV geoSetAccuracy(double acc)
46 {
47  if (acc <= 1.0) // make sure the accuracy is sane
48  {
49  geoAccuracy = acc;
50  return(GEO_OK);
51  }
52  return(GEO_ERROR);
53 }
54 
55 //-----------------------------------------------------------------
56 /*! \brief This routine will convert earth centered Cartesian
57 coordinates (E,F,G), into geodetic coordinates (latitude
58 \f$\phi\f$, longitude \f$\lambda\f$, and ellipsoid height \f$h\f$).
59 
60  \param int datum
61  \param double efg[] : EFG(xyz) in METERS
62  \param double *lat : Latitude in RADIANS
63  \param double *lon : Longitude in RADIANS
64  \param double *hgt : Height in METERS
65  \return nothing
66 
67 \note This routine will be exact only for WGS84 coordinates. All other
68 datums will be slightly off.
69 
70 */
71 
72 void DLL_CALLCONV geoEfg2Llh (int datum, double efg[],
73 double *lat, double *lon, double *hgt)
74 {
75  double p, u, u_prime, a, flat, b, e2, ee2,sign = 1.0,cu_prime,su,cu,t1;
76 
77  /* Get the ellipsoid parameters */
78  geoGetEllipsoid (&a, &b, &e2, &ee2, &flat, datum);
79  its=0;
80 
81  // Computes EFG to lat, lon & height
82  p = sqrt(sqr(efg[GEO_E]) + sqr(efg[GEO_F]));
83  u = atan((efg[GEO_G]/p) * (a/b));
84  // u = atan2(efg[GEO_G] * a , p * b);
85  su = sin(u);
86  cu = cos(u);
87 
88  *lat = atan((efg[GEO_G] + ee2 * b * cube(su) ) /
89  ( p - e2 * a * cube(cu) ) );
90 
91  u_prime = atan((1.0 - flat) * tan(*lat));
92  cu_prime = p - a * cos(u_prime);
93  if(cu_prime < 0.0) sign= -1.0; // determine sign
94 
95  // *hgt = p / cos(*lat) - ( a / (sqrt(1.0 - e2 * pow(sin(*lat),2.0)))); //same results
96  t1 = efg[GEO_G] - b * sin(u_prime);
97  *hgt = sign * sqrt( sqr(cu_prime) + sqr(t1) );
98  *lon = atan2(efg[GEO_F],efg[GEO_E]); // atan(f/e)
99 
100 }
101 
102 double DLL_CALLCONV geoEfg2Lat (int datum, double e, double f, double g)
103 {
104  double efg[3];
105  double lat, lon, hgt;
106  efg[GEO_E] = e;
107  efg[GEO_F] = f;
108  efg[GEO_G] = g;
109 
110  geoEfg2Llh (datum,efg,&lat,&lon,&hgt);
111  return(RAD_TO_DEG * lat);
112 }
113 
114 double DLL_CALLCONV geoEfg2Lon (int datum, double e, double f, double g)
115 {
116  double efg[3];
117  double lat, lon, hgt;
118  efg[GEO_E] = e;
119  efg[GEO_F] = f;
120  efg[GEO_G] = g;
121 
122  geoEfg2Llh (datum,efg,&lat,&lon,&hgt);
123  return(RAD_TO_DEG * lon);
124 }
125 double DLL_CALLCONV geoEfg2Hgt (int datum, double e, double f, double g)
126 {
127  double efg[3];
128  double lat, lon, hgt;
129  efg[GEO_E] = e;
130  efg[GEO_F] = f;
131  efg[GEO_G] = g;
132 
133  geoEfg2Llh (datum,efg,&lat,&lon,&hgt);
134  return(hgt);
135 }
136 
137 
138 /*! \brief Hirvonen & Moritz Method - This routine will convert earth centered Cartesian
139 coordinates (E,F,G), into geodetic coordinates (latitude
140 \f$\phi\f$, longitude \f$\lambda\f$, and ellipsoid height \f$h\f$).
141 
142  \param int datum
143  \param double efg[] : EFG(xyz) in METERS
144  \param double *lat : Latitude in RADIANS
145  \param double *lon : Longitude in RADIANS
146  \param double *hgt : Height in METERS
147  \return nothing
148 
149 \note This method is a successive approximation
150 
151 */
152 
153 void DLL_CALLCONV geoEfg2Llh_hm (int datum, double efg[],
154 double *lat, double *lon, double *hgt)
155 {
156  double a, flat, b, e2, ee2;
157  double lat1,tmp,n=0.0,delta,p,slat;
158  // int done=1;
159  int i;
160 
161  if(efg[GEO_G] == 0.0) efg[GEO_G] = 0.00000001; //cheat to avoid bad numbers
162 
163  /* Get the ellipsoid parameters */
164  geoGetEllipsoid (&a, &b, &e2, &ee2, &flat, datum);
165  p = sqrt(sqr(efg[GEO_E])+ sqr(efg[GEO_F]));
166  *lon = atan2(efg[GEO_F],efg[GEO_E]); // atan(f/e)
167  its=0;
168 
169  // Prepare the initial value of lat
170  tmp = efg[GEO_G] / p;
171  lat1 = atan((1.0 / (1.0 - e2))*tmp);
172  *lat=lat1;
173 
174  // iterate until delta lat is tiny
175  for(i=0;i<GEO_EFG2LLH_MAX_ITS;i++)
176  {
177  its++;
178  slat = sin(lat1);
179  n= a/sqrt(1.0 - e2 * sqr(slat));
180  *lat=atan(tmp*(1.0 + (e2*n*slat)/efg[GEO_G]));
181  delta = *lat - lat1;
182  if(fabs(delta) <= geoAccuracy) break;
183  // if(fabs(delta) == 0.0) break;
184  lat1=*lat;
185  }
186  *hgt = (p / cos(*lat)) - n;
187 
188 }
189 
190 double DLL_CALLCONV geoEfg2Lat_hm (int datum, double e, double f, double g)
191 {
192  double efg[3];
193  double lat, lon, hgt;
194  efg[GEO_E] = e;
195  efg[GEO_F] = f;
196  efg[GEO_G] = g;
197 
198  geoEfg2Llh_hm (datum,efg,&lat,&lon,&hgt);
199  return(RAD_TO_DEG * lat);
200 }
201 
202 double DLL_CALLCONV geoEfg2Lon_hm (int datum, double e, double f, double g)
203 {
204  double efg[3];
205  double lat, lon, hgt;
206  efg[GEO_E] = e;
207  efg[GEO_F] = f;
208  efg[GEO_G] = g;
209 
210  geoEfg2Llh_hm (datum,efg,&lat,&lon,&hgt);
211  return(RAD_TO_DEG * lon);
212 }
213 double DLL_CALLCONV geoEfg2Hgt_hm (int datum, double e, double f, double g)
214 {
215  double efg[3];
216  double lat, lon, hgt;
217  efg[GEO_E] = e;
218  efg[GEO_F] = f;
219  efg[GEO_G] = g;
220 
221  geoEfg2Llh_hm (datum,efg,&lat,&lon,&hgt);
222  return(hgt);
223 }
224 
225 double DLL_CALLCONV geoEfg2Hgt_hm_its (int datum, double e, double f, double g)
226 {
227  double efg[3];
228  double lat, lon, hgt;
229  efg[GEO_E] = e;
230  efg[GEO_F] = f;
231  efg[GEO_G] = g;
232 
233  geoEfg2Llh_hm (datum,efg,&lat,&lon,&hgt);
234 
235  return((double)geoEfg2LlhGetIts());
236 }
237 
238 /*! \brief Torge Method - This routine will convert earth centered Cartesian
239 coordinates (E,F,G), into geodetic coordinates (latitude
240 \f$\phi\f$, longitude \f$\lambda\f$, and ellipsoid height \f$h\f$).
241 
242  \param int datum
243  \param double efg[] : EFG(xyz) in METERS
244  \param double *lat : Latitude in RADIANS
245  \param double *lon : Longitude in RADIANS
246  \param double *hgt : Height in METERS
247  \return nothing
248 
249 \note This method is a successive approximation
250 
251 */
252 
253 void DLL_CALLCONV geoEfg2Llh_torge (int datum, double efg[],
254 double *lat, double *lon, double *hgt)
255 {
256  double a, flat, b, e2, ee2;
257  double lat1,hgt1,tmp,n=0.0,deltal,deltah,p,slat;
258  //,clat;
259  // int done=1;
260  int i;
261 
262  /* Get the ellipsoid parameters */
263  geoGetEllipsoid (&a, &b, &e2, &ee2, &flat, datum);
264  p = sqrt(sqr(efg[GEO_E])+ sqr(efg[GEO_F]));
265  *lon = atan2(efg[GEO_F],efg[GEO_E]); // atan(f/e)
266  its=0;
267 
268  // Prepare the initial value of lat
269  tmp = efg[GEO_G] / p;
270  // lat1 = atan(tmp/(1.0 - e2)); //assume hgt = 0 to start
271  lat1 = atan(tmp); //assume hgt = 0 to start
272  *lat=lat1;
273  *hgt=0.0;
274 
275 
276  // iterate until delta lat is tiny
277  for(i=0;i<GEO_EFG2LLH_MAX_ITS;i++)
278  {
279  its++;
280  lat1=*lat;
281  hgt1 = *hgt;
282 
283  slat = sin(lat1);
284  //clat = cos(lat1);
285 
286  // radius of curvature estimate
287  n= a/sqrt(1.0 - e2 * sqr(slat));
288 
289  // Next: estimates of lat/hgt
290  *lat=atan2(tmp , (1.0 - (e2 * ( n / (n + *hgt)) ) ) );
291  *hgt = (p / cos(*lat)) - n;
292  deltal = *lat - lat1;
293  deltah = *hgt - hgt1;
294  if((fabs(deltah) <= .00001) && (fabs(deltal) <= geoAccuracy)) break;
295  }
296 
297 }
298 
299 double DLL_CALLCONV geoEfg2Lat_torge (int datum, double e, double f, double g)
300 {
301  double efg[3];
302  double lat, lon, hgt;
303  efg[GEO_E] = e;
304  efg[GEO_F] = f;
305  efg[GEO_G] = g;
306 
307  geoEfg2Llh_torge (datum,efg,&lat,&lon,&hgt);
308  return(RAD_TO_DEG * lat);
309 }
310 
311 double DLL_CALLCONV geoEfg2Lon_torge (int datum, double e, double f, double g)
312 {
313  double efg[3];
314  double lat, lon, hgt;
315  efg[GEO_E] = e;
316  efg[GEO_F] = f;
317  efg[GEO_G] = g;
318 
319  geoEfg2Llh_torge (datum,efg,&lat,&lon,&hgt);
320  return(RAD_TO_DEG * lon);
321 }
322 double DLL_CALLCONV geoEfg2Hgt_torge (int datum, double e, double f, double g)
323 {
324  double efg[3];
325  double lat, lon, hgt;
326  efg[GEO_E] = e;
327  efg[GEO_F] = f;
328  efg[GEO_G] = g;
329 
330  geoEfg2Llh_torge (datum,efg,&lat,&lon,&hgt);
331  return(hgt);
332 }
333 
334 double DLL_CALLCONV geoEfg2Hgt_torge_its (int datum, double e, double f, double g)
335 {
336  double efg[3];
337  double lat, lon, hgt;
338  efg[GEO_E] = e;
339  efg[GEO_F] = f;
340  efg[GEO_G] = g;
341 
342  geoEfg2Llh_torge (datum,efg,&lat,&lon,&hgt);
343 
344  return((double)geoEfg2LlhGetIts());
345 }
346 
347 //-----------------------------------------------------------------
348 /*! \brief Bowring Method - This routine will convert earth centered Cartesian
349 coordinates (E,F,G), into geodetic coordinates (latitude
350 \f$\phi\f$, longitude \f$\lambda\f$, and ellipsoid height \f$h\f$).
351 
352  \param int datum
353  \param double efg[] : EFG(xyz) in METERS
354  \param double *lat : Latitude in RADIANS
355  \param double *lon : Longitude in RADIANS
356  \param double *hgt : Height in METERS
357  \return nothing
358 
359 \note This method is a successive approximation
360 
361 */
362 void DLL_CALLCONV geoEfg2Llh_bowring (int datum, double efg[],
363 double *lat, double *lon, double *hgt)
364 {
365  double a, flat, b, e2, ee2;
366  double tmp,n,b0,b1,p,slat,clat,ee2b,ae2,boa;
367  // int done=1;
368  int i;
369 
370  /* Get the ellipsoid parameters */
371  geoGetEllipsoid (&a, &b, &e2, &ee2, &flat, datum);
372  *lon = atan2(efg[GEO_F],efg[GEO_E]); // atan(f/e)
373  p = sqrt(sqr(efg[GEO_E])+ sqr(efg[GEO_F]));
374  tmp = efg[GEO_G] / p;
375  its=0;
376 
377  // temps pulled out of the loop
378  ae2 = a * e2;
379  ee2b = ee2 * b;
380  boa = b/a;
381 
382  // Prepare the initial value of lat
383  b0 = atan((a/b)*tmp);
384 
385  // iterate until delta lat is tiny
386  for(i=0;i<GEO_EFG2LLH_MAX_ITS;i++)
387  {
388  its++;
389  b1=b0;
390  slat = sin(b1);
391  clat = cos(b1);
392  *lat = atan((efg[GEO_G]+ ee2b * cube(slat) ) /
393  (p - ae2 * cube(clat) ) );
394  b0 = atan((boa) * tan(*lat));
395  if(fabs(b0-b1) <= geoAccuracy) break;
396  // if(fabs(b1-b0)== 0.0) break;
397  // if(b1-b0== 0.0) break;
398  }
399  slat = sin(*lat);
400  n= a/sqrt(1.0-e2*sqr(slat));
401 
402  *hgt = (p / cos(*lat)) - n;
403 
404 }
405 
406 double DLL_CALLCONV geoEfg2Lat_bowring (int datum, double e, double f, double g)
407 {
408  double efg[3];
409  double lat, lon, hgt;
410  efg[GEO_E] = e;
411  efg[GEO_F] = f;
412  efg[GEO_G] = g;
413 
414  geoEfg2Llh_bowring (datum,efg,&lat,&lon,&hgt);
415  return(RAD_TO_DEG * lat);
416 }
417 
418 double DLL_CALLCONV geoEfg2Lon_bowring (int datum, double e, double f, double g)
419 {
420  double efg[3];
421  double lat, lon, hgt;
422  efg[GEO_E] = e;
423  efg[GEO_F] = f;
424  efg[GEO_G] = g;
425 
426  geoEfg2Llh_bowring (datum,efg,&lat,&lon,&hgt);
427  return(RAD_TO_DEG * lon);
428 }
429 double DLL_CALLCONV geoEfg2Hgt_bowring (int datum, double e, double f, double g)
430 {
431  double efg[3];
432  double lat, lon, hgt;
433  efg[GEO_E] = e;
434  efg[GEO_F] = f;
435  efg[GEO_G] = g;
436 
437  geoEfg2Llh_bowring (datum,efg,&lat,&lon,&hgt);
438  return(hgt);
439 }
440 
441 double DLL_CALLCONV geoEfg2Hgt_bowring_its (int datum, double e, double f, double g)
442 {
443  double efg[3];
444  double lat, lon, hgt;
445  efg[GEO_E] = e;
446  efg[GEO_F] = f;
447  efg[GEO_G] = g;
448 
449  geoEfg2Llh_bowring (datum,efg,&lat,&lon,&hgt);
450  return((double)geoEfg2LlhGetIts());
451 }
452 
453 //-----------------------------------------------------------------
454 /*! \brief Astronomical Almanac 2002 Method - This routine will convert earth centered Cartesian
455 coordinates (E,F,G), into geodetic coordinates (latitude
456 \f$\phi\f$, longitude \f$\lambda\f$, and ellipsoid height \f$h\f$).
457 
458  \param int datum
459  \param double efg[] : EFG(xyz) in METERS
460  \param double *lat : Latitude in RADIANS
461  \param double *lon : Longitude in RADIANS
462  \param double *hgt : Height in METERS
463  \return nothing
464 
465 \note This method is a successive approximation from Astronomical Almanac 2002 K11-12
466 
467 */
468 void DLL_CALLCONV geoEfg2Llh_aa (int datum, double efg[],
469 double *lat, double *lon, double *hgt)
470 {
471  double phi,p,n=1.0,phi1;
472  int i;
473 
474  double a, flat, b, e2, ee2;
475  double slat;
476 
477  /* Get the ellipsoid parameters */
478  geoGetEllipsoid (&a, &b, &e2, &ee2, &flat, datum);
479 
480  // Computes EFG to lat, lon & height
481  // from the Astronomical Almanac 2002 - K11-12
482  *lon = atan2(efg[GEO_F],efg[GEO_E]); // atan(f/e)
483  p = sqrt(sqr(efg[GEO_E]) + sqr(efg[GEO_F]));
484 
485  phi = atan(efg[GEO_G]/p); // first approximation
486  // phi = atan(efg[GEO_G]/((1.0-e2)*p)); // faster
487 
488  its=0;
489 
490  // compute until the error in phi is no more
491  for(i=0;i<GEO_EFG2LLH_MAX_ITS;i++)
492  {
493  its++;
494  phi1 = phi;
495  slat = sin(phi1);
496  // c = 1.0/sqrt(1.0 - e2 * sqr(slat));
497  // phi = atan((efg[GEO_G] + a * c * e2 * slat)/r);
498  n = a/sqrt(1.0 - e2 * sqr(slat));
499  phi = atan((efg[GEO_G] + n * e2 * slat)/p);
500  if( fabs(phi1-phi) <= geoAccuracy) break;
501  }
502  *lat = phi;
503  // *hgt = r / cos(phi) - a * c;
504  *hgt = p / cos(phi) - n;
505 
506 }
507 
508 double DLL_CALLCONV geoEfg2Lat_aa (int datum, double e, double f, double g)
509 {
510  double efg[3];
511  double lat, lon, hgt;
512  efg[GEO_E] = e;
513  efg[GEO_F] = f;
514  efg[GEO_G] = g;
515 
516  geoEfg2Llh_aa (datum,efg,&lat,&lon,&hgt);
517  return(RAD_TO_DEG * lat);
518 }
519 
520 double DLL_CALLCONV geoEfg2Lon_aa (int datum, double e, double f, double g)
521 {
522  double efg[3];
523  double lat, lon, hgt;
524  efg[GEO_E] = e;
525  efg[GEO_F] = f;
526  efg[GEO_G] = g;
527 
528  geoEfg2Llh_aa (datum,efg,&lat,&lon,&hgt);
529  return(RAD_TO_DEG * lon);
530 }
531 double DLL_CALLCONV geoEfg2Hgt_aa (int datum, double e, double f, double g)
532 {
533  double efg[3];
534  double lat, lon, hgt;
535  efg[GEO_E] = e;
536  efg[GEO_F] = f;
537  efg[GEO_G] = g;
538 
539  geoEfg2Llh_aa (datum,efg,&lat,&lon,&hgt);
540  return(hgt);
541 }
542 
543 double DLL_CALLCONV geoEfg2Hgt_aa_its (int datum, double e, double f, double g)
544 {
545  double efg[3];
546  double lat, lon, hgt;
547  efg[GEO_E] = e;
548  efg[GEO_F] = f;
549  efg[GEO_G] = g;
550 
551  geoEfg2Llh_aa (datum,efg,&lat,&lon,&hgt);
552  return((double)geoEfg2LlhGetIts());
553 }
554 
555 
556 //-----------------------------------------------------------------
557 /*! \brief Borkowski Method - This routine will convert earth centered Cartesian
558 coordinates (E,F,G), into geodetic coordinates (latitude
559 \f$\phi\f$, longitude \f$\lambda\f$, and ellipsoid height \f$h\f$).
560 
561  \param int datum
562  \param double efg[] : EFG(xyz) in METERS
563  \param double *lat : Latitude in RADIANS
564  \param double *lon : Longitude in RADIANS
565  \param double *hgt : Height in METERS
566  \return nothing
567 
568 \note This method is closed exact method from Bull. Geod., 63 (1989), pp. 50-56
569 
570 */
571 void DLL_CALLCONV geoEfg2Llh_borkowski (int datum, double efg[],
572 double *lat, double *lon, double *hgt)
573 {
574  double a, flat, b, e2, ee2;
575  double r,e,f,p,q,d,v,g,t,sign,bg,samb,sqrtd,sqre;
576 
577 
578  /* Get the ellipsoid parameters */
579  geoGetEllipsoid (&a, &b, &e2, &ee2, &flat, datum);
580  its=0;
581 
582  // Computes EFG to lat, lon & height
583  // from Bull. Geod., 63 (1989), pp. 50-56
584  // Accurate Algorithms to Transform Geocentric to Geodetic Coordinates
585  // By K.M. Borkowski
586 
587  *lon = atan2(efg[GEO_F],efg[GEO_E]); // atan(f/e)
588 
589  //use the sign for Z
590  if(efg[GEO_G] < 0.0) sign = -1.0;
591  else sign = 1.0;
592  b = fabs(b) * sign;
593  bg = b * efg[GEO_G]; //temp
594  samb = sqr(a)-sqr(b) ;//temp
595  r = efg[GEO_E] / cos(*lon);
596  e = (bg - samb) / (a*r);
597  sqre = sqr(e);
598  f = (bg + samb) / (a*r);
599  p = 4.0/3.0*(e*f+1.0);
600  q = 2.0 * (sqre-sqr(f));
601  d = cube(p) + sqr(q);
602  sqrtd = sqrt(d);
603  if(d < 0.0)
604  v = 2.0 * sqrt(-p) * cos(acos(q/pow(-p,3.0/2.0))/3.0);
605  else
606  v = pow((sqrtd - q),1.0/3.0) - pow((sqrtd + q),1.0/3.0);
607  g = (sqrt(sqre+v) + e ) / 2.0;
608  t = sqrt(sqr(g)+((f-v*g)/(2.0*g-e))) - g;
609  *lat = atan(a*(1.0-sqr(t))/(2.0*b*t));
610 
611  *hgt = (r - a*t) * cos(*lat) + (efg[GEO_G]-b)*sin(*lat);
612 
613 }
614 
615 double DLL_CALLCONV geoEfg2Lat_borkowski (int datum, double e, double f, double g)
616 {
617  double efg[3];
618  double lat, lon, hgt;
619  efg[GEO_E] = e;
620  efg[GEO_F] = f;
621  efg[GEO_G] = g;
622 
623  geoEfg2Llh_borkowski (datum,efg,&lat,&lon,&hgt);
624  return(RAD_TO_DEG * lat);
625 }
626 
627 double DLL_CALLCONV geoEfg2Lon_borkowski (int datum, double e, double f, double g)
628 {
629  double efg[3];
630  double lat, lon, hgt;
631  efg[GEO_E] = e;
632  efg[GEO_F] = f;
633  efg[GEO_G] = g;
634 
635  geoEfg2Llh_borkowski (datum,efg,&lat,&lon,&hgt);
636  return(RAD_TO_DEG * lon);
637 }
638 double DLL_CALLCONV geoEfg2Hgt_borkowski (int datum, double e, double f, double g)
639 {
640  double efg[3];
641  double lat, lon, hgt;
642  efg[GEO_E] = e;
643  efg[GEO_F] = f;
644  efg[GEO_G] = g;
645 
646  geoEfg2Llh_borkowski (datum,efg,&lat,&lon,&hgt);
647  return(hgt);
648 }
649 
650 
651 //-----------------------------------------------------------------
652 /*! \brief Vermeille Method - This routine will convert earth centered Cartesian
653 coordinates (E,F,G), into geodetic coordinates (latitude
654 \f$\phi\f$, longitude \f$\lambda\f$, and ellipsoid height \f$h\f$).
655 
656  \param int datum
657  \param double efg[] : EFG(xyz) in METERS
658  \param double *lat : Latitude in RADIANS
659  \param double *lon : Longitude in RADIANS
660  \param double *hgt : Height in METERS
661  \return nothing
662 
663 \note This method is closed exact method from ....
664 
665 */
666 void DLL_CALLCONV geoEfg2Llh_vermeille (int datum, double efg[],
667 double *lat, double *lon, double *hgt)
668 {
669  double a, flat, b, e2, ee2;
670  double a2,x2,y2,z2,d2,e4,p,q,r,s,t,u,v,w,k,d,sqrtx2py2,sqrtd2pz2;
671 
672 
673  /* Get the ellipsoid parameters */
674  geoGetEllipsoid (&a, &b, &e2, &ee2, &flat, datum);
675  its=0;
676 
677  // precompute variables
678  x2 = sqr(efg[GEO_E]);
679  y2 = sqr(efg[GEO_F]);
680  z2 = sqr(efg[GEO_G]);
681  a2 = sqr(a);
682  e4 = sqr(e2);
683  sqrtx2py2 = sqrt(x2+y2);
684 
685  // main algorithm
686  p = (x2 + y2)/ a2;
687  q = ((1.0-e2)/ a2)*z2;
688  r = (p+q-e4)/6.0;
689  s = e4 * (p*q)/(4.0*pow(r,3.0));
690  t = pow(1.0+s+sqrt(s*(2.0+s)),1.0/3.0);
691  u = r*(1.0+t+(1.0/t));
692  v = sqrt(sqr(u)+e4*q);
693  w = e2*(u+v-q)/(2.0*v);
694  k = sqrt(u+v+sqr(w)) - w;
695  d = k*sqrtx2py2/(k+e2);
696  d2 = sqr(d);
697  sqrtd2pz2 = sqrt(d2+z2);
698 
699  *lon = 2 * atan(efg[GEO_F]/(efg[GEO_E] + sqrtx2py2));
700 
701  *lat = 2 * atan(efg[GEO_G]/(d + sqrtd2pz2));
702 
703  *hgt = ((k + e2 - 1.0) / k) * sqrtd2pz2;
704 
705 }
706 
707 double DLL_CALLCONV geoEfg2Lat_vermeille (int datum, double e, double f, double g)
708 {
709  double efg[3];
710  double lat, lon, hgt;
711  efg[GEO_E] = e;
712  efg[GEO_F] = f;
713  efg[GEO_G] = g;
714 
715  geoEfg2Llh_vermeille (datum,efg,&lat,&lon,&hgt);
716  return(RAD_TO_DEG * lat);
717 }
718 
719 double DLL_CALLCONV geoEfg2Lon_vermeille (int datum, double e, double f, double g)
720 {
721  double efg[3];
722  double lat, lon, hgt;
723  efg[GEO_E] = e;
724  efg[GEO_F] = f;
725  efg[GEO_G] = g;
726 
727  geoEfg2Llh_vermeille (datum,efg,&lat,&lon,&hgt);
728  return(RAD_TO_DEG * lon);
729 }
730 double DLL_CALLCONV geoEfg2Hgt_vermeille (int datum, double e, double f, double g)
731 {
732  double efg[3];
733  double lat, lon, hgt;
734  efg[GEO_E] = e;
735  efg[GEO_F] = f;
736  efg[GEO_G] = g;
737 
738  geoEfg2Llh_vermeille (datum,efg,&lat,&lon,&hgt);
739  return(hgt);
740 }
741 
742 
743 //-----------------------------------------------------------------
744 /*! \brief heikkinen Method - This routine will convert earth centered Cartesian
745 coordinates (E,F,G), into geodetic coordinates (latitude
746 \f$\phi\f$, longitude \f$\lambda\f$, and ellipsoid height \f$h\f$).
747 
748  \param int datum
749  \param double efg[] : EFG(xyz) in METERS
750  \param double *lat : Latitude in RADIANS
751  \param double *lon : Longitude in RADIANS
752  \param double *hgt : Height in METERS
753  \return nothing
754 
755 \note This method is closed exact method from ....
756 
757 */
758 void DLL_CALLCONV geoEfg2Llh_heikkinen (int datum, double efg[],
759 double *lat, double *lon, double *hgt)
760 {
761  double a, flat, b, e2, ee2;
762  double a2,b2,x2,y2,z2,d,d2,e4,m1e2,m1e2z2,
763  W,W2,F,G,p,q,q1,s,sp,W0,wew0,U,V,Z0;
764 
765 
766  /* Get the ellipsoid parameters */
767  geoGetEllipsoid (&a, &b, &e2, &ee2, &flat, datum);
768  its=0;
769 
770  // precompute variables
771  x2 = sqr(efg[GEO_E]);
772  y2 = sqr(efg[GEO_F]);
773  z2 = sqr(efg[GEO_G]);
774  a2 = sqr(a);
775  b2 = sqr(b);
776  e4 = sqr(e2);
777  m1e2 = 1.0 - e2;
778  m1e2z2 = m1e2*z2;
779 
780  // main algorithm
781  W = sqrt(x2+y2);
782  W2= W*W;
783 
784  F = 54.0*b2*z2;
785 
786  G = W2 + m1e2z2 - e2*(a2-b2);
787  d = e4*F*W2/(G*G*G);
788  d2=d*d;
789 
790  s = pow(1.0+d+sqrt(d2+2*d),1.0/3.0);
791  sp=s+1.0/s+1.0;
792  p = F/(3.0*sqr(sp)*G*G);
793  q = sqrt(1.0+2.0*e4*p);
794  q1 = q+1.0;
795 
796  W0= (-p*e2*W / q1) + sqrt(0.5*a2*(1.0+1.0/q)-(p*m1e2z2/(q*q1))-(p*W2/2.0));
797  wew0 = W-e2*W0;
798  U = sqrt(sqr(wew0)+z2);
799  V = sqrt(sqr(wew0)+m1e2z2);
800  Z0= b2*efg[GEO_G]/(a*V);
801  *hgt = U*(1.0-(b2/(a*V)));
802  *lat = atan2((efg[GEO_G]+ee2*Z0),W);
803  *lon = atan2(efg[GEO_F],efg[GEO_E]);
804 
805 
806 }
807 
808 double DLL_CALLCONV geoEfg2Lat_heikkinen (int datum, double e, double f, double g)
809 {
810  double efg[3];
811  double lat, lon, hgt;
812  efg[GEO_E] = e;
813  efg[GEO_F] = f;
814  efg[GEO_G] = g;
815 
816  geoEfg2Llh_heikkinen (datum,efg,&lat,&lon,&hgt);
817  return(RAD_TO_DEG * lat);
818 }
819 
820 double DLL_CALLCONV geoEfg2Lon_heikkinen (int datum, double e, double f, double g)
821 {
822  double efg[3];
823  double lat, lon, hgt;
824  efg[GEO_E] = e;
825  efg[GEO_F] = f;
826  efg[GEO_G] = g;
827 
828  geoEfg2Llh_heikkinen (datum,efg,&lat,&lon,&hgt);
829  return(RAD_TO_DEG * lon);
830 }
831 double DLL_CALLCONV geoEfg2Hgt_heikkinen (int datum, double e, double f, double g)
832 {
833  double efg[3];
834  double lat, lon, hgt;
835  efg[GEO_E] = e;
836  efg[GEO_F] = f;
837  efg[GEO_G] = g;
838 
839  geoEfg2Llh_heikkinen (datum,efg,&lat,&lon,&hgt);
840  return(hgt);
841 }
842 
843 //-----------------------------------------------------------------
844 /*! \brief toms Method - This routine will convert earth centered Cartesian
845 coordinates (E,F,G), into geodetic coordinates (latitude
846 \f$\phi\f$, longitude \f$\lambda\f$, and ellipsoid height \f$h\f$).
847 
848  \param int datum
849  \param double efg[] : EFG(xyz) in METERS
850  \param double *lat : Latitude in RADIANS
851  \param double *lon : Longitude in RADIANS
852  \param double *hgt : Height in METERS
853  \return nothing
854 
855 \note This method is closed exact method from ....
856 
857 */
858 #define aDc 1.0026000
859 
860 void DLL_CALLCONV geoEfg2Llh_toms (int datum, double efg[],
861 double *lat, double *lon, double *hgt)
862 {
863  double a, flat, b, e2, ee2;
864  double x2,y2,W,W2,t0,s0,sbo,cbo,t1,wae,s1,sin1,cos1,n;
865 
866 
867  /* Get the ellipsoid parameters */
868  geoGetEllipsoid (&a, &b, &e2, &ee2, &flat, datum);
869  its=0;
870 
871  // precompute variables
872  x2 = sqr(efg[GEO_E]);
873  y2 = sqr(efg[GEO_F]);
874  // z2 = sqr(efg[GEO_G]);
875 
876  // main algorithm
877  W = sqrt(x2+y2);
878  W2= W*W;
879 
880  t0 = efg[GEO_G] * aDc;
881  s0 = sqrt(sqr(t0)+W2);
882  sbo= t0/s0;
883  cbo= W/s0;
884 
885  t1= efg[GEO_G] + b * ee2 * cube(sbo);
886  wae = W - a * e2 * cube(cbo);
887  s1 = sqrt(sqr(t1) + sqr(wae));
888  sin1 = t1/s1;
889  cos1 = wae/s1;
890  n = a/sqrt(1.0-e2*sqr(sin1));
891  if(sin1 >= 0.3826834323650897717284599840304)
892  {
893  sin1 = sqrt(sqr(sin1));
894  *hgt = efg[GEO_G]/sin1 + n*(e2-1.0);
895  }
896  else
897  {
898  *hgt = W/cos1 - n;
899  }
900 
901  *lat = atan2(sin1,cos1);
902  *lon = atan2(efg[GEO_F],efg[GEO_E]);
903 
904 
905 }
906 
907 double DLL_CALLCONV geoEfg2Lat_toms (int datum, double e, double f, double g)
908 {
909  double lat, lon, hgt;
910  double efg[] = {e,f,g};
911 // efg[GEO_E] = e;
912 // efg[GEO_F] = f;
913 // efg[GEO_G] = g;
914 
915  geoEfg2Llh_toms (datum,efg,&lat,&lon,&hgt);
916  return(RAD_TO_DEG * lat);
917 }
918 
919 double DLL_CALLCONV geoEfg2Lon_toms (int datum, double e, double f, double g)
920 {
921  double efg[3];
922  double lat, lon, hgt;
923  efg[GEO_E] = e;
924  efg[GEO_F] = f;
925  efg[GEO_G] = g;
926 
927  geoEfg2Llh_toms (datum,efg,&lat,&lon,&hgt);
928  return(RAD_TO_DEG * lon);
929 }
930 double DLL_CALLCONV geoEfg2Hgt_toms (int datum, double e, double f, double g)
931 {
932  double efg[3];
933  double lat, lon, hgt;
934  efg[GEO_E] = e;
935  efg[GEO_F] = f;
936  efg[GEO_G] = g;
937 
938  geoEfg2Llh_toms (datum,efg,&lat,&lon,&hgt);
939  return(hgt);
940 }
941 
942 
943 
  • Generated by doxygen 1.8.4