docs: def distribute_q(network, allocation='p_nom'):
      
      """ Function that distributes reactive power at bus to all installed
      generators and storages. 
  
      Parameters
      ----------
      network : :class:`pypsa.Network
          Overall container of PyPSA
      allocation: str
          Choose key to distribute reactive power:
          'p_nom' to dirstribute via p_nom
          'p' to distribute via p_set
  
      Returns
      -------
  
          
      """
      network.allocation = allocation
      if allocation == 'p':
          p_sum = network.generators_t['p'].\
              groupby(network.generators.bus, axis=1).sum().\
              add(network.storage_units_t['p'].abs().groupby(
                  network.storage_units.bus, axis=1).sum(), fill_value=0)
          q_sum = network.generators_t['q'].\
              groupby(network.generators.bus, axis=1).sum()
  
          q_distributed = network.generators_t.p / \
              p_sum[network.generators.bus.sort_index()].values * \
              q_sum[network.generators.bus.sort_index()].values
  
          q_storages = network.storage_units_t.p / \
              p_sum[network.storage_units.bus.sort_index()].values *\
              q_sum[network.storage_units.bus.sort_index()].values
  
      if allocation == 'p_nom':
  
          q_bus = network.generators_t['q'].\
              groupby(network.generators.bus, axis=1).sum().add(
                  network.storage_units_t.q.groupby(
                  network.storage_units.bus, axis = 1).sum(), fill_value=0)
  
          p_nom_dist = network.generators.p_nom_opt.sort_index()
          p_nom_dist[p_nom_dist.index.isin(network.generators.index
                                           [network.generators.carrier ==
                                            'load shedding'])] = 0
  
          q_distributed = q_bus[
              network.generators.bus].multiply(p_nom_dist.values) /\
              (network.generators.p_nom_opt[network.generators.carrier !=
                                            'load shedding'].groupby(
                  network.generators.bus).sum().add(
                  network.storage_units.p_nom_opt.groupby
                  (network.storage_units.bus).sum(), fill_value=0))[
              network.generators.bus.sort_index()].values
  
          q_distributed.columns = network.generators.index
  
          q_storages = q_bus[network.storage_units.bus]\
              .multiply(network.storage_units.p_nom_opt.values) / \
              ((network.generators.p_nom_opt[network.generators.carrier !=
                                            'load shedding'].groupby(
                  network.generators.bus).sum().add(
                  network.storage_units.p_nom_opt.
                  groupby(network.storage_units.bus).sum(), fill_value=0))[
              network.storage_units.bus].values)
  
          q_storages.columns = network.storage_units.index
  
      q_distributed[q_distributed.isnull()] = 0
      q_distributed[q_distributed.abs() == np.inf] = 0
      q_storages[q_storages.isnull()] = 0
      q_storages[q_storages.abs() == np.inf] = 0
      network.generators_t.q = q_distributed
      network.storage_units_t.q = q_storages
  
      return network