How to use the starry.Map function in starry

def test_bad_case():
    """
    Test pathological wrong case identification.

    """
    map = starry.Map(reflected=True)

    # These values lead to a (very) wrong flux
    theta0 = -0.0409517311212404
    b0 = -0.83208413089546
    bo0 = 12.073565287605442
    ro = 12.155639360414618

    # Perturb theta in the vicinity of theta0
    delta = np.linspace(0, 1e-6, 100)
    theta = np.concatenate((theta0 - delta[::-1], theta0 + delta))

    # Compute the flux
    b = b0 * np.ones_like(theta)
    bo = bo0 * np.ones_like(theta)
    sT, *_ = map.ops._sT.func(b, theta, bo, ro, 0.0)
    flux = sT[:, 0]

"""Test the stability of the limb darkening calculations near singular points."""
import starry
import numpy as np
import pytest

# We'll run all tests on tenth degree maps
map = starry.Map(10)
map[0, 0] = 1
map[:] = 1
map_multi = starry.Map(10, multi=True)
map_multi[0, 0] = 1
map_multi[:] = 1


def test_b_near_zero(ftol=1e-11, gtol=1e-11):
    b = 10.0 ** np.linspace(-18, -6, 100)
    b[0] = 0
    f = map.flux(xo=b, yo=0.0, ro=0.1)
    f_multi = map_multi.flux(xo=b, yo=0.0, ro=0.1)
    assert(np.max(np.abs(f - f_multi)) < ftol)
    _, g = map.flux(xo=b, yo=0.0, ro=0.1, gradient=True)
    _, g_multi = map_multi.flux(xo=b, yo=0.0, ro=0.1, gradient=True)
    for key in g.keys():

def test_flux_quad_ld(abs_tol=1e-5, rel_tol=1e-5, eps=1e-7):
    with change_flags(compute_test_value="off"):
        map = starry.Map(udeg=2)
        xo = np.linspace(-1.5, 1.5, 10)
        yo = np.ones_like(xo) * 0.3
        zo = 1.0 * np.ones_like(xo)
        ro = 0.1
        np.random.seed(14)
        u = np.array([-1.0] + list(np.random.randn(2)))
        func = lambda *args: map.ops.flux(*args)

        verify_grad(
            func,
            (xo, yo, zo, ro, u),
            abs_tol=abs_tol,
            rel_tol=rel_tol,
            eps=eps,
            n_tests=1,
        )

