docs: def grid_powerlaw(exp, insert=0, skip=0, end=-1, step=0, axis=None, **kwargs):
    r'''
  Draw a power-law grid: a grid that respects a certain power-law exponent. The grid-lines start from
  the positions of the ticks.
  
  :arguments:
  
    **exp** (``float``)
      The power-law exponent.
  
  :options:
  
    **insert** (``<int>``)
      Insert extra lines in between the default lines set by the tick positions.
  
    **skip, end, step** (``<int>``)
      Select from the lines based on ``coor = coor[skip:end:step]``.
  
    **axis** ([``plt.gca()``] | ...)
      Specify the axis to which to apply the limits.
  
    ...
      Any ``plt.plot(...)`` option.
  
  :returns:
  
    The handle of the ``plt.plot(...)`` command.
    '''
  
    # default axis
    if axis is None: axis = plt.gca()
  
    # default plot settings
    kwargs.setdefault('color'    , 'k' )
    kwargs.setdefault('linestyle', '--')
    kwargs.setdefault('linewidth',  1  )
  
    # check
    if axis.get_xscale() != 'log' or axis.get_yscale() != 'log':
      raise IOError('This function only works on a log-log scale, where the power-law is a straight line')
  
    # zero-exponent: draw horizontal lines
    if exp == 0:
  
      # y-coordinate of the start positions
      starty = abs2rel_y(axis.get_yticks(), axis=axis)
  
      # insert extra coordinates
      if insert > 0:
  
        n  = len(starty)
        x  = np.linspace(0,1,n+(n-1)*int(insert))
        xp = np.linspace(0,1,n)
  
        starty = np.interp(x, xp, starty)
  
      # skip coordinates
      starty = starty[int(skip):int(end):int(1+step)]
  
      # set remaining coordinates
      endy   = starty
      startx = np.zeros((len(starty)))
      endx   = np.ones ((len(starty)))
  
    # all other exponents
    else:
  
      # get the axis' size in real coordinates
      # - get the limits
      xmin, xmax = axis.get_xlim()
      ymin, ymax = axis.get_ylim()
      # - compute the size in both directions
      deltax = np.log10(xmax) - np.log10(xmin)
      deltay = np.log10(ymax) - np.log10(ymin)
  
      # convert the exponent in real coordinates to an exponent in relative coordinates
      b = np.abs(exp) * deltax / deltay
  
      # x-coordinate of the start positions
      startx = abs2rel_x(axis.get_xticks(), axis=axis)
  
      # compute how many labels need to be prepended
      Dx   = startx[1] - startx[0]
      nneg = int(np.floor(1./(b*Dx))) - 1
  
      # add extra to be sure
      if insert > 0:
        nneg += 1
  
      # prepend
      if nneg > 0:
        startx = np.hstack(( startx[0]+np.cumsum(-Dx * np.ones((nneg)))[::-1], startx ))
  
      # insert extra coordinates
      if insert > 0:
  
        n  = len(startx)
        x  = np.linspace(0,1,n+(n-1)*int(insert))
        xp = np.linspace(0,1,n)
  
        startx = np.interp(x, xp, startx)
  
      # skip coordinates
      if step > 0:
  
        startx = startx[int(skip)::int(1+step)]
  
      # x-coordinate of the end of the lines
      endx = startx + 1/b
  
      # y-coordinate of the start and the end of the lines
      if exp > 0:
        starty = np.zeros((len(startx)))
        endy   = np.ones ((len(startx)))
      else:
        starty = np.ones ((len(startx)))
        endy   = np.zeros((len(startx)))
  
    # convert to real coordinates
    startx = rel2abs_x(startx, axis)
    endx   = rel2abs_x(endx  , axis)
    starty = rel2abs_y(starty, axis)
    endy   = rel2abs_y(endy  , axis)
  
    # plot
    lines = axis.plot(np.vstack(( startx, endx )), np.vstack(( starty, endy )), **kwargs)
  
    # remove access in labels
    plt.setp(lines[1:], label="_")
  
    # return handles
    return lines