How to use the exoplanet.orbits function in exoplanet

def test_ylm():
    map = starry.Map(ydeg=1)
    map[1, :] = [0.3, 0.2, 0.1]
    orbit = exo.orbits.KeplerianOrbit(period=1.0, m_star=1.0, r_star=1.0)
    t = np.linspace(-0.25, 0.25, 100)
    theta = 30.

    # Compute the whole light curve with Theano
    f1 = map.flux(t=t, orbit=orbit, ro=0.1, theta=theta, use_in_transit=False).eval()

    # Compute just the transit with Theano
    f2 = map.flux(t=t, orbit=orbit, ro=0.1, theta=theta, use_in_transit=True).eval()

    # Compute the whole light curve without Theano
    coords = orbit.get_relative_position(t)
    xo = (coords[0] / orbit.r_star).eval()
    yo = (coords[1] / orbit.r_star).eval()
    zo = -(coords[2] / orbit.r_star).eval()
    f3 = map.flux(xo=xo, yo=yo, zo=zo, ro=0.1, theta=theta)

def test_ylm_occ():
    texp = 0.05
    map = starry.Map(ydeg=2)
    np.random.seed(11)
    map[1:, :] = 0.1 * np.random.randn(8)
    orbit = exo.orbits.KeplerianOrbit(period=1.0, m_star=1.0, r_star=1.0)
    t = np.linspace(-0.2, 0.2, 10000)
    flux = map.flux(t=t, orbit=orbit, ro=0.1).eval()
    xo = orbit.get_relative_position(t)[0].eval()
    yo = orbit.get_relative_position(t)[1].eval()
    flux = map.flux(xo=xo, yo=yo, ro=0.1)
    fluence_mavg = moving_average(flux, int(texp / (t[1] - t[0])))
    fluence_starry = map.flux(t=t, orbit=orbit, ro=0.1, 
                              texp=texp, oversample=30).eval()
    fluence_starry_vec = map.flux(t=t, orbit=orbit, ro=0.1, 
                              texp=np.ones_like(t) * texp, oversample=30).eval()
    assert np.allclose(fluence_mavg, fluence_starry, fluence_starry_vec)

w=60,
        length_unit=u.Rsun,
        mass_unit=u.Msun,
        angle_unit=u.degree,
        time_unit=u.day,
    )

    # Define the system
    sys = starry.System(A, b)

    # Compute with starry
    time = np.linspace(-0.5, 0.5, 1000)
    rv1 = sys.rv(time, keplerian=True)

    # Compute with exoplanet
    orbit = exoplanet.orbits.KeplerianOrbit(
        period=1.0,
        t0=0.0,
        incl=86.0 * np.pi / 180,
        ecc=0.3,
        omega=60 * np.pi / 180,
        m_planet=0.01,
        m_star=1.0,
        r_star=1.0,
    )
    rv2 = orbit.get_radial_velocity(time).eval()

    assert np.allclose(rv1, rv2)

def test_ld():
    texp = 0.05
    map = starry.Map(udeg=2)
    map[1:] = [0.4, 0.26]
    orbit = exo.orbits.KeplerianOrbit(period=1.0, m_star=1.0, r_star=1.0)
    t = np.linspace(-0.2, 0.2, 10000)
    flux = map.flux(t=t, orbit=orbit, ro=0.1).eval()
    fluence_mavg = moving_average(flux, int(texp / (t[1] - t[0])))
    fluence_starry = map.flux(t=t, orbit=orbit, ro=0.1, 
                              texp=texp, oversample=30).eval()
    fluence_starry_vec = map.flux(t=t, orbit=orbit, ro=0.1, 
                              texp=np.ones_like(t) * texp, oversample=30).eval()
    assert np.allclose(fluence_mavg, fluence_starry, fluence_starry_vec)

def test_ld():
    map = starry.Map(udeg=2)
    map[1:] = [0.4, 0.26]
    orbit = exo.orbits.KeplerianOrbit(period=1.0, m_star=1.0, r_star=1.0)
    t = np.linspace(-0.25, 0.25, 100)

    # Compute the whole light curve with Theano
    f1 = map.flux(t=t, orbit=orbit, ro=0.1, use_in_transit=False).eval()

    # Compute just the transit with Theano
    f2 = map.flux(t=t, orbit=orbit, ro=0.1, use_in_transit=True).eval()

    # Compute the whole light curve without Theano
    coords = orbit.get_relative_position(t)
    xo = (coords[0] / orbit.r_star).eval()
    yo = (coords[1] / orbit.r_star).eval()
    b = np.sqrt(xo * xo + yo * yo)
    zo = -(coords[2] / orbit.r_star).eval()
    f3 = map.flux(b=b, zo=zo, ro=0.1)

import exoplanet
from packaging import version
import theano.tensor as tt
import numpy as np
from ..ops import autocompile


# NOTE: In version 0.1.7, DFM changed the coordinates
# so that the z-axis points TOWARD the observer!
if version.parse(exoplanet.__version__) > version.parse('0.1.7.dev0'):
    z_sign = 1
else:
    z_sign = -1


