def plot_fermi_surface(data, structure, cbm, energy_levels=[],
                       multiple_figure=True,
                       mlab_figure=None, kpoints_dict={}, color=(0, 0, 1),
                       transparency_factor=[], labels_scale_factor=0.05,
                       points_scale_factor=0.02, interative=True):
    """
    Plot the Fermi surface at specific energy value.

    Args:
        data: energy values in a 3D grid from a CUBE file 
              via read_cube_file function, or from a 
              BoltztrapAnalyzer.fermi_surface_data
        structure: structure object of the material
        energy_levels: list of energy value of the fermi surface. 
                       By default 0 eV correspond to the VBM, as in 
                       the plot of band structure along symmetry line.
                       Default: max energy value + 0.01 eV
        cbm: Boolean value to specify if the considered band is 
             a conduction band or not
        multiple_figure: if True a figure for each energy level will be shown.
                         If False all the surfaces will be shown in the same figure.
                         In this las case, tune the transparency factor.
        mlab_figure: provide a previous figure to plot a new surface on it.
        kpoints_dict: dictionary of kpoints to show in the plot.
                      example: {"K":[0.5,0.0,0.5]}, 
                      where the coords are fractional.
        color: tuple (r,g,b) of integers to define the color of the surface.
        transparency_factor: list of values in the range [0,1] to tune
                             the opacity of the surfaces.
        labels_scale_factor: factor to tune the size of the kpoint labels
        points_scale_factor: factor to tune the size of the kpoint points
        interative: if True an interactive figure will be shown.
                    If False a non interactive figure will be shown, but
                    it is possible to plot other surfaces on the same figure.
                    To make it interactive, run mlab.show().
        
    Returns:
        a Mayavi figure and a mlab module to control the plot.

    Note: Experimental. 
          Please, double check the surface shown by using some 
          other software and report issues.
    """

    try:
        from mayavi import mlab
    except ImportError:
        raise BoltztrapError(
            "Mayavi package should be installed to use this function")

    bz = structure.lattice.reciprocal_lattice.get_wigner_seitz_cell()
    cell = structure.lattice.reciprocal_lattice.matrix

    fact = 1 if cbm == False else -1
    en_min = np.min(fact * data.ravel())
    en_max = np.max(fact * data.ravel())

    if energy_levels == []:
        energy_levels = [en_min + 0.01] if cbm == True else \
            [en_max - 0.01]
        print("Energy level set to: " + str(energy_levels[0]) + " eV")

    else:
        for e in energy_levels:
            if e > en_max or e < en_min:
                raise BoltztrapError("energy level " + str(e) +
                                     " not in the range of possible energies: [" +
                                     str(en_min) + ", " + str(en_max) + "]")

    if transparency_factor == []:
        transparency_factor = [1] * len(energy_levels)

    if mlab_figure:
        fig = mlab_figure

    if mlab_figure == None and not multiple_figure:
        fig = mlab.figure(size=(1024, 768), bgcolor=(1, 1, 1))
        for iface in range(len(bz)):
            for line in itertools.combinations(bz[iface], 2):
                for jface in range(len(bz)):
                    if iface < jface and any(np.all(line[0] == x)
                                             for x in bz[jface]) and \
                            any(np.all(line[1] == x)
                                for x in bz[jface]):
                        mlab.plot3d(*zip(line[0], line[1]), color=(0, 0, 0),
                                    tube_radius=None, figure=fig)
        for label, coords in kpoints_dict.items():
            label_coords = structure.lattice.reciprocal_lattice \
                .get_cartesian_coords(coords)
            mlab.points3d(*label_coords, scale_factor=points_scale_factor,
                          color=(0, 0, 0), figure=fig)
            mlab.text3d(*label_coords, text=label, scale=labels_scale_factor,
                        color=(0, 0, 0), figure=fig)

    for isolevel, alpha in zip(energy_levels, transparency_factor):
        if multiple_figure:
            fig = mlab.figure(size=(1024, 768), bgcolor=(1, 1, 1))

            for iface in range(len(bz)):
                for line in itertools.combinations(bz[iface], 2):
                    for jface in range(len(bz)):
                        if iface < jface and any(np.all(line[0] == x)
                                                 for x in bz[jface]) and \
                                any(np.all(line[1] == x)
                                    for x in bz[jface]):
                            mlab.plot3d(*zip(line[0], line[1]), color=(0, 0, 0),
                                        tube_radius=None, figure=fig)

            for label, coords in kpoints_dict.items():
                label_coords = structure.lattice.reciprocal_lattice \
                    .get_cartesian_coords(coords)
                mlab.points3d(*label_coords, scale_factor=points_scale_factor,
                              color=(0, 0, 0), figure=fig)
                mlab.text3d(*label_coords, text=label,
                            scale=labels_scale_factor, color=(0, 0, 0),
                            figure=fig)

        cp = mlab.contour3d(fact * data, contours=[isolevel], transparent=True,
                            colormap='hot', color=color, opacity=alpha,
                            figure=fig)

        polydata = cp.actor.actors[0].mapper.input
        pts = np.array(polydata.points)  # - 1
        polydata.points = np.dot(pts,
                                 cell / np.array(data.shape)[:, np.newaxis])

        cx, cy, cz = [np.mean(np.array(polydata.points)[:, i])
                      for i in range(3)]

        polydata.points = (np.array(polydata.points) - [cx, cy, cz]) * 2

        #mlab.view(distance='auto')
        fig.scene.isometric_view() 
    if interative == True:
        mlab.show()

    return fig, mlab