Adding experimental custom stellar spectrum feature

[bmorris3]

For @mhooton!

Here's a demo of the new syntax for specifying a non-blackbody stellar spectrum:

import numpy as np
from scipy.stats import binned_statistic
import astropy.units as u

from kelp import Model, Planet, Filter, StellarSpectrum

p = Planet.from_name('MASCARA-1')
filt_cheops = Filter.from_name('CHEOPS')
filt_irac2 = Filter.from_name('IRAC 2')

# Load the spectrum from Daniel Kitzmann
mascara1_spectrum = np.loadtxt('MASCARA-1_spectrum.dat', skiprows=1)

# Bin down the stellar spectrum to a sensible resolution
bs = binned_statistic(
    mascara1_spectrum[:, 0], mascara1_spectrum[:, 1], 
    bins=1000, statistic=np.median
)
bin_wavelength = 0.5 * (bs.bin_edges[1:] + bs.bin_edges[:-1]) * u.um

# construct a StellarSpectrum object to use when computing the planet-to-star flux ratio (phase curve)
spectrum = StellarSpectrum(
    bin_wavelength, bs.statistic * u.W / u.m**2 / u.um / np.pi
)

# Here's the old syntax which assumes a blackbody for the stellar SED:
m_cheops = Model(
    0, 0.575, 4.5, 0, [[0], [0, 0.3, 0]], 1, 
    planet=p, filt=filt_cheops
)
m_irac2 = Model(
    0, 0.575, 4.5, 0, [[0], [0, 0.3, 0]], 1, 
    planet=p, filt=filt_irac2
)

# Here's the new syntax which interpolates the stellar spectrum:
m_cheops_stellarspec = Model(
    0, 0.575, 4.5, 0, [[0], [0, 0.3, 0]], 1, 
    planet=p, filt=filt_cheops, stellar_spectrum=spectrum
)
m_irac2_stellarspec = Model(
    0, 0.575, 4.5, 0, [[0], [0, 0.3, 0]], 1, 
    planet=p, filt=filt_irac2, stellar_spectrum=spectrum
)

# Compute phase curves with these models
xi = np.linspace(-np.pi, np.pi, 100)
pc_cheops = m_cheops.thermal_phase_curve(xi)
pc_irac2 = m_irac2.thermal_phase_curve(xi)

pc_cheops_ss = m_cheops_stellarspec.thermal_phase_curve(xi)
pc_irac2_ss = m_irac2_stellarspec.thermal_phase_curve(xi)

pc_cheops.plot(color='C0', lw=3, label='CHEOPS BB')
pc_irac2.plot(color='C0', lw=3, ls='--', label='IRAC-2 BB')

# Plot the results:
pc_cheops_ss.plot(color='C1', lw=3, label='CHEOPS spectrum')
pc_irac2_ss.plot(color='C1', lw=3, ls='--', label='IRAC-2 spectrum')
plt.legend(loc='upper left')
plt.ylabel('Fp/Fs [ppm]')
plt.xlabel('$\\xi$')

which produces the following figure: