Skip to content

Commit 3e4908f

Browse files
adamshapiro0scivision
authored andcommitted
Don't skip Beta calculation if some values would divide-by-zero.
If one or more values would divide by zero and needs to use a fallback Beta value, it's still necessary to calculate the correct Beta for any other values that would have succeeded. Previously, this used the same fallback value for all positions if even just one failed.
1 parent e847c32 commit 3e4908f

File tree

1 file changed

+16
-8
lines changed

1 file changed

+16
-8
lines changed

src/pymap3d/ecef.py

Lines changed: 16 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
import warnings
66

77
try:
8-
from numpy import asarray, finfo, where
8+
from numpy import asarray, finfo, full_like, logical_and, logical_or, ndarray, where
99

1010
from .eci import ecef2eci, eci2ecef
1111
except ImportError:
@@ -144,16 +144,24 @@ def ecef2geodetic(
144144

145145
# eqn. 4b
146146
try:
147+
calculated_idx = full_like(z, False)
148+
Beta = ndarray(z.shape)
147149
with warnings.catch_warnings(record=True):
148150
warnings.simplefilter("error")
149-
Beta = atan(huE / u * z / hypot(x, y))
151+
# Preempt any possible divide-by-zero errors by only calculating values that will succeed. Then, below, use
152+
# the fallback Beta values for any values that would have failed.
153+
xy_hypot = hypot(x, y)
154+
calculated_idx = ~logical_or(isclose(u, 0), isclose(xy_hypot, 0))
155+
Beta[calculated_idx] = atan(huE[calculated_idx] / u[calculated_idx] * z[calculated_idx] /
156+
xy_hypot[calculated_idx])
150157
except (ArithmeticError, RuntimeWarning):
151-
if any(isclose(z, 0)):
152-
Beta = 0
153-
elif z > 0:
154-
Beta = pi / 2
155-
else:
156-
Beta = -pi / 2
158+
pass
159+
160+
if not all(calculated_idx):
161+
not_calculated_idx = ~calculated_idx
162+
Beta[logical_and(not_calculated_idx, z > 0)] = pi / 2
163+
Beta[logical_and(not_calculated_idx, z < 0)] = -pi / 2
164+
Beta[logical_and(not_calculated_idx, isclose(z, 0))] = 0
157165

158166
# eqn. 13
159167
dBeta = ((ell.semiminor_axis * u - ell.semimajor_axis * huE + E**2) * sin(Beta)) / (

0 commit comments

Comments
 (0)