### Install Specutils from Source

Source: https://github.com/astropy/specutils/blob/main/docs/installation.rst

Commands to obtain the source code and perform a local installation.

```console
$ git clone git://github.com/astropy/specutils
```

```console
$ curl  -OL https://github.com/astropy/specutils/tarball/main
```

```console
$ python setup.py install
```

--------------------------------

### Install Specutils for Local Development

Source: https://github.com/astropy/specutils/blob/main/docs/contributing.rst

Install Specutils in an editable mode within a virtual environment, including test dependencies. This setup is recommended for local development.

```shell
conda create --name specutils-dev python=3.11 pip
conda activate specutils-dev
cd specutils/
pip install -e .'[test]'
```

--------------------------------

### Initialize Spectral Fitting Environment

Source: https://github.com/astropy/specutils/blob/main/docs/fitting.rst

Setup imports required for spectral fitting examples.

```python
import numpy as np
import matplotlib.pyplot as plt

from astropy.modeling import models
from astropy import units as u

from specutils.spectra import Spectrum
```

--------------------------------

### Spectrum1D Instantiation Patterns

Source: https://github.com/astropy/specutils/wiki/28Feb2012-Discussion

Examples of proposed initialization methods for the Spectrum1D object.

```python
s = Spectrum1D(flux=<flux array>, dispersion=<dispersion array>)
```

```python
s = Spectrum1D.from_whatever()
```

--------------------------------

### Setup Figure for Doctests

Source: https://github.com/astropy/specutils/blob/main/docs/index.rst

Ensures proper figure initialization for quantity support during testing.

```python
>>> fig = plt.figure()  # necessary because otherwise the doctests fail due to quantity_support and the flux units being different from the last figure
```

--------------------------------

### Install Specutils via Package Managers

Source: https://github.com/astropy/specutils/blob/main/docs/installation.rst

Use these commands to install the stable release of Specutils using conda or pip.

```console
$ conda install -c conda-forge specutils
```

```console
$ pip install specutils
```

--------------------------------

### Spectrum1D Instantiation Example

Source: https://github.com/astropy/specutils/wiki/Spectrum-1D-concept-&-issues

This shows a potential way to instantiate a Spectrum1D object using flux and dispersion arrays. It is suggested that I/O operations should not be part of the Spectrum1D class itself.

```python
s = Spectrum1D(flux=<flux array>, dispersion=<dispersion array>)
```

--------------------------------

### Create a Noisy Gaussian Spectrum

Source: https://github.com/astropy/specutils/blob/main/docs/analysis.rst

Sets up a synthetic spectrum with Gaussian profile and noise for use in analysis examples.

```python
>>> import numpy as np
>>> from astropy import units as u
>>> from astropy.nddata import StdDevUncertainty
>>> from astropy.modeling import models
>>> from specutils import Spectrum, SpectralRegion
>>> np.random.seed(42)
>>> spectral_axis = np.linspace(11., 1., 200) * u.GHz
>>> spectral_model = models.Gaussian1D(amplitude=5*(2*np.pi*0.8**2)**-0.5*u.Jy, mean=5*u.GHz, stddev=0.8*u.GHz)
>>> flux = spectral_model(spectral_axis)
>>> flux += np.random.normal(0., 0.05, spectral_axis.shape) * u.Jy
>>> uncertainty = StdDevUncertainty(0.2*np.ones(flux.shape)*u.Jy)
>>> noisy_gaussian = Spectrum(spectral_axis=spectral_axis, flux=flux, uncertainty=uncertainty)
>>> import matplotlib.pyplot as plt #doctest:+SKIP
>>> plt.step(noisy_gaussian.spectral_axis, noisy_gaussian.flux) #doctest:+SKIP
```

--------------------------------

### Read Spectrum using Astroquery

Source: https://github.com/astropy/specutils/blob/main/docs/spectrum.rst

Load a Spectrum object using astroquery to download data, then pass the downloaded object to Spectrum.read(). This requires the astroquery package to be installed.

```python
from astroquery.sdss import SDSS  # doctest: +REMOTE_DATA
specs = SDSS.get_spectra(plate=751, mjd=52251, fiberID=160, data_release=14)  # doctest: +REMOTE_DATA
Spectrum.read(specs[0], format="SDSS-III/IV spec")  # doctest: +REMOTE_DATA
```

--------------------------------

### Create SpectrumCollection from a list of Spectra

Source: https://github.com/astropy/specutils/blob/main/docs/spectrum_collection.rst

Use this method when you have a list of individual Spectrum objects. This is useful for combining spectra that may have been processed or loaded separately. Warnings related to spectral axis differences are suppressed in this example.

```python
import warnings
import numpy as np
from astropy import units as u
from specutils import Spectrum, SpectrumCollection
spec = Spectrum(spectral_axis=np.linspace(0, 50, 50) * u.AA,
flux=np.random.randn(50) * u.Jy,
uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
with warnings.catch_warnings():  # Ignore warnings
    warnings.simplefilter('ignore')
    spec1 = Spectrum(spectral_axis=np.linspace(20, 60, 50) * u.AA,
                       flux=np.random.randn(50) * u.Jy,
                       uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
    spec_coll = SpectrumCollection.from_spectra([spec, spec1])

spec_coll.shape
spec_coll.flux.unit
spec_coll.spectral_axis.shape
spec_coll.spectral_axis.unit
```

--------------------------------

### Initialize Specutils Environment

Source: https://github.com/astropy/specutils/blob/main/docs/index.rst

Imports necessary libraries and enables unit support for matplotlib axes.

```python
from astropy.io import fits
from astropy import units as u
import numpy as np
from matplotlib import pyplot as plt
from astropy.visualization import quantity_support
quantity_support()  # for getting units on the axes below
```

--------------------------------

### Query loaders by extension

Source: https://github.com/astropy/specutils/blob/main/docs/custom_loading.rst

Retrieve a list of available loaders for a specific file extension.

```python
from specutils.io import get_loaders_by_extension

loaders = get_loaders_by_extension('fits')
```

--------------------------------

### Create SpectralRegion from Line List

Source: https://github.com/astropy/specutils/blob/main/docs/spectral_regions.rst

Demonstrates creating a SpectralRegion from a QTable returned by line finding functions.

```python
>>> from astropy import units as u
>>> import numpy as np
>>> from specutils import Spectrum, SpectralRegion
>>> from astropy.modeling.models import Gaussian1D
>>> from specutils.fitting import find_lines_derivative

>>> g1 = Gaussian1D(1, 4.6, 0.2)
>>> g2 = Gaussian1D(2.5, 5.5, 0.1)
>>> g3 = Gaussian1D(-1.7, 8.2, 0.1)

>>> x = np.linspace(0, 10, 200)
>>> y = g1(x) + g2(x) + g3(x)
>>> spectrum = Spectrum(flux=y * u.Jy, spectral_axis=x * u.um)

>>> lines = find_lines_derivative(spectrum, flux_threshold=0.01)
>>> spec_reg = SpectralRegion.from_line_list(lines)
>>> spec_reg
Spectral Region, 3 sub-regions:
  (4.072864321608041 um, 5.072864321608041 um)
  (4.977386934673367 um, 5.977386934673367 um)
  (7.690954773869347 um, 8.690954773869347 um)
```

--------------------------------

### Apply Convenience Smoothing Functions

Source: https://github.com/astropy/specutils/blob/main/docs/manipulation.rst

Demonstrates the use of box, gaussian, and trapezoid smoothing convenience functions on a Spectrum object.

```python
>>> from specutils import Spectrum
>>> import astropy.units as u
>>> import numpy as np
>>> from specutils.manipulation import box_smooth, gaussian_smooth, trapezoid_smooth

>>> spec1 = Spectrum(spectral_axis=np.arange(1, 50) * u.nm,
...                    flux=np.random.default_rng(12345).random(49)*u.Jy)
>>> spec1_bsmooth = box_smooth(spec1, width=3)
>>> spec1_gsmooth = gaussian_smooth(spec1, stddev=3)
>>> spec1_tsmooth = trapezoid_smooth(spec1, width=3)
>>> gaussian_smooth(spec1, stddev=3) # doctest: +FLOAT_CMP
<Spectrum(flux=[0.25860917267578276 ... 0.15868783272966752] Jy (shape=(49,), mean=0.48009 Jy); spectral_axis=<SpectralAxis [ 1.  2.  3. ... 47. 48. 49.] nm> (length=49))>
```

--------------------------------

### Identify SDSS MaNGA Cube Format

Source: https://github.com/astropy/specutils/blob/main/docs/identify.rst

Use `identify_spectrum_format` to determine the format of a downloaded SDSS MaNGA data cube file. Requires `astropy.utils.data.download_file` for fetching the example file.

```python
>>> from astropy.utils.data import download_file
>>> from specutils.io.registers import identify_spectrum_format
>>> 
>>> url = 'https://dr17.sdss.org/sas/dr17/manga/spectro/redux/v3_1_1/8485/stack/manga-8485-1901-LOGCUBE.fits.gz'
>>> dd = download_file(url)  # doctest: +REMOTE_DATA
>>> identify_spectrum_format(dd)  # doctest: +REMOTE_DATA
'MaNGA cube'
```

--------------------------------

### Identify JWST x1d Spectrum Format

Source: https://github.com/astropy/specutils/blob/main/docs/identify.rst

Identify the format of a JWST extracted 1D spectral file using `identify_spectrum_format`. This example is marked with `doctest: +SKIP` as it assumes a local file path.

```python
>>> from specutils.io.registers import identify_spectrum_format
>>> path = '/data/jwst/jw00626-o030_s00000_nirspec_f170lp-g235m_x1d.fits'
>>> identify_spectrum_format(path)  # doctest: +SKIP
'JWST x1d'
```

--------------------------------

### Clone Specutils Repository

Source: https://github.com/astropy/specutils/blob/main/docs/contributing.rst

Clone your forked repository locally to begin development. Ensure you replace 'your_name_here' with your GitHub username.

```shell
git clone git@github.com:your_name_here/specutils.git
```

--------------------------------

### Load and Prepare Spectrum Data

Source: https://github.com/astropy/specutils/blob/main/docs/manipulation.rst

Imports necessary libraries and prepares a Spectrum object from FITS data.

```python
>>> from astropy.io import fits
>>> from astropy import units as u
>>> import numpy as np
>>> from matplotlib import pyplot as plt
>>> from astropy.visualization import quantity_support
>>> quantity_support()  # for getting units on the axes below  # doctest: +IGNORE_OUTPUT

>>> filename = 'https://data.sdss.org/sas/dr16/sdss/spectro/redux/26/spectra/1323/spec-1323-52797-0012.fits'
>>> # The spectrum is in the second HDU of this file.
>>> with fits.open(filename) as f:  # doctest: +IGNORE_OUTPUT +REMOTE_DATA
...     specdata = f[1].data[1020:1250]  # doctest: +REMOTE_DATA

>>> from specutils import Spectrum
>>> lamb = 10**specdata['loglam'] * u.AA # doctest: +REMOTE_DATA
>>> flux = specdata['flux'] * 10**-17 * u.Unit('erg cm-2 s-1 AA-1') # doctest: +REMOTE_DATA
>>> input_spec = Spectrum(spectral_axis=lamb, flux=flux) # doctest: +REMOTE_DATA

>>> f, ax = plt.subplots()  # doctest: +IGNORE_OUTPUT +REMOTE_DATA
>>> ax.step(input_spec.spectral_axis, input_spec.flux) # doctest: +IGNORE_OUTPUT +REMOTE_DATA
```

