Retrieving Model Spectra¶

Often, we want to know the spectrum of light emanating from the surface of a star, a planet, or some other object. For chromatic, we added a wrapper around the Husser et al. (2013) library of high-resolution stellar spectra to make it easier to retrieve model spectra at a variety of resolutions. The tools used to achieve this is the get_phoenix_photons function, as demonstrated below.

Let's imagine you want to plot the optical and near-infrared spectrum of a star with a particular effective temperature, surface gravity, and metallicity. The following commands show how to use chromatic to retrieve the model flux of a star like the Sun, making use of the get_phoenix_photons function.

In [1]:

from chromatic import get_phoenix_photons, get_planck_photons, version
import matplotlib.pyplot as plt, numpy as np
import astropy.units as u, astropy.constants as con

from chromatic import get_phoenix_photons, get_planck_photons, version import matplotlib.pyplot as plt, numpy as np import astropy.units as u, astropy.constants as con

In [2]:

version()

version()

Out[2]:

'0.4.12'

In [3]:

# retrieve a wavelength array and a flux array
wavelength, surface_flux = get_phoenix_photons(
    temperature=5780, logg=4.43, metallicity=0.0, R=1000
)

# plot the wavelength array
plt.figure(figsize=(8, 3), dpi=300)
plt.plot(wavelength, surface_flux, label="PHOENIX model")

# for comparison, plot a Planck spectrum
w_planck, f_planck = get_planck_photons(temperature=5780)
plt.plot(w_planck, f_planck, label="Planck function")

# add some labels
plt.xlabel(f"Wavelength ({wavelength.unit.to_string('latex_inline')})")
plt.ylabel(f"Surface Flux ({surface_flux.unit.to_string('latex_inline')})")
plt.legend(frameon=False);

# retrieve a wavelength array and a flux array wavelength, surface_flux = get_phoenix_photons( temperature=5780, logg=4.43, metallicity=0.0, R=1000 ) # plot the wavelength array plt.figure(figsize=(8, 3), dpi=300) plt.plot(wavelength, surface_flux, label="PHOENIX model") # for comparison, plot a Planck spectrum w_planck, f_planck = get_planck_photons(temperature=5780) plt.plot(w_planck, f_planck, label="Planck function") # add some labels plt.xlabel(f"Wavelength ({wavelength.unit.to_string('latex_inline')})") plt.ylabel(f"Surface Flux ({surface_flux.unit.to_string('latex_inline')})") plt.legend(frameon=False);

The solar spectrum looks approximately like the Planck thermal emission spectrum, but with some absorption features and redistribution of flux due radiative transfer through the stellar atmosphere.

How do we set the spectral resolution?¶

If you specify a spectral resolution $R = \lambda/d\lambda$, the full extent of the model wavelength range (approximately $0.05-5 \mu m$) will be returned, binned to that resolution. Different resolutions are stored in different files that are loaded as needed. Your code will run a lot faster if you load the lowest resolution that will meet your needs. For a single metallicity, files range from a very manageable 400 kilobytes for $R=10$ to an annoyingling large 3 gigabytes for $R=100000$. Trying to calculate an extremely high resolution model pre-existing metallicity grid point will try to load in all metallicities, which can cause slow-downs because of memory issues.

In [4]:

plt.figure(figsize=(8, 3), dpi=300)
for R in [10000, 1000, 100, 10]:
    w, f = get_phoenix_photons(R=R)
    plt.plot(w, f, alpha=0.5, label=f"R={R}")
plt.xlabel(f"Wavelength ({wavelength.unit.to_string('latex_inline')})")
plt.ylabel(f"Surface Flux ({surface_flux.unit.to_string('latex_inline')})")
plt.legend(frameon=False);

plt.figure(figsize=(8, 3), dpi=300) for R in [10000, 1000, 100, 10]: w, f = get_phoenix_photons(R=R) plt.plot(w, f, alpha=0.5, label=f"R={R}") plt.xlabel(f"Wavelength ({wavelength.unit.to_string('latex_inline')})") plt.ylabel(f"Surface Flux ({surface_flux.unit.to_string('latex_inline')})") plt.legend(frameon=False);

Should we specify custom wavelengths?¶

For many applications, you may need to generate lots of spectra over a very narrow wavelength range. In those cases, it'd be inefficient to interpolate the a very large spectral range only to trim down to a few specific wavelengths. If you know exactly the wavelengths you want and you're generating more than one spectrum on the same wavelength grid, your code will run likely faster if you specify the wavelengths you need directly to get_phoenix_photons, as shown below.

