# Learning Goals¶

This Notebook is designed to walk the user (you) through: Working with the COS Line Spread Function (LSF) to simulate, fit, or model COS observations.

1. Understanding the Line Spread Function

- 1.1. What is an LSF?

2. Taking a look at the LSF kernels

- 2.1. Reading in an LSF

- 2.2. Plotting an LSF kernel

3. Convolving an LSF

# 0. Introduction¶

The Cosmic Origins Spectrograph (COS) is an ultraviolet spectrograph on-board the Hubble Space Telescope(HST) with capabilities in the near ultraviolet (NUV) and far ultraviolet (FUV).

This tutorial aims to prepare you to work with the COS data of your choice by walking you through convolving a template or high-resolution spectrum with the COS LSF - both for FUV and NUV data. This is an important step for many spectral applications, notably fitting spectral lines.

## We will import the following packages:¶

• numpy to handle arrays and functions
• astropy.io fits and astropy.table Table for accessing FITS files
• astropy.modeling functional_models to generate synthetic spectral line shapes
• astropy.convolution convolve for a convolution of two discretized functions
• astroquery.mast Mast and Observations for finding and downloading data from the MAST archive
• scipy.interpolate interp1d for interpolating two discretized functions to the same sampling
• matplotlib.pyplot for plotting
• urllib for downloading files stored at a web address online
• tarfile for unzipping some compressed STIS data
• pathlib Path for managing system paths

If you have an existing astroconda environment, it may or may not already have the necessary packages to run this Notebook. To create a Python environment capable of running all the data analyses in these COS Notebooks, please see Section 1 of our Notebook tutorial on setting up an environment.

In [1]:
# Import for: array manipulation
import numpy as np

# Import for: reading fits files
from astropy.table import Table
from astropy.io import fits

# Import for: line profile functions
from astropy.modeling import functional_models

# Import for: convolutions
from astropy.convolution import convolve

from astroquery.mast import Observations

#Import for: interpolating a discretized function
from scipy.interpolate import interp1d

#Import for: plotting
from matplotlib import pyplot as plt

import urllib

#Import for: unzipping a compressed file of STIS data
import tarfile

#Import for: working with system paths
from pathlib import Path


## We will also define a few directories we will need:¶

In [2]:
# These will be important directories for the Notebook

cwd = Path(".")
outputdir = cwd / "output"
plotsdir = cwd / "output" / "plots"

# Make the directories if they don't exist

Out[2]:
(None, None, None)

## And we will need to download the data we wish to filter and analyze¶

We choose the exposures with the association obs_id: LCRS51010, because we know there is also STIS data on the source, AV75. For more information on downloading COS data, see our notebook tutorial on downloading COS data.

In [3]:
# Search for obs_id of files, generate data product list
pl = Observations.get_product_list(
Observations.query_criteria(obs_id="LCRS51010")
)

)


Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HST/product/lcrs51010_x1dsum.fits to data/mastDownload/HST/lcrs51010/lcrs51010_x1dsum.fits ... [Done]


# 1. Understanding the Line Spread Function¶

Each COS observation is taken with a specific grating and a specific central wavelength setting (cenwave). Each such configuration (grating and cenwave) and - thus each COS dataset - has a set of corresponding Line Spread Functions (LSFs).

## 1.1. What is an LSF?¶

A Line Spread Function (LSF) is a model of a spectrograph's response to a monochromatic light source: it explains how an infinitely thin spectral line (a delta function) input into the spectrograph would be output at the focal plane. A spectrum with infinite spectral resolution convolved with the COS LSF will reproduce the COS spectral line profiles:

$$True\ Input\ Spectrum \ast COS\ LSF = True\ Output\ Spectrum$$

Of course, we don't have access to a "true", infinite resolution input spectrum, nor can we know the infinite resolution LSF. The best we can do is use a model spectrum, or a spectrum from a higher resolution spectrograph (STIS often works well) and convolve it with a kernel of our LSF model. Convolving these yields our model output spectrum:

$$Model\ Input\ Spectrum \ast COS\ LSF\ Model = Model\ Output\ Spectrum$$

Note that there is a corresponding spread function in the cross-dispersion direction, known as the cross-dispersion spread function (CDSF). These functions are also modelled by the COS team and hosted on their website listed below. Working with the CDSFs is out of the scope of this Notebook.

How does COS handle the LSF?

The COS LSFs are generated using an optical model of the spectrograph in the program Code V, and are then validated using real spectral data obtained with the instrument. The LSF kernels are sampled at regular intervals over the wavelength range of each COS configuration. Each of these wavelengths thus corresponds to an LSF kernel. The kernel size (in COS "pixels") of each LSF varies for the near ultraviolet (NUV) and far ultraviolet (FUV) modes, as does the sample rate (how many Angstroms between sampled kernels). These values are shown below:

### Table 1.1: Line Spread Function (LSF) kernel parameters for the COS instrument¶

LSF kernels sampled every... (Å) Actual kernel size (pixels) Approximate kernel size in wavelength space(Å)
COS/NUV 100 Å 101 pixels -
NUV M Gratings 100 Å 101 pixels 3.7 Å
NUV L Gratings 100 Å 101 pixels 39 Å
COS/FUV 5 Å 321 pixels -
FUV M Gratings 5 Å 321 pixels 3.5 Å
FUV L Gratings 5 Å 321 pixels 26 Å

