How to use the harmonica.get_ellipsoid function in harmonica

>>> print(", ".join("{:.4f}".format(i) for i in spherical))
    0.0000, 90.0000, 6356752.3142
    >>> print("{:.4f}".format(hm.get_ellipsoid().semiminor_axis))
    6356752.3142

    In  the equator, it should be the semi-major axis:

    >>> spherical = hm.geodetic_to_spherical(longitude=0, latitude=0, height=0)
    >>> print(", ".join("{:.4f}".format(i) for i in spherical))
    0.0000, 0.0000, 6378137.0000
    >>> print("{:.4f}".format(hm.get_ellipsoid().semimajor_axis))
    6378137.0000

    """
    # Get ellipsoid
    ellipsoid = get_ellipsoid()
    # Convert latitude to radians
    latitude_rad = np.radians(latitude)
    prime_vertical_radius = ellipsoid.semimajor_axis / np.sqrt(
        1 - ellipsoid.first_eccentricity ** 2 * np.sin(latitude_rad) ** 2
    )
    # Instead of computing X and Y, we only comupute the projection on the XY
    # plane: xy_projection = sqrt( X**2 + Y**2 )
    xy_projection = (height + prime_vertical_radius) * np.cos(latitude_rad)
    z_cartesian = (
        height + (1 - ellipsoid.first_eccentricity ** 2) * prime_vertical_radius
    ) * np.sin(latitude_rad)
    radius = np.sqrt(xy_projection ** 2 + z_cartesian ** 2)
    spherical_latitude = np.degrees(np.arcsin(z_cartesian / radius))
    return longitude, spherical_latitude, radius

def _spherical_to_geodetic_parameters(spherical_latitude, radius):
    "Compute parameters for spherical to geodetic coordinates conversion"
    # Get ellipsoid
    ellipsoid = get_ellipsoid()
    # Convert latitude to radians
    spherical_latitude_rad = np.radians(spherical_latitude)
    big_z = radius * np.sin(spherical_latitude_rad)
    p_0 = (
        radius ** 2
        * np.cos(spherical_latitude_rad) ** 2
        / ellipsoid.semimajor_axis ** 2
    )
    q_0 = (
        (1 - ellipsoid.first_eccentricity ** 2)
        / ellipsoid.semimajor_axis ** 2
        * big_z ** 2
    )
    r_0 = (p_0 + q_0 - ellipsoid.first_eccentricity ** 4) / 6
    s_0 = ellipsoid.first_eccentricity ** 4 * p_0 * q_0 / 4 / r_0 ** 3
    t_0 = np.cbrt(1 + s_0 + np.sqrt(2 * s_0 + s_0 ** 2))

...     spherical_latitude=0,
    ...     radius=hm.get_ellipsoid().semimajor_axis,
    ... )
    >>> print(", ".join("{:.1f}".format(i) for i in geodetic))
    0.0, 0.0, 0.0
    >>> geodetic = hm.spherical_to_geodetic(
    ...     longitude=0,
    ...     spherical_latitude=-90,
    ...     radius=hm.get_ellipsoid().semiminor_axis + 2
    ... )
    >>> print(", ".join("{:.1f}".format(i) for i in geodetic))
    0.0, -90.0, 2.0

    """
    # Get ellipsoid
    ellipsoid = get_ellipsoid()
    k, big_d, big_z = _spherical_to_geodetic_parameters(spherical_latitude, radius)
    latitude = np.degrees(
        2 * np.arctan(big_z / (big_d + np.sqrt(big_d ** 2 + big_z ** 2)))
    )
    height = (
        (k + ellipsoid.first_eccentricity ** 2 - 1)
        / k
        * np.sqrt(big_d ** 2 + big_z ** 2)
    )
    return longitude, latitude, height