def test_check_kwargs_map(caplog):
    """Test that we capture bad keyword arguments."""
    starry.config.quiet = False
    caplog.clear()
    map = starry.Map(giraffe=10)
    assert len(caplog.records) >= 1
    assert any(
        [
            "Invalid keyword `giraffe`" in str(rec.message)
            for rec in caplog.records
        ]

assert np.allclose(map.p, np.array([0, 0, 0,
                                        np.sqrt(3 / (4 * np.pi))]))

    map.axis = [0, 0, 1]
    map.rotate(-90)
    assert np.allclose(map.y, np.array([0, 0, 0, norm]))
    assert np.allclose(map.p, np.array([0,
                                        np.sqrt(3 / (4 * np.pi)), 0, 0]))
    map.axis = [0, 1, 0]
    map.rotate(-90)
    assert np.allclose(map.y, np.array([0, 0, norm, 0]))
    assert np.allclose(map.p, np.array([0, 0,
                                        np.sqrt(3 / (4 * np.pi)), 0]))

    # A more complex rotation
    map = starry.Map(5, multi=multi)
    map[:, :] = 1
    map.axis = [1, 2, 3]
    map.rotate(60)
    benchmark = np.array([1.,  1.39148148,  0.91140212,  0.48283069,
                          1.46560344, 0.68477955,  0.30300625,  1.33817773,
                          -0.70749636, 0.66533701, 1.5250326,  0.09725931,
                          -0.13909678,  1.06812054, 0.81540106, -1.54823596,
                          -0.50475248,  1.90009363, 0.68002942, -0.10159448,
                          -0.48406777,  0.59834505, 1.22007458, -0.27302899,
                          -1.58323797, -1.37266583, 1.44638769,  1.36239055,
                          0.22257365, -0.24387785, -0.62003044,  0.03888137,
                          1.05768142,  0.87317586, -1.46092763, -0.81070502])
    assert np.allclose(map.y, benchmark)

    # Test rotation caching
    map = starry.Map(2, multi=multi)

def test_approximation(plot=False):
    """
    Test our polynomial approximation to the Oren-Nayar intensity.

    """
    # Simple map
    map = starry.Map(reflected=True)

    # Approximate and exact intensities
    map.roughness = 90
    img_approx = map.render(xs=1, ys=2, zs=3)
    img_exact = map.render(xs=1, ys=2, zs=3, on94_exact=True)
    img_diff = img_exact - img_approx
    diff = img_diff.reshape(-1)

    mu = np.nanmean(diff)
    std = np.nanstd(diff)
    maxabs = np.nanmax(np.abs(diff))

    if plot:
        fig, ax = plt.subplots(1, 3, figsize=(14, 2.5))
        im = ax[0].imshow(
            img_exact,

def test_high_order_ld():
    """Test transit light curve generation for 8th order limb darkening."""
    # Input params
    u = [0.4, 0.26, 0.3, 0.5, -0.2, 0.5, -0.7, 0.3]
    npts = 25
    r = 0.1
    b = np.linspace(0, 1 + r + 0.1, npts)

    # Compute the starry flux
    map = Map(len(u))
    map[0, 0] = 1
    map[:] = u
    sF = map.flux(xo=b, yo=0, ro=r)

    # Numerical flux
    nF = np.zeros_like(b)
    for i in range(npts):
        nF[i] = NumericalFlux(b[i], r, u)

    # Compute the error, check that it's better than 1 ppb
    error = np.max((np.abs(nF - sF)))
    assert error < 1e-9, "Error in flux exceeds 1 ppb (= %.3e)." % error

def test_spectral_no_grad():
    map = starry.Map(ydeg=ydeg, udeg=0, nw=nw)
    flux = map.flux(theta=theta, xo=xo, yo=yo, ro=ro, gradient=False)
    assert(flux.shape == (npts, nw))
    intensity = map(theta=theta, x=xo, y=yo)
    assert(intensity.shape == (npts, nw))

"""
import starry
import numpy as np
from scipy.linalg import cho_solve
from scipy.stats import multivariate_normal
import pytest
import itertools

# Parameter combinations we'll test
vals = ["scalar", "vector", "matrix", "cholesky"]
woodbury = [False, True]
solve_inputs = itertools.product(vals, vals)
lnlike_inputs = itertools.product(vals, vals, woodbury)

# Instantiate a dipole map
map = starry.Map(ydeg=1, reflected=True)
amp_true = 0.75
inc_true = 60
y_true = np.array([1, 0.1, 0.2, 0.3])
map.amp = amp_true
map[1, :] = y_true[1:]
map.inc = inc_true

# Generate a synthetic light curve with just a little noise
theta = np.linspace(0, 360, 100)
phi = 3.5 * theta
xs = np.cos(phi * np.pi / 180)
ys = 0.1 * np.cos(phi * np.pi / 180)
zs = np.sin(phi * np.pi / 180)
kwargs = dict(theta=theta, xs=xs, ys=ys, zs=zs)
flux = map.flux(**kwargs).eval()
sigma = 1e-5

# Standard map
    ylm = Map(l)
    ylm._r_M = np.inf
    ylm._r_quartic = np.inf
    ylm._b_J = 0

    # Taylor expand `M` (odd nu)
    # and approximate occultor limb as parabola (even nu)
    ylmtaylor = Map(l)
    ylmtaylor._r_M = 0
    ylmtaylor._r_quartic = 0
    ylmtaylor._b_J = 0

    # Multiprecision (~exact)
    ylm128 = Map(l)
    ylm128.use_mp = True

    # Set up
    fig, ax = pl.subplots(2, 2, figsize=(12, 5))
    fig.suptitle("Degree $l = %d$" % l, fontsize=14)

    #  Loop over the orders
    noise = np.zeros((2 * l + 1, npts)) * np.nan
    error = np.zeros((2 * l + 1, npts)) * np.nan
    noisetaylor = np.zeros((2 * l + 1, npts)) * np.nan
    errortaylor = np.zeros((2 * l + 1, npts)) * np.nan
    for j, m in tqdm(enumerate(range(-l, l + 1)), total=2 * l + 1):
        # Set this coefficient
        parity = (l + m) % 2
        ylm.reset()
        ylm[l, m] = 1