How are the LSF files structured?

In short, the LSF files are structured as a list of LSF profile kernels. It begins with a space-separated list of central sample wavelengths and the input monochromatic lines used by Code V. These are not to be confused with COS central wavelength settings, or cenwaves. To avoid confusing the two concepts, we'll call the central wavelengths of the LSF kernels the LSF wavelengths.

## 1.2. Getting the files you need¶

Which LSF files will you need?

This depends on your data's parameters, specifically those listed in Table 1.2.

Parameter Corresponding Header Keyword Example used in this Notebook
Wavelength range DETECTOR FUV
Grating OPT_ELEM G130M
Central wavelength setting CENWAVE 1291

If your data was taken in COS' NUV configuration, there's only one NUV LSF file which contains all the LSF profiles for the entire NUV. FUV data is more complicated: you will need to choose an LSF file based on the COS Lifetime position, grating, and central wavelength setting used to capture your data.

You'll also need to get the DISPTAB file associated with your data. This Dispersion Coefficient Table, or DISPTAB, file contains polynomial coefficients to convert from pixel number to wavelength. You can find the relevant DISPTAB in your data's primary fits header. You can figure out the relevant DISPTAB file by examining your data's header in DS9 or other software, but we'll briefly do this programmatically below.

Note that you should always ensure that your LSF kernels are normalized such that they integrate to 1. The files we use in this Notebook are normalized to this contraint.

In [4]:
fuvHeader0 = fits.getheader(fuvFile, ext=0)  # Select the primary header
print(f"For the file {fuvFile}, the relevant parameters are: ")
param_dict = {}  # Make a dict to store what you find here

