def METHOD_NAME(x, y, z=None, I=None, dip=0, strike=0, var=1):
    """Geostatistical covariance matrix for given points.
    Parameters
    ----------
    **kwargs
        I : float or list of floats
            correlation lengths (range) in individual directions
        dip : float
            dip angle (in degrees) of major axis (I[0])
        strike : float
            strike angle (for 3D)
    Returns
    -------
    Cm : np.array (square matrix of size cellCount/nodeCount)
        covariance matrix
    """
    if I is None:
        I = [1, 1, 1]
    elif isinstance(I, (float, int)):
        I = [I, I, I]
    elif len(I) < 3:
        I.append(I[-1])
    if z is None:
        z = np.zeros_like(x)
    hx = x - x[:, np.newaxis]
    hy = y - y[:, np.newaxis]
    hz = z - z[:, np.newaxis]
    alpha = -dip * pi / 180  # rotation of operator
    beta = -strike * pi / 180
    # compute lags, normalized by correlation lengths
    Hx = (hx*cos(alpha)-hy*sin(alpha))*cos(beta) / I[0]
    Hy = (hx*sin(alpha)+hy*cos(alpha))*cos(beta) / I[1]
    Hz = hz * sin(beta) / I[2]
    CM = var*np.exp(-np.sqrt(Hx**2+Hy**2+Hz**2))  # Covariance matrix
    return CM