@@ -12,6 +12,64 @@ static double kFirstEccentricitySquared = 6.69437999014 * 0.001;
1212static double kSecondEccentricitySquared = 6.73949674228 * 0.001 ;
1313static double kFlattening = 1 / 298.257223563 ;
1414
15+
16+ inline double rad2Deg (const double radians)
17+ {
18+ return (radians / M_PI) * 180.0 ;
19+ }
20+
21+ inline double deg2Rad (const double degrees)
22+ {
23+ return (degrees / 180.0 ) * M_PI;
24+ }
25+
26+ inline void Geodetic2Ecef (const double latitude, const double longitude, const double altitude,
27+ double * x, double * y, double * z)
28+ {
29+ // Convert geodetic coordinates to ECEF.
30+ // http://code.google.com/p/pysatel/source/browse/trunk/coord.py?r=22
31+ const double lat_rad = deg2Rad (latitude);
32+ const double lon_rad = deg2Rad (longitude);
33+ const double sLat = sin (lat_rad);
34+ const double sLon = sin (lon_rad);
35+ const double cLat = cos (lat_rad);
36+ const double cLon = cos (lon_rad);
37+
38+ double xi = sqrt (1 - kFirstEccentricitySquared * sLat * sLat );
39+ *x = (kSemimajorAxis / xi + altitude) * cLat * cLon;
40+ *y = (kSemimajorAxis / xi + altitude) * cLat * sLon ;
41+ *z = (kSemimajorAxis / xi * (1 - kFirstEccentricitySquared ) + altitude) * sLat ;
42+ }
43+
44+ inline void Ecef2Geodetic (const double x, const double y, const double z,
45+ double * latitude, double * longitude, double * altitude)
46+ {
47+ // Convert ECEF coordinates to geodetic coordinates.
48+ // J. Zhu, "Conversion of Earth-centered Earth-fixed coordinates
49+ // to geodetic coordinates," IEEE Transactions on Aerospace and
50+ // Electronic Systems, vol. 30, pp. 957-961, 1994.
51+
52+ double r = sqrt (x * x + y * y);
53+ double Esq = kSemimajorAxis * kSemimajorAxis - kSemiminorAxis * kSemiminorAxis ;
54+ double F = 54 * kSemiminorAxis * kSemiminorAxis * z * z;
55+ double G = r * r + (1 - kFirstEccentricitySquared ) * z * z - kFirstEccentricitySquared * Esq;
56+ double C = (kFirstEccentricitySquared * kFirstEccentricitySquared * F * r * r) / pow (G, 3 );
57+ double S = cbrt (1 + C + sqrt (C * C + 2 * C));
58+ double P = F / (3 * pow ((S + 1 / S + 1 ), 2 ) * G * G);
59+ double Q = sqrt (1 + 2 * kFirstEccentricitySquared * kFirstEccentricitySquared * P);
60+ double r_0 = -(P * kFirstEccentricitySquared * r) / (1 + Q)
61+ + sqrt (
62+ 0.5 * kSemimajorAxis * kSemimajorAxis * (1 + 1.0 / Q)
63+ - P * (1 - kFirstEccentricitySquared ) * z * z / (Q * (1 + Q)) - 0.5 * P * r * r);
64+ double U = sqrt (pow ((r - kFirstEccentricitySquared * r_0), 2 ) + z * z);
65+ double V = sqrt (
66+ pow ((r - kFirstEccentricitySquared * r_0), 2 ) + (1 - kFirstEccentricitySquared ) * z * z);
67+ double Z_0 = kSemiminorAxis * kSemiminorAxis * z / (kSemimajorAxis * V);
68+ *altitude = U * (1 - kSemiminorAxis * kSemiminorAxis / (kSemimajorAxis * V));
69+ *latitude = rad2Deg (atan ((z + kSecondEccentricitySquared * Z_0) / r));
70+ *longitude = rad2Deg (atan2 (y, x));
71+ }
72+
1573class GeodeticConverter
1674{
1775 public:
@@ -69,55 +127,13 @@ class GeodeticConverter
69127 Geodetic2Ecef (latitude, longitude, altitude, x, y, z);
70128 }
71129
72- static void Geodetic2Ecef (const double latitude, const double longitude, const double altitude,
73- double * x, double * y, double * z)
74- {
75- // Convert geodetic coordinates to ECEF.
76- // http://code.google.com/p/pysatel/source/browse/trunk/coord.py?r=22
77- double lat_rad = deg2Rad (latitude);
78- double lon_rad = deg2Rad (longitude);
79- double xi = sqrt (1 - kFirstEccentricitySquared * sin (lat_rad) * sin (lat_rad));
80- *x = (kSemimajorAxis / xi + altitude) * cos (lat_rad) * cos (lon_rad);
81- *y = (kSemimajorAxis / xi + altitude) * cos (lat_rad) * sin (lon_rad);
82- *z = (kSemimajorAxis / xi * (1 - kFirstEccentricitySquared ) + altitude) * sin (lat_rad);
83- }
84-
85130 // added just to keep the old API.
86131 inline void ecef2Geodetic (const double x, const double y, const double z,
87132 double * latitude, double * longitude, double * altitude)
88133 {
89134 Ecef2Geodetic (x, y, z, latitude, longitude, altitude);
90135 }
91136
92- static void Ecef2Geodetic (const double x, const double y, const double z,
93- double * latitude, double * longitude, double * altitude)
94- {
95- // Convert ECEF coordinates to geodetic coordinates.
96- // J. Zhu, "Conversion of Earth-centered Earth-fixed coordinates
97- // to geodetic coordinates," IEEE Transactions on Aerospace and
98- // Electronic Systems, vol. 30, pp. 957-961, 1994.
99-
100- double r = sqrt (x * x + y * y);
101- double Esq = kSemimajorAxis * kSemimajorAxis - kSemiminorAxis * kSemiminorAxis ;
102- double F = 54 * kSemiminorAxis * kSemiminorAxis * z * z;
103- double G = r * r + (1 - kFirstEccentricitySquared ) * z * z - kFirstEccentricitySquared * Esq;
104- double C = (kFirstEccentricitySquared * kFirstEccentricitySquared * F * r * r) / pow (G, 3 );
105- double S = cbrt (1 + C + sqrt (C * C + 2 * C));
106- double P = F / (3 * pow ((S + 1 / S + 1 ), 2 ) * G * G);
107- double Q = sqrt (1 + 2 * kFirstEccentricitySquared * kFirstEccentricitySquared * P);
108- double r_0 = -(P * kFirstEccentricitySquared * r) / (1 + Q)
109- + sqrt (
110- 0.5 * kSemimajorAxis * kSemimajorAxis * (1 + 1.0 / Q)
111- - P * (1 - kFirstEccentricitySquared ) * z * z / (Q * (1 + Q)) - 0.5 * P * r * r);
112- double U = sqrt (pow ((r - kFirstEccentricitySquared * r_0), 2 ) + z * z);
113- double V = sqrt (
114- pow ((r - kFirstEccentricitySquared * r_0), 2 ) + (1 - kFirstEccentricitySquared ) * z * z);
115- double Z_0 = kSemiminorAxis * kSemiminorAxis * z / (kSemimajorAxis * V);
116- *altitude = U * (1 - kSemiminorAxis * kSemiminorAxis / (kSemimajorAxis * V));
117- *latitude = rad2Deg (atan ((z + kSecondEccentricitySquared * Z_0) / r));
118- *longitude = rad2Deg (atan2 (y, x));
119- }
120-
121137 void ecef2Ned (const double x, const double y, const double z,
122138 double * north, double * east, double * down)
123139 {
@@ -240,18 +256,6 @@ class GeodeticConverter
240256 return ret;
241257 }
242258
243- inline static
244- double rad2Deg (const double radians)
245- {
246- return (radians / M_PI) * 180.0 ;
247- }
248-
249- inline static
250- double deg2Rad (const double degrees)
251- {
252- return (degrees / 180.0 ) * M_PI;
253- }
254-
255259 double initial_latitude_;
256260 double initial_longitude_;
257261 double initial_altitude_;
0 commit comments