How to use the poliastro.core.elements.rv2coe function in poliastro

def test_convert_between_coe_and_rv_is_transitive(classical):
    k = Earth.k.to(u.km ** 3 / u.s ** 2).value  # u.km**3 / u.s**2
    res = rv2coe(k, *coe2rv(k, *classical))
    assert_allclose(res, classical)

def test_convert_coe_and_rv_circular_equatorial(circular_equatorial):
    k, expected_res = circular_equatorial
    res = rv2coe(k, *coe2rv(k, *expected_res))
    assert_allclose(res, expected_res, atol=1e-8)

def test_convert_coe_and_rv_hyperbolic(hyperbolic):
    k, expected_res = hyperbolic
    res = rv2coe(k, *coe2rv(k, *expected_res))
    assert_allclose(res, expected_res, atol=1e-8)

orbit = Orbit.from_vectors(Earth, r0 * u.km, v0 * u.km / u.s)

    tofs = [48.0] * u.h
    rr, vv = cowell(
        Earth.k,
        orbit.r,
        orbit.v,
        tofs,
        ad=J2_perturbation,
        J2=Earth.J2.value,
        R=Earth.R.to(u.km).value,
    )

    k = Earth.k.to(u.km ** 3 / u.s ** 2).value

    _, _, _, raan0, argp0, _ = rv2coe(k, r0, v0)
    _, _, _, raan, argp, _ = rv2coe(k, rr[0].to(u.km).value, vv[0].to(u.km / u.s).value)

    raan_variation_rate = (raan - raan0) / tofs[0].to(u.s).value  # type: ignore
    argp_variation_rate = (argp - argp0) / tofs[0].to(u.s).value  # type: ignore

    raan_variation_rate = (raan_variation_rate * u.rad / u.s).to(u.deg / u.h)
    argp_variation_rate = (argp_variation_rate * u.rad / u.s).to(u.deg / u.h)

    assert_quantity_allclose(raan_variation_rate, -0.172 * u.deg / u.h, rtol=1e-2)
    assert_quantity_allclose(argp_variation_rate, 0.282 * u.deg / u.h, rtol=1e-2)

