How to use the starry.Planet function in starry

nF = np.zeros_like(time)
    for i in range(npts):
        nF[i] = NumericalFlux(b[i], rplanet, u1, u2)
    nF /= np.nanmax(nF)
    den = (1 - nF)
    den[den == 0] = 1e-10

    # Compute the starry flux
    # Instantiate a second-order map and a third-order map with u(3) = 0
    # The second-order map is optimized for speed and uses different
    # equations, but they should yield identical results.
    for lmax in [2, 3]:
        star = Star(lmax)
        star.map[1] = u1
        star.map[2] = u2
        planet = Planet(r=rplanet, a=a, inc=inc, porb=P, lambda0=90)
        system = System([star, planet])
        system.compute(time)
        sF = np.array(star.flux)
        sF /= sF[0]

        # Compute the error, check that it's better than 1 ppb
        error = np.max((np.abs(nF - sF) / den)[np.where(nF < 1)])
        assert error < 1e-9

spider_params.p_u2 = 0.0            # Planetary limb darkening parameter

    # SPIDERMAN spherical harmonics parameters
    ratio = 1.0e-3                                  # Planet-star flux ratio
    spider_params.sph = [ratio, ratio / 2, 0, 0]      # vector of Ylm coeffs
    spider_params.degree = 2
    spider_params.la0 = 0.0
    spider_params.lo0 = 0.0

    # Define starry model parameters to match SPIDERMAN system

    # Define star
    star = starry.Star()

    # Define planet
    planet = starry.Planet(lmax=2,
                           lambda0=90.0,
                           w=spider_params.w,
                           r=spider_params.rp,
                           # Factor of pi to match SPIDERMAN normalization
                           L=1.0e-3 * np.pi,
                           inc=spider_params.inc,
                           a=spider_params.a,
                           porb=spider_params.per,
                           tref=spider_params.t0,
                           prot=spider_params.per,
                           ecc=spider_params.ecc)

    # Define spherical harmonic coefficients
    planet.map[0, 0] = 1.0
    planet.map[1, -1] = 0.0
    planet.map[1, 0] = 0.0

"""Test the luminosity units."""
import matplotlib.pyplot as pl
import starry
import numpy as np

star = starry.Star()
planet = starry.Planet(L=1e-2, r=0.1, a=3)
sys = starry.System([star, planet])
time = np.linspace(0.25, 1.25, 10000)
sys.compute(time)
pl.plot(time, sys.flux, label="Uniform")

star[1] = 0.4
star[2] = 0.26
sys.compute(time)
pl.plot(time, sys.flux, label="Limb-Darkened")

pl.legend()
pl.show()

"""Testing light travel time delay."""
import numpy as np
import matplotlib.pyplot as pl
import starry

star = starry.Star()
planet = starry.Planet(L=1e-1, prot=1, porb=1, r=1, Omega=1)
sys = starry.System([star, planet])
planet.a /= 50
time = np.linspace(0.4999, 0.5001, 100000)

for R in [0, 1, 5, 10]:
    sys.scale = R
    sys.compute(time)
    a = planet.a * R * 0.00464913034
    pl.plot(time, planet.flux, label="a=%.3f AU" % a)

    print("Secondary eclipse time: ", time[np.argmin(planet.flux)])
    print("Analytic solution: ", 0.5 + 2 * planet.a * R * 6.957e8 /
          299792458 / 86400)

pl.legend()
pl.show()

npts = 500
time = np.linspace(-0.25, 0.25, npts)

# Compute the semi-major axis from Kepler's third law in units of rstar
a = ((P * 86400) ** 2 * (1.32712440018e20 * mstar) /
     (4 * np.pi ** 2)) ** (1. / 3.) / (6.957e8 * rstar)

# Get the inclination in degrees
inc = np.arccos(b0 / a) * 180 / np.pi

# Compute and plot the starry flux
star = Star()
star.map[1] = u1
star.map[2] = u2

planet = Planet(r=rplanet, inc=inc, porb=P, a=a, lambda0=90)
system = System([star, planet])
system.compute(time)
sF = np.array(star.flux)
sF /= sF[0]
ax[0].plot(time, sF, '-', color='C0', label='starry')

# Compute and plot the flux from numerical integration
print("Computing numerical flux...")
f = 2 * np.pi / P * time
b = a * np.sqrt(1 - np.sin(np.pi / 2. + f) ** 2
                * np.sin(inc * np.pi / 180) ** 2)
nF = np.zeros_like(time)
for i in tqdm(range(npts)):
    nF[i] = NumericalFlux(b[i], rplanet, u1, u2)
nF /= np.nanmax(nF)
ax[0].plot(time[::20], nF[::20], 'o', color='C4', label='numerical')

sF /= sF[0]
ax[0].plot(time, sF, '-', color='C0')


# System #2: A three-planet system with planet-planet occultations
# ----------------------------------------------------------------

# Instantiate the star
star = Star()

# Give the star a quadratic limb darkening profile
star.map[1] = 0.4
star.map[2] = 0.26

# Instantiate planet b
b = Planet(r=0.04584,
           L=1e-3,
           inc=89.0,
           porb=2.1,
           prot=2.1,
           a=6.901084,
           lambda0=90,
           tref=0.5)

# Give the planet a wonky map
b.map[1, 0] = 0.5
b.map[2, 1] = 0.5

# Instantiate planet c
c = Planet(r=0.07334,
           L=1e-3,
           inc=90.0,

ax[1].set_xlim(0, 20)
time = np.linspace(0, 20, 10000)


# System #1: A one-planet system with a hotspot offset
# ----------------------------------------------------

# Instantiate the star
star = Star()

# Give the star a quadratic limb darkening profile
star.map[1] = 0.4
star.map[2] = 0.26

# Instantiate planet b
b = Planet(r=0.091679,
           L=5e-3,
           inc=90,
           porb=4.3,
           prot=4.3,
           a=11.127991,
           lambda0=90,
           tref=2,
           axis=[0, 1, 0])

# Give the planet a simple dipole map
b.map[1, 0] = 0.5

# Rotate the planet map to produce a hotspot offset of 15 degrees
b.map.rotate(theta=15)

# Compute and plot the starry flux