def calculate_surface_energy(deltamux, deltamuy, x_energy, y_energy,
    xexcess, yexcess, normalised_bulk):
    """Calculates the surface for a given chemical potential of
    species x and species y for a single phase.

    .. math::
        \\gamma_{Surf} = \\frac{1}{2S} \\Bigg( E_{MO}^{slab} -
        \\frac{nCat_{Slab}}{nCat_{Bulk}} E_{MO}^{Bulk} \\Bigg) -
         \\Gamma_O \\mu_O - \\Gamma_{H_2O} \\mu_{H_2O} -
         \\Gamma_O  \\mu_O (T) - \\Gamma_{H_2O} \\mu_{H_2O} (T)

    where S is the surface area, :math:`E_{MO}^{slab}` is the DFT energy of
    the stoichiometric slab, :math:`nCat_{Slab}` is the number of cations
    in the slab, :math:`nCat_{Slab}` is the number of cations in the bulk
    unit cell, :math:`E_{MO}^{Bulk}` is the DFT energy of the bulk unit cell,
    :math:`\\Gamma_O`  :math:`\\Gamma_{H_2O}` is the excess oxygen / water at
    the surface and :math:`\\mu_O` :math:`\\mu_{H_2O}` is the oxygen /
    water chemcial potential.

    Parameters
    ----------
    deltamux : array like
        Chemical potential of species x
    deltamuy : array like
        Chemical potential of species y
    x_energy : float
        DFT energy or temperature corrected DFT energy
    y_energy : float
        DFT energy or temperature corrected DFT energy
    xexcess : float
        Surface excess of species x
    yexcess : float
        Surface excess of species y
    normalised_bulk : float
        Slab energy normalised to the bulk value.

    Returns
    -------
    array like:
        2D array of surface energies as a function of
        chemical potential of x and y
    """
    return (normalised_bulk - deltamux * xexcess - deltamuy * yexcess - 
        x_energy * xexcess - y_energy * yexcess)
