gollum’s documentation

The goal of this repo is to provide a Python Application Programming Interface (API) to several different synthetic spectral models. gollum will be built on the astropy affiliated package specutils , and will be inspired by the API design of lightkurve . This project is loosely related to the parallel muler framework that is built on specutils and focuses on data. This project is all about models. The code itself and will have some overlap with functionality in Starfish , and this project could one day become a microservice to Starfish, rather than duplicate code.

gollum makes it easy to chain routine operations in sequence:

Installing the development version

Note

Conda installation is not yet available.

Pip installation is available

pip install gollum

Currently we recommend installing the developer version from source to help us with beta testing. You can install gollum into its own isolated conda python environment using the instructions below. Line 3 will create a conda environment titled gollum_dev that contains all the dependencies required to install gollum, and even all the affiliated code to make the documentation and run unit tests from scratch. Very useful! You can activate this conda environment using the code on line 4 conda activate gollum_dev. You could alternatively activate some other existing conda environment on your computer, but it is not guaranteed to have all of the extra packages. Our sibling package muler offers an identical conda environment to this one, so if you made one for muler already, you can do: conda activate muler_dev

git clone https://github.com/BrownDwarf/gollum.git
cd gollum
conda env create -f environment.yml
conda activate gollum_dev # or conda activate muler_dev if you have that already
python setup.py develop

Eventually you can run the tests in the tests/ directory to double-check that everything installed correctly. Currently we are evaluating the best way to specify paths to voluminous libraries to make testing robust and quick across machines. Stay tuned!

py.test -vs test_precomputed.py  # should work for everyone
py.test -vs test_phoenix.py # requires local PHOENIX models
py.test -vs test_sonora.py # requires local Sonora models

Requirements

The project appears to work on modern Mac OS, Linux, and Windows operating systems, and has been tested for Python 3.7 and above. It may work on Python 3.6, and will not work on Python 2.

Downloading model grids

Sonora

Currently we only support the 2018 Sonora Model grids. In a web browser, navigate to the Sonora 2018 Zenodo website. Click on the Download button next to spectra.tar     1.2 GB. That will download the spectra.tar file to your ~/Downloads/ directory or equivalent default location. Untar that directory. The pattern should look something like below, depending on what your operating system defaults are, so copy-ing and past-ing this code is not likely to work— try each step at a time.

tar -xzvf spectra.tar
mv spectra ~/libraries/raw/Sonora

The end result should look something like this:

$ ls ~/libraries/raw/Sonora/
            [...]
sp_t2100g3160nc_m0.0.gz  sp_t450g316nc_m0.0.gz    sp_t950g31nc_m0.0.gz
sp_t2100g316nc_m0.0.gz   sp_t450g31nc_m0.0.gz     sp_t950g562nc_m0.0.gz
sp_t2100g31nc_m0.0.gz    sp_t450g562nc_m0.0.gz    sp_t950g56nc_m0.0.gz
            [...]

PHOENIX

Note

Downloading all the PHOENIX models can take hours or days! Start the downloading early.

The PHOENIX models total over 100 GB, and generally download at a relatively slow bandwidth. The files are arranged into sub-directories for metallicity and alpha-element abundance.

You can get all of the PHOENIX models in one-fell-swoop from the command line if you have wget

cd ~/Downloads
wget -r -l 0 ftp://phoenix.astro.physik.uni-goettingen.de/HiResFITS/

If you don’t have wget on your computer, please help the gollum grow by filing a GitHub Issue with how you resolved the problem, or what problems you are encountering.

As noted, this process will take a while as each individual file is painstakingly downloaded from a single German computer. The commandline script, as written, will preserve the directory structure— that’s good! The gollum code demands that the directory structure is preserved. Once it’s all downloaded it should look like this

