docs: def csep_periodic(ra, rb, L):
      """Return separation vectors between each pair of the two sets of points.
  
      Parameters
      ----------
      ra, rb: float array-like, shape (n, d) and (m, d) in d dimensions.
          Two sets of points.
      L: float array, shape (d,)
          System lengths.
  
      Returns
      -------
      csep: float array-like, shape (n, m, d)
          csep[i, j] is the separation vector from point j to point i.
          Note the un-intuitive vector direction.
      """
      seps = ra[:, np.newaxis, :] - rb[np.newaxis, :, :]
      for i_dim in range(ra.shape[1]):
          seps_dim = seps[:, :, i_dim]
          seps_dim[seps_dim > L[i_dim] / 2.0] -= L[i_dim]
          seps_dim[seps_dim < -L[i_dim] / 2.0] += L[i_dim]
      return seps