Comparing photometric fluxes and model spectra#

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.

Initiating species#

We start by importing species.

[1]:
import species

And initiating the workflow with the SpeciesInit class. This will create the configuration file and the HDF5 database.

[2]:
species.SpeciesInit()
Initiating species v0.5.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]
[2]:
<species.core.init.SpeciesInit at 0x14e25ee00>

Adding model spectra#

We create now a Database object which is used for importing various data into the database.

[3]:
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 \(T_\mathrm{eff}\) range into the database.

[4]:
database.add_model('atmo', teff_range=(2500., 3000.))
Downloading ATMO model spectra (425 MB)... [DONE]
Unpacking ATMO model spectra (425 MB)... [DONE]
Please cite Phillips et al. (2020) when using ATMO in a publication
Reference URL: https://ui.adsabs.harvard.edu/abs/2020A%26A...637A..38P/abstract
Wavelength range (um) = 0.4 - 6000
Spectral resolution = 10000
Teff range (K) = 2500.0 - 3000.0
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]

Adding companion data#

Next, we add the parallax 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.

[5]:
database.add_companion('PZ Tel B', verbose=False)
Downloading Vega spectrum (270 kB)... [DONE]
Adding Vega spectrum... [DONE]
Adding object: PZ Tel B [DONE]

Alternatively, the add_object method of Database can be used for manually adding magnitudes and spectra of an individual object. Before continuing, let’s check the content of the database.

[6]:
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 "<f8">
                        - det_type: energy
                - NIRI.H2S1v2-1-G0220: <HDF5 dataset "NIRI.H2S1v2-1-G0220": shape (129, 2), type "<f8">
                        - det_type: energy
        - Paranal: <HDF5 group "/filters/Paranal" (12 members)>
                - NACO.H: <HDF5 dataset "NACO.H": shape (23, 2), type "<f8">
                        - det_type: energy
                - NACO.J: <HDF5 dataset "NACO.J": shape (20, 2), type "<f8">
                        - det_type: energy
                - NACO.Ks: <HDF5 dataset "NACO.Ks": shape (27, 2), type "<f8">
                        - det_type: energy
                - NACO.Lp: <HDF5 dataset "NACO.Lp": shape (31, 2), type "<f8">
                        - det_type: energy
                - NACO.Mp: <HDF5 dataset "NACO.Mp": shape (18, 2), type "<f8">
                        - det_type: energy
                - NACO.NB405: <HDF5 dataset "NACO.NB405": shape (67, 2), type "<f8">
                        - det_type: energy
                - SPHERE.IRDIS_D_H23_2: <HDF5 dataset "SPHERE.IRDIS_D_H23_2": shape (113, 2), type "<f8">
                        - det_type: energy
                - SPHERE.IRDIS_D_H23_3: <HDF5 dataset "SPHERE.IRDIS_D_H23_3": shape (180, 2), type "<f8">
                        - det_type: energy
                - SPHERE.IRDIS_D_K12_1: <HDF5 dataset "SPHERE.IRDIS_D_K12_1": shape (175, 2), type "<f8">
                        - det_type: energy
                - SPHERE.IRDIS_D_K12_2: <HDF5 dataset "SPHERE.IRDIS_D_K12_2": shape (191, 2), type "<f8">
                        - det_type: energy
                - SPHERE.ZIMPOL_I_PRIM: <HDF5 dataset "SPHERE.ZIMPOL_I_PRIM": shape (189, 2), type "<f8">
                        - det_type: energy
                - SPHERE.ZIMPOL_R_PRIM: <HDF5 dataset "SPHERE.ZIMPOL_R_PRIM": shape (169, 2), type "<f8">
                        - det_type: energy
- models: <HDF5 group "/models" (1 members)>
        - atmo: <HDF5 group "/models/atmo" (4 members)>
                - flux: <HDF5 dataset "flux": shape (6, 7, 97055), type "<f8">
                - logg: <HDF5 dataset "logg": shape (7,), type "<f8">
                - teff: <HDF5 dataset "teff": shape (6,), type "<f8">
                - wavelength: <HDF5 dataset "wavelength": shape (97055,), 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
                - parallax: <HDF5 dataset "parallax": 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.

Reading model spectra#

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.

[7]:
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.

[8]:
readmodel.get_parameters()
[8]:
['teff', 'logg']

And also the parameter boundaries of the grid that is stored in the database.

[9]:
readmodel.get_bounds()
[9]:
{'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 (\(R_\mathrm{J}\)) 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.

[10]:
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\).

[11]:
modelbox = readmodel.get_model(model_param, spec_res=100., smooth=True)

Reading companion data#

The photometric data of PZ Tel B are also read from the database and stored in an ObjectBox.

