@@ -12,6 +12,48 @@ static double kFirstEccentricitySquared = 6.69437999014 * 0.001;
12
12
static double kSecondEccentricitySquared = 6.73949674228 * 0.001 ;
13
13
static double kFlattening = 1 / 298.257223563 ;
14
14
15
+ inline void Geodetic2Ecef (const double latitude, const double longitude, const double altitude,
16
+ double * x, double * y, double * z)
17
+ {
18
+ // Convert geodetic coordinates to ECEF.
19
+ // http://code.google.com/p/pysatel/source/browse/trunk/coord.py?r=22
20
+ double lat_rad = deg2Rad (latitude);
21
+ double lon_rad = deg2Rad (longitude);
22
+ double xi = sqrt (1 - kFirstEccentricitySquared * sin (lat_rad) * sin (lat_rad));
23
+ *x = (kSemimajorAxis / xi + altitude) * cos (lat_rad) * cos (lon_rad);
24
+ *y = (kSemimajorAxis / xi + altitude) * cos (lat_rad) * sin (lon_rad);
25
+ *z = (kSemimajorAxis / xi * (1 - kFirstEccentricitySquared ) + altitude) * sin (lat_rad);
26
+ }
27
+
28
+ void Ecef2Geodetic (const double x, const double y, const double z,
29
+ double * latitude, double * longitude, double * altitude)
30
+ {
31
+ // Convert ECEF coordinates to geodetic coordinates.
32
+ // J. Zhu, "Conversion of Earth-centered Earth-fixed coordinates
33
+ // to geodetic coordinates," IEEE Transactions on Aerospace and
34
+ // Electronic Systems, vol. 30, pp. 957-961, 1994.
35
+
36
+ double r = sqrt (x * x + y * y);
37
+ double Esq = kSemimajorAxis * kSemimajorAxis - kSemiminorAxis * kSemiminorAxis ;
38
+ double F = 54 * kSemiminorAxis * kSemiminorAxis * z * z;
39
+ double G = r * r + (1 - kFirstEccentricitySquared ) * z * z - kFirstEccentricitySquared * Esq;
40
+ double C = (kFirstEccentricitySquared * kFirstEccentricitySquared * F * r * r) / pow (G, 3 );
41
+ double S = cbrt (1 + C + sqrt (C * C + 2 * C));
42
+ double P = F / (3 * pow ((S + 1 / S + 1 ), 2 ) * G * G);
43
+ double Q = sqrt (1 + 2 * kFirstEccentricitySquared * kFirstEccentricitySquared * P);
44
+ double r_0 = -(P * kFirstEccentricitySquared * r) / (1 + Q)
45
+ + sqrt (
46
+ 0.5 * kSemimajorAxis * kSemimajorAxis * (1 + 1.0 / Q)
47
+ - P * (1 - kFirstEccentricitySquared ) * z * z / (Q * (1 + Q)) - 0.5 * P * r * r);
48
+ double U = sqrt (pow ((r - kFirstEccentricitySquared * r_0), 2 ) + z * z);
49
+ double V = sqrt (
50
+ pow ((r - kFirstEccentricitySquared * r_0), 2 ) + (1 - kFirstEccentricitySquared ) * z * z);
51
+ double Z_0 = kSemiminorAxis * kSemiminorAxis * z / (kSemimajorAxis * V);
52
+ *altitude = U * (1 - kSemiminorAxis * kSemiminorAxis / (kSemimajorAxis * V));
53
+ *latitude = rad2Deg (atan ((z + kSecondEccentricitySquared * Z_0) / r));
54
+ *longitude = rad2Deg (atan2 (y, x));
55
+ }
56
+
15
57
class GeodeticConverter
16
58
{
17
59
public:
@@ -69,55 +111,13 @@ class GeodeticConverter
69
111
Geodetic2Ecef (latitude, longitude, altitude, x, y, z);
70
112
}
71
113
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
-
85
114
// added just to keep the old API.
86
115
inline void ecef2Geodetic (const double x, const double y, const double z,
87
116
double * latitude, double * longitude, double * altitude)
88
117
{
89
118
Ecef2Geodetic (x, y, z, latitude, longitude, altitude);
90
119
}
91
120
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
-
121
121
void ecef2Ned (const double x, const double y, const double z,
122
122
double * north, double * east, double * down)
123
123
{
0 commit comments