In [5]:

plt.figure(figsize=(8, 3), dpi=300)
w, f = get_phoenix_photons(wavelength=np.linspace(0.4, 0.7, 1000) * u.micron)
plt.plot(w, f)
plt.xlabel(f"Wavelength ({wavelength.unit.to_string('latex_inline')})")
plt.ylabel(f"Surface Flux ({surface_flux.unit.to_string('latex_inline')})");

plt.figure(figsize=(8, 3), dpi=300) w, f = get_phoenix_photons(wavelength=np.linspace(0.4, 0.7, 1000) * u.micron) plt.plot(w, f) plt.xlabel(f"Wavelength ({wavelength.unit.to_string('latex_inline')})") plt.ylabel(f"Surface Flux ({surface_flux.unit.to_string('latex_inline')})");

What are the units?¶

The variables wavelength and surface_flux returned by get_phoenix_photons have astropy units attached to them, to reduce the risk of unit mistakes in subsequent calculations.

In [6]:

wavelength.unit

wavelength.unit

Out[6]:

$\mathrm{\mu m}$

In [7]:

surface_flux.unit

surface_flux.unit

Out[7]:

$\mathrm{\frac{ph}{nm\,s\,m^{2}}}$

Because chromatic primarily uses these model spectra for photon-noise calculations, the flux is provided in photons/second rather than energy units like joules/second (= watts) or ergs/second.

What does "surface flux" mean?¶

In general, you could characterize the brightness of a star at lots of different locations: the surface of the star, at a fixed orbital distance from the star, or at some even larger astronomical distance from the star. The flux returned by get_phoenix_photons is a surface flux, in the sense that it represents the rate of photons flowing out of a star's photosphere per unit area. To show how this might be used, the following example calculation integrates a model spectrum to estimate the bolometric luminosity of the Sun.

In [8]:

# retrieve a wavelength array and a flux array
wavelength, surface_flux = get_phoenix_photons(
    temperature=5780, logg=4.43, metallicity=0.0, R=1000
)

# calculate the energy per photon for each wavelength [J/photon]
energy_per_photon = (con.h * con.c / wavelength) / u.photon

# calculate the surface flux in power units [W/(m**2 nm)]
surface_flux_power = energy_per_photon * surface_flux

# integrate over wavelength to get bolometric surface flux [W/m**2]
ok = np.isfinite(surface_flux_power)
bolometric_surface_flux_power = np.trapz(surface_flux_power[ok], wavelength[ok])

# calculate the surface area of the Sun [m**2]
surface_area = 4 * np.pi * (1 * u.Rsun) ** 2

# calculate the luminosity as bolometric flux * area
luminosity = bolometric_surface_flux_power * surface_area
luminosity.to(u.Lsun)

# retrieve a wavelength array and a flux array wavelength, surface_flux = get_phoenix_photons( temperature=5780, logg=4.43, metallicity=0.0, R=1000 ) # calculate the energy per photon for each wavelength [J/photon] energy_per_photon = (con.h * con.c / wavelength) / u.photon # calculate the surface flux in power units [W/(m**2 nm)] surface_flux_power = energy_per_photon * surface_flux # integrate over wavelength to get bolometric surface flux [W/m**2] ok = np.isfinite(surface_flux_power) bolometric_surface_flux_power = np.trapz(surface_flux_power[ok], wavelength[ok]) # calculate the surface area of the Sun [m**2] surface_area = 4 * np.pi * (1 * u.Rsun) ** 2 # calculate the luminosity as bolometric flux * area luminosity = bolometric_surface_flux_power * surface_area luminosity.to(u.Lsun)

Out[8]:

$0.99032735 \; \mathrm{L_{\odot}}$

To review, we converted from photons to energy, integrated over wavelength, and multiplied by the surface area of the Sun. Our estimate for the luminosity of the Sun came out pretty close, especially considering we didn't actually integrate the wavelengths from $\lambda=0$ to $\lambda=\infty$.

What kind of interpolation is happening?¶