for hdrKeyword in [
"DETECTOR",
"OPT_ELEM",
"CENWAVE",
"DISPTAB",
]:  # Print out the relevant values
try:  # For DISPTAB
value = fuvHeader0[hdrKeyword].split("$")[1] # Save the key/value pairs to the dictionary param_dict[hdrKeyword] = value # DISPTAB needs the split here except: # For other params value = fuvHeader0[hdrKeyword] # Save the key/value pairs to the dictionary param_dict[hdrKeyword] = value print(f"{hdrKeyword} = {value}") # Print the key/value pairs  For the file data/mastDownload/HST/lcrs51010/lcrs51010_x1dsum.fits, the relevant parameters are: DETECTOR = FUV OPT_ELEM = G130M LIFE_ADJ = 3 CENWAVE = 1291 DISPTAB = 5b919205l_disp.fits  Now that you know which LSF files you need, we'll go gather them from the COS team's website: Where can you find the LSF files? The COS team maintains up-to-date LSF files on the COS Spectral Resolution page. Opening up this link leads to a page like that shown in Fig. 1.1, where the LSF files are discussed in detail. The bottom part of this page has links to all the relavent files. The links at the top of the page will take you to the relevant section. In Fig. 1.1, we have circled in black the link to the section pertaining to our data: FUV at the Lifetime Position: 3. ### Fig 1.1: Screenshot of the COS Spectral Resolution Site¶ Clicking on the circled link takes us to the table of hyperlinks to all the files perataining to data taken with the FUV, Lifetime Postition 3 configutation, shown in Fig. 1.2: ### Fig 1.2: Screenshot of the COS Spectral Resolution Site - Focus on LP-POS 3¶ Circled in solid red is the button to download the LSF file we need for our data with CENWAVE = 1291. Circled in dashed black is the corresponding CDSF. You can click on the solid red circled button to download the LSF file and move it to this current working directory or right-click and select "Download Linked File As..." and download directly to this directory. • This step may look somewhat different depending on your browser • If you're unsure of the current working directory, un-comment the print statement in the next cell, and run the cell. In [5]: # equivalent of !pwd # print(cwd.resolve())  These files can also be downloaded programmatically from with python, with the function we define below: In [6]: def fetch_files(det, grating, lpPos, cenwave, disptab): """ Given all the inputs: (detector, grating, LP-POS, cenwave, dispersion table,) this will download both the LSF file and Disptab file you should use in the convolution and return their paths. Returns: LSF_file_name (str): filename of the new downloaded LSF file disptab_path (str): path to the new downloaded disptab file """ COS_site_rootname = ( "https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/cos/" "performance/spectral-resolution/_documents/" ) # Link to where all the files live - split into 2 lines if det == "NUV": # Only one file for NUV LSF_file_name = "nuv_model_lsf.dat" elif det == "FUV": # FUV files follow a naming pattern LSF_file_name = f"aa_LSFTable_{grating}_{cenwave}_LP{lpPos}_cn.dat" LSF_file_webpath = COS_site_rootname + LSF_file_name # Where to find file online urllib.request.urlretrieve( LSF_file_webpath, str(datadir / LSF_file_name) ) # Where to save file to locally print(f"Downloaded LSF file to {str(datadir/ LSF_file_name)}") # And we'll need to get the DISPTAB file as well disptab_path = str(datadir / disptab) urllib.request.urlretrieve( f"https://hst-crds.stsci.edu/unchecked_get/references/hst/{disptab}", disptab_path, ) print(f"Downloaded DISPTAB file to {disptab_path}") return LSF_file_name, disptab_path  We'll run this function with the parameters we saved earlier to gather the proper LSF file: In [7]: # We'll pass that fetch function the parameters we determined earlier # This "unpacked argument" phrasing works because of the order in which we added the params to the dict LSF_file_name, disptab_path = fetch_files(*list(param_dict.values()))  Downloaded LSF file to data/aa_LSFTable_G130M_1291_LP3_cn.dat Downloaded DISPTAB file to data/5b919205l_disp.fits  # 2. Taking a look at the LSF kernels¶ ## 2.1. Reading in an LSF¶ Below, we create a simple function to read in an LSF file as an astropy table: In [8]: def read_lsf(filename): # This is the table of all the LSFs: called "lsf" # The first column is a list of the wavelengths corresponding to the line profile, so we set our header accordingly if "nuv_" in filename: # If its an NUV file, header starts 1 line later ftype = "nuv" else: # assume its an FUV file ftype = "fuv" hs = 0 lsf = Table.read(filename, format="ascii", header_start=hs) # This is the range of each LSF in pixels (for FUV from -160 to +160, inclusive) # middle pixel of the lsf is considered zero ; center is relative zero pix = np.arange(len(lsf)) - len(lsf) // 2 # integer division to yield whole pixels # the column names returned as integers. lsf_wvlns = np.array([int(float(k)) for k in lsf.keys()]) return lsf, pix, lsf_wvlns  Now let's read the LSF file we downloaded and display the first 5 columns: Because this is an FUV file, the columns are 321 values long, corresponding to the 321 pixel size of an LSF kernal in the FUV. In [9]: lsf, pix, lsf_wvlns = read_lsf(str(datadir / LSF_file_name)) lsf[lsf.colnames[:5]]  Out[9]: Table length=321 11341139114411491154 float64float64float64float64float64 4.885797114010963e-064.787755984607906e-064.713489842513227e-064.735199692705489e-064.7784485923017495e-06 4.916023102906068e-064.873309105308387e-064.752957668230419e-064.749528365900715e-064.8460096205578055e-06 4.941950149310319e-064.954122351536572e-064.829737664959035e-064.819638005845053e-064.959435869169078e-06 4.975031355401172e-065.043803066313383e-064.943558000555235e-064.944166202963326e-065.105983074318574e-06 5.024945714942217e-065.132957330887238e-065.08655216831887e-065.110203690134147e-065.257588174065823e-06 5.1062647374008795e-065.226843223712312e-065.249334040421822e-065.292011369183484e-065.4004601352656915e-06 5.22334211429479e-065.329101494747061e-065.421819580311945e-065.476126509683101e-065.537769117220119e-06 5.363562444717191e-065.438680038755283e-065.589977083978292e-065.661455266614987e-065.679786309281309e-06 5.503054802394619e-065.554937875874203e-065.742187109577016e-065.841583683466683e-065.828667384063064e-06 5.624258813784741e-065.675603233615175e-065.872209611591034e-065.9970180710284226e-065.968283703978237e-06 ............... 5.538829078895442e-065.61452382644541e-065.785171777030175e-065.938726800824539e-066.006545611677446e-06 5.4710939230511685e-065.533244491994195e-065.726741058629339e-065.822288317661374e-065.894784543374904e-06 5.387839481809153e-065.458103038764225e-065.6470377934756154e-065.691283605758636e-065.750873288909732e-06 5.2992064517062305e-065.387313076900991e-065.540020267903346e-065.551977719200976e-065.5868670831158e-06 5.221160697526395e-065.320366681451691e-065.415092210912367e-065.402179516688138e-065.4169404397396914e-06 5.164281884145832e-065.2577636408949785e-065.289233084543136e-065.246140365509941e-065.258296979675213e-06 5.127032374917626e-065.1935627777544765e-065.173617813851988e-065.106882557229389e-065.115096193521994e-06 5.0992742707624055e-065.125375786306556e-065.076003036824709e-065.01082452583297e-064.988808887440265e-06 5.061823427754464e-065.04372509224167e-064.984863851741219e-064.952955359048715e-064.883922635184996e-06 0.00.00.00.00.0 ## 2.2. Plotting an LSF kernel¶ We can plot the LSF kernels themselves and take a look. Here's a very simple plot of the first LSF kernel sampled at 1134 Å: In [10]: # Set up the figure plt.figure(figsize=(10, 8)) # Fill the figure with a plot of the data plt.plot(pix, lsf["1134"]) # Give figure the title and labels plt.title("Figure 2.1\nCOS FUV LSF kernel sampled at 1134 Å", size=20) plt.xlabel("Pixel", size=20) plt.ylabel("Flux [$normalized,unitless$]", size=20) # format and save the figure plt.tight_layout() plt.savefig(str(plotsdir / "oneKernel.png"), bbox_inches="tight")  We'll now make a more complex plot showing all of the lines contained in the LSF file, distributed in wavelength space. Note that while the centers of the lines are correctly mapped to the LSF wavelength at which they were sampled, their kernel width in pixels has only been roughly translated to a wavelength range in Å. In short, the x-axis is not strictly to-scale. This will be rectified by the remapping shown in Fig 3.1. In [11]: fig, (ax0, ax1) = plt.subplots(2, 1, figsize=(10, 8), dpi=100) # Loop through the lsf kernels for i, col in enumerate(lsf.colnames): line_wvln = int(col) # central position contents = lsf[col].data # actual shape if line_wvln < 1291: # split into the FUVB segment xrange = 0.00997 * pix + line_wvln # ROUGHLY convert pix to wvln ax1.plot(xrange, contents, linewidth=0.5) # Plot the kernel ax1.axvline(line_wvln, linewidth=0.1) # plot the LSF wvln as a faded line elif line_wvln > 1291: # split into the FUVA segment xrange = 0.00997 * pix + line_wvln # ROUGHLY convert pix to wvln ax0.plot(xrange, contents, linewidth=0.5) # Plot the kernel ax0.axvline(line_wvln, linewidth=0.1) # plot the LSF wvln as a faded line # Format the figure ax0.set_xlim(1290, 1435) ax1.set_xlim(1125, 1295) # Add labels, titles and text fig.suptitle(f"Fig 2.2\nAll the LSF kernels in {LSF_file_name}", fontsize=20) ax1.set_title("FUVB Detector") ax0.set_title("FUVA Detector") ax1.set_xlabel("Approximate Wavelength [$\AA$]", size=14) fig.text(s="Flux [$normalized,unitless$]", x=-0.01, y=0.35, rotation="vertical", size=14) plt.tight_layout() # Save the figure plt.savefig(str(plotsdir / "allKernels.png"), bbox_inches="tight")  ### We will also demonstrate all of the above for a NUV dataset:¶ First, downloading the NUV data and getting the relevant header parameters: In [12]: #### NUV: download the data pl = Observations.get_product_list( Observations.query_criteria(obs_id="LBY606010") ) # search for file, generate data product list download = Observations.download_products( pl[pl["productSubGroupDescription"] == "X1DSUM"], # filter and download searched files download_dir=str(datadir), ) nuvFile = download["Local Path"][0] # Give the program the path to your downloaded data #### NUV: get relavent params nuvHeader0 = fits.getheader(nuvFile, ext=0) # Select the primary header print(f"For the file {nuvFile}, the relevant parameters are: ") nuv_param_dict = {} # Make a dict to store what you find here for hdrKeyword in [ "DETECTOR", "OPT_ELEM", "LIFE_ADJ", "CENWAVE", "DISPTAB", ]: # print out the relevant values try: # For DISPTAB value = nuvHeader0[hdrKeyword].split("$")[1]  # Save the key/value pairs to the dictionary,
nuv_param_dict[hdrKeyword] = value            # DISPTAB needs the split here
except:  # For other params
value = nuvHeader0[hdrKeyword]  # Save the key/value pairs to the dictionary
nuv_param_dict[hdrKeyword] = value
print(f"{hdrKeyword} = {value}")  # Print the key/value pairs

Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HST/product/lby606010_x1dsum.fits to data/mastDownload/HST/lby606010/lby606010_x1dsum.fits ... [Done]
DETECTOR = NUV
OPT_ELEM = G185M
CENWAVE = 1921
DISPTAB = 63p1559jl_disp.fits