# From inside this directory: ~/Downloads/phoenix.astro.physik.uni-goettingen.de/
$ tree
.
└── HiResFITS
    ├── PHOENIX-ACES-AGSS-COND-2011
    │   └── Z-0.0
    │       ├── lte02300-0.00-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits
    │       ├── lte02300-0.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits
    │       ├── lte02300-1.00-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits
    │       ├── lte02300-1.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits
                                    [ ... ]
    │       ├── lte04000-2.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits
    │       ├── lte04000-3.00-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits
    │       └── lte04000-3.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits
    │                                [ ... ]
    │    ├── Z-0.0
    │    ├── Z+0.5
    │    ├── Z-0.5
    │    ├── Z+1.0
    │    ├── Z-1.0
    │    ├── Z-1.5
    │    ├── Z-2.0
    │    ├── Z-3.0
    │    └── Z-4.0
    └── WAVE_PHOENIX-ACES-AGSS-COND-2011.fits

While you could leave these directories in say ~/Downloads/phoenix.astro.physik.uni-goettingen.de/, I recommend making a more permanent and recognizable home for these models. In particular gollum attempts to search a single default path for models: ~/libraries/raw/, where the tilde ~/ denotes your home directory

mv ~/Downloads/phoenix.astro.physik.uni-goettingen.de/HiResFITS/PHOENIX-ACES-AGSS-COND-2011 ~/libraries/raw/PHOENIX/
ls ~/libraries/raw/PHOENIX/
Z-0.0  Z+0.5  Z-0.5  Z+1.0  Z-1.0  Z-1.5  Z-2.0  Z-3.0  Z-4.0

Finally, you must copy the wavelength file into this directory as well. Notice that this placement breaks the native directory structure, so you must complete this step in order for gollum to work.

mv ~/Downloads/phoenix.astro.physik.uni-goettingen.de/HiResFITSWAVE_PHOENIX-ACES-AGSS-COND-2011.fits ~/libraries/raw/PHOENIX/
ls ~/libraries/raw/PHOENIX/
WAVE_PHOENIX-ACES-AGSS-COND-2011.fits  Z-0.0  Z+0.5  Z-0.5  Z+1.0  Z-1.0  Z-1.5  Z-2.0  Z-3.0  Z-4.0

gollum Quickstart

%config InlineBackend.figure_format='retina'

from gollum.phoenix import PHOENIXSpectrum

Simply provide the \(T_{\text{eff}}\) and \(\log{g}\) values you desire:

spec = PHOENIXSpectrum(teff=5000, logg=4)

Normalize the spectrum by the median:

normalized_spectrum = spec.normalize()

The spectrum has wavelength, with units:

normalized_spectrum.wavelength

$[8038.01,~8038.02,~8038.03,~\dots,~12848.94,~12848.96,~12848.98] \; \mathrm{\mathring{A}}$

The flux vector is the same length as wavelength:

normalized_spectrum.flux.dtype, normalized_spectrum.shape

(dtype('float64'), (338649,))

ax = normalized_spectrum.plot()
ax.set_ylim(0, 2)

(0.0, 2.0)

Right now we truncate the spectrum to the near-IR by default, you can change that with keyword arguments from the beginning:

spec = PHOENIXSpectrum(teff=5000, logg=4, wl_lo=10_800, wl_hi=10_860)

ax = spec.normalize().plot()
ax.set_ylim(0)

(0.0, 1.2)

Neat! gollum is still under active development. Help us develop the tool by engaging with us on our GitHub Issues page. You can suggest a feature, or help us brainstorm how to build this project more. Thanks!

gollum tutorials

Working with the PHOENIX models

The PHOENIX models are under active development. Check back soon!

from gollum.phoenix import PHOENIXSpectrum
import numpy as np
%config InlineBackend.figure_format='retina'

spec = PHOENIXSpectrum(teff=2800, logg=5.0)

spec.wavelength # "Angstroms"

$[8038.01,~8038.02,~8038.03,~\dots,~12848.94,~12848.96,~12848.98] \; \mathrm{\mathring{A}}$

spec.flux.unit # "egs / s / cm^2 / cm"

$\mathrm{\frac{erg}{s\,cm^{3}}}$

We can normalize the spectrum, which divides the spectrum by the median value, rendering the flux units dimensionless

spectrum = spec.normalize()
np.median(spectrum.flux)

$1 \; \mathrm{}$

ax = spectrum.plot()
ax.set_xlim(8160, 8220)
ax.set_ylim(0, 1.5)

