How to use the poliastro.bodies.Earth.k.to function in poliastro

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)
            incs.append(angles[0].value)
            raans.append(angles[1].value)
            argps.append(angles[2].value)

        # averaging over 5 last values in the way Curtis does
        inc_f, raan_f, argp_f = (
            np.mean(incs[-5:]),
            np.mean(raans[-5:]),
            np.mean(argps[-5:]),
        )

        assert_quantity_allclose(
            [
                (raan_f * u.rad).to(u.deg) - test_params["orbit"][3],

def test_propagation_parabolic(propagator):
    # example from Howard Curtis (3rd edition), section 3.5, problem 3.15
    # TODO: add parabolic solver in some parabolic propagators, refer to #417
    if propagator in [mikkola, gooding]:
        pytest.skip()

    p = 2.0 * 6600 * u.km
    _a = 0.0 * u.deg
    orbit = Orbit.parabolic(Earth, p, _a, _a, _a, _a)
    orbit = orbit.propagate(0.8897 / 2.0 * u.h, method=propagator)

    _, _, _, _, _, nu0 = rv2coe(
        Earth.k.to(u.km ** 3 / u.s ** 2).value,
        orbit.r.to(u.km).value,
        orbit.v.to(u.km / u.s).value,
    )
    assert_quantity_allclose(nu0, np.deg2rad(90.0), rtol=1e-4)

    orbit = Orbit.parabolic(Earth, p, _a, _a, _a, _a)
    orbit = orbit.propagate(36.0 * u.h, method=propagator)
    assert_quantity_allclose(norm(orbit.r), 304700.0 * u.km, rtol=1e-4)

def test_atmospheric_drag_exponential():
    # http://farside.ph.utexas.edu/teaching/celestial/Celestialhtml/node94.html#sair (10.148)
    # given the expression for \dot{r} / r, aproximate \Delta r \approx F_r * \Delta t

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

    # parameters of a circular orbit with h = 250 km (any value would do, but not too small)
    orbit = Orbit.circular(Earth, 250 * u.km)
    r0, _ = orbit.rv()
    r0 = r0.to(u.km).value

    # 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
    B = C_D * A_over_m

    # parameters of the atmosphere
    rho0 = rho0_earth.to(u.kg / u.km ** 3).value  # kg/km^3
    H0 = H0_earth.to(u.km).value  # km

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 custom(initial,choice,state,time,f_time):                                          # requires modification and validation
    R = Earth.R.to(u.km).value
    k = Earth.k.to(u.km ** 3 / u.s ** 2).value
    #s0 = Orbit.from_classical(Earth, x[0] * u.km, x[1] * u.one, x[4] * u.deg, * u.deg, x[3] * u.deg, 0 * u.deg, epoch=Time(0, format='jd', scale='tdb'))

    #orbit = Orbit.circular(
    #    Earth, 250 * u.km, epoch=Time(0.0, format="jd", scale="tdb"))

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

    # parameters of the atmosphere
    rho0 = Earth.rho0.to(u.kg / u.km ** 3).value  # kg/km^3
    H0 = Earth.H0.to(u.km).value

    #J2,J3 parameters

def atmd_earth(initial,time,f_time):                                               # perturbation for atmospheric drag
    R = Earth.R.to(u.km).value
    k = Earth.k.to(u.km ** 3 / u.s ** 2).value
    #s0 = Orbit.from_classical(Earth, x[0] * u.km, x[1] * u.one, x[4] * u.deg, * u.deg, x[3] * u.deg, 0 * u.deg, epoch=Time(0, format='jd', scale='tdb'))

    #orbit = Orbit.circular(
    #    Earth, 250 * u.km, epoch=Time(0.0, format="jd", scale="tdb"))

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

    # parameters of the atmosphere
    rho0 = Earth.rho0.to(u.kg / u.km ** 3).value  # kg/km^3
    H0 = Earth.H0.to(u.km).value
    
    print("Use default constants or you want to customize?\n1.Default.\n2.Custom.")

def j3_earth(initial,time,f_time):                                                   # perturbation for oblateness
    R = Earth.R.to(u.km).value
    k = Earth.k.to(u.km ** 3 / u.s ** 2).value
    #s0 = Orbit.from_classical(Earth, x[0] * u.km, x[1] * u.one, x[4] * u.deg, * u.deg, x[3] * u.deg, 0 * u.deg, epoch=Time(0, format='jd', scale='tdb'))

    #orbit = Orbit.circular(
    #    Earth, 250 * u.km, epoch=Time(0.0, format="jd", scale="tdb"))

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

    J3=Earth.J3.value
    print("Use default constants or you want to customize?\n1.Default.\n2.Custom.")
    check=input()
    if(check=='1'):
        pass

def j2_earth(initial,time,f_time):                                                  # perturbation for oblateness
    R = Earth.R.to(u.km).value
    k = Earth.k.to(u.km ** 3 / u.s ** 2).value
    #s0 = Orbit.from_classical(Earth, x[0] * u.km, x[1] * u.one, x[4] * u.deg, * u.deg, x[3] * u.deg, 0 * u.deg, epoch=Time(0, format='jd', scale='tdb'))

    #orbit = Orbit.circular(
    #    Earth, 250 * u.km, epoch=Time(0.0, format="jd", scale="tdb"))

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

    J2=Earth.J2.value
    print("Use default constants or you want to customize?\n1.Default.\n2.Custom.")
    check=input()
    if(check=='1'):
        pass