In [13]:
#### NUV: Fetch disptab + LSF files
# We'll pass that fetch function the parameters we determined earlier
# This phrasing works because of the order in which we added the params to the dict
nuv_LSF_file_name, nuv_disptab_path = fetch_files(*list(nuv_param_dict.values()))

Downloaded LSF file to data/nuv_model_lsf.dat


Read in the NUV LSF kernels and plot them:

In [14]:
#### Read in the NUV LSF kernels:

In [15]:
# Set up the figure
plt.figure(figsize=(10, 8))

# Fill the figure with a plot of the data
plt.plot(pix_nuv, lsf_nuv["1800"])

# Give fig the title and labels
plt.title("Figure 2.3\nCOS NUV LSF kernel sampled at 1800 Å", size=20)
plt.xlabel("Pixel", size=20)
plt.ylabel("Flux [$normalized,unitless$]", size=20)

# format and save the figure
plt.tight_layout()
plt.savefig(str(plotsdir / "oneKernel_nuv.png"), bbox_inches="tight")

In [16]:
### Get the approximate NUV dispersion to convert pixel --> wavelength
approx_disp = nuv_dtab[
np.where(
(nuv_dtab["CENWAVE"] == 1921)
& (nuv_dtab["SEGMENT"] == "NUVB")
& (nuv_dtab["APERTURE"] == "PSA")
)
]["COEFF"][0][1]

### Now plot
fig, (ax0) = plt.subplots(1, 1, figsize=(10, 5), dpi=100)

# Loop through the lsf kernels
for i, col in enumerate(lsf_nuv.colnames):
line_wvln = int(col)  # central position
contents = lsf_nuv[col].data  # actual data on line shape

xrange = approx_disp * pix_nuv + line_wvln  # Roughly convert pix to wvln
ax0.plot(xrange, contents, linewidth=0.5)  # Plot the kernel
ax0.axvline(line_wvln, linewidth=0.1)  # plot the LSF wvln as a faded line

# Add labels, titles and text
fig.suptitle(f"Fig 2.4\nAll the NUV LSF kernels in {nuv_LSF_file_name}", fontsize=20)
ax0.set_xlabel("Wavelength [$\AA$]", size=14)
ax0.set_ylabel("Flux [$normalized,unitless$]", size=14)

# format and save the figure
plt.tight_layout()
plt.savefig(str(plotsdir / "allKernels_nuv.png"), bbox_inches="tight")


