docs: def propagate(self, date):
          """Compute state of orbit at a given date, past or future
  
          Args:
              date (Date)
          Return:
              Orbit:
          """
  
          i0, Ω0, e0, ω0, M0, n0 = self.tle
          n0 *= 60  # conversion to min⁻¹
          if isinstance(date, Date):
              t0 = self.tle.date.datetime
              tdiff = (date.datetime - t0).total_seconds() / 60.
          elif isinstance(date, timedelta):
              tdiff = date.total_seconds() / 60.
              date = self.tle.date + date
          else:
              raise TypeError("Unhandled type for 'date': %s" % type(date))
  
          bstar = self.tle.complements['bstar']
          µ = self.gravity.µ_e
          r_e = self.gravity.r_e
          k_e = self.gravity.k_e
  
          # retrieve initialized variables
          _i = self._init
          n0 = _i.n0
  
          Mdf = M0 + _i.Mdot * n0 * tdiff
          ωdf = ω0 + _i.ωdot * n0 * tdiff
          Ωdf = Ω0 + _i.Ωdot * n0 * tdiff
  
          delta_ω = bstar * _i.C3 * cos(ω0) * tdiff
          delta_M = 0.
          if e0 > 1e-4:
              delta_M = - 2 / 3 * (_i.q0 - _i.s) ** 4 * bstar * _i.ξ ** 4 / (e0 * _i.η) * ((1 + _i.η * cos(Mdf)) ** 3 - (1 + _i.η * cos(M0)) ** 3)
  
          Mp = (Mdf + delta_ω + delta_M) % (2 * np.pi)
          ω = ωdf - delta_ω - delta_M
          Ω = Ωdf - 21 * n0 * _i.k2 * _i.θ / (2 * _i.a0 ** 2 * _i.β_0 ** 2) * _i.C1 * tdiff ** 2
          e = e0 - bstar * _i.C4 * tdiff - bstar * _i.C5 * (sin(Mp) - sin(M0))
  
          if e < 1e-6:
              e = 1e-6
  
          a = _i.a0 * (1 - _i.C1 * tdiff - _i.D2 * tdiff ** 2 - _i.D3 * tdiff ** 3 - _i.D4 * tdiff ** 4) ** 2
  
          L = Mp + ω + Ω + n0 * (3 / 2 * _i.C1 * tdiff ** 2 + (_i.D2 + 2 * _i.C1 ** 2) * tdiff ** 3 + 1 / 4 * (3 * _i.D3 + 12 * _i.C1 * _i.D2 + 10 * _i.C1 ** 3) * tdiff ** 4 + 1 / 5 * (3 * _i.D4 + 12 * _i.C1 * _i.D3 + 6 * _i.D2 ** 2 + 30 * _i.C1 ** 2 * _i.D2 + 15 * _i.C1 ** 4) * tdiff ** 5)
  
          β = sqrt(1 - e ** 2)
          n = µ / (a ** (3 / 2))
  
          # Long-period terms
          axN = e * cos(ω)
          ayNL = _i.A30 * sin(i0) / (4 * _i.k2 * a * β ** 2)
          tmp = (1 + _i.θ) if (1 + _i.θ) > 1.5e-12 else 1.5e-12
          L_L = ayNL / 2 * axN * ((3 + 5 * _i.θ) / tmp)
  
          L_T = L + L_L
          ayN = e * sin(ω) + ayNL
  
          # Resolving of kepler equation
          U = (L_T - Ω) % (2 * np.pi)
          Epω = U
          for xxx in range(10):
              delta_Epω = (U - ayN * cos(Epω) + axN * sin(Epω) - Epω) / (1 - ayN * sin(Epω) - axN * cos(Epω))
              if abs(delta_Epω) < 1e-12:
                  break
              Epω = Epω + delta_Epω
  
          # Short-period terms
          ecosE = axN * cos(Epω) + ayN * sin(Epω)
          esinE = axN * sin(Epω) - ayN * cos(Epω)
          e_L = sqrt(axN ** 2 + ayN ** 2)
          p_L = a * (1 - e_L ** 2)
          r = a * (1 - ecosE)
          rdot = sqrt(a) / r * esinE
          rfdot = sqrt(p_L) / r
  
          cosu = a / r * (cos(Epω) - axN + ayN * esinE / (1 + sqrt(1 - e_L ** 2)))
          sinu = a / r * (sin(Epω) - ayN - axN * esinE / (1 + sqrt(1 - e_L ** 2)))
          u = arctan2(sinu, cosu)
  
          Delta_r = _i.k2 / (2 * p_L) * (1 - _i.θ ** 2) * cos(2 * u)
          Delta_u = - _i.k2 / (4 * p_L ** 2) * (7 * _i.θ ** 2 - 1) * sin(2 * u)
          Delta_Ω = 3 * _i.k2 * _i.θ / (2 * p_L ** 2) * sin(2 * u)
          Delta_i = 3 * _i.k2 * _i.θ / (2 * p_L ** 2) * sin(i0) * cos(2 * u)
          Delta_rdot = - n * _i.k2 * (1 - _i.θ ** 2) * sin(2 * u) / (p_L * µ)
          Delta_rfdot = _i.k2 * n * ((1 - _i.θ ** 2) * cos(2 * u) - 3 / 2 * (1 - 3 * _i.θ ** 2)) / (p_L * µ)
  
          rk = r * (1 - 3 / 2 * _i.k2 * sqrt(1 - e_L ** 2) / (p_L ** 2) * (3 * _i.θ ** 2 - 1)) + Delta_r
          uk = u + Delta_u
          Ωk = Ω + Delta_Ω
          ik = i0 + Delta_i
          rdotk = rdot + Delta_rdot
          rfdotk = rfdot + Delta_rfdot
  
          # Vectors
          vM = np.array([- sin(Ωk) * cos(ik), cos(Ωk) * cos(ik), sin(ik)])
          vN = np.array([cos(Ωk), sin(Ωk), 0])
  
          vU = vM * sin(uk) + vN * cos(uk)
          vV = vM * cos(uk) - vN * sin(uk)
  
          vR = rk * vU * r_e
          vRdot = (rdotk * vU + rfdotk * vV) * (r_e * k_e / 60.)
  
          vector = np.concatenate((vR, vRdot)) * 1000  # conversion to meters
  
          return self.tle.__class__(date, vector, 'cartesian', 'TEME', self.__class__(), **self.tle.complements)