def vrh(kclay, kqtz, vclay):
    """
    Voigt-Reuss-Hill average to find Kmatrix from clay and qtz components.

    From Smith et al, Geophysics 68(2), 2003.

    Works for any two components.

    Returns Kvrh, AKA Kmatrix.
    """
    vqtz = 1 - vclay
    kreuss = 1.0 / (vclay / kclay + vqtz / kqtz)
    kvoigt = vclay * kclay + vqtz * kqtz
    kvrh = 0.5 * (kreuss + kvoigt)
    return kvrh