# 3. Convolving an LSF¶

Now we come to the central goal of the Notebook, convolving a spectrum with the COS LSF. We'll begin by defining some functions we'll use for the convolution.

## 3.1. Defining some functions for LSF convolution¶

We've already defined a function to read in the LSF file, but we'll need the main function to take a spectrum and convolve it with the LSF kernel. We'll call this function:convolve_lsf. This function will, in turn, call two short functions: get_disp_params to and redefine_lsf:

### 3.1.1. Getting the dispersion relationship parameters¶

First, we'll define a function, get_disp_params, to interpret the DISPTAB to find the dispersion relationship. If provided with a range of pixel values, it will apply the dispersion relationship to those values and return the equivalent wavelength values as well as the dispersion coefficients.

In [17]:
def get_disp_params(disptab, cenwave, segment, x=[]):
"""
Helper function to redefine_lsf(). Reads through a DISPTAB file and gives relevant\
dispersion relationship/wavelength solution over input pixels.
Parameters:
disptab (str): Path to your DISPTAB file.
cenwave (str): Cenwave for calculation of dispersion relationship.
segment (str): FUVA or FUVB?
x (list): Range in pixels over which to calculate wvln with dispersion relationship (optional).
Returns:
disp_coeff (list): Coefficients of the relevant polynomial dispersion relationship
wavelength (list; if applicable): Wavelengths corresponding to input x pixels
"""
with fits.open(disptab) as d:
wh_disp = np.where(
(d[1].data["cenwave"] == cenwave)
& (d[1].data["segment"] == segment)
& (d[1].data["aperture"] == "PSA")
)[0]
disp_coeff = d[1].data[wh_disp]["COEFF"][0] # 0 is needed as this returns nested list [[arr]]
d_tv03 = d[1].data[wh_disp]["D_TV03"]  # Offset from WCA to PSA in Thermal Vac. 2003 data
d_orbit = d[1].data[wh_disp]["D"]  # Current offset from WCA to PSA

delta_d = d_tv03 - d_orbit

if len(x):  # If given a pixel range, build up a polynomial wvln solution pix -> λ
wavelength = np.polyval(p=disp_coeff[::-1], x=np.arange(16384))
return disp_coeff, wavelength
else:  # If x is empty:
return disp_coeff


### 3.1.2. Redefining the LSF in terms of wavelength¶

Now, we'll define a function, redefine_lsf, to apply the dispersion relationship to the LSF kernels. This effectively converts an LSF kernel from a function of pixel to a function of wavelength, making it compatible with a spectrum.

In [18]:
def redefine_lsf(lsf_file, cenwave, disptab, detector="FUV"):
"""
Helper function to convolve_lsf(). Converts the LSF kernels in the LSF file from a fn(pixel) -> fn(λ)\
which can then be used by convolve_lsf() and re-bins the kernels.
Parameters:
lsf_file (str): path to your LSF file
cenwave (str): Cenwave for calculation of dispersion relationship
disptab (str): path to your DISPTAB file
detector (str): FUV or NUV?
Returns:
new_lsf (numpy.ndarray): Remapped LSF kernels.
new_w (numpy.ndarray): New LSF kernel's LSF wavelengths.
step (float): first order coefficient of the FUVA dispersion relationship; proxy for Δλ/Δpixel.
"""

if detector == "FUV":
xfull = np.arange(16384)

# Read in the dispersion relationship here for the segments
### FUVA is simple
disp_coeff_a, wavelength_a = get_disp_params(disptab, cenwave, "FUVA", x=xfull)
### FUVB isn't taken for cenwave 1105, nor 800:
if (cenwave != 1105) & (cenwave != 800):
disp_coeff_b, wavelength_b = get_disp_params(
disptab, cenwave, "FUVB", x=xfull)
elif cenwave == 1105:
# 1105 doesn't have an FUVB so set it to something arbitrary and clearly not real:
wavelength_b = [-99.0, 0.0]

# Get the step size info from the FUVA 1st order dispersion coefficient
step = disp_coeff_a[1]

# Read in the lsf file

# take median spacing between original LSF kernels
deltaw = np.median(np.diff(w))

lsf_array = [np.array(lsf[key]) for key in lsf.keys()]
if (deltaw < len(pix) * step * 2):  # resamples if the spacing of the original LSF wvlns is too narrow
# this is all a set up of the bins we want to use
# The wvln difference between kernels of the new LSF should be about twice their width
new_deltaw = round(len(pix) * step * 2.0)
new_nw = (int(round((max(w) - min(w)) / new_deltaw)) + 1)  # nw = number of LSF wavelengths
new_w = min(w) + np.arange(new_nw) * new_deltaw  # new version of lsf_wvlns

# populating the lsf with the proper bins
new_lsf = np.zeros((len(pix), new_nw))  # empty 2-D array to populate
for i, current_w in enumerate(new_w):
dist = abs(current_w - w)  # Find closest original LSF wavelength to new LSF wavelength
lsf_index = np.argmin(dist)
orig_lsf_wvln_key = lsf.keys()[lsf_index]  # column name corresponding to closest orig LSF wvln
new_lsf[:, i] = np.array(lsf[orig_lsf_wvln_key])  # assign new LSF wvln the kernel of the closest original lsf wvln
else:
new_lsf = lsf
new_w = w
return new_lsf, new_w, step

