Model spectra and synthetic photometry#

In this tutorial, we will have a look at some spectra of the DRIFT-PHOENIX atmospheric model. The spectra are first downloaded and added to the database. Then we will use the functionalities of ReadModel to extract a spectrum and calculate a photometric flux.

Getting started#

We start by importing the required Python modules.

[1]:
from species import SpeciesInit
from species.data.database import Database
from species.read.read_model import ReadModel
from species.plot.plot_spectrum import plot_spectrum

Then we initialize species with SpeciesInit, which creates a default configuration file and the HDF5 database.

[2]:
SpeciesInit()
==============
species v0.7.4
==============

Working folder: /Users/tomasstolker/applications/species/docs/tutorials
Creating species_config.ini... [DONE]
Configuration settings:
   - Database: /Users/tomasstolker/applications/species/docs/tutorials/species_database.hdf5
   - Data folder: /Users/tomasstolker/applications/species/docs/tutorials/data
   - Interpolation method: linear
   - Magnitude of Vega: 0.03
Creating species_database.hdf5... [DONE]
Creating data folder... [DONE]
[2]:
<species.core.species_init.SpeciesInit at 0x1112e8510>

Adding model spectra to the database#

To store the spectra, we first create an instance of Database.

[3]:
database = Database()

Let’s check which atmospheric models are available by running the available_models method of the Database object.

