How to use the starry.Secondary function in starry

def test_bodies():
    pri = starry.Primary(starry.Map())
    sec = starry.Secondary(starry.Map(ydeg=1), porb=1.0)
    sys = starry.System(pri, sec)
    assert sys.primary == pri
    assert sys.secondaries[0] == sec

def test_compare_to_map_rv():
    """Ensure we get the same result by calling `map.rv()` and `sys.rv()`.
    """
    # Define the map
    map = starry.Map(ydeg=1, udeg=2, rv=True, amp=1, veq=1, alpha=0)
    map[1, 0] = 0.5

    # Define the star
    A = starry.Primary(
        map, r=1.0, m=1.0, prot=0, length_unit=u.Rsun, mass_unit=u.Msun
    )

    # Define the planet
    b = starry.Secondary(
        starry.Map(rv=True, amp=1, veq=0),
        r=0.1,
        porb=1.0,
        m=0.01,
        t0=0.0,
        inc=86.0,
        ecc=0.3,
        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)

def test_default_system_units():
    pri = starry.Primary(starry.Map())
    sec = starry.Secondary(starry.Map(), porb=1.0)
    sys = starry.System(pri, sec)
    assert sys.time_unit == u.day

# Compute the time of eclipse
    E = np.arctan2(
        np.sqrt(1 - ecc ** 2) * np.sin(f + np.pi), ecc + np.cos(f + np.pi)
    )
    M = E - ecc * np.sin(E)
    t_ecl = (M - M0) * porb / (2 * np.pi)

    # This is the required phase offset such that the map coefficients
    # correspond to what the observer sees at secondary eclipse
    theta0 = -(t_ecl - t0) * 360

    # Version 1
    pri = starry.Primary(starry.Map(udeg=2))
    pri.map[1:] = u
    sec = starry.Secondary(
        starry.Map(ydeg=ydeg, amp=amp),
        porb=porb,
        r=r,
        m=m,
        inc=90,
        Omega=0,
        ecc=ecc,
        w=w,
        t0=t0,
        theta0=theta0,
    )
    sec.map[1:, :] = y
    sys = starry.System(pri, sec)
    flux = sys.flux(t)

    # Compare

def test_system_rv_show(mp4=False):
    pri = starry.Primary(starry.Map(rv=True))
    sec = starry.Secondary(starry.Map(rv=True), porb=1.0)
    sys = starry.System(pri, sec)
    sys.show(0.1, file="tmp.pdf")
    os.remove("tmp.pdf")
    if mp4:
        sys.show([0.1, 0.2], file="tmp.mp4")
        os.remove("tmp.mp4")

woodbury = [False, True]
solve_inputs = itertools.product(vals, vals)
lnlike_inputs = itertools.product(vals, vals, woodbury)

# Instantiate a star with a dipole map
A = starry.Primary(starry.Map(ydeg=1), prot=0.0)
amp_true = 0.75
y_true = np.array([1, 0.1, 0.2, 0.3])
inc_true = 60
A.map.amp = amp_true
A.map[1, :] = y_true[1:]
A.map.inc = inc_true

# Instantiate two transiting planets with different longitudes of
# ascending node. This ensures there's no null space!
b = starry.Secondary(starry.Map(amp=0), porb=1.0, r=0.1, t0=-0.05, Omega=30.0)
c = starry.Secondary(starry.Map(amp=0), porb=1.0, r=0.1, t0=0.05, Omega=-30.0)
sys = starry.System(A, b, c)

# Generate a synthetic light curve with just a little noise
t = np.linspace(-0.1, 0.1, 100)
flux = sys.flux(t)
sigma = 1e-5
np.random.seed(1)
flux += np.random.randn(len(t)) * sigma


@pytest.mark.parametrize("L,C", solve_inputs)
def test_solve(L, C):
    # Place a generous prior on the map coefficients
    if L == "scalar":
        A.map.set_prior(L=1)

def test_normalization(d=10, r=1):
    """Test the normalization of the flux."""
    # Instantiate a system. Planet has radius `r` and is at
    # distance `d` from a point illumination source.
    planet = starry.Secondary(starry.Map(reflected=True), a=d, r=r)
    star = starry.Primary(starry.Map(), r=0)
    sys = starry.System(star, planet)

    # Get the star & planet flux when it's at full phase
    t_full = 0.5 * sys._get_periods()[0]
    f_star, f_planet = sys.flux(t=t_full, total=False)

    # Star should have unit flux
    assert np.allclose(f_star, 1.0)

    # Planet should have flux equal to (2 / 3) r^2 / d^2
    assert np.allclose(f_planet, (2.0 / 3.0) * r ** 2 / d ** 2)

def test_bad_sys_data():
    pri = starry.Primary(starry.Map(ydeg=1))
    sec = starry.Secondary(starry.Map(ydeg=1), porb=1.0)
    sys = starry.System(pri, sec)

    # User didn't provide the covariance
    with pytest.raises(ValueError) as e:
        sys.set_data([0.0])

    # User didn't provide a dataset
    with pytest.raises(ValueError) as e:
        sys.solve(t=[0.0])

    # Provide a dummy dataset
    sys.set_data([0.0], C=1.0)

    # User didn't provide a prior for the primary
    with pytest.raises(ValueError) as e:
        sys.solve(t=[0.0])

def test_system_show(mp4=False):
    pri = starry.Primary(starry.Map())
    sec = starry.Secondary(starry.Map(), porb=1.0)
    sys = starry.System(pri, sec)
    sys.show(0.1, file="tmp.pdf")
    os.remove("tmp.pdf")
    if mp4:
        sys.show([0.1, 0.2], file="tmp.mp4")
        os.remove("tmp.mp4")
    sys.show([0.1, 0.2], file="tmp.gif")
    os.remove("tmp.gif")

def test_period_semi():
    # Check that an error is raised if neither a nor porb is given
    with pytest.raises(ValueError) as e:
        body = starry.Secondary(starry.Map())
    assert "Must provide a value for either `porb` or `a`" in str(e.value)

    # Check that the semi --> period conversion works
    pri = starry.Primary(starry.Map(), m=1.0, mass_unit=u.Msun)
    sec = starry.Secondary(
        starry.Map(), a=10.0, m=1.0, length_unit=u.AU, mass_unit=u.Mearth
    )
    sys = starry.System(pri, sec)
    period = sys._get_periods()[0]
    true_period = (
        (
            (2 * np.pi)
            * (sec.a * sec.length_unit) ** (3 / 2)
            / (np.sqrt(G * (pri.m * pri.mass_unit + sec.m * sec.mass_unit)))
        )
        .to(u.day)