elif detector == "NUV":
xfull = np.arange(1024)
# Read in the dispersion relationship here for the segments
disp_coeff_a, wavelength_a = get_disp_params(disptab, cenwave, "NUVA", x=xfull)
disp_coeff_b, wavelength_b = get_disp_params(disptab, cenwave, "NUVB", x=xfull)
disp_coeff_c, wavelength_c = get_disp_params(disptab, cenwave, "NUVC", x=xfull)

# Get the step size info from the NUVB 1st order dispersion coefficient
step = disp_coeff_b[1]

# Read in the lsf file

# take median spacing between original LSF kernels
deltaw = np.median(np.diff(w))

lsf_array = [np.array(lsf[key]) for key in lsf.keys()]

# this section is a set up of the new bins we want to use:
new_deltaw = round(len(pix) * step * 2.0)  # The wvln difference between kernels of the new LSF should be about twice their width
new_nw = (int(round((max(w) - min(w)) / new_deltaw)) + 1)  # nw = number of LSF wavelengths
new_w = min(w) + np.arange(new_nw) * new_deltaw  # new version of lsf_wvlns

# populating the lsf with the proper bins
new_lsf = np.zeros((len(pix), new_nw))  # empty 2-D array to populate
for i, current_w in enumerate(new_w):
dist = abs(current_w - w)  # Find closest original LSF wavelength to new LSF wavelength
lsf_index = np.argmin(dist)
orig_lsf_wvln_key = lsf.keys()[lsf_index]  # column name corresponding to closest orig LSF wvln
new_lsf[:, i] = np.array(lsf[orig_lsf_wvln_key])  # assign new LSF wvln the kernel of the closest original lsf wvln
return new_lsf, new_w, step


Let's make a version of the Fig. 2.2 with these newly remapped LSF kernels:

In [19]:
# Generate the redefined lsf for the plot
new_lsf, new_w, step = redefine_lsf(
param_dict["CENWAVE"],
)

fig, (ax0, ax1) = plt.subplots(2, 1, figsize=(10, 8), dpi=100)

# Loop through the new remapped lsf kernels
for i, col in enumerate(new_w):
line_wvln = int(col)  # central position
contents = new_lsf[:, i]  # actual shape
if line_wvln < 1291:  # split into the FUVB segment
xrange = new_w[i] + pix * step  # ROUGHLY convert pix to wvln
ax1.plot(xrange, contents, linewidth=0.5)  # Plot the kernel
ax1.axvline(line_wvln, linewidth=0.1)  # plot the LSF wvln as a faded line
elif line_wvln > 1291:  # split into the FUVA segment
xrange = new_w[i] + pix * step  # ROUGHLY convert pix to wvln
ax0.plot(xrange, contents, linewidth=0.5)  # Plot the kernel
ax0.axvline(line_wvln, linewidth=0.1)  # plot the LSF wvln as a faded line

# Add labels, titles and text
fig.suptitle(f"Fig 3.1\nAll the $remapped$ LSF kernels in {LSF_file_name}\n")
ax1.set_title("FUVB Detector")
ax0.set_title("FUVA Detector")
ax0.set_xlim(1290, 1435)
ax1.set_xlim(1125, 1295)

ax1.set_xlabel("Wavelength [$\AA$]", size=14)
fig.text(
s="Flux [$normalized,unitless$]", x=-0.01, y=0.35, rotation="vertical", size=14
)

# format and save the figure
plt.tight_layout()
plt.savefig(str(plotsdir / "allKernels_new.png"), bbox_inches="tight")


### 3.1.3. Applying the convolution¶

Finally, we'll define the main function, convolve_lsf, to convolve the template spectrum with the redefined LSF:

In [20]:
def convolve_lsf(wavelength, spec, cenwave, lsf_file, disptab, detector="FUV"):
"""
Main function; Convolves an input spectrum - i.e. template or STIS spectrum - with the COS LSF.
Parameters:
wavelength (list or array): Wavelengths of the spectrum to convolve.
spec (list or array): Fluxes or intensities of the spectrum to convolve.
cenwave (str): Cenwave for calculation of dispersion relationship
lsf_file (str): Path to your LSF file
disptab (str): Path to your DISPTAB file
detector (str) : Assumes an FUV detector, but you may specify 'NUV'.
Returns:
wave_cos (numpy.ndarray): Wavelengths of convolved spectrum.!Different length from input wvln
final_spec (numpy.ndarray): New LSF kernel's LSF wavelengths.!Different length from input spec
"""
# First calls redefine to get right format of LSF kernels
new_lsf, new_w, step = redefine_lsf(lsf_file, cenwave, disptab, detector=detector)

# sets up new wavelength scale used in the convolution
nstep = round((max(wavelength) - min(wavelength)) / step) - 1
wave_cos = min(wavelength) + np.arange(nstep) * step

# resampling onto the input spectrum's wavelength scale
interp_func = interp1d(wavelength, spec)  # builds up interpolated function from input spectrum
spec_cos = interp_func(wave_cos)  # builds interpolated initial spectrum at COS' wavelength scale for convolution
final_spec = interp_func(wave_cos)  # Initializes final spectrum to the interpolated input spectrum