Earth.k,
            initial.r,
            initial.v,
            np.linspace(0, (tof).to(u.s).value, 4000) * u.s,
            rtol=1e-8,
            ad=radiation_pressure,
            R=Earth.R.to(u.km).value,
            C_R=2.0,
            A_over_m=2e-4 / 100,
            Wdivc_s=Wdivc_sun.value,
            star=sun_normalized,
        )

        delta_eccs, delta_incs, delta_raans, delta_argps = [], [], [], []
        for ri, vi in zip(rr.to(u.km).value, vv.to(u.km / u.s).value):
            orbit_params = rv2coe(Earth.k.to(u.km ** 3 / u.s ** 2).value, ri, vi)
            delta_eccs.append(orbit_params[1] - initial.ecc.value)
            delta_incs.append(
                (orbit_params[2] * u.rad).to(u.deg).value - initial.inc.value
            )
            delta_raans.append(
                (orbit_params[3] * u.rad).to(u.deg).value - initial.raan.value
            )
            delta_argps.append(
                (orbit_params[4] * u.rad).to(u.deg).value - initial.argp.value
            )

        # averaging over 5 last values in the way Curtis does
        index = int(
            1.0 * t_days / tof.to(u.day).value * 4000  # type: ignore
        )
        delta_ecc, delta_inc, delta_raan, delta_argp = (

Otherwise, if time is provided, propagate this `Orbit` some `time` and return the result.

        Parameters
        ----------
        value : Multiple options
            True anomaly values or time values. If given an angle, it will always propagate forward.
        rtol : float, optional
            Relative tolerance for the propagation algorithm, default to 1e-10.
        method : function, optional
            Method used for propagation
        **kwargs
            parameters used in perturbation models

        """
        if hasattr(value, "unit") and value.unit in ('rad', 'deg'):
            p, ecc, inc, raan, argp, _ = rv2coe(self.attractor.k.to(u.km ** 3 / u.s ** 2).value,
                                                self.r.to(u.km).value,
                                                self.v.to(u.km / u.s).value)

            # Compute time of flight for correct epoch
            M = nu_to_M(self.nu, self.ecc)
            new_M = nu_to_M(value, self.ecc)
            time_of_flight = Angle(new_M - M).wrap_at(360 * u.deg) / self.n

            return self.from_classical(self.attractor, p / (1.0 - ecc ** 2) * u.km,
                                       ecc * u.one, inc * u.rad, raan * u.rad,
                                       argp * u.rad, value,
                                       epoch=self.epoch + time_of_flight, plane=self._plane)
        else:
            if isinstance(value, time.Time) and not isinstance(value, time.TimeDelta):
                time_of_flight = value - self.epoch
            else:

rtol: float
        This method does not require of tolerance since it is non iterative.

    Returns
    -------
    rr : ~astropy.units.Quantity
        Propagated position vectors.
    vv : ~astropy.units.Quantity

    Note
    ----
    Original paper: https://doi.org/10.1007/BF01235850
    """

    # Solving for the classical elements
    p, ecc, inc, raan, argp, nu = rv2coe(k, r0, v0)
    M0 = nu_to_M(nu, ecc, delta=0)
    a = p / (1 - ecc ** 2)
    n = np.sqrt(k / np.abs(a) ** 3)
    M = M0 + n * tof

    # Solve for specific geometrical case
    if ecc < 1.0:
        # Equation (9a)
        alpha = (1 - ecc) / (4 * ecc + 1 / 2)
    else:
        alpha = (ecc - 1) / (4 * ecc + 1 / 2)

    beta = M / 2 / (4 * ecc + 1 / 2)

    # Equation (9b)
    if beta >= 0:

Relative error for accuracy of the method.

    Returns
    -------
    rr : 1x3 vector
        Propagated position vectors.
    vv : 1x3 vector

    Note
    ----
    This algorithm was developed by Danby in his paper *The solution of Kepler
    Equation* with DOI: https://doi.org/10.1007/BF01686811
    """

    # Solve first for eccentricity and mean anomaly
    p, ecc, inc, raan, argp, nu = rv2coe(k, r0, v0)
    M0 = nu_to_M(nu, ecc, delta=0)
    semi_axis_a = p / (1 - ecc ** 2)
    n = np.sqrt(k / np.abs(semi_axis_a) ** 3)
    M = M0 + n * tof

    # Range mean anomaly
    xma = M - 2 * np.pi * np.floor(M / 2 / np.pi)

    if ecc == 0:
        # Solving for circular orbit
        nu = xma
        return coe2rv(k, p, ecc, inc, raan, argp, nu)

    elif ecc < 1.0:
        # For elliptical orbit
        E = xma + 0.85 * np.sign(np.sin(xma)) * ecc

Time of flight.

    Returns
    -------
    rf: array
        Final position vector
    vf: array
        Final velocity vector

    Note
    ----
    The following algorithm was taken from http://dx.doi.org/10.1007/BF00691917.
    """

    # Solve first for eccentricity and mean anomaly
    p, ecc, inc, raan, argp, nu = rv2coe(k, r0, v0)

    M0 = nu_to_M(nu, ecc, delta=0)
    a = p / (1 - ecc ** 2)
    n = np.sqrt(k / a ** 3)
    M = M0 + n * tof

    # Range between -pi and pi
    M = M % (2 * np.pi)
    if M > np.pi:
        M = -(2 * np.pi - M)

    # Equation (20)
    alpha = (3 * np.pi ** 2 + 1.6 * (np.pi - np.abs(M)) / (1 + ecc)) / (np.pi ** 2 - 6)

    # Equation (5)
    d = 3 * (1 - ecc) + alpha * ecc