The PHOENIX model grid provides spectra at fixed intervals of temperature $T_{\rm eff}$ (100-200K), surface gravity $\log{g}$ (0.5 dex), and metallicity $[\mathrm{Fe/H}]$ (0.5 dex). If you request a model that exists in that grid, it will be returned exactly. If you request a model somewhere between these grid points, it will be interpolated from the closest available grid points. The interpolation is linear in the quantities $\log T_{\rm eff}$, $\log{g}$, and $[\mathrm{Fe/H}]$. We interpolate in the logarithm of temperature because per-wavelength fluxes tend to grow $\propto T_{\rm eff}^{x}$ where $x\ne 1$.

In [9]:

fi, ax = plt.subplots(1, 2, figsize=(8, 3), dpi=300, constrained_layout=True)

# loop over two different wavelengths
for a, i in zip(ax, [10, -10]):

    # point to the appropriate plotting axes
    plt.sca(a)

    # plot flux for the original grid points
    grid_temperatures = np.arange(4000, 4400, 100)
    grid_fluxes = []
    for T in grid_temperatures:
        w, f = get_phoenix_photons(temperature=T, logg=4.5)
        grid_fluxes.append(f)
    plt.plot(grid_temperatures, np.array(grid_fluxes)[:, i], marker="o", label="grid")

    # plot flux for interpolated temperatures
    interpolated_temperatures = np.arange(4000, 4300, 11)
    interpolated_fluxes = []
    for T in interpolated_temperatures:
        w, f = get_phoenix_photons(temperature=T, logg=4.5)
        interpolated_fluxes.append(f)
    plt.plot(
        interpolated_temperatures,
        np.array(interpolated_fluxes)[:, i],
        marker=".",
        label="interpolated",
    )

    # add some labels
    plt.title(f"$\lambda$={w[i]:.3f}")
    plt.xlabel("Temperature (K)")
plt.sca(ax[0])
plt.ylabel(f"Monochromatic Flux\n({f.unit.to_string('latex_inline')})")
plt.legend(frameon=False);

fi, ax = plt.subplots(1, 2, figsize=(8, 3), dpi=300, constrained_layout=True) # loop over two different wavelengths for a, i in zip(ax, [10, -10]): # point to the appropriate plotting axes plt.sca(a) # plot flux for the original grid points grid_temperatures = np.arange(4000, 4400, 100) grid_fluxes = [] for T in grid_temperatures: w, f = get_phoenix_photons(temperature=T, logg=4.5) grid_fluxes.append(f) plt.plot(grid_temperatures, np.array(grid_fluxes)[:, i], marker="o", label="grid") # plot flux for interpolated temperatures interpolated_temperatures = np.arange(4000, 4300, 11) interpolated_fluxes = [] for T in interpolated_temperatures: w, f = get_phoenix_photons(temperature=T, logg=4.5) interpolated_fluxes.append(f) plt.plot( interpolated_temperatures, np.array(interpolated_fluxes)[:, i], marker=".", label="interpolated", ) # add some labels plt.title(f"$\lambda$={w[i]:.3f}") plt.xlabel("Temperature (K)") plt.sca(ax[0]) plt.ylabel(f"Monochromatic Flux\n({f.unit.to_string('latex_inline')})") plt.legend(frameon=False);

<>:31: SyntaxWarning: invalid escape sequence '\l'
<>:31: SyntaxWarning: invalid escape sequence '\l'
/var/folders/sz/5zjhq1h17p7gsb7fr5jr1lch0000gp/T/ipykernel_52838/1642702033.py:31: SyntaxWarning: invalid escape sequence '\l'
  plt.title(f"$\lambda$={w[i]:.3f}")

Where are the spectra stored?¶

The get_phoenix_photons function automatically downloads the files it needs and caches them on your local computer. You can find out where those files are stored and how much space their taking up by interacting with the phoenix_library object (of which get_phoenix_photons is a method).

In [10]:

from chromatic import phoenix_library

from chromatic import phoenix_library

In [11]:

phoenix_library.get_cache_dir()

phoenix_library.get_cache_dir()

Out[11]:

'/Users/zach/.chromatic/cache'

In [12]:

phoenix_library.get_cache_size()

phoenix_library.get_cache_size()

Out[12]:

$41.779267 \; \mathrm{Gbyte}$

How long do model retrievals take?¶

We put some effort into speeding up the process of retrieving multiple different spectral models on similar wavelength grids. The following plot shows approximately how long different steps take on a 2020 MacBook Pro with M1 chip.

In [13]:

phoenix_library.plot_time_required();

phoenix_library.plot_time_required();