def K_separator_Watkins(x, rhol, rhog, horizontal=False, method='spline'):
    r'''Calculates the Sounders-Brown `K` factor as used in determining maximum 
    allowable gas velocity in a two-phase separator in either a horizontal or
    vertical orientation. This function approximates a graph published in [1]_ 
    to determine `K` as used in the following equation:
    
    .. math::
        v_{max} =  K_{SB}\sqrt{\frac{\rho_l-\rho_g}{\rho_g}}
    
    The graph has `K_{SB}` on its y-axis, and the following as its x-axis:
    
    .. math::
        \frac{m_l}{m_g}\sqrt{\rho_g/\rho_l}
        = \frac{(1-x)}{x}\sqrt{\rho_g/\rho_l}
    
    Cubic spline interpolation is the default method of retrieving a value
    from the graph, which was digitized with Engauge-Digitizer.
    
    Also supported are two published curve fits to the graph. The first is that
    of Blackwell (1984) [2]_, as follows:
    
    .. math::
        K_{SB} = \exp(-1.942936 -0.814894X -0.179390 X^2 -0.0123790 X^3
        + 0.000386235 X^4 + 0.000259550 X^5)

        X = \log\left[\frac{(1-x)}{x}\sqrt{\rho_g/\rho_l}\right]
    
    The second is that of Branan (1999), as follows:
    
    .. math::
        K_{SB} = \exp(-1.877478097 -0.81145804597X -0.1870744085 X^2 
        -0.0145228667 X^3 -0.00101148518 X^4)
    
        X = \log\left[\frac{(1-x)}{x}\sqrt{\rho_g/\rho_l}\right]

    Parameters
    ----------
    x : float
        Quality of fluid entering separator, [-]
    rhol : float
        Density of liquid phase [kg/m^3]
    rhog : float
        Density of gas phase [kg/m^3]
    horizontal : bool, optional
        Whether to use the vertical or horizontal value; horizontal is 1.25 
        higher
    method : str
        One of 'spline, 'blackwell', or 'branan'

    Returns
    -------
    K : float
        Sounders Brown horizontal or vertical `K` factor for two-phase 
        separator design only, [m/s]

    Notes
    -----
    Both the 'branan' and 'blackwell' models are used frequently. However,
    the spline is much more accurate.
    
    No limits checking is enforced. However, the x-axis spans only 0.006 to
    5.4, and the function should not be used outside those limits.

    Examples
    --------
    >>> K_separator_Watkins(0.88, 985.4, 1.3, horizontal=True)
    0.07951613600476297

    References
    ----------
    .. [1] Watkins (1967). Sizing Separators and Accumulators, Hydrocarbon 
       Processing, November 1967.
    .. [2] Blackwell, W. Wayne. Chemical Process Design on a Programmable 
       Calculator. New York: Mcgraw-Hill, 1984.
    .. [3] Branan, Carl R. Pocket Guide to Chemical Engineering. 1st edition. 
       Houston, Tex: Gulf Professional Publishing, 1999.
    '''
    factor = (1. - x)/x*(rhog/rhol)**0.5
    if method == 'spline':
        K = exp(float(splev(log(factor), tck_Watkins)))
    elif method == 'blackwell':
        X = log(factor)    
        A = -1.877478097
        B = -0.81145804597
        C = -0.1870744085
        D = -0.0145228667
        E = -0.00101148518
        K = exp(A + X*(B + X*(C + X*(D + E*X))))
    elif method == 'branan':
        X = log(factor)
        A = -1.942936
        B = -0.814894
        C = -0.179390
        D = -0.0123790
        E = 0.000386235
        F = 0.000259550
        K = exp(A + X*(B + X*(C + X*(D + X*(E + F*X)))))
    else:
        raise Exception("Only methods 'spline', 'branan', and 'blackwell' are supported.")
    K *= foot # Converts units of ft/s to m/s; the graph and all fits are in ft/s 
    if horizontal:
        K *= 1.25 # Watkins recommends a factor of 1.25 for horizontal separators over vertical separators
    return K