(0.0, 1.5)

Neat! There are many great uses for the PHOENIX models. We will be adding support for other model grids in the future.

Finding the best fit rotational broadening and radial velocity

Often in astronomy we have a data spectrum and we want to answer the question:

What is the best fit :math:`vsin{i}` and RV for this data?

In this demo we will show some simple ways to use gollum to find the model with the best fit rotational broadening and radial velocity, assuming a fixed template.

from gollum.phoenix import PHOENIXSpectrum, PHOENIXGrid
import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt
from tqdm.notebook import tqdm
%config InlineBackend.figure_format='retina'

For this demo we will need some real world example data

Let’s use data of an A0V star from HPF and our sibling package muler.

You can get free example data from the muler example data GitHub repository.

from muler.hpf import HPFSpectrumList, HPFSpectrum

local_file = '../../../muler_example_data/HPF/01_A0V_standards/Goldilocks_20210517T054403_v1.0_0060.spectra.fits'
raw_data = HPFSpectrumList.read(local_file)

def data_clean(data):
    """Clean the HPF data with standard post-processing techniques"""
    data = data.sky_subtract(method='vector')
    data = data.deblaze()
    data = data.trim_edges((4, 2042))
    data = data.normalize()
    data = data.stitch()
    return data

full_data = data_clean(raw_data)

As a final step, we will mask the telluric absorption lines. This step can sometimes benefit from hand-tuning.

data = full_data.mask_tellurics(threshold=0.999, dilation=13)

We will restrict our fits to the region with the highest density of H lines.

data = data[8500*u.Angstrom:8950*u.Angstrom]
data = data.normalize()

ax = data.plot(marker='.', linestyle='None');
ax.set_xlim(8500, 8950);

OK, that’s our data spectrum against which we will compare models. You can see large voids in the spectrum due to our telluric masking– that’s fine, the data need not be contiguous or evenly sampled to estimate a best fit model. We will resample the model to the data.

We can choose 3 dimensions in our grid: \(T_{\mathrm{eff}}, \log{g}, \left[\frac{\mathrm{Fe}}{\mathrm{H}}\right]\)

template = PHOENIXSpectrum(teff=9_600, logg=4.5, metallicity=0)

We then want to rotationally broaden and RV shift the spectrum. We’ll try some guesses first:

guess_vsini = 150
guess_RV = -50

spec = template.rotationally_broaden(guess_vsini)\
            .rv_shift(guess_RV)\
            .instrumental_broaden(resolving_power=55_000)\
            .resample(data)\
            .normalize()

ax = data.plot(label='Data')
spec.plot(ax=ax, label='Doctored Model')
ax.legend();

You can see that our guess RV is off. Let’s do a grid search for \(v\sin{i}\) and RV.

n_vsinis = 20
n_RVs = 20
vsinis = np.linspace(1, 150.0, n_vsinis)
RVs = np.linspace(-100.0, 100.0, n_RVs)

loss_values = np.zeros((n_vsinis, n_RVs))

We will compute the sum-of-the-squares-of-the-residuals (“chi-squared”) for each value of \(v\sin{i}\) and RV.

for i, vsini in tqdm(enumerate(vsinis), total=n_vsinis):
    for j, RV in enumerate(RVs):
            spec = template.rotationally_broaden(vsini)\
                .rv_shift(RV)\
                .instrumental_broaden(resolving_power=55_000)\
                .resample(data)\
                .normalize()
            residual = data.subtract(spec, handle_meta='ff')
            chi2_loss = np.sum(0.5*residual.flux.value**2/residual.uncertainty.array**2)

            loss_values[i,j] = chi2_loss

The best fit values minimizes the chi-squared.

best_i, best_j = np.unravel_index(np.argmin(loss_values), loss_values.shape)

best_vsini, best_RV = vsinis[best_i], RVs[best_j]
best_RV, best_vsini

(15.789473684210535, 102.94736842105263)

