Skip to content

Commit 5ab7858

Browse files
leouiedascivision
authored andcommitted
Don't recompute eccentricity
1 parent 4aff132 commit 5ab7858

File tree

1 file changed

+8
-18
lines changed

1 file changed

+8
-18
lines changed

src/pymap3d/spherical.py

Lines changed: 8 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -85,8 +85,6 @@ def geodetic2spherical(
8585
sinlat = sin(lat)
8686
coslat = sqrt(1 - sinlat**2)
8787

88-
first_eccentricity = _calculate_first_eccentricity(ell)
89-
9088
# radius of curvature of the prime vertical section
9189
N = ell.semimajor_axis**2 / hypot(
9290
ell.semimajor_axis * coslat, ell.semiminor_axis * sinlat,
@@ -95,7 +93,7 @@ def geodetic2spherical(
9593
# Instead of computing X and Y, we only compute the projection on the XY
9694
# plane: xy_projection = sqrt( X**2 + Y**2 )
9795
xy_projection = (alt + N) * coslat
98-
z_cartesian = (alt + (1 - first_eccentricity**2) * N) * sinlat
96+
z_cartesian = (alt + (1 - ell.eccentricity**2) * N) * sinlat
9997
radius = hypot(xy_projection, z_cartesian)
10098
slat = asin(z_cartesian / radius)
10199

@@ -152,32 +150,24 @@ def spherical2geodetic(
152150
sinlat = sin(lat)
153151
coslat = sqrt(1 - sinlat**2)
154152

155-
first_eccentricity = _calculate_first_eccentricity(ell)
156153
Z = radius * sinlat
157154
p_0 = pow(radius, 2) * coslat**2 / ell.semimajor_axis**2
158-
q_0 = (1 - first_eccentricity**2) / ell.semimajor_axis**2 * Z**2
159-
r_0 = (p_0 + q_0 - first_eccentricity**4) / 6
160-
s_0 = first_eccentricity**4 * p_0 * q_0 / 4 / r_0**3
155+
q_0 = (1 - ell.eccentricity**2) / ell.semimajor_axis**2 * Z**2
156+
r_0 = (p_0 + q_0 - ell.eccentricity**4) / 6
157+
s_0 = ell.eccentricity**4 * p_0 * q_0 / 4 / r_0**3
161158
t_0 = cbrt(1 + s_0 + sqrt(2 * s_0 + s_0**2))
162159
u_0 = r_0 * (1 + t_0 + 1 / t_0)
163-
v_0 = sqrt(u_0**2 + q_0 * first_eccentricity**4)
164-
w_0 = first_eccentricity**2 * (u_0 + v_0 - q_0) / 2 / v_0
160+
v_0 = sqrt(u_0**2 + q_0 * ell.eccentricity**4)
161+
w_0 = ell.eccentricity**2 * (u_0 + v_0 - q_0) / 2 / v_0
165162
k = sqrt(u_0 + v_0 + w_0**2) - w_0
166-
D = k * radius * coslat / (k + first_eccentricity**2)
163+
D = k * radius * coslat / (k + ell.eccentricity**2)
167164
hypotDZ = hypot(D, Z)
168165

169166
glat = 2 * atan2(Z, (D + hypotDZ))
170-
alt = (k + first_eccentricity**2 - 1) / k * hypotDZ
167+
alt = (k + ell.eccentricity**2 - 1) / k * hypotDZ
171168

172169
if deg:
173170
glat = degrees(glat)
174171
lon = degrees(lon)
175172

176173
return glat, lon, alt
177-
178-
179-
def _calculate_first_eccentricity(ell):
180-
"""
181-
Calculate the first eccentricity of an ellipsoid.
182-
"""
183-
return sqrt(ell.semimajor_axis**2 - ell.semiminor_axis**2) / ell.semimajor_axis

0 commit comments

Comments
 (0)