--------------------------------

### Perform basic spectrum arithmetic

Source: https://github.com/astropy/specutils/blob/main/docs/arithmetic.rst

Demonstrates addition of two Spectrum objects. Requires compatible WCS information between the two objects.

```python
>>> from specutils import Spectrum
>>> import astropy.units as u
>>> import numpy as np

>>> rng = np.random.default_rng(12345)
>>> spec1 = Spectrum(spectral_axis=np.arange(1, 50) * u.nm, flux=rng.random(49)*u.Jy)
>>> spec2 = Spectrum(spectral_axis=np.arange(1, 50) * u.nm, flux=rng.random(49)*u.Jy)
>>> spec3 = spec1 + spec2
>>> spec3  # doctest: +FLOAT_CMP
<Spectrum(flux=[0.8559405665668484 ... 0.9711264429515736] Jy (shape=(49,), mean=0.91592 Jy); spectral_axis=<SpectralAxis
   (observer to target:
      radial_velocity=0.0 km / s
      redshift=0.0)
  [ 1.  2.  3. ... 47. 48. 49.] nm> (length=49))>
```

--------------------------------

### Use custom loader

Source: https://github.com/astropy/specutils/blob/main/docs/custom_loading.rst

Load data using a custom format defined in the registry.

```python
from specutils import Spectrum

spec = Spectrum.read("path/to/data", format="my-format")
```

--------------------------------

### Fit with Different Units

Source: https://github.com/astropy/specutils/blob/main/docs/fitting.rst

Shows how the fitting process handles initial model guesses provided in units different from the spectrum's spectral axis.

```python
import numpy as np
import matplotlib.pyplot as plt

from astropy.modeling import models
from astropy import units as u

from specutils.spectra import Spectrum
from specutils.fitting import fit_lines

# Create a simple spectrum with a Gaussian.
np.random.seed(0)
x = np.linspace(0., 10., 200)
y = 3 * np.exp(-0.5 * (x- 6.3)**2 / 0.8**2)
y += np.random.normal(0., 0.2, x.shape)

# Create the spectrum
spectrum = Spectrum(flux=y*u.Jy, spectral_axis=x*u.um)

# Fit the spectrum
g_init = models.Gaussian1D(amplitude=3.*u.Jy, mean=61000*u.AA, stddev=10000.*u.AA)
g_fit = fit_lines(spectrum, g_init)
y_fit = g_fit(x*u.um)

plt.plot(x, y)
plt.plot(x, y_fit)
plt.title('Single fit peak, different model units')
plt.grid(True)
```

--------------------------------

### Create Spectrum with FITS WCS

Source: https://github.com/astropy/specutils/blob/main/docs/spectrum.rst

Initialize a Spectrum object with a FITS-style WCS header. Ensure necessary imports for WCS, units, and numpy are included. The spectral axis is automatically derived from the WCS.

```python
>>> from specutils.spectra import Spectrum
>>> import astropy.wcs as fitswcs
>>> import astropy.units as u
>>> import numpy as np
>>> my_wcs = fitswcs.WCS(header={
...     'CDELT1': 1, 'CRVAL1': 6562.8, 'CUNIT1': 'Angstrom', 'CTYPE1': 'WAVE',
...     'RESTFRQ': 1400000000, 'CRPIX1': 25})
>>> spec = Spectrum(flux=[5,6,7] * u.Jy, wcs=my_wcs)
>>> spec.spectral_axis  # doctest: +FLOAT_CMP
<SpectralAxis
   (observer to target:
      radial_velocity=0.0 km / s
      redshift=0.0
    doppler_rest=1400000000.0 Hz
    doppler_convention=None)
  [6.5388e-07, 6.5398e-07, 6.5408e-07] m>
>>> spec.wcs.pixel_to_world(np.arange(3))  # doctest: +FLOAT_CMP
<SpectralCoord [6.5388e-07, 6.5398e-07, 6.5408e-07] m>
```

--------------------------------

### Initialize Spectrum with Uncertainties

Source: https://github.com/astropy/specutils/blob/main/docs/spectrum.rst

Initializes a Spectrum with specific uncertainty classes to enable propagation during arithmetic operations.

```python
>>> from specutils import Spectrum
>>> from astropy.nddata import StdDevUncertainty

>>> spec = Spectrum(spectral_axis=np.arange(5000, 5010) * u.AA, flux=np.random.sample(10) * u.Jy, uncertainty=StdDevUncertainty(np.random.sample(10) * 0.1))
```

--------------------------------

### Create and Plot a Spectrum Object

Source: https://github.com/astropy/specutils/blob/main/docs/index.rst

Initializes a Spectrum object with spectral axis and flux units, then plots the data using matplotlib.

```python
>>> from specutils import Spectrum
>>> lamb = 10**specdata['loglam'] * u.AA # doctest: +REMOTE_DATA
>>> flux = specdata['flux'] * 10**-17 * u.Unit('erg cm-2 s-1 AA-1') # doctest: +REMOTE_DATA
>>> spec = Spectrum(spectral_axis=lamb, flux=flux) # doctest: +REMOTE_DATA

>>> f, ax = plt.subplots()  # doctest: +IGNORE_OUTPUT +REMOTE_DATA
>>> ax.step(spec.spectral_axis, spec.flux) # doctest: +IGNORE_OUTPUT +REMOTE_DATA
```

--------------------------------

### Import Median Smoothing

