How to use the pymap3d.vincenty.vdist function in pymap3d

def test_unit(lat, lon, lat1, lon1, srange, az):
    dist, az1 = vincenty.vdist(lat, lon, lat1, lon1)
    assert dist == approx(srange, rel=0.005)
    assert az1 == approx(az)

def test_identity(lat, lon, slantrange, az):
    lat1, lon1 = vincenty.vreckon(lat, lon, slantrange, az)

    dist, az1 = vincenty.vdist(lat, lon, lat1, lon1)

    assert dist == approx(slantrange)
    assert az1 == approx(az)

def bench_vdist(N: int) -> float:
    lat = np.random.random(N)
    lon = np.random.random(N)

    tic = time.monotonic()
    asr, aaz = vdist(*ll0, lat, lon)

    return time.monotonic() - tic

def test_vector():
    pytest.importorskip("numpy")
    asr, aaz = vincenty.vdist(10, 20, [10.02137267, 10.01917819], [20.0168471, 20.0193493])

    assert 3e3 == approx(asr)
    assert aaz == approx([38, 45])

print(exc, file=sys.stderr)
    eng = None


def matlab_func(lat1: float, lon1: float, lat2: float, lon2: float) -> typing.Tuple[float, float]:
    """ Using Matlab Engine to do same thing as Pymap3d """
    ell = eng.wgs84Ellipsoid()
    return eng.distance(lat1, lon1, lat2, lon2, ell, nargout=2)


dlast, alast = nan, nan
lon1, lon2 = 0.0, 1.0
for i in range(20):
    lat1 = lat2 = 10.0 ** (-i)
    try:
        dist_m, az_deg = vdist(lat1, lon1, lat2, lon2)
    except Exception as exc:
        print(exc, f"at latitudes {lat1} {lat2}", file=sys.stderr)
        continue
    assert dist_m != dlast
    assert az_deg != alast
    mat_match = True
    dist_matlab, az_matlab = matlab_func(lat1, lon1, lat2, lon2)
    if not isclose(dist_matlab, dist_m):
        mat_match = False
        print(f"MISMATCH: latitude {lat1} {lat2}: Python: {dist_m}  Matlab: {dist_matlab}", file=sys.stderr)
    if not isclose(az_matlab, az_deg):
        mat_match = False
        print(f"MISMATCH: latitude {lat1} {lat2}: Python: {az_matlab}  Matlab: {az_deg}", file=sys.stderr)
    if mat_match:
        print(f"latitudes {lat1} {lat2}: {dist_m} meters {az_deg} deg azimuth")

def main():
    p = ArgumentParser(description="vdist distance between WGS-84 coordinates")
    p.add_argument("lat1", help="latitude1 WGS-84 [degrees]", type=float)
    p.add_argument("lon1", help="longitude1 WGS-84 [degrees]", type=float)
    p.add_argument("lat2", help="latitude2 WGS-84 [degrees]", type=float)
    p.add_argument("lon2", help="longitude2 WGS-84 [degrees]", type=float)
    P = p.parse_args()

    dist_m = vdist(P.lat1, P.lon1, P.lat2, P.lon2)

    print("{:.3f} meters {:.3f} deg azimuth".format(*dist_m))

degrees input/output  (False: radians in/out)

    Returns
    -------
    radius: "ndarray"
        radius of sphere
    """
    if not deg:
        lat1, lon1, lat2, lon2 = degrees(lat1), degrees(lon1), degrees(lat2), degrees(lon2)
    if asarray is not None:
        lat1, lat2 = asarray(lat1), asarray(lat2)

    latmid = lat1 + (lat2 - lat1) / 2  # compute the midpoint

    # compute azimuth
    az = vdist(lat1, lon1, lat2, lon2, ell=ell)[1]

    #   compute meridional and transverse radii of curvature
    rho = rcurve_meridian(latmid, ell, deg=True)
    nu = rcurve_transverse(latmid, ell, deg=True)

    az = radians(az)
    den = rho * sin(az) ** 2 + nu * cos(az) ** 2

    #  compute radius of the arc from point 1 to point 2
    return rho * nu / den