def area_frac_vs_chempot_plot(self, ref_delu, chempot_range, delu_dict=None,
                                  delu_default=0, increments=10, no_clean=False, no_doped=False):
        """
        1D plot. Plots the change in the area contribution
        of each facet as a function of chemical potential.

        Args:
            ref_delu (sympy Symbol): The free variable chempot with the format:
                Symbol("delu_el") where el is the name of the element.
            chempot_range (list): Min/max range of chemical potential to plot along
            delu_dict (Dict): Dictionary of the chemical potentials to be set as
                constant. Note the key should be a sympy Symbol object of the
                format: Symbol("delu_el") where el is the name of the element.
            delu_default (float): Default value for all unset chemical potentials
            increments (int): Number of data points between min/max or point
                of intersection. Defaults to 10 points.

        Returns:
            (Pylab): Plot of area frac on the Wulff shape
                for each facet vs chemical potential.
        """

        delu_dict = delu_dict if delu_dict else {}
        chempot_range = sorted(chempot_range)
        all_chempots = np.linspace(min(chempot_range), max(chempot_range),
                                   increments)

        # initialize a dictionary of lists of fractional areas for each hkl
        hkl_area_dict = {}
        for hkl in self.all_slab_entries.keys():
            hkl_area_dict[hkl] = []

        # Get plot points for each Miller index
        for u in all_chempots:
            delu_dict[ref_delu] = u
            wulffshape = self.wulff_from_chempot(delu_dict=delu_dict, no_clean=no_clean,
                                                 no_doped=no_doped, delu_default=delu_default)

            for hkl in wulffshape.area_fraction_dict.keys():
                hkl_area_dict[hkl].append(wulffshape.area_fraction_dict[hkl])

        # Plot the area fraction vs chemical potential for each facet
        plt = pretty_plot(width=8, height=7)
        axes = plt.gca()

        for hkl in self.all_slab_entries.keys():
            clean_entry = list(self.all_slab_entries[hkl].keys())[0]
            # Ignore any facets that never show up on the
            # Wulff shape regardless of chemical potential
            if all([a == 0 for a in hkl_area_dict[hkl]]):
                continue
            else:
                plt.plot(all_chempots, hkl_area_dict[hkl],
                         '--', color=self.color_dict[clean_entry],
                         label=str(hkl))

        # Make the figure look nice
        plt.ylabel(r"Fractional area $A^{Wulff}_{hkl}/A^{Wulff}$")
        self.chempot_plot_addons(plt, chempot_range, str(ref_delu).split("_")[1],
                                 axes, rect=[-0.0, 0, 0.95, 1], pad=5, ylim=[0, 1])

        return plt