Refined HPF sky subtraction

The science and sky fibers in HPF have slightly different throughputs. In this notebook we show how muler can automatically apply a calibration correction based on twilight sky flats. We believe this change can reduce sky residuals by up to \(14\times\) over the naive subtraction of the sky fiber.

%config InlineBackend.figure_format='retina'
from muler.hpf import HPFSpectrum

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 2
      1 get_ipython().run_line_magic('config', "InlineBackend.figure_format='retina'")
----> 2 from muler.hpf import HPFSpectrum

File ~/checkouts/readthedocs.org/user_builds/muler/envs/latest/lib/python3.8/site-packages/muler/hpf.py:14
     12 import warnings
     13 import logging
---> 14 from muler.echelle import EchelleSpectrum, EchelleSpectrumList
     15 from muler.utilities import resample_list
     16 import numpy as np

File ~/checkouts/readthedocs.org/user_builds/muler/envs/latest/lib/python3.8/site-packages/muler/echelle.py:31
     29 from scipy.ndimage import median_filter, gaussian_filter1d
     30 import specutils
---> 31 from muler.utilities import apply_numpy_mask, is_list, resample_list
     34 # from barycorrpy import get_BC_vel
     35 from astropy.coordinates import SkyCoord, EarthLocation

File ~/checkouts/readthedocs.org/user_builds/muler/envs/latest/lib/python3.8/site-packages/muler/utilities.py:15
     13 from astropy.convolution import convolve, Gaussian1DKernel
     14 from scipy.ndimage import binary_dilation
---> 15 from astroquery.simbad import Simbad
     16 Simbad.add_votable_fields('V', 'B', 'J', 'H', 'K', 'parallax')
     17 LinInterpResampler = LinearInterpolatedResampler()

ModuleNotFoundError: No module named 'astroquery'

example_file = '../../data/HPF/Goldilocks_20210517T054403_v1.0_0060.spectra.fits'
spectrum = HPFSpectrum(file=example_file, order=19)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[2], line 2
      1 example_file = '../../data/HPF/Goldilocks_20210517T054403_v1.0_0060.spectra.fits'
----> 2 spectrum = HPFSpectrum(file=example_file, order=19)

NameError: name 'HPFSpectrum' is not defined

spectrum = spectrum.normalize()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 1
----> 1 spectrum = spectrum.normalize()

NameError: name 'spectrum' is not defined

We see the conspicuous blaze pattern in each spectra order. Let’s apply a correction to one of those orders:

ax= spectrum.sky.plot(ylo=0, yhi=0.1);

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[4], line 1
----> 1 ax= spectrum.sky.plot(ylo=0, yhi=0.1);

NameError: name 'spectrum' is not defined

These are the skylines we wish to remove. Since the sky emission is additive it should always be the first step before dealing with multiplicative terms, such as the blaze or telluric absorption.

ax = spectrum.plot(ylo=0.6, yhi=1.3, label='Observed')
spectrum.sky_subtract(method='scalar').plot(ax=ax, label='Sky Subtracted')
ax.legend();

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[5], line 1
----> 1 ax = spectrum.plot(ylo=0.6, yhi=1.3, label='Observed')
      2 spectrum.sky_subtract(method='scalar').plot(ax=ax, label='Sky Subtracted')
      3 ax.legend();

NameError: name 'spectrum' is not defined

You can see that that the sky subtraction successfully removes the emission lines. We can now proceed with deblazing:

spectrum.sky_subtract(method='scalar').deblaze().normalize().plot();

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[6], line 1
----> 1 spectrum.sky_subtract(method='scalar').deblaze().normalize().plot();

NameError: name 'spectrum' is not defined

Awesome! Let’s see if we can spot the differences among the three sky subtraction methods:

ax = spectrum.deblaze().normalize().plot(label='observed');
for method in ['naive', 'scalar', 'vector']:
    spectrum.sky_subtract(method=method).deblaze().normalize().plot(ax=ax, label=method);

ax.set_xlim(10_830, 10_840)
ax.set_ylim(0.98, 1.09)
ax.legend();

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[7], line 1
----> 1 ax = spectrum.deblaze().normalize().plot(label='observed');
      2 for method in ['naive', 'scalar', 'vector']:
      3     spectrum.sky_subtract(method=method).deblaze().normalize().plot(ax=ax, label=method);

NameError: name 'spectrum' is not defined

We see that the naive method oversubtracts the sky line, as we expected. The scalar and vector methods appear to perform equally well for this source in this wavelength range.