[4]:
_ = database.available_models()
Available model grids:
   - AMES-Cond:
      - Label = ames-cond
      - Model parameters: ['teff', 'logg']
      - Teff range (K): [100, 6600]
      - Wavelength range (um): [0.5, 40]
      - Resolution lambda/Dlambda: 4000
      - File size: 150 MB

   - AMES-Dusty:
      - Label = ames-dusty
      - Model parameters: ['teff', 'logg']
      - Teff range (K): [500, 4000]
      - Wavelength range (um): [0.5, 40]
      - Resolution lambda/Dlambda: 4000
      - File size: 58 MB

   - ATMO:
      - Label = atmo
      - Model parameters: ['teff', 'logg']
      - Teff range (K): [200, 3000]
      - Wavelength range (um): [0.4, 6000]
      - Resolution lambda/Dlambda: 10000
      - File size: 425 MB
      - Reference: Phillips et al. (2020)
      - URL: https://ui.adsabs.harvard.edu/abs/2020A%26A...637A..38P/abstract

   - ATMO CEQ:
      - Label = atmo-ceq
      - Model parameters: ['teff', 'logg']
      - Teff range (K): [200, 3000]
      - Wavelength range (um): [0.2, 6000]
      - Resolution lambda/Dlambda: 10000
      - File size: 455 MB
      - Reference: Phillips et al. (2020)
      - URL: https://ui.adsabs.harvard.edu/abs/2020A%26A...637A..38P/abstract

   - ATMO NEQ weak:
      - Label = atmo-neq-weak
      - Model parameters: ['teff', 'logg']
      - Teff range (K): [200, 1800]
      - Wavelength range (um): [0.2, 6000]
      - Resolution lambda/Dlambda: 10000
      - File size: 276 MB
      - Reference: Phillips et al. (2020)
      - URL: https://ui.adsabs.harvard.edu/abs/2020A%26A...637A..38P/abstract

   - ATMO NEQ strong:
      - Label = atmo-neq-strong
      - Model parameters: ['teff', 'logg']
      - Teff range (K): [200, 1800]
      - Wavelength range (um): [0.2, 6000]
      - Resolution lambda/Dlambda: 10000
      - File size: 274 MB
      - Reference: Phillips et al. (2020)
      - URL: https://ui.adsabs.harvard.edu/abs/2020A%26A...637A..38P/abstract

   - ATMO (Petrus et al. 2023):
      - Label = atmo-petrus2023
      - Model parameters: ['teff', 'logg', 'feh', 'c_o_ratio', 'ad_index']
      - Teff range (K): [800, 3000]
      - Wavelength range (um): [0.2, 30]
      - Resolution lambda/Dlambda: 3000
      - File size: 2.1 GB
      - Reference: Petrus et al. (2023)
      - URL: https://ui.adsabs.harvard.edu/abs/2023A%26A...670L...9P/abstract

   - blackbody:
      - Label = blackbody
      - Model parameters: ['teff']
      - Teff range (K): [10, 20000]
      - Wavelength range (um): [0.1, 5000]
      - Resolution lambda/Dlambda: 1000
      - File size: 56 MB

   - BT-Cond:
      - Label = bt-cond
      - Model parameters: ['teff', 'logg']
      - Teff range (K): [800, 4000]
      - Wavelength range (um): [0.1, 100]
      - Resolution lambda/Dlambda: 10000
      - File size: 162 MB

   - BT-Dusty:
      - Label = bt-dusty
      - Model parameters: ['teff', 'logg']
      - Teff range (K): [1000, 4000]
      - Wavelength range (um): [0.1, 100]
      - Resolution lambda/Dlambda: 10000
      - File size: 205 MB

   - BT-Cond [Fe/H]:
      - Label = bt-cond-feh
      - Model parameters: ['teff', 'logg', 'feh']
      - Teff range (K): [2600, 4000]
      - Wavelength range (um): [0.1, 100]
      - Resolution lambda/Dlambda: 10000
      - File size: 390 MB

   - BT-NextGen:
      - Label = bt-nextgen
      - Model parameters: ['teff', 'logg', 'feh']
      - Teff range (K): [2600, 30000]
      - Wavelength range (um): [0.1, 50]
      - Resolution lambda/Dlambda: 4000
      - File size: 368 MB

   - BT-NextGen subsolar:
      - Label = bt-nextgen-subsolar
      - Model parameters: ['teff', 'logg', 'feh']
      - Teff range (K): [2600, 30000]
      - Wavelength range (um): [0.1, 50]
      - Resolution lambda/Dlambda: 4000
      - Extra details: [alpha/Fe] = 0.4
      - File size: 1.5 GB

   - BT-Settl:
      - Label = bt-settl
      - Model parameters: ['teff', 'logg']
      - Teff range (K): [400, 4000]
      - Wavelength range (um): [0.1, 100]
      - Resolution lambda/Dlambda: 10000
      - File size: 227 MB

   - BT-Settl CIFIST:
      - Label = bt-settl-cifist
      - Model parameters: ['teff', 'logg']
      - Teff range (K): [1200, 7000]
      - Wavelength range (um): [0.1, 5000]
      - Resolution lambda/Dlambda: 20000
      - File size: 1.4 GB

   - DRIFT-PHOENIX:
      - Label = drift-phoenix
      - Model parameters: ['teff', 'logg', 'feh']
      - Teff range (K): [1000, 3000]
      - Wavelength range (um): [0.1, 50]
      - Resolution lambda/Dlambda: 4000
      - File size: 229 MB
      - Reference: Helling et al. (2008)
      - URL: https://ui.adsabs.harvard.edu/abs/2008ApJ...675L.105H/abstract

   - Exo-REM:
      - Label = exo-rem
      - Model parameters: ['teff', 'logg', 'feh', 'c_o_ratio']
      - Teff range (K): [400, 2000]
      - Wavelength range (um): [0.35, 250.0]
      - Resolution lambda/Dlambda: 500
      - File size: 530 MB
      - Reference: Charnay et al. (2018)
      - URL: https://ui.adsabs.harvard.edu/abs/2018ApJ...854..172C/abstract

   - Exo-REM:
      - Label = exo-rem-highres
      - Model parameters: ['teff', 'logg', 'feh', 'c_o_ratio']
      - Teff range (K): [400, 2000]
      - Wavelength range (um): [0.67, 250.0]
      - Resolution lambda/Dlambda: 20000
      - File size: 16 GB
      - Reference: Charnay et al. (2018)
      - URL: https://ui.adsabs.harvard.edu/abs/2018ApJ...854..172C/abstract

   - Koester WD:
      - Label = koester-wd
      - Model parameters: ['teff', 'logg']
      - Teff range (K): [5000, 80000]
      - Wavelength range (um): [0.09, 3]
      - Resolution lambda/Dlambda: 2000
      - File size: 291 MB
      - Reference: Koester (2010)
      - URL: https://ui.adsabs.harvard.edu/abs/2010MmSAI..81..921K/abstract

   - Morley et al. (2012):
      - Label = morley-2012
      - Model parameters: ['teff', 'logg', 'fsed']
      - Teff range (K): [400, 1300]
      - Wavelength range (um): [0.6, 30]
      - Resolution lambda/Dlambda: 10000
      - File size: 141 MB
      - Reference: Morley et al. (2012)
      - URL: https://ui.adsabs.harvard.edu/abs/2012ApJ...756..172M/abstract

   - petitCODE cool clear:
      - Label = petitcode-cool-clear
      - Model parameters: ['teff', 'logg', 'feh']
      - Teff range (K): [500, 1000]
      - Wavelength range (um): [0.1, 250]
      - Resolution lambda/Dlambda: 1000
      - File size: 86 MB
      - Reference: Mollière et al. (2015)
      - URL: https://ui.adsabs.harvard.edu/abs/2015ApJ...813...47M/abstract

   - petitCODE cool cloudy:
      - Label = petitcode-cool-cloudy
      - Model parameters: ['teff', 'logg', 'feh', 'fsed']
      - Teff range (K): [500, 850]
      - Wavelength range (um): [0.1, 250]
      - Resolution lambda/Dlambda: 1000
      - Extra details: log(Kzz) = 7.5
      - File size: 224 MB
      - Reference: Mollière et al. (2015)
      - URL: https://ui.adsabs.harvard.edu/abs/2015ApJ...813...47M/abstract

   - petitCODE hot clear:
      - Label = petitcode-hot-clear
      - Model parameters: ['teff', 'logg', 'feh', 'c_o_ratio']
      - Teff range (K): [1000, 2000]
      - Wavelength range (um): [0.1, 250]
      - Resolution lambda/Dlambda: 1000
      - File size: 42 MB
      - Reference: Mollière et al. (2015)
      - URL: https://ui.adsabs.harvard.edu/abs/2015ApJ...813...47M/abstract

   - petitCODE hot cloudy:
      - Label = petitcode-hot-cloudy
      - Model parameters: ['teff', 'logg', 'feh', 'c_o_ratio', 'fsed']
      - Teff range (K): [1000, 2000]
      - Wavelength range (um): [0.1, 250]
      - Resolution lambda/Dlambda: 1000
      - File size: 126 MB
      - Reference: Mollière et al. (2015)
      - URL: https://ui.adsabs.harvard.edu/abs/2015ApJ...813...47M/abstract

   - petitCODE clear (Linder et al. 2019):
      - Label = petitcode-linder2019-clear
      - Model parameters: ['teff', 'logg', 'feh']
      - Teff range (K): [150, 1000]
      - Wavelength range (um): [0.1, 250]
      - Resolution lambda/Dlambda: 1000
      - File size: 225 MB
      - Reference: Linder et al. (2019)
      - URL: https://ui.adsabs.harvard.edu/abs/2019A%26A...623A..85L/abstract

   - petitCODE cloudy (Linder et al. 2019):
      - Label = petitcode-linder2019-cloudy
      - Model parameters: ['teff', 'logg', 'feh', 'fsed']
      - Teff range (K): [150, 1000]
      - Wavelength range (um): [0.1, 250]
      - Resolution lambda/Dlambda: 1000
      - File size: 1.4 GB
      - Reference: Linder et al. (2019)
      - URL: https://ui.adsabs.harvard.edu/abs/2019A%26A...623A..85L/abstract

   - Saumon & Marley (2008) clear:
      - Label = saumon2008-clear
      - Model parameters: ['teff', 'logg']
      - Teff range (K): [500, 2400]
      - Wavelength range (um): [0.4, 50]
      - Resolution lambda/Dlambda: 5000
      - File size: 197 Mb

   - Saumon & Marley (2008) cloudy:
      - Label = saumon2008-cloudy
      - Model parameters: ['teff', 'logg']
      - Teff range (K): [500, 2400]
      - Wavelength range (um): [0.4, 50]
      - Resolution lambda/Dlambda: 5000
      - File size: 122 Mb

   - Sonora Cholla:
      - Label = sonora-cholla
      - Model parameters: ['teff', 'logg', 'log_kzz']
      - Teff range (K): [500, 1300]
      - Wavelength range (um): [1, 250]
      - Resolution lambda/Dlambda: 10000
      - Extra details: [Fe/H] = 0.0
      - File size: 497 MB
      - Reference: Karalidi et al. (2021)
      - URL: https://zenodo.org/record/4450269

   - Sonora Bobcat:
      - Label = sonora-bobcat
      - Model parameters: ['teff', 'logg', 'feh']
      - Teff range (K): [200, 2400]
      - Wavelength range (um): [0.4, 50]
      - Resolution lambda/Dlambda: 10000
      - File size: 1.1 GB
      - Reference: Marley et al. (2021)
      - URL: https://zenodo.org/record/5063476

   - Sonora Bobcat C/O:
      - Label = sonora-bobcat-co
      - Model parameters: ['teff', 'c_o_ratio']
      - Teff range (K): [200, 2400]
      - Wavelength range (um): [0.4, 50]
      - Resolution lambda/Dlambda: 10000
      - Extra details: [Fe/H] = 0.0, log(g) = 5.0
      - File size: 112 MB
      - Reference: Marley et al. (2021)
      - URL: https://zenodo.org/record/5063476

   - SPHINX:
      - Label = sphinx
      - Model parameters: ['teff', 'logg', 'feh', 'c_o_ratio']
      - Teff range (K): [2000, 4000]
      - Wavelength range (um): [0.1, 20]
      - Resolution lambda/Dlambda: 250
      - File size: 143 Mb
      - Reference: Aishwarya et al. (2022)
      - URL: https://zenodo.org/record/7416042

Next, we will import the model spectra with the add_model method of Database. This step will automatically download the DRIFT-PHOENIX spectra (R = 2000) to the data_folder. The dowloaded data will then be unpacked and added to the database. We restrict the temperature range to 1300-1700 K, so not all spectra are added to the databse.

[5]:
database.add_model(model='drift-phoenix', teff_range=(1300., 1700.))
Downloading data from 'https://home.strw.leidenuniv.nl/~stolker/species/drift-phoenix.tgz' to file '/Users/tomasstolker/applications/species/docs/tutorials/data/drift-phoenix.tgz'.
100%|████████████████████████████████████████| 240M/240M [00:00<00:00, 106GB/s]
SHA256 hash of downloaded file: ba71a5e4d3d399a6f8ae249590c2e174e90ec2b55e712d350dad8ca1ae83a907
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.
Unpacking DRIFT-PHOENIX model spectra (229 MB)... [DONE]
Please cite Helling et al. (2008) when using DRIFT-PHOENIX in a publication
Reference URL: https://ui.adsabs.harvard.edu/abs/2008ApJ...675L.105H/abstract
Wavelength range (um) = 0.1 - 50
Spectral resolution = 4000
Teff range (K) = 1300.0 - 1700.0
Adding DRIFT-PHOENIX model spectra... [DONE]
Grid points stored in the database:
   - Teff = [1300. 1400. 1500. 1600. 1700.]
   - log(g) = [3.  3.5 4.  4.5 5.  5.5]
   - [Fe/H] = [-0.6 -0.3 -0.   0.3]
Number of grid points per parameter:
   - teff: 5
   - logg: 6
   - feh: 4
Fix missing grid points with a linear interpolation:
   - teff = 1600.0, logg = 3.0, feh = 0.3
   - teff = 1600.0, logg = 5.5, feh = 0.3
Number of stored grid points: 120
Number of interpolated grid points: 2
Number of missing grid points: 0
/Users/tomasstolker/applications/species/species/util/data_util.py:379: RuntimeWarning: divide by zero encountered in log10
  flux = np.log10(flux)

Two of the grid points were missing in the original data and have been added with a linear, multidimensional interpolation.

Interpolating the model grid#

We will read the spectra from the database by creating an instance of ReadModel.

[6]:
read_model = ReadModel(model='drift-phoenix', wavel_range=(0.5, 10.))

Let’s see what the grid boundaries are from the spectra that are stored in the database.

[7]:
read_model.get_bounds()
[7]:
{'teff': (1300.0, 1700.0), 'logg': (3.0, 5.5), 'feh': (-0.6, 0.3)}

We will now interpolate the grid in the (\(T_\mathrm{eff}\), \(\log(g)\), \([\mathrm{Fe}/\mathrm{H}]\)) space at some specific parameter values, which need to be provided in a dictionary. The radius and distance are optional, otherwise the emitted flux is given at the planet surface.

[8]:
model_param = {'teff':1510., 'logg':4.1, 'feh':0.1, 'radius': 1., 'distance': 100.}

To interpolate a spectrum, we use the get_model method and provide the parameter dictionary, and also an optional spectral resolution. Together with smooth=True, the spectrum will be smoothed (but not resampeld) to the given spectral resolution.

[9]:
model_box = read_model.get_model(model_param=model_param, spec_res=200., smooth=True)

The data is stored in a ModelBox. Let’s have a look at its content.

[10]:
model_box.open_box()
Opening ModelBox...
model = drift-phoenix
type = None
wavelength = [ 0.49989727  0.50002105  0.50014486 ...  9.99710369  9.99957902
 10.00205496]
flux = [9.87533496e-20 9.80916961e-20 9.73738842e-20 ... 1.60596377e-18
 1.60964903e-18 1.61359501e-18]
parameters = {'teff': 1510.0, 'logg': 4.1, 'feh': 0.1, 'radius': 1.0, 'distance': 100.0, 'luminosity': 4.946207755957789e-05, 'mass': 5.079072973539381}
quantity = flux
contribution = None
bol_flux = None

We will now use the same atmospheric parameters but we will add some ISM extinction with the relation from Cardelli et al. (1989). Therefore, we add the V band extinction and reddening parameters to the dictionary.

[11]:
model_param = {'teff':1510., 'logg':4.1, 'feh':0.1, 'radius': 1., 'distance': 100., 'ism_ext': 5., 'ism_red': 3.}

We use again the get_model method and store the result in a different ModelBox.

[12]:
model_ext = read_model.get_model(model_param=model_param, spec_res=200., smooth=True)

The two boxes with the model spectra are provided to the plot_spectrum. We also include some filter profiles to indicate where the telluric windows are.

[13]:
fig = plot_spectrum(boxes=[model_box, model_ext],
                    filters=['MKO/NSFCam.J', 'MKO/NSFCam.H', 'MKO/NSFCam.K', 'MKO/NSFCam.Lp', 'MKO/NSFCam.Mp'],
                    offset=(-0.15, -0.04),
                    xlim=(0.8, 5.),
                    ylim=(0., 5.5e-17),
                    legend={'loc': 'lower right', 'frameon': False, 'fontsize': 10.},
                    figsize=(7., 3.),
                    output=None)
Adding filter: MKO/NSFCam.J... [DONE]
Adding filter: MKO/NSFCam.H... [DONE]
Adding filter: MKO/NSFCam.K... [DONE]
Adding filter: MKO/NSFCam.Lp... [DONE]
Adding filter: MKO/NSFCam.Mp... [DONE]
Plotting spectrum...
../_images/tutorials_model_spectra_31_1.png
 [DONE]

The plot_spectrum function returned the Figure object of the plot. The functionalities of Matplotlib can be used for further customization of the plot. For example, the axes of the plot are stored at the axes attribute of Figure.

[14]:
fig.axes
[14]:
[<Axes: xlabel='Wavelength (μm)', ylabel='$F_\\lambda$ (10$^{-17}$ W m$^{-2}$ μm$^{-1}$)'>,
 <Axes: ylabel='$T_\\lambda$'>]

Extracting a spectrum at a grid point#

It is also possible to extract a spectrum at one of the grid points, which doesn’t require any interpolation. Let’s check with the get_points method what parameter values are stored in the database.

[15]:
read_model.get_points()
[15]:
{'teff': array([1300., 1400., 1500., 1600., 1700.]),
 'logg': array([3. , 3.5, 4. , 4.5, 5. , 5.5]),
 'feh': array([-0.6, -0.3, -0. ,  0.3])}

We create a dictionary with values at one of the grid points.

[16]:
model_param = {'teff':1500., 'logg':4., 'feh':0.}

And now use the get_data method to extract a spectrum.

[17]:
model_full = read_model.get_data(model_param)

Let’s make another plot with plot_spectrum. The spectrum is now shown at the spectral resolution as stored in the database (\(R = 2000\)).

[18]:
fig = plot_spectrum(boxes=[model_full],
                    filters=None,
                    offset=(-0.12, -0.05),
                    xlim=(0.8, 5.),
                    ylim=(0., 1e5),
                    legend={'loc': 'upper right', 'fontsize': 10.},
                    figsize=(7., 3.),
                    output=None)
Plotting spectrum...
../_images/tutorials_model_spectra_42_1.png
 [DONE]

Calculating synthetic photometry#

The ReadModel class can also be used for calculating photometric fluxes and magnitudes. To do so, we create a new instance and set the filter_name argument to the VLT/NACO M’ filter. This will automatically download and add the filter profile.

[19]:
read_model = ReadModel(model='drift-phoenix', filter_name='Paranal/NACO.Mp')
Adding filter: Paranal/NACO.Mp... [DONE]

We create again a dictionary with the parameters but now run the get_flux method, which returns the flux in W m-2 um-1.

[20]:
model_param = {'teff':1510., 'logg':4.1, 'feh':0.1, 'radius': 1., 'distance': 100.}
flux = read_model.get_flux(model_param)
print(f'Flux (W m-2 um-1) = {flux[0]:.2e}')
Downloading data from 'https://archive.stsci.edu/hlsps/reference-atlases/cdbs/current_calspec/alpha_lyr_stis_011.fits' to file '/Users/tomasstolker/applications/species/docs/tutorials/data/alpha_lyr_stis_011.fits'.
100%|████████████████████████████████████████| 288k/288k [00:00<00:00, 135MB/s]
Adding Vega spectrum... [DONE]
Reference: Bohlin et al. 2014, PASP, 126
URL: https://ui.adsabs.harvard.edu/abs/2014PASP..126..711B/abstract
Flux (W m-2 um-1) = 1.39e-17

Since we provided a radius and distance, the emitted flux at the planet surface has been scaled by \((\mathrm{radius}/\mathrm{distance})^2\).

Similarly, we can use the get_magnitude method to calculate the magnitude for the NACO M’ filter. Note that the returned absolute magnitude is set to None if the parameter dictionary does not contain a radius and distance.

[21]:
app_mag, abs_mag = read_model.get_magnitude(model_param)
print(f'Apparent magnitude = {app_mag:.2f}')
print(f'Absolute magnitude = {abs_mag:.2f}')
Apparent magnitude = 15.46
Absolute magnitude = 10.46

As expected, at a distance of 100 pc, the difference between the apparent and absolute magnitude is 5.