Source: https://github.com/astropy/specutils/blob/main/docs/manipulation.rst

Imports the median_smooth function for applying median filtering to spectral flux.

```python
>>> from specutils.manipulation import median_smooth
```

--------------------------------

### Run Full Test Suite with Tox

Source: https://github.com/astropy/specutils/blob/main/docs/contributing.rst

Execute the full test suite using tox, which runs tests across multiple Python versions. This is recommended before opening a pull request.

```shell
tox
```

--------------------------------

### Initialize Spectrum with Multi-dimensional Flux

Source: https://github.com/astropy/specutils/blob/main/docs/spectrum.rst

Create a Spectrum object with a multi-dimensional flux array, where one dimension matches the spectral axis length. This allows for efficient operations on collections of spectra sharing the same spectral axis.

```python
>>> from specutils import Spectrum

>>> spec = Spectrum(spectral_axis=np.arange(5000, 5010) * u.AA,
...                 flux=np.random.default_rng(12345).random((5, 10)) * u.Jy)
>>> spec_slice = spec[0]
>>> spec_slice.spectral_axis
<SpectralAxis [5000., 5001., 5002., 5003., 5004., 5005., 5006., 5007., 5008., 5009.] Angstrom>
>>> spec_slice.flux
<Quantity [0.22733602, 0.31675834, 0.79736546, 0.67625467, 0.39110955,
```

--------------------------------

### Fit a Single Gaussian Peak

Source: https://github.com/astropy/specutils/blob/main/docs/fitting.rst

Demonstrates fitting a Gaussian model to a generated spectrum and plotting the results.

```python
from specutils.fitting import fit_lines

# Create a simple spectrum with a Gaussian.
np.random.seed(0)
x = np.linspace(0., 10., 200)
y = 3 * np.exp(-0.5 * (x- 6.3)**2 / 0.8**2)
y += np.random.normal(0., 0.2, x.shape)
spectrum = Spectrum(flux=y*u.Jy, spectral_axis=x*u.um)

# Fit the spectrum and calculate the fitted flux values (``y_fit``)
g_init = models.Gaussian1D(amplitude=3.*u.Jy, mean=6.1*u.um, stddev=1.*u.um)
g_fit = fit_lines(spectrum, g_init)
y_fit = g_fit(x*u.um)

# Plot the original spectrum and the fitted.
plt.plot(x, y, label="Original spectrum")
plt.plot(x, y_fit, label="Fit result")
plt.title('Single fit peak')
plt.grid(True)
plt.legend()
```

--------------------------------

### Apply Dust Extinction to Spectrum

Source: https://github.com/astropy/specutils/blob/main/docs/analysis.rst

Applies dust extinction to a Spectrum instance using the dust_extinction package. Requires defining the wavelength, flux, and extinction model.

```python
from astropy.modeling.blackbody import blackbody_lambda
from dust_extinction.parameter_averages import F99

wave = np.logspace(np.log10(1000), np.log10(3e4), num=10) * u.AA
flux = blackbody_lambda(wave, 10000 * u.K)
spec = Spectrum(spectral_axis=wave, flux=flux)

# define the model
ext = F99(Rv=3.1)

# extinguish (redden) the spectrum
flux_ext = spec.flux * ext.extinguish(spec.spectral_axis, Ebv=0.5)
spec_ext = Spectrum(spectral_axis=wave, flux=flux_ext)
```

--------------------------------

### Initialize Spectrum with Redshift or Radial Velocity

Source: https://github.com/astropy/specutils/blob/main/docs/spectrum.rst

Sets initial redshift or radial velocity values for a Spectrum object upon creation.

```python
>>> spec1 = Spectrum(spectral_axis=np.arange(5000, 5010) * u.AA, flux=np.random.sample(10) * u.Jy, redshift = 0.15)
>>> spec2 = Spectrum(spectral_axis=np.arange(5000, 5010) * u.AA, flux=np.random.sample(10) * u.Jy, radial_velocity = 1000 * u.Unit("km/s"))
```

--------------------------------

### Fit Continuum Over Specific Windows

Source: https://github.com/astropy/specutils/blob/main/docs/fitting.rst

Fits a continuum model to a spectrum restricted to specific wavelength intervals defined by the window parameter.

```python
>>> import warnings
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> import astropy.units as u

>>> from specutils.spectra.spectrum import Spectrum
>>> from specutils.fitting.continuum import fit_continuum

>>> np.random.seed(0)
>>> x = np.linspace(0., 10., 200)
>>> y = 3 * np.exp(-0.5 * (x - 6.3) ** 2 / 0.1 ** 2)
>>> y += np.random.normal(0., 0.2, x.shape)
>>> y += 3.2 * np.exp(-0.5 * (x - 5.6) ** 2 / 4.8 ** 2)

>>> spectrum = Spectrum(flux=y * u.Jy, spectral_axis=x * u.um)
>>> region = [(1 * u.um, 5 * u.um), (7 * u.um, 10 * u.um)]
>>> with warnings.catch_warnings():  # Ignore warnings
...     warnings.simplefilter('ignore')
...     fitted_continuum = fit_continuum(spectrum, window=region)
>>> y_fit = fitted_continuum(x*u.um)

>>> f, ax = plt.subplots()  # doctest: +IGNORE_OUTPUT
```

--------------------------------

### Define Custom Parameter Estimators

Source: https://github.com/astropy/specutils/blob/main/docs/fitting.rst

Demonstrates how to define custom estimators for models lacking predefined ones by assigning lambda functions to the model's parameter estimator attributes.

