<jupyter_start><jupyter_text>Masking and flattening a spectrumOften in high resolution spectroscopy we wish to designate the *continuum*---regions in the spectrum devoid of spectral lines.  The continuum is conspicuous in hot stars, but becomes ill-defined in cool M-dwarfs.  In this tutorial we show how to define the continuum using `muler` and `specutils`.<jupyter_code>from muler.hpf import HPFSpectrum
import numpy as np
import matplotlib.pyplot as plt
%config InlineBackend.figure_format='retina'
plt.rcParams['axes.facecolor']='white'
plt.rcParams['savefig.facecolor']='white'
filename = '../../../muler_example_data/HPF/01_A0V_standards/Goldilocks_20210801T083618_v1.0_0036.spectra.fits'
raw_spectrum = HPFSpectrum(file = filename, order=6)
spectrum = raw_spectrum.sky_subtract(method='vector').remove_nans().deblaze().normalize()
spectrum.plot(ylo=0);<jupyter_output><empty_output><jupyter_text>Look at those conspicuous sharp spectral lines and the smooth dip near 10940 Angstroms.  Let's divide the spectrum by a tuned model.<jupyter_code>from gollum.phoenix import PHOENIXSpectrum
native_template = PHOENIXSpectrum(teff=9600, logg=5.0).normalize()
doctored_template = native_template\
                        .rotationally_broaden(105.0)\
                        .rv_shift(-40.0)\
                        .resample(spectrum)\
                        .normalize()
ax = spectrum.plot(ylo=0, label='A0V Star');
doctored_template.plot(ax=ax, label='Doctored A0V Template'); ax.legend();
telluric_response = spectrum / doctored_template<jupyter_output><empty_output><jupyter_text>Normalize by the top 10th percentile, rather than the median...<jupyter_code>spectrum_peak = np.percentile(telluric_response.flux, 65)
telluric_response = telluric_response.divide(spectrum_peak)
ax = telluric_response.plot(ylo=0);
ax.axhline(1.0, color='k', linestyle='dotted');<jupyter_output><empty_output><jupyter_text>The continuum should all hover near $1.0$.  You can see that the spectrum has a slight downward curvature to it, with the blue-most values about $5\%$ lower than $1.0$, and some red-end values exceeding $1.0$ by 2 percentage points.  Let's mask all the lines, fit a smooth function to the continuum, and then flatten by that function, evaluated at all wavelength points.  Masking the lines and fitting the continuum<jupyter_code>import astropy.units as u
from specutils import SpectralRegion
from specutils.manipulation import extract_region
from specutils.fitting import find_lines_derivative, fit_continuum
continuum_subtracted = telluric_response-1.0<jupyter_output><empty_output><jupyter_text>Here we find all the absorption lines above some threshold.<jupyter_code>lines = find_lines_derivative(continuum_subtracted, flux_threshold=0.05)
lines = lines[lines['line_type'] == 'absorption']
len(lines)
ax = telluric_response.plot(ylo=0);
ax.axhline(1.0, color='k', linestyle='dashed');

for line in lines:
    ax.axvline(line['line_center'].value, linestyle='dotted', color='k', alpha=0.5)<jupyter_output><empty_output><jupyter_text>Neat, we see that we have found the most conspicuous lines.  Let's mask these lines with a $\pm1\;$A width.<jupyter_code>line_regions = SpectralRegion.from_line_list(lines, width=2.0*u.Angstrom)<jupyter_output><empty_output><jupyter_text>Fit the continuum while excluding the spectral lines.<jupyter_code>%%capture
continuum_model = fit_continuum(telluric_response, exclude_regions=line_regions)
continuum_fit = continuum_model(telluric_response.wavelength)
continuum_spectrum = HPFSpectrum(flux=continuum_fit, 
                                 spectral_axis = telluric_response.wavelength)
ax=telluric_response.plot(ylo=0, label='Telluric response')
ax.axhline(1.0, color='k', linestyle='dashed');
continuum_spectrum.plot(ax=ax, label='Continuum fit');
ax.legend();<jupyter_output><empty_output><jupyter_text>Nice!  That looks OK.  You can see that we are fitting to weak telluric lines at 10910 A, which pulls down the fit from 10890.  So it's not perfect.<jupyter_code>flattened_telluric = telluric_response/continuum_spectrum
ax=flattened_telluric.plot(ylo=0)
ax.axhline(1.0, color='k', linestyle='dashed');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1029868578255675, 'fraction_numerical': 0.021505376344086023, 'mean_word_length': 5.91900826446281, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '244426', 'license': 'MIT', 'n_tokens_mistral': 1456, 'n_tokens_neox': 1318, 'n_words': 329, 'path': 'docs/tutorials/Masking_and_flattening.ipynb', 'repo_name': 'OttoStruve/muler', 'size': '4185'}