for i, w in enumerate(new_w):  # Loop through the redefined LSF kernels
# First need to find the boundaries of each kernel's "jurisdiction": where it applies
# The first and last elements need to be treated separately
if i == 0:  # First kernel
diff_wave_left = 500
diff_wave_right = (new_w[i + 1] - w) / 2.0
elif i == len(new_w) - 1:  # Last kernel
diff_wave_right = 500
diff_wave_left = (w - new_w[i - 1]) / 2.0
else:  # All other kernels
diff_wave_left = (w - new_w[i - 1]) / 2.0
diff_wave_right = (new_w[i + 1] - w) / 2.0

# splitting up the spectrum into slices around the redefined LSF kernel wvlns
# will apply the kernel corresponding to that chunk to that region of the spectrum - its "jurisdiction"
chunk = np.where(
(wave_cos < w + diff_wave_right) & (wave_cos >= w - diff_wave_left)
)[0]
if len(chunk) == 0:
# off the edge, go to the next chunk
continue

current_lsf = new_lsf[:, i]  # selects the current kernel

if len(chunk) >= len(
current_lsf
):  # Makes sure that the kernel is smaller than the chunk
final_spec[chunk] = convolve(
spec_cos[chunk],
current_lsf,  # Applies the actual convolution
boundary="extend",
normalize_kernel=True,
)

return wave_cos, final_spec  # Remember, not the same length as input spectrum data!


## 3.2. Convolving simple line profiles¶

Let's demonstrate the convolution with a quick initial plot.

The cell below first creates a simple synthetic spectrum (with a wavelength range from 1147-1153Å and a Voigt profile emission line at 1150Å), then convolves it with the COS LSF. The next cell then plots the two spectra together.

In [21]:
##### Generate data:
# Define a model spectral line with a Voigt profile
voigt_shape = functional_models.Voigt1D(
x_0=1150, amplitude_L=1, fwhm_G=0.05, fwhm_L=0.05
)
# Make a simple spectrum with just that line at 1150 Å

# generate x data - Minimum size of Δ6Å for the kernel to apply here.
wvln_in = np.linspace(1147, 1153, int(1e5))

spec_in = voigt_shape(wvln_in)  # Generate the y data
spec_in /= max(spec_in)  # Normalize the y data to a max of 1

##### Run the convolution
wvln_out, spec_out = convolve_lsf(
wvln_in,
spec_in,
1291,
)

In [22]:
##### Make a plot from the data generated above:
plt.figure(figsize=(10, 8))  # Set up figure
# Plot the two spectra
plt.plot(wvln_in, spec_in, label="Un-Convolved Voigt line spectrum")
plt.plot(
wvln_out,spec_out,
linestyle="--",linewidth=2,
c="k",label="Convolved Voigt line spectrum",
)

# Format and give fig the title and labels
plt.title("Figure 3.2\nConvolution applied to a single synthetic line in the FUV", size=20)
plt.xlabel("Wavelength [$\AA$]", size=20)
plt.ylabel("Flux [$normalized,unitless$]", size=20)
plt.legend(fontsize=12, loc="upper right")

# Save the figure
plt.tight_layout()
plt.savefig(str(plotsdir / "applyConv1.png"), bbox_inches="tight")


As shown above, the LSF pushes power from the line's peak into the wings of the line profile.

And, as we confirm below, the total flux has not changed, but has just been spread out:

In [23]:
integral1, integral2 = np.trapz(x=wvln_in, y=spec_in), np.trapz(x=wvln_out, y=spec_out)
print(f"The integrated fluxes are within {100* (integral1-integral2)/integral2:.2f} % of eachother")

The integrated fluxes are within 0.06 % of eachother


We can now examine how the convolution differs across the spectrum. We do this by creating a "picket fence"/"frequency comb"-like spectrum of evenly-spaced, identical, normalized Voigt profiles. Each of these lines is normalized to a maximum of 1.

We'll make our synthetic spectrum slightly more realistic by converting from emission lines to absorption lines on a flat spectrum. To do this, we define a simple function, emit2abs(), which defines a flat continuum of 1, then multiplies this continuum by $e^{-emission\ spectrum}$

We then use convolve_lsf to convolve each of these synthetic spectral lines with the LSF kernel under whose "jurisdiction" they fall. For each line, we plot pre- and post- convolution spectrum in the 8 smaller panels. In the bottom, larger panel, we also plot the entire spectrum of all the lines.

In [24]:
def emit2abs(emitspec, continuum=-1):
"""
Takes an emission spectrum and "reverses" it, that is it:
Makes a default flat continuum of 1 and subtracts exp(line)
Inputs:
emitspec (1D arraylike): the emission spectrum
continuum (array or -1) : if -1, default of continuum of 1, \
otherwise must be same length as emitspec
"""
if type(continuum) == int:
if continuum == -1:
continuum = np.ones(len(emitspec))
abs_spec = continuum * np.exp(-1.0 * emitspec)
return abs_spec

In [25]:
########## Set up the DATA ##############

# Build up a synthetic wavelength range comparable to a real COS spectrum
wvln = np.linspace(1120, 1442, 16384 * 2)

# Initialize the spectrum to continuum of zero, we'll build it up with a series of Voigt profiles
flux = np.zeros(wvln.shape)
combo_flux = flux  # A copy of flux that we'll add each new line's flux to, in turn

