import numpy as np
import matplotlib.pyplot as plt
from astropy.modeling import models
from astropy import units as u
from specutils.spectra import Spectrum1D, SpectralRegion
from specutils.fitting import fit_lines
#
np.random.seed(42)
g1 = models.Gaussian1D(1, 4.6, 0.2)
g2 = models.Gaussian1D(2.5, 5.5, 0.1)
x = np.linspace(0, 10, 200)
y = g1(x) + g2(x) + np.random.normal(0., 0.2, x.shape)
#
spectrum = Spectrum1D(flux=y*u.Jy, spectral_axis=x*u.um)
#
gl_init = models.Gaussian1D(amplitude=1.*u.Jy, mean=4.8*u.um, stddev=0.2*u.um)
gl_fit = fit_lines(spectrum, gl_init, exclude_regions=[SpectralRegion(5.2*u.um, 5.8*u.um)])
yl_fit = gl_fit(x*u.um)
#
f, ax = plt.subplots()  # doctest: +IGNORE_OUTPUT
ax.plot(x, y)  # doctest: +IGNORE_OUTPUT
ax.plot(x, yl_fit)  # doctest: +IGNORE_OUTPUT
ax.set_title("Double Peak - Single Models and Exclude Region")  # doctest: +IGNORE_OUTPUT
ax.grid(True)  # doctest: +IGNORE_OUTPUT