docs: def diff_prof(step):
      """Diffusion.
  
      Args:
          step (:class:`~stagpy.stagyydata._Step`): a step of a StagyyData
              instance.
      Returns:
          tuple of :class:`numpy.array`: the diffusion and the radial position
          at which it is evaluated.
      """
      rbot, rtop = misc.get_rbounds(step)
      rad = step.rprof['r'].values + rbot
      tprof = step.rprof['Tmean'].values
      diff = (tprof[:-1] - tprof[1:]) / (rad[1:] - rad[:-1])
      # assume tbot = 1
      diff = np.insert(diff, 0, (1 - tprof[0]) / (rad[0] - rbot))
      # assume ttop = 0
      diff = np.append(diff, tprof[-1] / (rtop - rad[-1]))
      # actually computed at r_edges...
      return diff, np.append(rad, rtop)