Make it easier to include models in `.plot()` function.

[zkbt]

Right now, the .plot() function (which by default quietly calls .plot_lightcurves()) can make plots like this:

from chromatic import *

r = (SimulatedRainbow(wlim=[1,5]*u.micron)
     .inject_transit()
     .inject_noise())

r.bin(dw=1*u.micron).plot()

We need to add some abilities to make it just as easy to include models along with these too, to plot residuals, and to plot their histograms. There are a couple things already in the works: @mone0982 is drafting a "plot a histogram of the residuals for this wavelength" function in #115 and @zkbt recently added a RainbowWithModel object that does automatic calculation of a residual array based on an attached model. The remaining things that we should might include:

• [x] confirm that .plot can plot model quantities just as easily
• [x] make sure that it's possible over/underplot a model along with the flux
• [ ] write a little wrapper that would expand out an up-to-three-panel plot with flux, residuals, and histograms, including points for the data and lines for the model

@erinmmay , to help with this, a really helpful first step would be to describe exactly what you'd like for this plot and/or point to some examples you like? I remember the MOPPS plots being really nice! I have a little bit of work to do to tidy up the model management going on in the background, but then I think we can start iterating on versions of this?

[zkbt]

@erinmmay, here's what the some options for the plot_lightcurves_and_residuals_with_models look like. (These are all the same function, just with different keyword arguments supplied.)

def test_plot_lightcurves_and_residuals_with_models():
    s = SimulatedRainbow(R=3, dt=3 * u.minute)
    r = s.inject_transit(
        limb_dark="quadratic",
        u=np.transpose(
            [np.linspace(1.0, 0.0, s.nwave), np.linspace(0.5, 0.0, s.nwave)]
        ),
    ).inject_noise(signal_to_noise=1000)
    for i, options in enumerate(
        [
            dict(),
            dict(errorbar=True),
            dict(cmap=one2another("skyblue", "sienna")),
            dict(data_plotkw=dict(alpha=0.5), cmap="magma_r"),
            dict(figsize=(8, 4), cmap=one2another("orchid", "indigo")),
        ]
    ):
        r.plot_lightcurves_and_residuals_with_models(**options)
        plt.savefig(
            os.path.join(
                test_directory, f"rainbow-of-lightcurves-and-residuals-example{i}.png"
            )
        )

Thoughts on these? They should be easy to produce as long as we have wavelength, time, flux, uncertainty, model defined.

[erinmmay]

These look great! I'm partial to including the errorbars :) I think what we'll want for the paper version is the ability to grab every X light curves from that wavelength/time/flux/uncertainty/model rainbow so that we don't have to plot them all.