Comparison of data and model¶
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.4.0... [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 0x14f37f190>
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 Teff range into the database.
[4]:
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
Adding companion data¶
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.
[5]:
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.
[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 "<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, 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
- 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.
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 (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.
[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
Opening Box objects¶
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.
[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
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.5841804, 'Gemini/NIRI.H2S1v2-1-G0220': 2.2447143, 'Paranal/NACO.H': 1.6588092, 'Paranal/NACO.J': 1.2650998, 'Paranal/NACO.Ks': 2.1449544, 'Paranal/NACO.Lp': 3.8050284, 'Paranal/NACO.Mp': 4.780971, 'Paranal/NACO.NB405': 4.0558634, 'Paranal/SPHERE.IRDIS_D_H23_2': 1.5863509, 'Paranal/SPHERE.IRDIS_D_H23_3': 1.6661441, 'Paranal/SPHERE.IRDIS_D_K12_1': 2.1038554, 'Paranal/SPHERE.IRDIS_D_K12_2': 2.2551725, 'Paranal/SPHERE.ZIMPOL_I_PRIM': 0.78439957, 'Paranal/SPHERE.ZIMPOL_R_PRIM': 0.62781125}
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:
[17]:
modelbox.wavelength
[17]:
array([ 0.49998877, 0.50003831, 0.50008785, ..., 9.99834109,
9.99933175, 10.0003225 ])
Plotting model spectrum and photometry¶
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.
[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'),
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.
[19]:
from IPython.display import Image
Image('spectrum.png')
[19]:
