How to use the poliastro.twobody.propagation.cowell function in poliastro

def test_cowell_propagation_with_zero_acceleration_equals_kepler():
    # Data from Vallado, example 2.4

    r0 = np.array([1131.340, -2282.343, 6672.423]) * u.km
    v0 = np.array([-5.64305, 4.30333, 2.42879]) * u.km / u.s
    tofs = [40 * 60.0] * u.s

    orbit = Orbit.from_vectors(Earth, r0, v0)

    expected_r = np.array([-4219.7527, 4363.0292, -3958.7666]) * u.km
    expected_v = np.array([3.689866, -1.916735, -6.112511]) * u.km / u.s

    r, v = cowell(Earth.k, orbit.r, orbit.v, tofs, ad=None)

    assert_quantity_allclose(r[0], expected_r, rtol=1e-5)
    assert_quantity_allclose(v[0], expected_v, rtol=1e-4)

def test_cowell_propagation_circle_to_circle():
    # From [Edelbaum, 1961]
    accel = 1e-7

    def constant_accel(t0, u, k):
        v = u[3:]
        norm_v = (v[0] ** 2 + v[1] ** 2 + v[2] ** 2) ** 0.5
        return accel * v / norm_v

    ss = Orbit.circular(Earth, 500 * u.km)
    tofs = [20] * ss.period

    r, v = cowell(Earth.k, ss.r, ss.v, tofs, ad=constant_accel)

    ss_final = Orbit.from_vectors(Earth, r[0], v[0])

    da_a0 = (ss_final.a - ss.a) / ss.a
    dv_v0 = abs(norm(ss_final.v) - norm(ss.v)) / norm(ss.v)
    assert_quantity_allclose(da_a0, 2 * dv_v0, rtol=1e-2)

    dv = abs(norm(ss_final.v) - norm(ss.v))
    accel_dt = accel * u.km / u.s ** 2 * tofs[0]
    assert_quantity_allclose(dv, accel_dt, rtol=1e-2)

def test_J2_propagation_Earth():
    # from Curtis example 12.2:
    r0 = np.array([-2384.46, 5729.01, 3050.46])  # km
    v0 = np.array([-7.36138, -2.98997, 1.64354])  # km/s

    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

t_decay = 7.17 * u.d

    # parameters of a body
    C_D = 2.2  # dimentionless (any value would do)
    A_over_m = ((np.pi / 4.0) * (u.m ** 2) / (100 * u.kg)).to_value(
        u.km ** 2 / u.kg
    )  # km^2/kg

    tofs = [365] * u.d

    lithobrake_event = LithobrakeEvent(R)
    events = [lithobrake_event]

    coesa76 = COESA76()

    rr, _ = cowell(
        Earth.k,
        orbit.r,
        orbit.v,
        tofs,
        ad=atmospheric_drag_model,
        R=R,
        C_D=C_D,
        A_over_m=A_over_m,
        model=coesa76,
        events=events,
    )

    assert_quantity_allclose(norm(rr[0].to(u.km).value), R, atol=1)  # below 1km

    assert_quantity_allclose(lithobrake_event.last_t, t_decay, rtol=1e-2)

def test_elliptic_near_parabolic(ecc, propagator):
    # 'kepler fails if really close to parabolic'. Refer to issue #714.
    if propagator in [vallado] and ecc > 0.99:
        pytest.xfail()

    _a = 0.0 * u.rad
    tof = 1.0 * u.min
    ss0 = Orbit.from_classical(
        Earth, 10000 * u.km, ecc * u.one, _a, _a, _a, 1.0 * u.rad
    )

    ss_cowell = ss0.propagate(tof, method=cowell)
    ss_propagator = ss0.propagate(tof, method=propagator)

    assert_quantity_allclose(ss_propagator.r, ss_cowell.r)
    assert_quantity_allclose(ss_propagator.v, ss_cowell.v)

