docs: def interpolate(self, date, method=None, order=None):
          """Interpolate data at a given date
  
          Args:
              date (Date):
              method (str): Method of interpolation to use
              order (int): In case of ``LAGRANGE`` method is used
          Return:
              Orbit:
          """
  
          if not self.start <= date <= self.stop:
              raise ValueError("Date '%s' not in range" % date)
  
          prev_idx = 0
          ephem = self
  
          # Binary search of the orbit step just before the desired date
          while True:
              idx = len(ephem)
              if idx == 1:
                  break
              k = idx // 2
  
              if date > ephem[k].date:
                  prev_idx += k
                  ephem = ephem[k:]
              else:
                  ephem = ephem[:k]
  
          method = method if method is not None else self.method
          order = order if order is not None else self.order
  
          if method == self.LINEAR:
  
              y0 = self[prev_idx]
              y1 = self[prev_idx + 1]
  
              result = y0[:] + (y1[:] - y0[:]) * (date.mjd - y0.date.mjd) / (y1.date.mjd - y0.date.mjd)
  
          elif method == self.LAGRANGE:
  
              stop = prev_idx + 1 + order // 2 + order % 2
              start = prev_idx - order // 2 + 1
              if stop >= len(self):
                  start -= stop - len(self)
              elif start < 0:
                  stop -= start
                  start = 0
  
              # selection of the subset of data, of length 'order' around the desired value
              subset = self[start:stop]
              date_subset = np.array([x.date.mjd for x in subset])
  
              result = np.zeros(6)
  
              # Everything is on wikipedia
              #        k
              # L(x) = Σ y_j * l_j(x)
              #        j=0
              #
              # l_j(x) = Π (x - x_m) / (x_j - x_m)
              #     0 <= m <= k
              #        m != j
              for j in range(order):
                  # This mask is here to enforce the m != j in the lagrange polynomials
                  mask = date_subset != date_subset[j]
                  l_j = (date.mjd - date_subset[mask]) / (date_subset[j] - date_subset[mask])
                  result = result + l_j.prod() * subset[j]
  
          else:
              raise ValueError("Unkown interpolation method", method)
  
          orb = ephem[0]
  
          return orb.__class__(date, result, orb.form, orb.frame, orb.propagator)