[Bug]: Discontinuity in ecef2geodetic() outputs

[ghost]

What happened?

If I use the ecef2geodetic function with x >= 521860, y = 0, and z = 0, I get correct outputs of lat = 0, lon =0, and the appropriate alt.

However, if I use x <= 521850, y = 0, z = 0, I get a lat of 90 (incorrect), lon of 0, and an incorrect altitude.

The difference seems to occur in the calculation of u on line 136 in ecef.py, which is:

u = sqrt(0.5 * (r**2 - E**2) + 0.5 * sqrt((r**2 - E**2) ** 2 + 4 * E**2 * z**2))

If r < E, then u = 0.

This can be reproduced by running the following:

import pymap3d

from numpy import linspace

Xs = linspace(521000, 522000, 101)

for x in Xs:

    lat, lon, alt = pymap3d.ecef2geodetic(x, 0, 0)
    
    print(x)
    print(lat, lon, alt)
    print(" ")

You will see that the discontinuity occurs from x=521850 and x=521860
I don't think this should happen, but I'm not sure what the problem is. The equations seem to match the referenced paper.

Any ideas? Thanks!

Relevant log output

No response

[ghost]

Almost forgot (not sure if it matters or not):

In the calculation of Beta, if there is a divide by zero, I get "RuntimeWarning: divide by zero encountered in double_scalars" not a "ZeroDivisionError" and therefore it is not caught by the except statement on line 145.

To fix, I did

import warnings
warnings.filterwarnings('error')

and then changed the except statement to

except (ZeroDivisionError, RuntimeWarning):

[scivision]

Thank you!