```python
>>> from specutils import SpectralRegion
>>> from specutils.fitting import estimate_line_parameters
>>> from specutils.manipulation import extract_region
>>> from specutils.analysis import centroid, fwhm

>>> sub_region = SpectralRegion(4*u.um, 5*u.um)
>>> sub_spectrum = extract_region(spectrum, sub_region)

>>> ricker = models.RickerWavelet1D()
>>> ricker.amplitude.estimator = lambda s: max(s.flux)
>>> ricker.x_0.estimator = lambda *args: centroid(args[0], region=None)
>>> ricker.sigma.estimator = lambda *args: fwhm(args[0])

>>> estimate_line_parameters(spectrum, ricker)  # doctest: +FLOAT_CMP +IGNORE_WARNINGS
<RickerWavelet1D(amplitude=2.4220684 Jy, x_0=3.60454769 um, sigma=0.24416769 um)>
```

--------------------------------

### Load SDSS Spectral Data

Source: https://github.com/astropy/specutils/blob/main/docs/index.rst

Downloads and opens a FITS file from the SDSS archive to access spectral data.

```python
filename = 'https://data.sdss.org/sas/dr16/sdss/spectro/redux/26/spectra/1323/spec-1323-52797-0012.fits'
# The spectrum is in the second HDU of this file.
with fits.open(filename) as f:
    specdata = f[1].data
```

--------------------------------

### Fit Multiple Peaks

Source: https://github.com/astropy/specutils/blob/main/docs/fitting.rst

Demonstrates fitting a compound model to a spectrum containing two distinct peaks.

```python
import numpy as np
import matplotlib.pyplot as plt

from astropy.modeling import models
from astropy import units as u

from specutils.spectra import Spectrum
from specutils.fitting import fit_lines

# Create a simple spectrum with a Gaussian.
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)

# Create the spectrum to fit
spectrum = Spectrum(flux=y*u.Jy, spectral_axis=x*u.um)

# Fit the spectrum
g1_init = models.Gaussian1D(amplitude=2.3*u.Jy, mean=5.6*u.um, stddev=0.1*u.um)
g2_init = models.Gaussian1D(amplitude=1.*u.Jy, mean=4.4*u.um, stddev=0.1*u.um)
```

--------------------------------

### Fit and Normalize Spectrum Continuum

Source: https://github.com/astropy/specutils/blob/main/docs/fitting.rst

Shows how to fit a generic continuum to a spectrum and normalize the original spectrum by dividing it by the fitted model.

```python
>>> import warnings
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from astropy.modeling import models
>>> from astropy import units as u
>>> from specutils.spectra import Spectrum, SpectralRegion
>>> from specutils.fitting import fit_generic_continuum

>>> np.random.seed(0)
>>> x = np.linspace(0., 10., 200)
>>> y = 3 * np.exp(-0.5 * (x - 6.3)**2 / 0.1**2)
>>> y += np.random.normal(0., 0.2, x.shape)

>>> y_continuum = 3.2 * np.exp(-0.5 * (x - 5.6)**2 / 4.8**2)
>>> y += y_continuum

>>> spectrum = Spectrum(flux=y*u.Jy, spectral_axis=x*u.um)

>>> with warnings.catch_warnings():  # Ignore warnings
...     warnings.simplefilter('ignore')
...     g1_fit = fit_generic_continuum(spectrum)

>>> y_continuum_fitted = g1_fit(x*u.um)

>>> f, ax = plt.subplots()  # doctest: +IGNORE_OUTPUT
>>> ax.plot(x, y)  # doctest: +IGNORE_OUTPUT
>>> ax.plot(x, y_continuum_fitted)  # doctest: +IGNORE_OUTPUT
>>> ax.set_title("Continuum Fitting")  # doctest: +IGNORE_OUTPUT
>>> ax.grid(True)  # doctest: +IGNORE_OUTPUT
```

```python
>>> spec_normalized = spectrum / y_continuum_fitted

>>> f, ax = plt.subplots()  # doctest: +IGNORE_OUTPUT
>>> ax.plot(spec_normalized.spectral_axis, spec_normalized.flux)  # doctest: +IGNORE_OUTPUT
>>> ax.set_title("Continuum normalized spectrum")  # doctest: +IGNORE_OUTPUT
>>> ax.grid(True)  # doctest: +IGNORE_OUTPUT
```

--------------------------------

### Fit Spectral Lines with Excluded Regions

Source: https://github.com/astropy/specutils/blob/main/docs/fitting.rst

Demonstrates fitting a Gaussian model to a spectrum while ignoring a specific wavelength range using the SpectralRegion class.

```python
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from astropy.modeling import models
>>> from astropy import units as u
>>> from specutils.spectra import Spectrum, SpectralRegion
>>> from specutils.fitting import fit_lines

>>> # Create a simple spectrum with a Gaussian.
>>> 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)

>>> # Create the spectrum to fit
>>> spectrum = Spectrum(flux=y*u.Jy, spectral_axis=x*u.um)

>>> # Fit each peak
>>> 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
```

--------------------------------

### Create Synthetic Spectrum

Source: https://github.com/astropy/specutils/blob/main/docs/fitting.rst

Generates a synthetic spectrum with Gaussian emission and absorption lines, including noise. Requires numpy, astropy.modeling, astropy.units, and specutils.

```python
import numpy as np
from astropy.modeling import models
import astropy.units as u
from specutils import Spectrum, SpectralRegion

np.random.seed(42)
g1 = models.Gaussian1D(1, 4.6, 0.2)
g2 = models.Gaussian1D(2.5, 5.5, 0.1)
g3 = models.Gaussian1D(-1.7, 8.2, 0.1)
x = np.linspace(0, 10, 200)
y = g1(x) + g2(x) + g3(x) + np.random.normal(0., 0.2, x.shape)
spectrum = Spectrum(flux=y*u.Jy, spectral_axis=x*u.um)
```