# Set up figure
fig = plt.figure(figsize=(20, 15))
gs = fig.add_gridspec(nrows=10, ncols=4)  # Using gridspec to control panel sizes and locations

# Loop through the input spectrum's wavelength range and place a synthetic Voigt line every 40 Å
for i, discrete_wvln in enumerate(np.arange(int(min(wvln) + 5), max(wvln) - 5, 40)):
voigt_shape = functional_models.Voigt1D(
x_0=discrete_wvln,  # Center a Voigt profile there
amplitude_L=1,
fwhm_G=0.05,
fwhm_L=0.05)

flux = voigt_shape(wvln)  # Evaluate flux from that Voigt profile function
flux /= max(flux)  # Normalize that line's flux
combo_flux = combo_flux + flux  # Add each line's flux to total summed flux

# "Reverse" emission spectrum to absorption spectrum
combo_flux = emit2abs(combo_flux)

# Apply the convolution to the combined (many-line) synthetic spectrum to create an lsf_convolved wvln and flux
lwvln, lsf_combo_flux = convolve_lsf(
wvln,
combo_flux,
1291,
)

########## Make the plots ##############

# Loop through again to build up the plots
for i, discrete_wvln in enumerate(
np.arange(int(min(wvln) + 5), max(wvln) - 5, 40)
):  # This will make 8 subplots
###### Build the small subplots for each line
gs[4 * int(i / 4) : 4 * int(i / 4) + 4, i % 4 : (i) % 4 + 1]
)  # Add a plot at the correct position on the grid

# First plot the original, unconvolved line at each position
ax.plot(
wvln, combo_flux, label="Un-convolved"
)
# Now plot the convolved line
ax.plot(
lwvln, lsf_combo_flux, c="k", linestyle="--", linewidth=2, label="LSF-convolved"
)
# Now add the peak wvln
ax.axvline(
discrete_wvln, c="C1", linestyle="dotted", alpha=0.8, label="Peak of Voigt"
)
# Some formatting
ax.set_xlim(discrete_wvln - 0.5, discrete_wvln + 0.5)
ax.ticklabel_format(axis="x", style="plain", useOffset=True, useMathText=True)
if i == 0:  # Add a legend to the first subplot
ax.legend(fontsize=12, loc="upper left")

# Build up lower plot of all the lines:
low_ax.plot(wvln, combo_flux, label="Un-convolved")
low_ax.plot(lwvln, lsf_combo_flux, c="k", linestyle="--", label="LSF-convolved")
low_ax.legend(fontsize=14, loc="upper right")

fig.suptitle(
'Fig 3.3\nComparison of convolutions across a "picket fence" of absorption lines in the FUV\n',
size=25,
)

fig.text(0.5, -0.01, "Wavelength [$\AA$]", ha="center", fontsize=20)
fig.text(
-0.01,0.5,
"Flux [$normalized,\ unitless$]",
va="center",rotation="vertical",
fontsize=20,
)

plt.tight_layout()
plt.savefig(str(plotsdir / "applyConv_picketFence.png"), bbox_inches="tight", dpi=200)


The figure above demonstrates the way in which the LSF changes significantly throughout the wavelength range of a spectrum.

It can be difficult to distinguish the LSFs plotted on their own; however, the results of their convolution with a synthetic line shows some are more sharply peaked, or contain more flux in the wings, etc.

We'll finish this section with a brief example of convolving a line in the NUV with the relevant LSF:

In [26]:
##### Generate data:
# Define a model spectral line with a Voigt profile
voigt_shape = functional_models.Voigt1D(
x_0=1803, amplitude_L=1, fwhm_G=0.05, fwhm_L=0.05
)
# Make a simple spectrum with just that line at 1150 Å

# generate x data - Minimum size of Δ6Å for the kernel to apply here.
wvln_in = np.linspace(1800, 1806, int(1e4))

spec_in = voigt_shape(wvln_in)  # Generate the y data
spec_in /= max(spec_in)  # Normalize the y data to a max of 1
### Run the convolution
wvln_out, spec_out = convolve_lsf(
wavelength=wvln_in,
spec=spec_in,
cenwave=1786,
disptab=nuv_disptab_path,
detector="NUV",
)
############################################################################################
##### Make a plot from that data:
plt.figure(figsize=(10, 8))  # Set up figure
# Plot the two spectra
plt.plot(wvln_in, spec_in, label="Un-Convolved Voigt line spectrum")
plt.plot(wvln_out,spec_out,
linestyle="--",linewidth=2,
c="k",label="Convolved Voigt line spectrum",
)
plt.legend(fontsize=12, loc="upper right")
# Give fig the title and labels
plt.title("Figure 3.4\nConvolution applied to a single synthetic line in the NUV", size=20)
plt.xlabel("Wavelength [$\AA$]", size=20)
plt.ylabel("Flux [$normalized,unitless$]", size=20)
# format and save the figure
plt.tight_layout()
plt.savefig(str(plotsdir / "applyConv_nuv1.png"), bbox_inches="tight")

# Give the user a heads-up that the integrated fluxes should agree:
integral1, integral2 = np.trapz(x=wvln_in, y=spec_in), np.trapz(x=wvln_out, y=spec_out)
print(f"The integrated fluxes are within {100* (integral1-integral2)/integral2:.2f} % of eachother")

The integrated fluxes are within 0.12 % of eachother