<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def kmeans(*args, **kwargs):
 """ Find centers of multi-modal clusters of data using a kmeans approach modified for spherical measurements. Parameters *args : 2 or 3 sequences of measurements By default, this will be expected to be ``strike`` & ``dip``, both array-like sequences representing poles to planes. (Rake measurements require three parameters, thus the variable number of arguments.) The ``measurement`` kwarg controls how these arguments are interpreted. num : int The number of clusters to find. Defaults to 2. bidirectional : bool Whether or not the measurements are bi-directional linear/planar features or directed vectors. Defaults to True. tolerance : float Iteration will continue until the centers have not changed by more than this amount. Defaults to 1e-5. measurement : string, optional Controls how the input arguments are interpreted. Defaults to ``"poles"``. May be one of the following: ``"poles"`` : strikes, dips Arguments are assumed to be sequences of strikes and dips of planes. Poles to these planes are used for analysis. ``"lines"`` : plunges, bearings Arguments are assumed to be sequences of plunges and bearings of linear features. ``"rakes"`` : strikes, dips, rakes Arguments are assumed to be sequences of strikes, dips, and rakes along the plane. ``"radians"`` : lon, lat Arguments are assumed to be "raw" longitudes and latitudes in the stereonet's underlying coordinate system. Returns ------- centers : An Nx2 array-like Longitude and latitude in radians of the centers of each cluster. """

    lon, lat = _convert_measurements(args, kwargs.get('measurement', 'poles'))
    num = kwargs.get('num', 2)
    bidirectional = kwargs.get('bidirectional', True)
    tolerance = kwargs.get('tolerance', 1e-5)

    points = lon, lat
    dist = lambda x: stereonet_math.angular_distance(x, points, bidirectional)

    center_lon = np.random.choice(lon, num)
    center_lat = np.random.choice(lat, num)
    centers = np.column_stack([center_lon, center_lat])

    while True:
        dists = np.array([dist(item) for item in centers]).T
        closest = dists.argmin(axis=1)

        new_centers = []
        for i in range(num):
            mask = mask = closest == i
            _, vecs = cov_eig(lon[mask], lat[mask], bidirectional)
            new_centers.append(stereonet_math.cart2sph(*vecs[:,-1]))

        if np.allclose(centers, new_centers, atol=tolerance):
            break
        else:
            centers = new_centers

    return centers