--------------------------------

### Load a Spectral Cube

Source: https://github.com/astropy/specutils/blob/main/docs/spectral_cube.rst

Load a MaNGA cube from a remote URL into a Spectrum object.

```python
>>> from astropy.utils.data import download_file
>>> from specutils.spectra import Spectrum
>>> filename = "https://stsci.box.com/shared/static/28a88k1qfipo4yxc4p4d40v4axtlal8y.fits"
>>> file = download_file(filename, cache=True) # doctest: +REMOTE_DATA
>>> sc = Spectrum.read(file, format='MaNGA cube') # doctest: +REMOTE_DATA
```

```python
>>> sc.shape # doctest: +REMOTE_DATA
(4563, 74, 74)
```

```python
>>> sc[2000:2003, 30:33, 30:33] # doctest: +REMOTE_DATA
<Spectrum(flux=[[[0.4892023205757141 ... 0.5994223356246948]]] 1e-17 erg / (Angstrom s spaxel cm2) (shape=(3, 3, 3), mean=0.54165 1e-17 erg / (Angstrom s spaxel cm2)); spectral_axis=<SpectralAxis
   (observer to target:
      radial_velocity=0.0 km / s
      redshift=0.0)
  [5.73984286e-07 5.74116466e-07 5.74248676e-07] m> (length=3); uncertainty=InverseVariance)>
```

--------------------------------

### Spectrum1D Instantiation with Dispersion and Flux

Source: https://github.com/astropy/specutils/wiki/28Feb2012-Discussion

This is a suggested signature for instantiating a Spectrum1D object with dispersion and flux arrays. WCS will be an optional argument within NDData.

```python
_init_(dispersion,flux,dispersion_units,...ndddata...)
```

--------------------------------

### specutils.utils.wcs_utils

Source: https://github.com/astropy/specutils/blob/main/docs/wcs_utils.rst

Utility functions for spectral WCS conversions, including air-to-vacuum and vacuum-to-air transformations.

```APIDOC
## specutils.utils.wcs_utils

### Description
Provides functions for converting spectral values between air and vacuum. Includes refractive index calculations for air. Default method is Morton (2000) for consistency with IAU standards.

### Notes
- Methods other than Greisen (2006) have mathematical singularities in the far UV and are only valid for wavelengths > 200 nm.
- Greisen (2006) assumes an air temperature of 0C, while other methods assume 15C.
```

--------------------------------

### Read FITS file with 1D WCS

Source: https://github.com/astropy/specutils/blob/main/docs/custom_loading.rst

Load a spectrum from a FITS file using the Spectrum.read method with the wcs1d-fits format.

```python
import os
from specutils import Spectrum

file_path = os.path.join('path/to/folder', 'file_with_1d_wcs.fits')

spec = Spectrum.read(file_path, format='wcs1d-fits')
```

--------------------------------

### Run Codestyle Checks

Source: https://github.com/astropy/specutils/blob/main/docs/contributing.rst

Execute the codestyle checks using tox to ensure your code adheres to project standards before committing.

```shell
tox -e codestyle
```

--------------------------------

### Create SpectrumCollection from multi-dimensional arrays

Source: https://github.com/astropy/specutils/blob/main/docs/spectrum_collection.rst

Use this method when you have flux, spectral_axis, wcs, uncertainty, mask, and meta as multi-dimensional arrays. Ensure all input arrays have compatible shapes.

```python
import numpy as np
import astropy.units as u
from astropy.nddata import StdDevUncertainty
from specutils import SpectrumCollection
from specutils.utils.wcs_utils import gwcs_from_array

flux = u.Quantity(np.random.sample((5, 10)), unit='Jy')
spectral_axis = u.Quantity(np.arange(50).reshape((5, 10)), unit='AA')
wcs = np.array([gwcs_from_array(x, x.shape) for x in spectral_axis])
uncertainty = StdDevUncertainty(np.random.sample((5, 10)), unit='Jy')
mask = np.ones((5, 10)).astype(bool)
meta = [{'test': 5, 'info': [1, 2, 3]} for i in range(5)]

spec_coll = SpectrumCollection(
flux=flux, spectral_axis=spectral_axis, wcs=wcs,
uncertainty=uncertainty, mask=mask, meta=meta)

spec_coll.shape
spec_coll.flux.unit
spec_coll.spectral_axis.shape
spec_coll.spectral_axis.unit
```

--------------------------------

### Commit and Push Changes

Source: https://github.com/astropy/specutils/blob/main/docs/contributing.rst

Commit your changes with a descriptive message and push the branch to your GitHub repository. This prepares your contribution for a pull request.

```shell
git add .
git commit -m "Your detailed description of your changes."
git push origin name-of-your-bugfix-or-feature
```

--------------------------------

### Template Matching with Unknown Redshift

Source: https://github.com/astropy/specutils/blob/main/docs/analysis.rst

Performs template matching to find the best redshift for an observed spectrum against a template. Requires setting up spectral axes, observed spectrum, template spectrum, and redshift grid.

```python
>>> from specutils.analysis import template_comparison
>>> spec_axis = np.linspace(0, 50, 50) * u.AA
>>> observed_redshift = 2.0
>>> min_redshift = 1.0
>>> max_redshift = 3.0
>>> delta_redshift = .25
>>> resample_method = "flux_conserving"
>>> rs_values = np.arange(min_redshift, max_redshift+delta_redshift, delta_redshift)

>>> observed_spectrum = Spectrum(spectral_axis=spec_axis*(1+observed_redshift), flux=np.random.randn(50) * u.Jy, uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
>>> spectral_template = Spectrum(spectral_axis=spec_axis, flux=np.random.randn(50) * u.Jy, uncertainty=StdDevUncertainty(np.random.sample(50), unit='Jy'))
>>> tm_result = template_comparison.template_match(observed_spectrum=observed_spectrum, spectral_templates=spectral_template, resample_method=resample_method, redshift=rs_values) # doctest:+FLOAT_CMP
```