class KeplerianOrbit(exoplanet.orbits.KeplerianOrbit):
    """
    A wrapper around `exoplanet.orbits.KeplerianOrbit` that
    plays nice with `starry`. Refer to the docs of that class
    for all accepted keywords. In addition to those, this class
    accepts the following keyword arguments:

    Args:
        r_planet: The radius of the planet in ``R_sun``. Default is 
            the radius of the Earth.
        rot_period: The period of rotation of the planet in days.
            Default ``1.0``. Set to ``None`` to disable rotation.
        theta0: The rotational phase in degrees at ``t=t0``.
            Default ``0.0``
        lazy: 

    """

# The time of a reference transit for each planet
    t0 = pm.Normal("t0", mu=t0_true, sd=1.0)
    
    # The log period; also tracking the period itself
    logP = pm.Normal("logP", mu=np.log(period_true), sd=0.1)
    period = pm.Deterministic("period", pm.math.exp(logP))
    
    # Normal distributions for the map coeffs
    y = pm.Normal("y", mu=y_true, sd=1.0, shape=len(y_true))

    # Normal distributions for r and b
    r = pm.Normal("r", mu=0.06, sd=0.001)
    b = pm.Normal("b", mu=0.4, sd=0.03)
    
    # Set up a Keplerian orbit for the planets
    orbit = xo.orbits.KeplerianOrbit(period=period, t0=t0, b=b)
    
    # Compute the model light curve using starry
    light_curve = op.get_light_curve(orbit=orbit, r=r, t=t, y=y) + mean
    
    # Here we track the value of the model light curve for plotting
    # purposes
    pm.Deterministic("light_curve", light_curve)
    
    # In this line, we simulate the dataset that we will fit
    flux = xo.eval_in_model(light_curve)
    flux += ferr * np.random.randn(len(flux))
    
    # The likelihood function assuming known Gaussian uncertainty
    pm.Normal("obs", mu=light_curve, sd=ferr, observed=flux)
    
    # Fit for the maximum a posteriori parameters given the simuated

names = ['t_peri', 'obs']
        to_remove = []
        for name in names:
            for i, par in enumerate(model.vars):
                if par.name == name:
                    to_remove.append(par)
                    del model.named_vars[name]
        [model.vars.remove(p) for p in to_remove]

        p = self.prior.pars
        with model:
            t_peri = pm.Deterministic('t_peri',
                                      p['P'] * p['M0'] / (2*np.pi))

            # Set up the orbit model
            orbit = xo.orbits.KeplerianOrbit(period=p['P'],
                                             ecc=p['e'],
                                             omega=p['omega'],
                                             t_periastron=t_peri)

        # design matrix
        M = get_trend_design_matrix(data, ids, self.prior.poly_trend)

        # deal with v0_offsets, trend here:
        _, offset_names = validate_n_offsets(self.prior.n_offsets)
        _, vtrend_names = validate_poly_trend(self.prior.poly_trend)

        with model:
            v_pars = ([p['v0']]
                      + [p[name] for name in offset_names]
                      + [p[name] for name in vtrend_names[1:]])  # skip v0
            v_trend_vec = tt.stack(v_pars, axis=0)

# We're not fitting for theta
    # @dfm: How do I prevent pymc3 from fitting for it?
    theta = np.ones_like(t)  # np.ones_like(t) * pm.Uniform("theta", 0, 1)
    
    # The map Ylm degree is zero, so there are no Ylms to fit
    y = np.empty(0)
    # y = tt.as_tensor_variable([], name='y')
    # y.name = 'y'

    # Vectorize the occultor radius
    rs = np.ones_like(t) * r
    rs.name = 'r'

    # Set up a Keplerian orbit for the planet
    orbit = xo.orbits.KeplerianOrbit(period=period, t0=t0, b=b)
    coords = orbit.get_relative_position(t)
    _x = coords[0]
    _y = coords[1]
    _z = -coords[2]

    # Compute the model light curve using starry
    rv_model = starry_op(
        y,
        u,
        inc,
        obl,
        veq, 
        alpha,
        theta,
        _x,
        _y,

"b":        0.25,
    "r":        0.1,
    "u":        [0.4, 0.26],
    "rv_err":   0.0005
}

# Generate a synthetic dataset
t = np.linspace(-0.15, 0.15, 1000)
udeg = len(truths["u"])
map = starry.DopplerMap(udeg=udeg)
map[1:] = truths["u"]
map.inc = truths["inc"]
map.obl = truths["obl"]
map.alpha = truths["alpha"]
map.veq = truths["veq"]
orbit = xo.orbits.KeplerianOrbit(period=truths["period"], 
                                 t0=truths["t0"], b=truths["b"])
coords = orbit.get_relative_position(t)
x = coords[0].eval()
y = coords[1].eval()
z = -coords[2].eval()
truths["rv"] = map.rv(xo=x, yo=y, zo=z, ro=truths["r"])

# Noise it
rv = truths["rv"] + truths["rv_err"] * np.random.randn(len(t))

# Plot it
plt.plot(t, rv, 'k.', alpha=0.3, ms=3)
plt.plot(t, truths["rv"])
plt.show()

# Sample it