docs: def diff(xi, yi, order=1) -> np.ndarray:
      """Take the numerical derivative of a 1D array.
  
      Output is mapped onto the original coordinates  using linear interpolation.
      Expects monotonic xi values.
  
      Parameters
      ----------
      xi : 1D array-like
          Coordinates.
      yi : 1D array-like
          Values.
      order : positive integer (optional)
          Order of differentiation.
  
      Returns
      -------
      1D numpy array
          Numerical derivative. Has the same shape as the input arrays.
      """
      yi = np.array(yi).copy()
      flip = False
      if xi[-1] < xi[0]:
          xi = np.flipud(xi.copy())
          yi = np.flipud(yi)
          flip = True
      midpoints = (xi[1:] + xi[:-1]) / 2
      for _ in range(order):
          d = np.diff(yi)
          d /= np.diff(xi)
          yi = np.interp(xi, midpoints, d)
      if flip:
          yi = np.flipud(yi)
      return yi