--------------------------------

### Extract using spectral_slab

Source: https://github.com/astropy/specutils/blob/main/docs/spectral_regions.rst

Uses spectral_slab as an alternative entry point for extracting a range from a spectrum.

```python
>>> from astropy import units as u
>>> import numpy as np
>>> from specutils import Spectrum, SpectralRegion
>>> from specutils.manipulation import spectral_slab

>>> spectrum = Spectrum(spectral_axis=np.arange(1, 50) * u.nm,
...                       flux=np.random.default_rng(12345).random(49)*u.Jy)
>>> sub_spectrum = spectral_slab(spectrum, 8*u.nm, 20*u.nm)
>>> sub_spectrum.spectral_axis
<SpectralAxis [ 8.,  9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20.] nm>
```

--------------------------------

### Estimate Gaussian Line Parameters

Source: https://github.com/astropy/specutils/blob/main/docs/fitting.rst

Extracts a spectral region and estimates parameters for a Gaussian1D model.

```python
>>> from specutils import SpectralRegion
>>> from specutils.fitting import estimate_line_parameters
>>> from specutils.manipulation import extract_region

>>> sub_region = SpectralRegion(4*u.um, 5*u.um)
>>> sub_spectrum = extract_region(spectrum, sub_region)
```

```python
>>> result = estimate_line_parameters(sub_spectrum, models.Gaussian1D())
>>> print(result.amplitude)  # doctest: +FLOAT_CMP
Parameter('amplitude', value=1.1845669151078486, unit=Jy)

>>> print(result.mean)  # doctest: +FLOAT_CMP
Parameter('mean', value=4.57517271067525, unit=um)

>>> print(result.stddev)  # doctest: +FLOAT_CMP
Parameter('stddev', value=0.19373372929165977, unit=um, bounds=(1.1754943508222875e-38, None))
```

--------------------------------

### Create a New Branch for Development

Source: https://github.com/astropy/specutils/blob/main/docs/contributing.rst

Create a new branch for your bug fixes or feature development. This helps in organizing changes and managing contributions.

```shell
git checkout -b name-of-your-bugfix-or-feature
```

--------------------------------

### Write Spectrum to FITS

Source: https://github.com/astropy/specutils/blob/main/docs/spectrum.rst

Saves a Spectrum object to a FITS file. Explicitly passing the format parameter is recommended for better control over the file type.

```python
>>> spec1d.write("/path/to/output.fits")  # doctest: +SKIP
```

--------------------------------

### Normalize Continuum and Calculate Equivalent Width

Source: https://github.com/astropy/specutils/blob/main/docs/index.rst

Normalizes the spectrum using a generic continuum fit and calculates the equivalent width for a specific spectral region.

```python
>>> import warnings
>>> from specutils.fitting import fit_generic_continuum
>>> with warnings.catch_warnings():  # Ignore warnings
...     warnings.simplefilter('ignore')
...     cont_norm_spec = spec / fit_generic_continuum(spec)(spec.spectral_axis) # doctest: +REMOTE_DATA

>>> f, ax = plt.subplots()  # doctest: +IGNORE_OUTPUT +REMOTE_DATA
>>> ax.step(cont_norm_spec.wavelength, cont_norm_spec.flux)  # doctest: +IGNORE_OUTPUT +REMOTE_DATA
>>> ax.set_xlim(654 * u.nm, 660 * u.nm)  # doctest: +IGNORE_OUTPUT +REMOTE_DATA

>>> from specutils import SpectralRegion
>>> from specutils.analysis import equivalent_width
>>> equivalent_width(cont_norm_spec, regions=SpectralRegion(6562 * u.AA, 6575 * u.AA)) # doctest: +REMOTE_DATA +FLOAT_CMP
<Quantity -14.82013888 Angstrom>
```

--------------------------------

### Read Spectrum from URL

Source: https://github.com/astropy/specutils/blob/main/docs/spectrum.rst

Load a Spectrum object directly from a URL pointing to a FITS file. The format string is required when reading from non-file objects like URL streams.

```python
from specutils import Spectrum
import urllib
spec = urllib.request.urlopen('https://data.sdss.org/sas/dr14/sdss/spectro/redux/26/spectra/0751/spec-0751-52251-0160.fits') # doctest: +REMOTE_DATA
Spectrum.read(spec, format="SDSS-III/IV spec") # doctest: +REMOTE_DATA
```

--------------------------------

### Estimate Spectral Widths with specutils.analysis

Source: https://github.com/astropy/specutils/blob/main/docs/analysis.rst

Provides functions to estimate spectral widths. gaussian_sigma_width and gaussian_fwhm assume a Gaussian shape. fwhm and fwzi do not assume a Gaussian shape. The analytic argument can be set to False to use Monte Carlo distributions for uncertainty.

```python
>>> from specutils.analysis import gaussian_sigma_width, gaussian_fwhm, fwhm, fwzi
>>> gaussian_sigma_width(noisy_gaussian) # doctest: +FLOAT_CMP
<Quantity 0.74075431 GHz>
```

```python
>>> gaussian_fwhm(noisy_gaussian) # doctest: +FLOAT_CMP
<Quantity 1.74434311 GHz>
```

```python
>>> fwhm(noisy_gaussian) # doctest: +FLOAT_CMP
<Quantity 1.86047666 GHz>
```

```python
>>> fwzi(noisy_gaussian) # doctest: +FLOAT_CMP
<Quantity 94.99997484 GHz>
```

