docs: def plot_line_loading(
          network,
          timesteps=range(1,2),
          filename=None,
          boundaries=[],
          arrows=False):
      """
      Plots line loading as a colored heatmap.
  
      Line loading is displayed as relative to nominal capacity in %.
      
      Parameters
      ----------
      network : PyPSA network container
          Holds topology of grid including results from powerflow analysis
      timesteps : range 
          Defines which timesteps are considered. If more than one, an 
          average line loading is calculated.
      filename : str
          Specify filename
          If not given, figure will be show directly
      boundaries : list
          If given, the colorbar is fixed to a given min and max value
      arrows : bool
          If True, the direction of the power flows is displayed as 
          arrows.
      """
      # TODO: replace p0 by max(p0,p1) and analogously for q0
      # TODO: implement for all given snapshots
  
      # calculate relative line loading as S/S_nom
      # with S = sqrt(P^2 + Q^2)
      cmap = plt.cm.jet
      array_line = [['Line'] * len(network.lines), network.lines.index]
      array_link = [['Link'] * len(network.links), network.links.index]
      
      if network.lines_t.q0.empty:
  
          loading_lines = pd.Series((network.lines_t.p0.mul(
              network.snapshot_weightings, axis=0).loc[network.snapshots[
              timesteps]].abs().sum() / (network.lines.s_nom_opt)).data,
              index=array_line)
  
      else:
          loading_lines = pd.Series(((network.lines_t.p0.mul(
                  network.snapshot_weightings, axis=0)\
                      .loc[network.snapshots[timesteps]].abs().sum() ** 2 +\
                      network.lines_t.q0.mul(
                              network.snapshot_weightings, axis=0)\
                      .loc[network.snapshots[timesteps]].abs().sum() ** 2).\
                      apply(sqrt) / (network.lines.s_nom_opt)).data, index =
                              array_line)
  
      # Aviod covering of bidirectional links
      network.links['linked_to'] = 0
      for i,  row in network.links.iterrows():
          if not (network.links.index[(network.links.bus0 == row['bus1']) &
                                    (network.links.bus1 == row['bus0']) &
                                    (network.links.length == row['length']
                                    )]).empty:
  
              l = network.links.index[(network.links.bus0 == row['bus1']) &
                                    (network.links.bus1 == row['bus0']) &
                                    (network.links.length == row['length'])]
  
              network.links.set_value(i, 'linked_to',l.values[0])
  
      network.links.linked_to = network.links.linked_to.astype(str)
      link_load = network.links_t.p0[network.links.index[
              network.links.linked_to == '0']]
  
      for i, row in network.links[network.links.linked_to != '0'].iterrows():
          load = pd.DataFrame(index = network.links_t.p0.index,
                              columns = ['to', 'from'])
          load['to'] = network.links_t.p0[row['linked_to']]
          load['from'] = network.links_t.p0[i]
          link_load[i] = load.abs().max(axis = 1)
  
      loading_links = pd.Series((link_load.mul(
              network.snapshot_weightings, axis=0).loc[network.snapshots[
              timesteps]].abs().sum()[network.links.index] / (
          network.links.p_nom_opt)).data, index=array_link).dropna()
  
      load_links_rel = (loading_links/  
                        network.snapshot_weightings\
                              [network.snapshots[timesteps]].sum())* 100
          
      load_lines_rel = (loading_lines / network.snapshot_weightings\
                              [network.snapshots[timesteps]].sum()) * 100
  
      loading = load_lines_rel.append(load_links_rel)
  
      ll = network.plot(line_colors=loading, line_cmap=cmap,
                        title="Line loading", line_widths=0.55)
      # add colorbar, note mappable sliced from ll by [1]
  
      if not boundaries:
          v = np.linspace(min(loading), max(loading), 101)
          boundaries = [min(loading), max(loading)]
          
      else:
          v = np.linspace(boundaries[0], boundaries[1], 101)
          
      cb = plt.colorbar(ll[1], boundaries=v,
                        ticks=v[0:101:10])
      cb_Link = plt.colorbar(ll[2], boundaries=v,
                        ticks=v[0:101:10])
  
      cb.set_clim(vmin=boundaries[0], vmax=boundaries[1])
      
      cb_Link.set_clim(vmin=boundaries[0], vmax=boundaries[1])
      
      cb_Link.remove()
      
      cb.set_label('Line loading in %')
  
      if arrows:
          ax = plt.axes()
          path = ll[1].get_segments()
          x_coords_lines = np.zeros([len(path)])
          cmap = cmap
          colors = cmap(ll[1].get_array() / 100)
          for i in range(0, len(path)):
              x_coords_lines[i] = network.buses.loc[str(
                  network.lines.iloc[i, 2]), 'x']
              color = colors[i]
              if (x_coords_lines[i] == path[i][0][0] and load_lines_rel[i] >= 0):
                  arrowprops = dict(arrowstyle="->", color=color)
              else:
                  arrowprops = dict(arrowstyle="<-", color=color)
              ax.annotate(
                  "",
                  xy=abs(
                      (path[i][0] - path[i][1]) * 0.51 - path[i][0]),
                  xytext=abs(
                      (path[i][0] - path[i][1]) * 0.49 - path[i][0]),
                  arrowprops=arrowprops,
                  size=10)
  
      if filename is None:
          plt.show()
      else:
          plt.savefig(filename)
          plt.close()