extent=[RVs.min(), RVs.max(), vsinis.min(), vsinis.max()]
plt.imshow(loss_values, extent=extent, aspect=1, origin='lower')
plt.scatter(best_RV,best_vsini,  marker='*', c='w', ec='k', s=200)
plt.xlabel('RV ($\mathrm{km}\;\mathrm{s}^{-1}$)')
plt.ylabel('$v\sin{i}$ ($\mathrm{km}\;\mathrm{s}^{-1}$)');

Awesome, we have found the RV and \(v\sin{i}\) with the closest match to the data.

best_spec_full = template.rotationally_broaden(best_vsini)\
            .rv_shift(best_RV)\
            .instrumental_broaden(resolving_power=55_000)
best_spec = best_spec_full[data.wavelength.min():data.wavelength.max()].normalize()

How does the best fit look by-eye?

ax=(data/data.flux.max()).plot(marker='.', linestyle='None')
(best_spec/best_spec.flux.max()).plot(ax=ax)
title1 = '$v\sin{i} = $'+'{:0.0f}'.format(best_vsini)+' km/s,   '
title2 = 'RV = '+'{:0.1f}'.format(best_RV)+' km/s'
title = title1 + title2
ax.set_title(title);

Excellent! That’s much better than our initial guess.

It’s still not perfect because the model template has imperfections.

Divide by a black body

Often in real-world astronomical applications, we want to see how a stellar spectrum varies in wavelength without the large and omnipresent crest of the black body curve.

In this tutorial we show how to remove the black body curve from a precomputed model spectrum.

from gollum.phoenix import PHOENIXSpectrum
from gollum.precomputed_spectrum import PrecomputedSpectrum
from astropy.modeling.physical_models import BlackBody
import astropy.units as u
import numpy as np
%config InlineBackend.figure_format='retina'

T_eff = 5_700

original_spectrum = PHOENIXSpectrum(teff=T_eff, logg=4.5, metallicity=0,
                                    wl_lo=0, wl_hi=1e10) # Get the whole spectrum

The PHOENIX spectra have units of \(\frac{\mathrm{erg}}{\mathrm{s\;cm^2 \;cm}}\)

original_spectrum.flux.unit

$\mathrm{\frac{erg}{s\,cm^{3}}}$

original_spectrum.flux.min()

$5.5922428 \times 10^{-15} \; \mathrm{\frac{erg}{s\,cm^{3}}}$

Make a black body spectrum with the same temperature

blackbody_model = BlackBody(temperature=T_eff*u.Kelvin)
blackbody_flux_per_sr = blackbody_model(original_spectrum.wavelength)

The Black Body spectra have units of \(\frac{\mathrm{erg}}{\mathrm{s\;cm^2 \;Hz \;sr}}\). To convert between the two conventions we have to do two things. First we have to pretend we are standing on the surface of the star, and multiply by \(\pi\) steradians. Second, we have to convert flux density in per frequency to per wavelength units. We can do that with astropy’s equivalencies kwarg.

blackbody_flux_per_Hz = blackbody_flux_per_sr * np.pi * u.steradian

blackbody_flux = blackbody_flux_per_Hz.to(original_spectrum.flux.unit,
                  equivalencies=u.spectral_density(original_spectrum.wavelength))

blackbody = PrecomputedSpectrum(flux=blackbody_flux,
                                spectral_axis=original_spectrum.wavelength)

ax = original_spectrum.plot()
blackbody.plot(ax=ax)
ax.set_ylim(1e0, 1e16)
ax.set_yscale('log')
ax.set_ylabel('Flux (erg/s/cm$^3$)');

The plot spans 16 orders of magnitude– a huge dynamic range! Notice that the spectra have similar—but not identical—broadband spectral shapes. They should have the identical area under the curve, by the definition of effective temperature. Let’s see if they do!

from scipy.integrate import trapezoid

original_flux = trapezoid(original_spectrum.flux, x=original_spectrum.wavelength.to(u.cm))
black_body_flux = trapezoid(blackbody.flux, x=original_spectrum.wavelength.to(u.cm))

original_flux/black_body_flux

$0.99048532 \; \mathrm{}$

The two spectral models have the same flux to within 1%, which is close enough to identical given that more spectrum resides outside the extent of the plot, and numerical artifacts limitations in the spectral modeling procedure. Let’s compute the ratio spectrum to see how flat the spectrum looks. We’ll first plot it over the same dynamic range as before to emphasize how much more compressed it is.

