def plot_shapefun(pos, out, mode):
    """
  Creates a simple matplotlib image to show the shapefunctions given it's positions in the unit cell, the atoms's vertices in `ut` and the plotting mode

  :author: Philipp Ruessmann

  :param pos: positions of the centers of the cells
  :param verts: array of vertices of the shapefunction (outlines of shapes)
  :param mode: 'all' or 'single' determines whether or not all shapes are combined in a single figure or plotted as individual figures

  :returns ax: return the axis in which the plot was done (useful to pass to 'change_zoom' and 'zoom_in' functions of this module
  """
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    from numpy import array, shape
    if mode == 'all':
        fig = plt.figure()
        ax = fig.add_subplot(111, projection='3d')
    Natm = len(pos)
    for j in range(1, Natm + 1):
        if mode == 'single':
            fig = plt.figure()
            ax = fig.add_subplot(111, projection='3d')
        for i in list(out.keys()):
            if i[0] == j:
                linecolor = 'r'
            if j in [2, 3, 6, 7]:
                linecolor = 'b'
            d = array(out[i])
            if len(shape(pos)) > 1 and len(pos) > 2 and mode == 'all':
                dx, dy, dz = pos[j - 1, 0], pos[j - 1, 1], pos[j - 1, 2]
            else:
                dx, dy, dz = 0, 0, 0
            ax.plot(d[:, (0)] + dx, d[:, (1)] + dy, d[:, (2)] + dz, linecolor)
    if mode == 'all':
        xy_scale = 1
        scale = 1
        ax.set_xlim(-scale * xy_scale, scale * xy_scale)
        ax.set_ylim(-scale * xy_scale, scale * xy_scale)
        ax.set_zlim(-scale, scale)
    return ax