--------------------------------

### Create a 1D Spectrum Object

Source: https://github.com/astropy/specutils/blob/main/docs/spectrum.rst

Create a Spectrum object from NumPy arrays or Astropy Quantities for flux and spectral axis. The spectral_axis can be ascending or descending but must be monotonic.

```python
import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt
from specutils import Spectrum
flux = np.random.randn(200)*u.Jy
wavelength = np.arange(5100, 5300)*u.AA
spec1d = Spectrum(spectral_axis=wavelength, flux=flux)
ax = plt.subplots()[1]  # doctest: +SKIP
ax.plot(spec1d.spectral_axis, spec1d.flux)  # doctest: +SKIP
ax.set_xlabel("Dispersion")  # doctest: +SKIP
ax.set_ylabel("Flux")  # doctest: +SKIP
```

--------------------------------

### Compute moment maps for wavelength ranges

Source: https://github.com/astropy/specutils/blob/main/docs/spectral_cube.rst

Extracts sub-cubes using spectral_slab and calculates moment maps for H-alpha line analysis.

```python
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.utils.data import download_file
from specutils import Spectrum, SpectralRegion
from specutils.analysis import moment
from specutils.manipulation import spectral_slab

filename = "https://stsci.box.com/shared/static/28a88k1qfipo4yxc4p4d40v4axtlal8y.fits"
fn = download_file(filename, cache=True)
spec1d = Spectrum.read(fn)

# Extract H-alpha sub-cube for moment maps using spectral_slab
subspec = spectral_slab(spec1d, 6745.*u.AA, 6765*u.AA)
ha_wave = subspec.spectral_axis

# Extract wider sub-cube covering H-alpha and [N II] using spectral_slab
subspec_wide = spectral_slab(spec1d, 6705.*u.AA, 6805*u.AA)
ha_wave_wide= subspec_wide.spectral_axis

# Convert flux density to microJy and correct negative flux offset for
# this particular dataset
ha_flux = (np.sum(subspec.flux.value, axis=(1,2)) + 0.0093) * 1.0E-6*u.Jy
ha_flux_wide = (np.sum(subspec_wide.flux.value, axis=(1,2)) + 0.0093) * 1.0E-6*u.Jy

# Compute moment maps for H-alpha line
moment0_halpha = moment(subspec, order=0)
moment1_halpha = moment(subspec, order=1)

# Convert moment1 from AA to velocity
# H-alpha is redshifted to 6755 AA for this galaxy
print(moment1_halpha[40,40])
vel_map = 3.0E5 * (moment1_halpha.value - 6.755E-7) / 6.755E-7

# Plot results in 3 panels (subspec_wide,  H-alpha line flux, H-alpha velocity map)
f,(ax1,ax2,ax3) = plt.subplots(1, 3, figsize=(15, 5))
ax1.plot(ha_wave_wide, (ha_flux_wide)*1000.)
ax1.set_xlabel('Angstrom', fontsize=14)
ax1.set_ylabel('uJy', fontsize=14)
ax1.tick_params(axis="both", which='major', labelsize=14, length=8, width=2, direction='in', top=True, right=True)
ax2.imshow(moment0_halpha.value, origin='lower')
ax2.set_title('moment = 0')
ax2.set_xlabel('x pixels', fontsize=14)
ax3.imshow(vel_map, vmin=-100., vmax=100., cmap='rainbow', origin='lower')
ax3.set_title('moment = 1')
ax3.set_xlabel('x pixels', fontsize=14)
```

--------------------------------

### Crop Spectrum using World Coordinates with WCS

Source: https://github.com/astropy/specutils/blob/main/docs/spectrum.rst

Crop a Spectrum object based on world coordinates, suitable for spectral cubes with spatial dimensions. Requires a WCS object and astropy.coordinates objects for the lower and upper bounds of the desired region.

```python
from astropy.coordinates import SpectralCoord, SkyCoord
from astropy import units as u
from astropy.wcs import WCS

w = WCS({'WCSAXES': 3, 'CRPIX1': 38.0, 'CRPIX2': 38.0, 'CRPIX3': 1.0,
         'CRVAL1': 205.4384, 'CRVAL2': 27.004754, 'CRVAL3': 4.890499866509344,
         'CTYPE1': 'RA---TAN', 'CTYPE2': 'DEC--TAN', 'CTYPE3': 'WAVE',
         'CUNIT1': 'deg', 'CUNIT2': 'deg', 'CUNIT3': 'um',
         'CDELT1': 3.61111097865634E-05, 'CDELT2': 3.61111097865634E-05, 'CDELT3': 0.001000000047497451,
         'PC1_1 ': -1.0, 'PC1_2 ': 0.0, 'PC1_3 ': 0,
         'PC2_1 ': 0.0, 'PC2_2 ': 1.0, 'PC2_3 ': 0,
         'PC3_1 ': 0, 'PC3_2 ': 0, 'PC3_3 ': 1,
         'DISPAXIS': 2, 'VELOSYS': -2538.02,
         'SPECSYS': 'BARYCENT', 'RADESYS': 'ICRS', 'EQUINOX': 2000.0,
         'LONPOLE': 180.0, 'LATPOLE': 27.004754})
spec = Spectrum(flux=np.random.default_rng(12345).random((20, 5, 10)) * u.Jy, wcs=w)  # doctest: +IGNORE_WARNINGS
lower = [SkyCoord(ra=205, dec=26, unit=u.deg), SpectralCoord(4.9, unit=u.um)]
upper = [SkyCoord(ra=205.5, dec=27.5, unit=u.deg), SpectralCoord(4.9, unit=u.um)]
spec.crop(lower, upper)  # doctest: +IGNORE_WARNINGS +FLOAT_CMP
```