ratio_spectrum = original_spectrum.divide(blackbody)

The resulting spectrum is a ratio of fluxes with the same units, so it is dimensionless.

ratio_spectrum.flux.unit == u.dimensionless_unscaled

True

ax = ratio_spectrum.plot();

ax.set_ylim(1e-12, 1e4)
ax.set_yscale('log')
ax.set_ylabel('Normalized flux');

ax.axhline(1.0, linestyle='dashed', color='k');

At this dramatic zoom level, the flux looks pretty flat (except for the extreme ultraviolet portion of the spectrum). Let’s zoom in on a region of interest from \(8000 - 13000\; Å\).

ax = ratio_spectrum.plot(ylo=0, yhi=1.15);
ax.set_ylabel('Normalized flux');

ax.axhline(1.0, linestyle='dashed', color='k');
ax.set_xlim(8_000, 13_000);

OK, looks good! We have successfully used the blackbody curve to coarsely flatten the spectrum!

Interactively fit Brown Dwarf Spectra with the gollum dashboard

In this tutorial we will see how the spectra of brown dwarfs vary as a function of their intrinsic properties. We will fit observed spectra of a particular brown dwarf with the gollum dashboard, a dashboard which fits models based on properties including effective temperature, surface gravity, metallicity, rotational broadening, and radial velocity. The fitting for this tutorial will be based on the Sonora-Bobcat 2021 models, which takes into account effective temperature, surface gravity, and metallicity as intrinsic values.

from gollum.sonora import SonoraGrid
from specutils import Spectrum1D
import pandas as pd
import astropy.units as u

from IPython.display import HTML
from IPython.display import Image

First, we will read in an example spectrum of this ultracool dwarf:

2MASS J05591914-1404488

We got its data from the Keck Telescope’s NIRSPEC spectrograph. A specific section of this data is displayed below.

df = pd.read_csv('../../data/2mass0559_59.dat',
                 delim_whitespace=True,
                 comment='#',
                 names=['wave', 'flux'])

df.head()

The unit for wavelength here is microns and the unit for flux is “counts”.

bdss_spectrum = Spectrum1D(spectral_axis=df.wave.values*u.micron,
                           flux=df.flux.values*u.ct)

wl_lo, wl_hi = (bdss_spectrum.wavelength.value.min(),
                bdss_spectrum.wavelength.value.max())

Next, we can read in the Sonora-Bobcat grid and show an interactive dashboard.

grid = SonoraGrid(wl_lo=wl_lo, wl_hi=wl_hi)

Awesome! Now you can hand-in a data spectrum to overlay it onto the grid and begin fitting using the interactive sliders.

grid.show_dashboard(data=bdss_spectrum, show_telluric=False)

The dashboard looks great!

Demo of the PHOENIX Interactive Dashboard for Stellar Spectra

In this tutorial we will see how the spectra of stars vary as a function of their intrinsic properties.

from gollum.phoenix import PHOENIXGrid
from specutils import Spectrum1D
import pandas as pd
import astropy.units as u
import numpy as np
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))

Fetch example IGRINS data

from muler.igrins import IGRINSSpectrum

path = 'https://github.com/OttoStruve/muler_example_data/raw/main/IGRINS/01_IGRINS_test_data/'
filename='SDCH_20201202_0059.spec_a0v.fits'
full_path = path + filename

spec = IGRINSSpectrum(file=full_path, order=12).normalize().remove_nans().trim_edges()
spec.plot(color=None, ylo=0.7, yhi=1.1);

Telling gollum where to find your local files

The online installation guide shows how to download and store the PHOENIX models. In order to run the code below, you will need to update the my_path variable to the correct path on your local machine.

In this demo, the PHOENIX/ folder is houses all of the voluminous models:

grid = PHOENIXGrid(teff_range=(2500, 7000), logg_range=(2, 5), metallicity_range=(0, 0.5), wl_lo= 16200, wl_hi= 16400)

grid.show_dashboard(data=spec)

API

The API is under active development. Feedback is welcomed.