epoch = Time(j_date, format="jd", scale="tdb")

        initial = Orbit.from_classical(
            Earth,
            10085.44 * u.km,
            0.025422 * u.one,
            88.3924 * u.deg,
            45.38124 * u.deg,
            227.493 * u.deg,
            343.4268 * u.deg,
            epoch=epoch,
        )
        # in Curtis, the mean distance to Sun is used. In order to validate against it, we have to do the same thing
        sun_normalized = functools.partial(normalize_to_Curtis, sun_r=sun_r)

        rr, vv = cowell(
            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)

Earth.k,
        orbit.r,
        orbit.v,
        tofs,
        ad=J2_perturbation,
        J2=Earth.J2.value,
        R=Earth.R.to(u.km).value,
        rtol=1e-8,
    )

    def a_J2J3(t0, u_, k_):
        j2 = J2_perturbation(t0, u_, k_, J2=Earth.J2.value, R=Earth.R.to(u.km).value)
        j3 = J3_perturbation(t0, u_, k_, J3=Earth.J3.value, R=Earth.R.to(u.km).value)
        return j2 + j3

    r_J3, v_J3 = cowell(Earth.k, orbit.r, orbit.v, tofs, ad=a_J2J3, rtol=1e-8)

    a_values_J2 = np.array(
        [
            rv2coe(k, ri, vi)[0] / (1.0 - rv2coe(k, ri, vi)[1] ** 2)
            for ri, vi in zip(r_J2.to(u.km).value, v_J2.to(u.km / u.s).value)
        ]
    )
    a_values_J3 = np.array(
        [
            rv2coe(k, ri, vi)[0] / (1.0 - rv2coe(k, ri, vi)[1] ** 2)
            for ri, vi in zip(r_J3.to(u.km).value, v_J3.to(u.km / u.s).value)
        ]
    )
    da_max = np.max(np.abs(a_values_J2 - a_values_J3))

    ecc_values_J2 = np.array(

def test_3rd_body_Curtis(test_params):
    # based on example 12.11 from Howard Curtis
    body = test_params["body"]
    with solar_system_ephemeris.set("builtin"):
        j_date = 2454283.0 * u.day
        tof = (test_params["tof"]).to(u.s).value
        body_r = build_ephem_interpolant(
            body,
            test_params["period"],
            (j_date, j_date + test_params["tof"]),
            rtol=1e-2,
        )

        epoch = Time(j_date, format="jd", scale="tdb")
        initial = Orbit.from_classical(Earth, *test_params["orbit"], epoch=epoch)
        rr, vv = cowell(
            Earth.k,
            initial.r,
            initial.v,
            np.linspace(0, tof, 400) * u.s,
            rtol=1e-10,
            ad=third_body,
            k_third=body.k.to(u.km ** 3 / u.s ** 2).value,
            perturbation_body=body_r,
        )

        incs, raans, argps = [], [], []
        for ri, vi in zip(rr.to(u.km).value, vv.to(u.km / u.s).value):
            angles = Angle(
                rv2coe(Earth.k.to(u.km ** 3 / u.s ** 2).value, ri, vi)[2:5] * u.rad
            )  # inc, raan, argp
            angles = angles.wrap_at(180 * u.deg)

return np.array([0, 0, 0])

        if gravity is EarthGravity.J2:
            perturbations[J2_perturbation] = {
                "J2": Earth.J2.value,
                "R": Earth.R.to(u.km).value,
            }
        if atmosphere is not None:
            perturbations[atmospheric_drag_model] = {
                "R": Earth.R.to(u.km).value,
                "C_D": self.spacecraft.C_D,
                "A_over_m": (self.spacecraft.A / self.spacecraft.m),
                "model": atmosphere,
            }
        ad_kwargs.update(perturbations=perturbations)
        new_orbit = self.orbit.propagate(value=tof, method=cowell, ad=ad, **ad_kwargs)
        return EarthSatellite(new_orbit, self.spacecraft)