How to use the exoplanet.eval_in_model function in exoplanet

# 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
    # dataset
    map_soln = pm.find_MAP(start=model.test_point)


np.random.seed(42)
sampler = xo.PyMC3Sampler(window=100, finish=200)
with model:
    burnin = sampler.tune(tune=2500, start=map_soln, step_kwargs=dict(target_accept=0.9))
    trace = sampler.sample(draws=3000)

# 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)
    y = y_true + yerr * np.random.randn(len(yerr))

    # Compute the prediction
    vrad_pred = orbit.get_radial_velocity(t, K=tt.exp(logK))
    if N_pl == 1:
        vrad_pred = vrad_pred[:, None]
    A_pred = np.vander(t - 0.5 * (x.min() + x.max()), 3)
    bkg_pred = tt.dot(A_pred, trend)
    rv_model_pred = tt.sum(vrad_pred, axis=-1) + bkg_pred

    # Likelihood
    err = yerr
    # err = tt.sqrt(yerr**2 + tt.exp(2*logs))
    pm.Normal("obs", mu=rv_model, sd=err, observed=y)

    # Optimize