[12]:
objectbox = database.get_object(object_name='PZ Tel B')
Getting object: PZ Tel B... [DONE]

Synthetic photometry for all filters#

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.

[13]:
synphotbox = species.multi_photometry(datatype='model',
                                      spectrum='atmo',
                                      filters=objectbox.filters,
                                      parameters=model_param)
Calculating synthetic photometry... [DONE]

Creating flux residuals#

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.

[14]:
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
Reduced chi2 = 10.05
Number of degrees of freedom = 11

Opening Box objects#

The open_box method can be used to view the content of any Box object. For example, the ModelBox contains several attributes, including the wavelengths and fluxes.

[15]:
modelbox.open_box()
Opening ModelBox...
model = atmo
type = None
wavelength = [ 0.49998877  0.50003831  0.50008785 ...  9.99834109  9.99933175
 10.0003225 ]
flux = [1.99712687e-15 2.00092575e-15 2.00491964e-15 ... 5.56173692e-17
 5.56061304e-17 5.55950840e-17]
parameters = {'teff': 2900.0, 'logg': 4.5, 'radius': 2.2, 'distance': 47.13, 'luminosity': 0.003114426265448587, 'mass': 59.04988128500266}
quantity = flux
contribution = None
bol_flux = None

Similarly, an ObjectBox contains a dictionary with the magnitudes and a dictionary with the fluxes.

[16]:
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']
mean_wavel = {'Gemini/NICI.ED286': 1.5841803431418238, 'Gemini/NIRI.H2S1v2-1-G0220': 2.2447142746110718, 'Paranal/NACO.H': 1.6588090664617747, 'Paranal/NACO.J': 1.265099894847529, 'Paranal/NACO.Ks': 2.144954491491888, 'Paranal/NACO.Lp': 3.8050282724280526, 'Paranal/NACO.Mp': 4.780970919324577, 'Paranal/NACO.NB405': 4.055862923806052, 'Paranal/SPHERE.IRDIS_D_H23_2': 1.5863509078883227, 'Paranal/SPHERE.IRDIS_D_H23_3': 1.6661442175885708, 'Paranal/SPHERE.IRDIS_D_K12_1': 2.1038552712775034, 'Paranal/SPHERE.IRDIS_D_K12_2': 2.255172356268582, 'Paranal/SPHERE.ZIMPOL_I_PRIM': 0.7843997176190827, 'Paranal/SPHERE.ZIMPOL_R_PRIM': 0.6278112553204571}
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.78256313e-14, 3.59792047e-15]), 'Gemini/NIRI.H2S1v2-1-G0220': array([1.05904381e-14, 1.36936891e-15]), 'Paranal/NACO.H': array([1.96875866e-14, 2.54565189e-15]), 'Paranal/NACO.J': array([3.11068448e-14, 5.76255332e-15]), 'Paranal/NACO.Ks': array([1.12416273e-14, 7.25276710e-16]), 'Paranal/NACO.Lp': array([2.02006237e-15, 4.12126891e-16]), 'Paranal/NACO.Mp': array([9.18778530e-16, 2.53900194e-17]), 'Paranal/NACO.NB405': array([1.67046283e-15, 1.07773346e-16]), 'Paranal/SPHERE.IRDIS_D_H23_2': array([2.54130005e-14, 4.46991835e-15]), 'Paranal/SPHERE.IRDIS_D_H23_3': array([2.42699565e-14, 4.26886718e-15]), 'Paranal/SPHERE.IRDIS_D_K12_1': array([1.15478090e-14, 9.58329872e-16]), 'Paranal/SPHERE.IRDIS_D_K12_2': array([1.14281211e-14, 1.05405764e-15]), 'Paranal/SPHERE.ZIMPOL_I_PRIM': array([1.08628801e-14, 1.20305572e-15]), 'Paranal/SPHERE.ZIMPOL_R_PRIM': array([1.82633135e-15, 5.28569079e-16])}
spectrum = None
parallax = [21.1621  0.0223]
distance = None

The attributes in a Box object can be extracted for further analysis or creating plots. For example, to extract the array with wavelengths from the ModelBox:

[17]:
modelbox.wavelength
[17]:
array([ 0.49998877,  0.50003831,  0.50008785, ...,  9.99834109,
        9.99933175, 10.0003225 ])

Plotting model spectrum and photometric fluxes#

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.

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.

[18]:
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'),
                      offset=(-0.45, -0.04),
                      legend={'loc': 'upper right', 'fontsize': 11.},
                      figsize=(7., 3.),
                      quantity='flux',
                      output=None)
Plotting spectrum... [DONE]
../_images/tutorials_data_model_47_1.png