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

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

def test_propagate_to_anomaly_in_the_past_fails_for_open_orbits():
    r0 = [Earth.R.to(u.km).value + 300, 0, 0] * u.km
    v0 = [0, 15, 0] * u.km / u.s
    orb = Orbit.from_vectors(Earth, r0, v0)

    with pytest.raises(ValueError, match="True anomaly -0.02 rad not reachable"):
        orb.propagate_to_anomaly(orb.nu - 1 * u.deg)

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

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

    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)

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

if perturbations:
                return np.sum(
                    [f(t0=t0, state=state, k=k, **p) for f, p in perturbations.items()],
                    axis=0,
                )
            else:
                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)

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'):