def get_default_guess_params(model, df, pressure_key, loading_key):
    """
    Get dictionary of default parameters for starting guesses in data fitting
    routine.

    The philosophy behind the default starting guess is that (1) the saturation
    loading is close to the highest loading observed in the data, and (2) the
    default assumption is a Langmuir isotherm.

    Reminder: pass your own guess via `param_guess` in instantiation if these
    default guesses do not lead to a converged set of parameters.

    :param model: String name of analytical model
    :param df: DataFrame adsorption isotherm data
    :param pressure_key: String key for pressure column in df
    :param loading_key: String key for loading column in df
    """
    # guess saturation loading to 10% more than highest loading
    saturation_loading = 1.1 * df[loading_key].max()
    # guess Langmuir K using the guess for saturation loading and lowest
    #   pressure point (but not zero)
    df_nonzero = df[df[loading_key] != 0.0]
    idx_min = df_nonzero[loading_key].idxmin()
    langmuir_k = df_nonzero[loading_key].loc[idx_min] / \
                 df_nonzero[pressure_key].loc[idx_min] / (
                     saturation_loading - df_nonzero[loading_key].loc[idx_min])

    if model == "Langmuir":
        return {"M": saturation_loading, "K": langmuir_k}

    if model == "Freundlich":
        return {"n_inv": 0.5, "K": 10}

    if model == "Quadratic":
        # Quadratic = Langmuir when Kb = Ka^2. This is our default assumption.
        # Also, M is half of the saturation loading in the Quadratic model.
        return {
            "M": saturation_loading / 2.0,
            "Ka": langmuir_k,
            "Kb": langmuir_k**2.0
        }

    if model == "BET":
        # BET = Langmuir when Kb = 0.0. This is our default assumption.
        return {
            "M": saturation_loading,
            "Ka": langmuir_k,
            "Kb": langmuir_k * 0.01
        }

    if model == "DSLangmuir":
        return {
            "M1": 0.5 * saturation_loading,
            "K1": 0.4 * langmuir_k,
            "M2": 0.5 * saturation_loading,
            "K2": 0.6 * langmuir_k
        }

    if model == "Henry":
        return {"KH": saturation_loading * langmuir_k}

    if model == "TemkinApprox":
        # equivalent to Langmuir model if theta = 0.0
        return {"M": saturation_loading, "K": langmuir_k, "theta": 0.0}