How to use the exoplanet.orbits.KeplerianOrbit function in exoplanet

def test_ylm_phase():
    texp = 0.05
    map = starry.Map(ydeg=2)
    np.random.seed(11)
    map[1:, :] = 0.1 * np.random.randn(8)
    theta = np.linspace(0, 360, 10000)
    t = np.linspace(-0.2, 0.2, 10000)
    orbit = exo.orbits.KeplerianOrbit(period=1.0)
    flux = map.flux(theta=theta)
    window = int(texp / (t[1] - t[0]))
    fluence_mavg = moving_average(flux, window)
    fluence_starry = map.flux(t=t, orbit=orbit, theta=theta, 
                              texp=texp, oversample=50).eval()

    # The error is primarily coming from our moving average
    # integrator, so let's be lenient
    f1 = fluence_mavg[window:-window]
    f2 = fluence_starry[window:-window]
    assert np.allclose(f1, f2, atol=1e-4, rtol=1e-4)

dt = np.linspace(-0.5, 0.5, oversample)
                stencil[1:-1:2] = 4
                stencil[2:-1:2] = 2
            else:
                raise ValueError("Parameter `order` must be <= 2")
            stencil /= np.sum(stencil)

            if texp.ndim == 0:
                dt = texp * dt
            else:
                dt = tt.shape_padright(texp) * dt
            t = tt.shape_padright(t) + dt
            t = tt.reshape(t, (-1,))

        # Compute the relative positions of all bodies
        orbit = exoplanet.orbits.KeplerianOrbit(
            period=sec_porb,
            t0=sec_t0,
            incl=sec_iorb,
            ecc=sec_ecc,
            omega=sec_w,
            Omega=sec_Omega,
            m_planet=sec_m,
            m_star=pri_m,
            r_star=pri_r,
        )
        try:
            x, y, z = orbit.get_relative_position(
                t, light_delay=self.light_delay
            )
        except TypeError:
            if self.light_delay:

# The amlitudes should be sorted
    pm.Potential("logK_order", tt.switch(logK[1:] > logK[:-1], -np.inf, 0.0))

    # We also want to keep period physical but this probably won't be hit
    pm.Potential("P_bound", tt.switch(P <= 0, -np.inf, 0.0))

    # Eccentricity & argument of periasteron
    ecc = pm.Uniform("ecc", lower=0, upper=0.99, shape=N_pl, testval=eccs)
    omega = xo.distributions.Angle("omega", shape=N_pl, testval=omegas)

    # Jitter & a quadratic RV trend
    # logs = pm.Normal("logs", mu=np.log(np.median(yerr)), sd=5.0)
    trend = pm.Normal("trend", mu=0, sd=10.0 ** -np.arange(3)[::-1], shape=3)

    # Set up the orbit
    orbit = xo.orbits.KeplerianOrbit(period=P, t0=t0, ecc=ecc, omega=omega)

    # Set up the RV model and save it as a deterministic
    # for plotting purposes later
    vrad = orbit.get_radial_velocity(x, K=tt.exp(logK))
    if N_pl == 1:
        vrad = vrad[:, None]

    # Define the background model
    A = np.vander(x - 0.5 * (x.min() + x.max()), 3)
    bkg = tt.dot(A, trend)

    # Sum over planets and add the background to get the full model
    rv_model = tt.sum(vrad, axis=-1) + bkg

    # Simulate the data
    y_true = xo.eval_in_model(rv_model)

sec_porb,
        sec_ecc,
        sec_w,
        sec_Omega,
        sec_iorb,
        sec_inc,
        sec_obl,
        sec_y,
        sec_u,
        sec_f,
        sec_alpha,
        sec_sigr,
    ):
        """Render all of the bodies in the system."""
        # Compute the relative positions of all bodies
        orbit = exoplanet.orbits.KeplerianOrbit(
            period=sec_porb,
            t0=sec_t0,
            incl=sec_iorb,
            ecc=sec_ecc,
            omega=sec_w,
            Omega=sec_Omega,
            m_planet=sec_m,
            m_star=pri_m,
            r_star=pri_r,
        )
        try:
            x, y, z = orbit.get_relative_position(
                t, light_delay=self.light_delay
            )
        except TypeError:
            if self.light_delay: