In this tutorial, we will compare the photometric fluxes of the brown dwarf companion PZ Tel B with a synthetic spectrum from the ATMO grid.
We start by importing species.
import species
And initiating the workflow with the SpeciesInit
class. This will create the configuration file and the HDF5 database.
species.SpeciesInit()
Initiating species v0.3.1... [DONE] Creating species_config.ini... [DONE] Database: /Users/tomasstolker/applications/species/docs/tutorials/species_database.hdf5 Data folder: /Users/tomasstolker/applications/species/docs/tutorials/data Working folder: /Users/tomasstolker/applications/species/docs/tutorials Creating species_database.hdf5... [DONE] Creating data folder... [DONE]
<species.core.setup.SpeciesInit at 0x144c1b1d0>
We create now a Database
object which is used for importing various data into the database.
database = species.Database()
The spectra of ATMO are downloaded and added to the database with the add_model
method of Database
. This requires sufficient disk storage in the data_folder that is set in the configuration file. The full ATMO grid is downloaded but the teff_range
parameter can be used to only import a certain Teff range into the database.
database.add_model('atmo', teff_range=(2500., 3000.))
Downloading ATMO model spectra (430 MB)... [DONE] Unpacking ATMO model spectra (430 MB)... [DONE] Adding ATMO model spectra... [DONE] Grid points stored in the database: - Teff = [2500. 2600. 2700. 2800. 2900. 3000.] - log(g) = [2.5 3. 3.5 4. 4.5 5. 5.5] Number of grid points per parameter: - teff: 6 - logg: 7 Fix missing grid points with a linear interpolation: Number of stored grid points: 42 Number of interpolated grid points: 0 Number of missing grid points: 0
Next, we add the distance and magnitudes of PZ Tel B to the database with the `add_companion`` method. This will automatically download the required filter profiles and a flux-calibrated spectrum of Vega. These are used to convert the magnitudes into fluxes.
database.add_companion('PZ Tel B')
Downloading Vega spectrum (270 kB)... [DONE] Adding Vega spectrum... [DONE] Adding filter: Paranal/SPHERE.ZIMPOL_R_PRIM... [DONE] Adding filter: Paranal/SPHERE.ZIMPOL_I_PRIM... [DONE] Adding filter: Paranal/SPHERE.IRDIS_D_H23_2... [DONE] Adding filter: Paranal/SPHERE.IRDIS_D_H23_3... [DONE] Adding filter: Paranal/SPHERE.IRDIS_D_K12_1... [DONE] Adding filter: Paranal/SPHERE.IRDIS_D_K12_2... [DONE] Adding filter: Paranal/NACO.J... [DONE] Adding filter: Paranal/NACO.H... [DONE] Adding filter: Paranal/NACO.Ks... [DONE] Adding filter: Paranal/NACO.Lp... [DONE] Adding filter: Paranal/NACO.NB405... [DONE] Adding filter: Paranal/NACO.Mp... [DONE] Adding filter: Gemini/NICI.ED286... [DONE] Adding filter: Gemini/NIRI.H2S1v2-1-G0220... [DONE] Adding object: PZ Tel B - Distance (pc) = 47.13 +/- 0.13 - Paranal/SPHERE.ZIMPOL_R_PRIM: - Apparent magnitude = 17.84 +/- 0.31 - Flux (W m-2 um-1) = 1.83e-15 +/- 5.29e-16 - Paranal/SPHERE.ZIMPOL_I_PRIM: - Apparent magnitude = 15.16 +/- 0.12 - Flux (W m-2 um-1) = 1.09e-14 +/- 1.20e-15 - Paranal/SPHERE.IRDIS_D_H23_2: - Apparent magnitude = 11.78 +/- 0.19 - Flux (W m-2 um-1) = 2.54e-14 +/- 4.47e-15 - Paranal/SPHERE.IRDIS_D_H23_3: - Apparent magnitude = 11.65 +/- 0.19 - Flux (W m-2 um-1) = 2.43e-14 +/- 4.27e-15 - Paranal/SPHERE.IRDIS_D_K12_1: - Apparent magnitude = 11.56 +/- 0.09 - Flux (W m-2 um-1) = 1.15e-14 +/- 9.58e-16 - Paranal/SPHERE.IRDIS_D_K12_2: - Apparent magnitude = 11.29 +/- 0.10 - Flux (W m-2 um-1) = 1.14e-14 +/- 1.05e-15 - Paranal/NACO.J: - Apparent magnitude = 12.47 +/- 0.20 - Flux (W m-2 um-1) = 3.11e-14 +/- 5.76e-15 - Paranal/NACO.H: - Apparent magnitude = 11.93 +/- 0.14 - Flux (W m-2 um-1) = 1.97e-14 +/- 2.55e-15 - Paranal/NACO.Ks: - Apparent magnitude = 11.53 +/- 0.07 - Flux (W m-2 um-1) = 1.12e-14 +/- 7.25e-16 - Paranal/NACO.Lp: - Apparent magnitude = 11.04 +/- 0.22 - Flux (W m-2 um-1) = 2.02e-15 +/- 4.12e-16 - Paranal/NACO.NB405: - Apparent magnitude = 10.94 +/- 0.07 - Flux (W m-2 um-1) = 1.67e-15 +/- 1.08e-16 - Paranal/NACO.Mp: - Apparent magnitude = 10.93 +/- 0.03 - Flux (W m-2 um-1) = 9.19e-16 +/- 2.54e-17 - Gemini/NICI.ED286: - Apparent magnitude = 11.68 +/- 0.14 - Flux (W m-2 um-1) = 2.78e-14 +/- 3.60e-15 - Gemini/NIRI.H2S1v2-1-G0220: - Apparent magnitude = 11.39 +/- 0.14 - Flux (W m-2 um-1) = 1.06e-14 +/- 1.37e-15
Alternatively, the add_object
method of Database
can be used for manually adding magnitudes and spectra of an individual object. Before coninuing, let's check the content of the database.
database.list_content()
Database content: - filters: <HDF5 group "/filters" (2 members)> - Gemini: <HDF5 group "/filters/Gemini" (2 members)> - NICI.ED286: <HDF5 dataset "NICI.ED286": shape (387, 2), type "<f4"> - det_type: energy - NIRI.H2S1v2-1-G0220: <HDF5 dataset "NIRI.H2S1v2-1-G0220": shape (129, 2), type "<f4"> - det_type: energy - Paranal: <HDF5 group "/filters/Paranal" (12 members)> - NACO.H: <HDF5 dataset "NACO.H": shape (23, 2), type "<f4"> - det_type: energy - NACO.J: <HDF5 dataset "NACO.J": shape (20, 2), type "<f4"> - det_type: energy - NACO.Ks: <HDF5 dataset "NACO.Ks": shape (27, 2), type "<f4"> - det_type: energy - NACO.Lp: <HDF5 dataset "NACO.Lp": shape (31, 2), type "<f4"> - det_type: energy - NACO.Mp: <HDF5 dataset "NACO.Mp": shape (18, 2), type "<f4"> - det_type: energy - NACO.NB405: <HDF5 dataset "NACO.NB405": shape (67, 2), type "<f4"> - det_type: energy - SPHERE.IRDIS_D_H23_2: <HDF5 dataset "SPHERE.IRDIS_D_H23_2": shape (113, 2), type "<f4"> - det_type: energy - SPHERE.IRDIS_D_H23_3: <HDF5 dataset "SPHERE.IRDIS_D_H23_3": shape (180, 2), type "<f4"> - det_type: energy - SPHERE.IRDIS_D_K12_1: <HDF5 dataset "SPHERE.IRDIS_D_K12_1": shape (175, 2), type "<f4"> - det_type: energy - SPHERE.IRDIS_D_K12_2: <HDF5 dataset "SPHERE.IRDIS_D_K12_2": shape (191, 2), type "<f4"> - det_type: energy - SPHERE.ZIMPOL_I_PRIM: <HDF5 dataset "SPHERE.ZIMPOL_I_PRIM": shape (189, 2), type "<f4"> - det_type: energy - SPHERE.ZIMPOL_R_PRIM: <HDF5 dataset "SPHERE.ZIMPOL_R_PRIM": shape (169, 2), type "<f4"> - det_type: energy - models: <HDF5 group "/models" (1 members)> - atmo: <HDF5 group "/models/atmo" (4 members)> - flux: <HDF5 dataset "flux": shape (6, 7, 98116), type "<f8"> - logg: <HDF5 dataset "logg": shape (7,), type "<f8"> - teff: <HDF5 dataset "teff": shape (6,), type "<f8"> - wavelength: <HDF5 dataset "wavelength": shape (98116,), type "<f8"> - objects: <HDF5 group "/objects" (1 members)> - PZ Tel B: <HDF5 group "/objects/PZ Tel B" (3 members)> - Gemini: <HDF5 group "/objects/PZ Tel B/Gemini" (2 members)> - NICI.ED286: <HDF5 dataset "NICI.ED286": shape (4,), type "<f8"> - n_phot: 1 - NIRI.H2S1v2-1-G0220: <HDF5 dataset "NIRI.H2S1v2-1-G0220": shape (4,), type "<f8"> - n_phot: 1 - Paranal: <HDF5 group "/objects/PZ Tel B/Paranal" (12 members)> - NACO.H: <HDF5 dataset "NACO.H": shape (4,), type "<f8"> - n_phot: 1 - NACO.J: <HDF5 dataset "NACO.J": shape (4,), type "<f8"> - n_phot: 1 - NACO.Ks: <HDF5 dataset "NACO.Ks": shape (4,), type "<f8"> - n_phot: 1 - NACO.Lp: <HDF5 dataset "NACO.Lp": shape (4,), type "<f8"> - n_phot: 1 - NACO.Mp: <HDF5 dataset "NACO.Mp": shape (4,), type "<f8"> - n_phot: 1 - NACO.NB405: <HDF5 dataset "NACO.NB405": shape (4,), type "<f8"> - n_phot: 1 - SPHERE.IRDIS_D_H23_2: <HDF5 dataset "SPHERE.IRDIS_D_H23_2": shape (4,), type "<f8"> - n_phot: 1 - SPHERE.IRDIS_D_H23_3: <HDF5 dataset "SPHERE.IRDIS_D_H23_3": shape (4,), type "<f8"> - n_phot: 1 - SPHERE.IRDIS_D_K12_1: <HDF5 dataset "SPHERE.IRDIS_D_K12_1": shape (4,), type "<f8"> - n_phot: 1 - SPHERE.IRDIS_D_K12_2: <HDF5 dataset "SPHERE.IRDIS_D_K12_2": shape (4,), type "<f8"> - n_phot: 1 - SPHERE.ZIMPOL_I_PRIM: <HDF5 dataset "SPHERE.ZIMPOL_I_PRIM": shape (4,), type "<f8"> - n_phot: 1 - SPHERE.ZIMPOL_R_PRIM: <HDF5 dataset "SPHERE.ZIMPOL_R_PRIM": shape (4,), type "<f8"> - n_phot: 1 - distance: <HDF5 dataset "distance": shape (2,), type "<f8"> - spectra: <HDF5 group "/spectra" (1 members)> - calibration: <HDF5 group "/spectra/calibration" (1 members)> - vega: <HDF5 dataset "vega": shape (3, 8827), type "<f8">
We see the various groups, subgroups, datasets, and attributes that are stored in the HDF5 database.
Model spectra are read from the database by first creating an instance of ReadModel
. The model name and optionally a wavelength range are provided as arguments.
readmodel = species.ReadModel('atmo', wavel_range=(0.5, 10.))
Before extracting a spectrum, let's check which parameters are required for the ATMO model spectra.
readmodel.get_parameters()
['teff', 'logg']
And also the parameter boundaries of the grid that is stored in the database.
readmodel.get_bounds()
{'teff': (2500.0, 3000.0), 'logg': (2.5, 5.5)}
The parameters are provided in a dictionary for which we have to make sure that chose values are within the grid boundaries. The radius (RJup) and distance (pc) will scale the emitted spectrum to the observer. Without these values, the spectrum fluxes are provided at the surface of the atmosphere.
model_param = {'teff': 2900., 'logg': 4.5, 'radius': 2.2, 'distance': 47.13}
We now use the get_model
method of ReadModel
to linearly interpolate the grid of spectra and store the extracted spectrum in a ModelBox
. The spectrum is smoothed to a spectral resolution of R = 100.
modelbox = readmodel.get_model(model_param, spec_res=100., smooth=True)
The photometric data of PZ Tel B are also read from the database and stored in an ObjectBox
.
objectbox = database.get_object(object_name='PZ Tel B')
Getting object: PZ Tel B... [DONE]
For comparison, we create synthetic photometry from the extracted ATMO spectrum for all filters of PZ Tel B. The synthetic fluxes are stored in a SynphotBox
.
synphotbox = species.multi_photometry(datatype='model',
spectrum='atmo',
filters=objectbox.filters,
parameters=model_param)
Calculating synthetic photometry... [DONE]
The get_residuals
function is now used to calculate the difference between the observed fluxes and the synthetic fluxes from the model spectrum. The residuals are stored in a ResidualsBox
.
res_box = species.get_residuals(datatype='model',
spectrum='atmo',
parameters=model_param,
objectbox=objectbox,
inc_phot=True,
inc_spec=False)
Calculating synthetic photometry... [DONE] Calculating residuals... [DONE] Residuals (sigma): - Gemini/NICI.ED286: 1.35 - Gemini/NIRI.H2S1v2-1-G0220: -0.02 - Paranal/NACO.H: -0.54 - Paranal/NACO.J: -0.19 - Paranal/NACO.Ks: 0.85 - Paranal/NACO.Lp: 0.02 - Paranal/NACO.Mp: 5.45 - Paranal/NACO.NB405: 0.31 - Paranal/SPHERE.IRDIS_D_H23_2: 0.54 - Paranal/SPHERE.IRDIS_D_H23_3: 0.00 - Paranal/SPHERE.IRDIS_D_K12_1: 1.00 - Paranal/SPHERE.IRDIS_D_K12_2: 1.00 - Paranal/SPHERE.ZIMPOL_I_PRIM: -7.73 - Paranal/SPHERE.ZIMPOL_R_PRIM: -3.96
The open_box
function can be used to view the content of a Box
object. For example, the ModelBox
contains several attributes, including the wavelengths and fluxes.
modelbox.open_box()
Opening ModelBox... model = atmo type = None wavelength = [ 0.4999641 0.50001364 0.50006318 ... 9.99858711 9.99957782 10.00056862] flux = [2.05385551e-15 2.05647449e-15 2.05928733e-15 ... 5.56118275e-17 5.56005834e-17 5.55895317e-17] parameters = {'teff': 2900.0, 'logg': 4.5, 'radius': 2.2, 'distance': 47.13, 'mass': 59.04988128500266, 'luminosity': 0.003114426265448587} quantity = flux
Similarly, an ObjectBox
contains a dictionary with the magnitudes and a dictionary with the fluxes.
objectbox.open_box()
Opening ObjectBox... name = PZ Tel B filters = ['Gemini/NICI.ED286', 'Gemini/NIRI.H2S1v2-1-G0220', 'Paranal/NACO.H', 'Paranal/NACO.J', 'Paranal/NACO.Ks', 'Paranal/NACO.Lp', 'Paranal/NACO.Mp', 'Paranal/NACO.NB405', 'Paranal/SPHERE.IRDIS_D_H23_2', 'Paranal/SPHERE.IRDIS_D_H23_3', 'Paranal/SPHERE.IRDIS_D_K12_1', 'Paranal/SPHERE.IRDIS_D_K12_2', 'Paranal/SPHERE.ZIMPOL_I_PRIM', 'Paranal/SPHERE.ZIMPOL_R_PRIM'] magnitude = {'Gemini/NICI.ED286': array([11.68, 0.14]), 'Gemini/NIRI.H2S1v2-1-G0220': array([11.39, 0.14]), 'Paranal/NACO.H': array([11.93, 0.14]), 'Paranal/NACO.J': array([12.47, 0.2 ]), 'Paranal/NACO.Ks': array([11.53, 0.07]), 'Paranal/NACO.Lp': array([11.04, 0.22]), 'Paranal/NACO.Mp': array([10.93, 0.03]), 'Paranal/NACO.NB405': array([10.94, 0.07]), 'Paranal/SPHERE.IRDIS_D_H23_2': array([11.78, 0.19]), 'Paranal/SPHERE.IRDIS_D_H23_3': array([11.65, 0.19]), 'Paranal/SPHERE.IRDIS_D_K12_1': array([11.56, 0.09]), 'Paranal/SPHERE.IRDIS_D_K12_2': array([11.29, 0.1 ]), 'Paranal/SPHERE.ZIMPOL_I_PRIM': array([15.16, 0.12]), 'Paranal/SPHERE.ZIMPOL_R_PRIM': array([17.84, 0.31])} flux = {'Gemini/NICI.ED286': array([2.78256301e-14, 3.59792032e-15]), 'Gemini/NIRI.H2S1v2-1-G0220': array([1.05904383e-14, 1.36936892e-15]), 'Paranal/NACO.H': array([1.96875856e-14, 2.54565177e-15]), 'Paranal/NACO.J': array([3.11068456e-14, 5.76255347e-15]), 'Paranal/NACO.Ks': array([1.12416276e-14, 7.25276730e-16]), 'Paranal/NACO.Lp': array([2.02006233e-15, 4.12126882e-16]), 'Paranal/NACO.Mp': array([9.18778503e-16, 2.53900187e-17]), 'Paranal/NACO.NB405': array([1.67046281e-15, 1.07773345e-16]), 'Paranal/SPHERE.IRDIS_D_H23_2': array([2.54130005e-14, 4.46991835e-15]), 'Paranal/SPHERE.IRDIS_D_H23_3': array([2.42699566e-14, 4.26886720e-15]), 'Paranal/SPHERE.IRDIS_D_K12_1': array([1.15478089e-14, 9.58329864e-16]), 'Paranal/SPHERE.IRDIS_D_K12_2': array([1.14281212e-14, 1.05405765e-15]), 'Paranal/SPHERE.ZIMPOL_I_PRIM': array([1.08628802e-14, 1.20305574e-15]), 'Paranal/SPHERE.ZIMPOL_R_PRIM': array([1.82633135e-15, 5.28569080e-16])} distance = [47.13 0.13] spectrum = None
The attributes in a Box
object can be extracted for further analyis or creating plots. For example, to extract the array with wavelengths from the ModelBox
:
modelbox.wavelength
array([ 0.4999641 , 0.50001364, 0.50006318, ..., 9.99858711, 9.99957782, 10.00056862])
Finally, we will combine the model spectrum and the photometric fluxes in a plot with plot_spectrum
. A list with Box
objects is provided as an argument of boxes
. These are interpreted accordingly by the plot_spectrum
function. Also a list with filter names can be provided as argument of filters
to show the filter profiles. The ResidualsBox
is provided as arguments of residuals
. Finally, the optional argument of plot_kwargs
contains a list with optional dictionaries to tune the visualization of the plotted data. The number of items in the list of plot_kwargs
should be equal to the number of items in the list of boxes
. For the SynphotBox
, we can set the item in plot_kwargs
to None
such that the marker design is based on the data from ObjectBox
.
species.plot_spectrum(boxes=[modelbox, objectbox, synphotbox],
filters=objectbox.filters,
residuals=res_box,
plot_kwargs=[{'ls': '-', 'lw': 1., 'color': 'black'},
{'Gemini/NICI.ED286': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},
'Gemini/NIRI.H2S1v2-1-G0220': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},
'Paranal/NACO.H': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},
'Paranal/NACO.J': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},
'Paranal/NACO.Ks': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},
'Paranal/NACO.Lp': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},
'Paranal/NACO.Mp': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},
'Paranal/NACO.NB405': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},
'Paranal/SPHERE.IRDIS_D_H23_2': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},
'Paranal/SPHERE.IRDIS_D_H23_3': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},
'Paranal/SPHERE.IRDIS_D_K12_1': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},
'Paranal/SPHERE.IRDIS_D_K12_2': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},
'Paranal/SPHERE.ZIMPOL_I_PRIM': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},
'Paranal/SPHERE.ZIMPOL_R_PRIM': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'}},
None],
xlim=(0.5, 5.5),
ylim=(-5e-15, 5.5e-14),
ylim_res=(-8, 8),
scale=('linear', 'linear'),
title='PZ Tel B and ATMO spectrum',
offset=(-0.45, -0.04),
legend={'loc': 'upper right', 'frameon': False, 'fontsize': 12},
figsize=(9., 4),
quantity='flux',
output='spectrum.png')
Plotting spectrum: spectrum.png... [DONE]
Let's have a look at the result. The plot is stored in the working folder. The blue squares are the photometric fluxes of PZ Tel B and the open squares are the synthetic photometry computed from the model spectrum. The residuals are shown relative to the uncertainties on the fluxes.
from IPython.display import Image
Image('spectrum.png')