Fitting data with a grid of model spectra

In this tutorial, we will fit spectra and photometric fluxes of beta Pic b with the synthetic spectra from a radiative-convective equilibrium atmosphere model. Here we will use DRIFT-PHOENIX but there are also several other models supported by species.

Getting started

We start by importing urllib and the species toolkit.

[1]:
import urllib.request
import species

We initiate species by running the SpeciesInit class. By doing so, both the configuration file and the HDF5 database are created in the working folder.

[2]:
species.SpeciesInit()
Initiating species v0.3.5... [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.setup.SpeciesInit at 0x10be68640>

Let’s now download the GRAVITY \(K\) band spectrum of beta Pic b that was published by Gravity Collaboration et al. 2020 and we will also make use of the GPI \(YJH\) band spectra from Chilcote et al. 2017.

[3]:
urllib.request.urlretrieve('https://home.strw.leidenuniv.nl/~stolker/species/BetaPictorisb_2018-09-22.fits',
                           'BetaPictorisb_2018-09-22.fits')
urllib.request.urlretrieve('https://home.strw.leidenuniv.nl/~stolker/species/betapicb_gpi_y.dat',
                           'betapicb_gpi_y.dat')
urllib.request.urlretrieve('https://home.strw.leidenuniv.nl/~stolker/species/betapicb_gpi_j.dat',
                           'betapicb_gpi_j.dat')
urllib.request.urlretrieve('https://home.strw.leidenuniv.nl/~stolker/species/betapicb_gpi_h.dat',
                           'betapicb_gpi_h.dat')
[3]:
('betapicb_gpi_h.dat', <http.client.HTTPMessage at 0x148840ca0>)

Adding model spectra to the database

Data will be added to the database by first creating an instance of Database.

[4]:
database = species.Database()

The Database object has a number of methods to read and write data. For adding model spectra, we use the add_model method which will download the grid of DRIFT-PHOENIX spectra to the data_folder (see configuration file) except if the data was already present. After downloading, it will import the spectra into the database. We will only require a limited range of T\(_\mathrm{eff}\) since we know the approximate temperature of beta Pic b.

[5]:
database.add_model(model='drift-phoenix', teff_range=(1500., 2000.))
Downloading DRIFT-PHOENIX model spectra (229 MB)... [DONE]
Unpacking DRIFT-PHOENIX model spectra (229 MB)... [DONE]
Adding DRIFT-PHOENIX model spectra... [DONE]
Grid points stored in the database:
   - Teff = [1500. 1600. 1700. 1800. 1900. 2000.]
   - 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: 6
   - 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
   - teff = 1900.0, logg = 4.5, feh = 0.3
   - teff = 1900.0, logg = 5.5, feh = 0.3
Number of stored grid points: 144
Number of interpolated grid points: 4
Number of missing grid points: 0
/Users/tomasstolker/applications/species/species/util/data_util.py:278: RuntimeWarning: divide by zero encountered in log10
  flux = np.log10(flux)

There were 4 grid points missing so these have been interpolated (in log-space) to make it an homogeneous grid of parameters.

Adding photometry and spectra

We will also add photometric data of beta Pic b that have been adopted in species from various references. Let’s first have a look from which directly imaged planets there is photometric data available.

[6]:
database.list_companions()
Object name = beta Pic b
Distance (pc) = 19.75 +/- 0.13
Magellan/VisAO.Ys (mag) = 15.53 +/- 0.34
Paranal/NACO.J (mag) = 14.11 +/- 0.21
Gemini/NICI.ED286 (mag) = 13.18 +/- 0.15
Paranal/NACO.H (mag) = 13.32 +/- 0.14
Paranal/NACO.Ks (mag) = 12.64 +/- 0.11
Paranal/NACO.NB374 (mag) = 11.25 +/- 0.23
Paranal/NACO.Lp (mag) = 11.3 +/- 0.06
Paranal/NACO.NB405 (mag) = 10.98 +/- 0.05
Paranal/NACO.Mp (mag) = 11.1 +/- 0.12

Object name = HIP 65426 b
Distance (pc) = 109.21 +/- 0.75
Paranal/SPHERE.IRDIS_D_H23_2 (mag) = 17.94 +/- 0.05
Paranal/SPHERE.IRDIS_D_H23_3 (mag) = 17.58 +/- 0.06
Paranal/SPHERE.IRDIS_D_K12_1 (mag) = 17.01 +/- 0.09
Paranal/SPHERE.IRDIS_D_K12_2 (mag) = 16.79 +/- 0.09
Paranal/NACO.Lp (mag) = 15.33 +/- 0.12
Paranal/NACO.NB405 (mag) = 15.23 +/- 0.22
Paranal/NACO.Mp (mag) = 14.65 +/- 0.29

Object name = 51 Eri b
Distance (pc) = 29.78 +/- 0.12
MKO/NSFCam.J (mag) = 19.04 +/- 0.4
MKO/NSFCam.H (mag) = 18.99 +/- 0.21
MKO/NSFCam.K (mag) = 18.67 +/- 0.19
Paranal/SPHERE.IRDIS_B_H (mag) = 19.45 +/- 0.29
Paranal/SPHERE.IRDIS_D_H23_2 (mag) = 18.41 +/- 0.26
Paranal/SPHERE.IRDIS_D_K12_1 (mag) = 17.55 +/- 0.14
Keck/NIRC2.Lp (mag) = 16.2 +/- 0.11
Keck/NIRC2.Ms (mag) = 16.1 +/- 0.5

Object name = HR 8799 b
Distance (pc) = 41.29 +/- 0.15
Subaru/CIAO.z (mag) = 21.22 +/- 0.29
Paranal/SPHERE.IRDIS_B_J (mag) = 19.78 +/- 0.09
Keck/NIRC2.H (mag) = 18.05 +/- 0.09
Paranal/SPHERE.IRDIS_D_H23_2 (mag) = 18.08 +/- 0.14
Paranal/SPHERE.IRDIS_D_H23_3 (mag) = 17.78 +/- 0.1
Keck/NIRC2.Ks (mag) = 17.03 +/- 0.08
Paranal/SPHERE.IRDIS_D_K12_1 (mag) = 17.15 +/- 0.06
Paranal/SPHERE.IRDIS_D_K12_2 (mag) = 16.97 +/- 0.09
Paranal/NACO.Lp (mag) = 15.52 +/- 0.1
Paranal/NACO.NB405 (mag) = 14.82 +/- 0.18
Keck/NIRC2.Ms (mag) = 16.05 +/- 0.3

Object name = HR 8799 c
Distance (pc) = 41.29 +/- 0.15
Paranal/SPHERE.IRDIS_B_J (mag) = 18.6 +/- 0.13
Keck/NIRC2.H (mag) = 17.06 +/- 0.13
Paranal/SPHERE.IRDIS_D_H23_2 (mag) = 17.09 +/- 0.12
Paranal/SPHERE.IRDIS_D_H23_3 (mag) = 16.78 +/- 0.1
Keck/NIRC2.Ks (mag) = 16.11 +/- 0.08
Paranal/SPHERE.IRDIS_D_K12_1 (mag) = 16.19 +/- 0.05
Paranal/SPHERE.IRDIS_D_K12_2 (mag) = 15.86 +/- 0.07
Paranal/NACO.Lp (mag) = 14.65 +/- 0.11
Paranal/NACO.NB405 (mag) = 13.97 +/- 0.11
Keck/NIRC2.Ms (mag) = 15.03 +/- 0.14

Object name = HR 8799 d
Distance (pc) = 41.29 +/- 0.15
Paranal/SPHERE.IRDIS_B_J (mag) = 18.59 +/- 0.37
Keck/NIRC2.H (mag) = 16.71 +/- 0.24
Paranal/SPHERE.IRDIS_D_H23_2 (mag) = 17.02 +/- 0.17
Paranal/SPHERE.IRDIS_D_H23_3 (mag) = 16.85 +/- 0.16
Keck/NIRC2.Ks (mag) = 16.09 +/- 0.12
Paranal/SPHERE.IRDIS_D_K12_1 (mag) = 16.2 +/- 0.07
Paranal/SPHERE.IRDIS_D_K12_2 (mag) = 15.84 +/- 0.1
Paranal/NACO.Lp (mag) = 14.55 +/- 0.14
Paranal/NACO.NB405 (mag) = 13.87 +/- 0.15
Keck/NIRC2.Ms (mag) = 14.65 +/- 0.35

Object name = HR 8799 e
Distance (pc) = 41.29 +/- 0.15
Paranal/SPHERE.IRDIS_B_J (mag) = 18.4 +/- 0.21
Paranal/SPHERE.IRDIS_D_H23_2 (mag) = 16.91 +/- 0.2
Paranal/SPHERE.IRDIS_D_H23_3 (mag) = 16.68 +/- 0.21
Keck/NIRC2.Ks (mag) = 15.91 +/- 0.22
Paranal/SPHERE.IRDIS_D_K12_1 (mag) = 16.12 +/- 0.1
Paranal/SPHERE.IRDIS_D_K12_2 (mag) = 15.82 +/- 0.11
Paranal/NACO.Lp (mag) = 14.49 +/- 0.21
Paranal/NACO.NB405 (mag) = 13.72 +/- 0.2

Object name = HD 95086 b
Distance (pc) = 86.44 +/- 0.24
Gemini/GPI.H (mag) = 20.51 +/- 0.25
Gemini/GPI.K1 (mag) = 18.99 +/- 0.2
Paranal/NACO.Lp (mag) = 16.27 +/- 0.19

Object name = PDS 70 b
Distance (pc) = 113.43 +/- 0.52
Paranal/SPHERE.IRDIS_D_H23_2 (mag) = 18.12 +/- 0.21
Paranal/SPHERE.IRDIS_D_H23_3 (mag) = 17.97 +/- 0.18
Paranal/SPHERE.IRDIS_D_K12_1 (mag) = 16.66 +/- 0.04
Paranal/SPHERE.IRDIS_D_K12_2 (mag) = 16.37 +/- 0.06
MKO/NSFCam.J (mag) = 20.04 +/- 0.09
MKO/NSFCam.H (mag) = 18.24 +/- 0.04
Paranal/NACO.Lp (mag) = 14.68 +/- 0.22
Paranal/NACO.NB405 (mag) = 14.68 +/- 0.27
Paranal/NACO.Mp (mag) = 13.8 +/- 0.27
Keck/NIRC2.Lp (mag) = 14.64 +/- 0.18

Object name = PDS 70 c
Distance (pc) = 113.43 +/- 0.52
Paranal/NACO.NB405 (mag) = 14.91 +/- 0.35
Keck/NIRC2.Lp (mag) = 15.5 +/- 0.46

Object name = 2M1207 b
Distance (pc) = 64.42 +/- 0.65
HST/NICMOS1.F090M (mag) = 22.58 +/- 0.35
HST/NICMOS1.F110M (mag) = 20.61 +/- 0.15
HST/NICMOS1.F145M (mag) = 19.05 +/- 0.03
HST/NICMOS1.F160W (mag) = 18.27 +/- 0.02
Paranal/NACO.J (mag) = 20.0 +/- 0.2
Paranal/NACO.H (mag) = 18.09 +/- 0.21
Paranal/NACO.Ks (mag) = 16.93 +/- 0.11
Paranal/NACO.Lp (mag) = 15.28 +/- 0.14

Object name = AB Pic B
Distance (pc) = 50.12 +/- 0.07
Paranal/NACO.J (mag) = 16.18 +/- 0.1
Paranal/NACO.H (mag) = 14.69 +/- 0.1
Paranal/NACO.Ks (mag) = 14.14 +/- 0.08

Object name = HD 206893 B
Distance (pc) = 40.81 +/- 0.11
Paranal/SPHERE.IRDIS_B_H (mag) = 16.79 +/- 0.06
Paranal/SPHERE.IRDIS_D_K12_1 (mag) = 15.2 +/- 0.1
Paranal/SPHERE.IRDIS_D_K12_2 (mag) = 14.88 +/- 0.09
Paranal/NACO.Lp (mag) = 13.79 +/- 0.31
Paranal/NACO.NB405 (mag) = 13.16 +/- 0.34
Paranal/NACO.Mp (mag) = 12.77 +/- 0.27

Object name = RZ Psc B
Distance (pc) = 195.86 +/- 4.03
Paranal/SPHERE.IRDIS_B_H (mag) = (13.71, 0.14) +/- (13.85, 0.26)
Paranal/SPHERE.IRDIS_B_Ks (mag) = 13.51 +/- 0.2

Object name = GQ Lup B
Distance (pc) = 151.82 +/- 1.1
HST/WFPC2-PC.F606W (mag) = 19.19 +/- 0.07
HST/WFPC2-PC.F814W (mag) = 17.67 +/- 0.05
HST/NICMOS2.F171M (mag) = 13.84 +/- 0.13
HST/NICMOS2.F190N (mag) = 14.08 +/- 0.2
HST/NICMOS2.F215N (mag) = 13.4 +/- 0.15
Magellan/VisAO.ip (mag) = 18.89 +/- 0.24
Magellan/VisAO.zp (mag) = 16.4 +/- 0.1
Magellan/VisAO.Ys (mag) = 15.88 +/- 0.1
Paranal/NACO.Ks (mag) = (13.474, 0.031) +/- (13.386, 0.032)
Subaru/CIAO.CH4s (mag) = 13.76 +/- 0.26
Subaru/CIAO.K (mag) = 13.37 +/- 0.12
Subaru/CIAO.Lp (mag) = 12.44 +/- 0.22

Object name = PZ Tel B
Distance (pc) = 47.13 +/- 0.13
Paranal/SPHERE.ZIMPOL_R_PRIM (mag) = 17.84 +/- 0.31
Paranal/SPHERE.ZIMPOL_I_PRIM (mag) = 15.16 +/- 0.12
Paranal/SPHERE.IRDIS_D_H23_2 (mag) = 11.78 +/- 0.19
Paranal/SPHERE.IRDIS_D_H23_3 (mag) = 11.65 +/- 0.19
Paranal/SPHERE.IRDIS_D_K12_1 (mag) = 11.56 +/- 0.09
Paranal/SPHERE.IRDIS_D_K12_2 (mag) = 11.29 +/- 0.1
Paranal/NACO.J (mag) = 12.47 +/- 0.2
Paranal/NACO.H (mag) = 11.93 +/- 0.14
Paranal/NACO.Ks (mag) = 11.53 +/- 0.07
Paranal/NACO.Lp (mag) = 11.04 +/- 0.22
Paranal/NACO.NB405 (mag) = 10.94 +/- 0.07
Paranal/NACO.Mp (mag) = 10.93 +/- 0.03
Gemini/NICI.ED286 (mag) = 11.68 +/- 0.14
Gemini/NIRI.H2S1v2-1-G0220 (mag) = 11.39 +/- 0.14

Object name = kappa And b
Distance (pc) = 50.06 +/- 0.87
Subaru/CIAO.J (mag) = 15.86 +/- 0.21
Subaru/CIAO.H (mag) = 14.95 +/- 0.13
Subaru/CIAO.Ks (mag) = 14.32 +/- 0.09
Keck/NIRC2.Lp (mag) = 13.12 +/- 0.1
LBT/LMIRCam.M_77K (mag) = 13.3 +/- 0.3

Object name = HD 1160 B
Distance (pc) = 125.9 +/- 1.2
MKO/NSFCam.J (mag) = 14.69 +/- 0.05
MKO/NSFCam.H (mag) = 14.21 +/- 0.02
MKO/NSFCam.Ks (mag) = 14.12 +/- 0.05
Paranal/NACO.Lp (mag) = 13.6 +/- 0.1
Keck/NIRC2.Ms (mag) = 13.81 +/- 0.24

Object name = ROXs 42 Bb
Distance (pc) = 144.16 +/- 1.53
Keck/NIRC2.J (mag) = 16.91 +/- 0.11
Keck/NIRC2.H (mag) = 15.88 +/- 0.05
Keck/NIRC2.Ks (mag) = 15.01 +/- 0.06
Keck/NIRC2.Lp (mag) = 13.97 +/- 0.06
Keck/NIRC2.Ms (mag) = 14.01 +/- 0.23

Object name = GJ 504 b
Distance (pc) = 17.54 +/- 0.08
Paranal/SPHERE.IRDIS_D_Y23_2 (mag) = 20.98 +/- 0.2
Paranal/SPHERE.IRDIS_D_Y23_3 (mag) = 20.14 +/- 0.09
Paranal/SPHERE.IRDIS_D_J23_3 (mag) = 19.01 +/- 0.17
Paranal/SPHERE.IRDIS_D_H23_2 (mag) = 18.95 +/- 0.3
Paranal/SPHERE.IRDIS_D_H23_3 (mag) = 21.81 +/- 0.35
Paranal/SPHERE.IRDIS_D_K12_1 (mag) = 18.77 +/- 0.2
Subaru/CIAO.J (mag) = 19.78 +/- 0.1
Subaru/CIAO.H (mag) = 20.01 +/- 0.14
Subaru/CIAO.Ks (mag) = 19.38 +/- 0.11
Subaru/CIAO.CH4s (mag) = 19.58 +/- 0.13
Subaru/IRCS.Lp (mag) = 16.7 +/- 0.17

Object name = GU Psc b
Distance (pc) = 47.61 +/- 0.16
Gemini/GMOS-S.z (mag) = 21.75 +/- 0.07
CFHT/Wircam.Y (mag) = 19.4 +/- 0.05
CFHT/Wircam.J (mag) = 18.12 +/- 0.03
CFHT/Wircam.H (mag) = 17.7 +/- 0.03
CFHT/Wircam.Ks (mag) = 17.4 +/- 0.03
WISE/WISE.W1 (mag) = 17.17 +/- 0.33
WISE/WISE.W2 (mag) = 15.41 +/- 0.22

Object name = 2M0103 ABb
Distance (pc) = 47.2 +/- 3.1
Paranal/NACO.J (mag) = 15.47 +/- 0.3
Paranal/NACO.H (mag) = 14.27 +/- 0.2
Paranal/NACO.Ks (mag) = 13.67 +/- 0.2
Paranal/NACO.Lp (mag) = 12.67 +/- 0.1

Object name = 1RXS 1609 B
Distance (pc) = 139.67 +/- 1.33
Gemini/NIRI.J-G0202w (mag) = 17.9 +/- 0.12
Gemini/NIRI.H-G0203w (mag) = 16.87 +/- 0.07
Gemini/NIRI.K-G0204w (mag) = 16.17 +/- 0.18
Gemini/NIRI.Lprime-G0207w (mag) = 14.8 +/- 0.3

Object name = GSC 06214 B
Distance (pc) = 108.84 +/- 0.51
MKO/NSFCam.J (mag) = 16.24 +/- 0.04
MKO/NSFCam.H (mag) = 15.55 +/- 0.04
MKO/NSFCam.Kp (mag) = 14.95 +/- 0.05
MKO/NSFCam.Lp (mag) = 13.75 +/- 0.07
LBT/LMIRCam.M_77K (mag) = 13.75 +/- 0.3

Object name = HD 72946 B
Distance (pc) = 25.87 +/- 0.03
Paranal/SPHERE.IRDIS_D_H23_2 (mag) = 14.56 +/- 0.07
Paranal/SPHERE.IRDIS_D_H23_3 (mag) = 14.4 +/- 0.07

Object name = HIP 64892 B
Distance (pc) = 125.2 +/- 1.42
Paranal/SPHERE.IRDIS_D_H23_2 (mag) = 14.21 +/- 0.17
Paranal/SPHERE.IRDIS_D_H23_3 (mag) = 13.94 +/- 0.17
Paranal/SPHERE.IRDIS_D_K12_1 (mag) = 13.77 +/- 0.17
Paranal/SPHERE.IRDIS_D_K12_2 (mag) = 13.45 +/- 0.19
Paranal/NACO.Lp (mag) = 13.09 +/- 0.17

Object name = TYC 8988 b
Distance (pc) = 94.6 +/- 0.3
Paranal/SPHERE.IRDIS_D_Y23_2 (mag) = 17.03 +/- 0.21
Paranal/SPHERE.IRDIS_D_Y23_3 (mag) = 16.67 +/- 0.16
Paranal/SPHERE.IRDIS_D_J23_2 (mag) = 16.27 +/- 0.08
Paranal/SPHERE.IRDIS_D_J23_3 (mag) = 15.73 +/- 0.07
Paranal/SPHERE.IRDIS_D_H23_2 (mag) = 15.11 +/- 0.08
Paranal/SPHERE.IRDIS_D_H23_3 (mag) = 14.78 +/- 0.07
Paranal/SPHERE.IRDIS_D_K12_1 (mag) = 14.44 +/- 0.04
Paranal/SPHERE.IRDIS_D_K12_2 (mag) = 14.07 +/- 0.04
Paranal/SPHERE.IRDIS_B_J (mag) = 15.73 +/- 0.38
Paranal/SPHERE.IRDIS_B_H (mag) = 15.87 +/- 0.38
Paranal/SPHERE.IRDIS_B_Ks (mag) = 14.7 +/- 0.14
Paranal/NACO.Lp (mag) = 13.3 +/- 0.08
Paranal/NACO.Mp (mag) = 13.08 +/- 0.2

Object name = TYC 8988 c
Distance (pc) = 94.6 +/- 0.3
Paranal/SPHERE.IRDIS_D_Y23_3 (mag) = 22.37 +/- 0.31
Paranal/SPHERE.IRDIS_D_J23_2 (mag) = 21.81 +/- 0.22
Paranal/SPHERE.IRDIS_D_J23_3 (mag) = 21.17 +/- 0.15
Paranal/SPHERE.IRDIS_D_H23_2 (mag) = 19.78 +/- 0.08
Paranal/SPHERE.IRDIS_D_H23_3 (mag) = 19.32 +/- 0.06
Paranal/SPHERE.IRDIS_D_K12_1 (mag) = 18.34 +/- 0.04
Paranal/SPHERE.IRDIS_D_K12_2 (mag) = 17.85 +/- 0.09
Paranal/SPHERE.IRDIS_B_H (mag) = 19.69 +/- 0.23
Paranal/NACO.Lp (mag) = 16.29 +/- 0.21

To add the photometric data and distance of beta Pic b to the database, we use the add_companion method, which will also convert the magnitudes into fluxes with a flux-calibrated spectrum of Vega and the filter profiles (downloaded from the SVO Filter Profile Service). It is also possible to manually add photometric data with the app_mag and flux_density arguments of add_object.

[7]:
database.add_companion(name='beta Pic b')
Downloading Vega spectrum (270 kB)... [DONE]
Adding Vega spectrum... [DONE]
Adding filter: Magellan/VisAO.Ys... [DONE]
Adding filter: Paranal/NACO.J... [DONE]
Adding filter: Gemini/NICI.ED286... [DONE]
Adding filter: Paranal/NACO.H... [DONE]
Adding filter: Paranal/NACO.Ks... [DONE]
Adding filter: Paranal/NACO.NB374... [DONE]
Adding filter: Paranal/NACO.Lp... [DONE]
Adding filter: Paranal/NACO.NB405... [DONE]
Adding filter: Paranal/NACO.Mp... [DONE]
Adding object: beta Pic b
   - Distance (pc) = 19.75 +/- 0.13
   - Magellan/VisAO.Ys:
      - Apparent magnitude = 15.53 +/- 0.34
      - Flux (W m-2 um-1) = 4.27e-15 +/- 1.36e-15
   - Paranal/NACO.J:
      - Apparent magnitude = 14.11 +/- 0.21
      - Flux (W m-2 um-1) = 6.87e-15 +/- 1.34e-15
   - Gemini/NICI.ED286:
      - Apparent magnitude = 13.18 +/- 0.15
      - Flux (W m-2 um-1) = 6.99e-15 +/- 9.69e-16
   - Paranal/NACO.H:
      - Apparent magnitude = 13.32 +/- 0.14
      - Flux (W m-2 um-1) = 5.47e-15 +/- 7.08e-16
   - Paranal/NACO.Ks:
      - Apparent magnitude = 12.64 +/- 0.11
      - Flux (W m-2 um-1) = 4.04e-15 +/- 4.10e-16
   - Paranal/NACO.NB374:
      - Apparent magnitude = 11.25 +/- 0.23
      - Flux (W m-2 um-1) = 1.69e-15 +/- 3.61e-16
   - Paranal/NACO.Lp:
      - Apparent magnitude = 11.30 +/- 0.06
      - Flux (W m-2 um-1) = 1.59e-15 +/- 8.79e-17
   - Paranal/NACO.NB405:
      - Apparent magnitude = 10.98 +/- 0.05
      - Flux (W m-2 um-1) = 1.61e-15 +/- 7.42e-17
   - Paranal/NACO.Mp:
      - Apparent magnitude = 11.10 +/- 0.12
      - Flux (W m-2 um-1) = 7.86e-16 +/- 8.70e-17

Next, we will add the spectra with the add_object method. Make sure to use the same database name as argument of object_name (i.e. beta Pic b in this case) such that the spectra are added to the object name for which previously photometric data had been stored with add_companion.

The argument of spectrum is a dictionary with the names of the spectra as key and a tuple with three elements as value. The first element is a text or FITS file with the spectrum, the second element is the optional covariance matrix, and the third value is the mandatory spectral resolution. The latter is used to smooth the model spectra to the instrument resolution.

For this example, we have adopted the VLTI/GRAVITY \(K\) band spectrum from Gravity Collaboration et al. 2020 and the Gemini/GPI \(YJH\) band spectra from Chilcote et al. 2017. Data that have been stored in the FITS format by the ExoGRAVITY collaboration is automatically recognized so the same FITS filename can be provided as first and second element in the dictionary.

[8]:
database.add_object('beta Pic b',
                    distance=None,
                    app_mag=None,
                    flux_density=None,
                    spectrum={'GRAVITY': ('BetaPictorisb_2018-09-22.fits', 'BetaPictorisb_2018-09-22.fits', 500.),
                              'GPI\\_Y': ('betapicb_gpi_y.dat', None, 40.),
                              'GPI\\_J': ('betapicb_gpi_j.dat', None, 40.),
                              'GPI\\_H': ('betapicb_gpi_h.dat', None, 40.)},
                              deredden=None)
Adding object: beta Pic b
   - GRAVITY spectrum:
      - Object: Unknown
      - Database tag: GRAVITY
      - Filename: BetaPictorisb_2018-09-22.fits
      - Data shape: (237, 3)
      - Wavelength range (um): 1.97 - 2.49
      - Mean flux (W m-2 um-1): 4.65e-15
      - Mean error (W m-2 um-1): 1.00e-16
   - Spectrum:
      - Database tag: GPI\_Y
      - Filename: betapicb_gpi_y.dat
      - Data shape: (29, 3)
      - Wavelength range (um): 0.98 - 1.13
      - Mean flux (W m-2 um-1): 5.40e-15
      - Mean error (W m-2 um-1): 1.87e-15
   - Spectrum:
      - Database tag: GPI\_J
      - Filename: betapicb_gpi_j.dat
      - Data shape: (32, 3)
      - Wavelength range (um): 1.13 - 1.34
      - Mean flux (W m-2 um-1): 6.84e-15
      - Mean error (W m-2 um-1): 4.79e-16
   - Spectrum:
      - Database tag: GPI\_H
      - Filename: betapicb_gpi_h.dat
      - Data shape: (34, 3)
      - Wavelength range (um): 1.51 - 1.79
      - Mean flux (W m-2 um-1): 5.63e-15
      - Mean error (W m-2 um-1): 3.86e-16
   - GRAVITY covariance matrix:
      - Object: Unknown
      - Database tag: GRAVITY
      - Filename: BetaPictorisb_2018-09-22.fits
      - Data shape: (237, 237)
   - Spectral resolution:
      - GRAVITY: 500.0
      - GPI\_Y: 40.0
      - GPI\_J: 40.0
      - GPI\_H: 40.0

The use of add_object is incremental so when rerunning the method it will add new data in case a filter or spectrum name hadn’t been used before, or it will overwrite existing data in case a filter or spectrum name does already exist. To remove existing data of an object, it is easiest to remove the group of data entirely from the database with the delete_data method of Database (e.g. delete_data('objects/beta Pic b') and rerun add_companion and/or add_object afterwards.

Bayesian inference with nested sampling

To fit the data with the grid of model spectra (i.e. grid retrieval), we start by creating an instance of FitModel, which requires the database tag of the object (beta Pic b) and the atmosphere model (drift-phoenix).

The dictionary of the bounds parameter contains the prior boundaries for the various parameters. By default, the parameters of the atmosphere model are set to the full range of the available spectra (e.g. \(\log(g)\) and \(\mathrm{[Fe/H]}\) in the example below).

In addition to the parameters of the atmosphere model, several optional parameters can be fitted, for example to account for calibration systematics, extinction, and excess emission from a disk. More details on the various parameters can be found in the API documentation of FitModel. As example, we fit a scaling parameter for the GPI \(H\) band spectrum.

The arguments of inc_phot and inc_spec are either a boolean or a list of photometry and spectra that were previously stored in the database for the object beta Pic b. We will use all the spectra but only the photometric data in the \(L\) and \(M\) bands.

To account for correlated noise, it is possible to adopt the approach by Wang et al. (2020) to estimate the covariances by using a Gaussian process. This will add two additional parameters per spectrum (correlation length and amplitude) but is not mandatory. For the example, we only select the GPI \(H\) band spectrum, to limit the computation time.

For simplicity we create a list with the 4 filter names that we want to use.

[9]:
inc_phot = ['Paranal/NACO.Lp', 'Paranal/NACO.NB374', 'Paranal/NACO.NB405', 'Paranal/NACO.Mp']

We create now the instance of FitModel. By doing so, the grid of model spectra will be interpolated to the wavelengths of the data.

[10]:
fit = species.FitModel(object_name='beta Pic b',
                       model='drift-phoenix',
                       bounds={'teff': (1500., 2000.),
                               'radius': (0.5, 2.),
                               'GPI\\_H': ((0.5, 1.5), None)},
                       inc_phot=inc_phot,
                       inc_spec=True,
                       fit_corr=['GPI\\_H'],
                       weights=None)
Getting object: beta Pic b... [DONE]
Interpolating Paranal/NACO.Lp... [DONE]
Interpolating Paranal/NACO.NB374... [DONE]
Interpolating Paranal/NACO.NB405... [DONE]
Interpolating Paranal/NACO.Mp... [DONE]
Interpolating GPI\_H... [DONE]
Interpolating GPI\_J... [DONE]
Interpolating GPI\_Y... [DONE]
Interpolating GRAVITY... [DONE]
Fitting 7 parameters:
   - teff
   - logg
   - feh
   - radius
   - corr_len_GPI\_H
   - corr_amp_GPI\_H
   - scaling_GPI\_H
Prior boundaries:
   - teff = (1500.0, 2000.0)
   - radius = (0.5, 2.0)
   - logg = (3.0, 5.5)
   - feh = (-0.6, 0.3)
   - corr_len_GPI\_H = (-3.0, 0.0)
   - corr_amp_GPI\_H = (0.0, 1.0)
   - scaling_GPI\_H = (0.5, 1.5)
Weights for the log-likelihood function:
   - GPI\_H = 1.00e+00
   - GPI\_J = 1.00e+00
   - GPI\_Y = 1.00e+00
   - GRAVITY = 1.00e+00
   - Paranal/NACO.Lp = 1.00e+00
   - Paranal/NACO.NB374 = 1.00e+00
   - Paranal/NACO.NB405 = 1.00e+00
   - Paranal/NACO.Mp = 1.00e+00

We are now ready to estimate the posterior distributions by using either MultiNest, UltraNest, or emcee. It is recommended to use MultiNest (see run_multinest method) or UltraNest (see run_ultranest method) since these nested sampling algorithms are more powerful in sampling multi-modal distributions. They will also estimate the marginalized likelihood (i.e. model evidence) which enables pair-wise model comparison through the Bayes factor.

In this example, we will use the run_ultranest method since UltraNest is somewhat more straightforward to install than MultiNest. We specify the database tag were the posterior samples will be stored, the minimum number of live points (see UltraNest documentation), the output folder where UltraNest will write its data, and we adopt the mass of beta Pic b from Gravity Collaboration et al. 2020 and use it as a Gaussian prior on the mass (which propagates into \(\log(g)\)).

Depending on the number of parameters, the posterior sampling can be computationally expensive. The runtime of the nested sampling part of this tutorial (i.e. the run_ultranest method below) was about 15 min on a MacBook Pro from 2020. For a first attempt, one could also reduce the number of parameters by setting fit_corr and the scaling parameter to None. Using MultiNest with a small number of live points (e.g. n_live_points=500) might also be faster.

To speed up the computation, it is possible to run the nested sampling in parallel (e.g. with mpirun) to benefit from the multiprocessing support by UltraNest and MultiNest. In that case it is important that any functions of species that will write to the Database will be commented out since simultaneous writing to the HDF5 database by different processes is not possible. It is therefore recommended to first add all the required data to the database and then only run SpeciesInit, FitModel, and the sampler (run_multinest or run_ultranest) in parallel with MPI.

[11]:
fit.run_ultranest(tag='betapic',
                  min_num_live_points=400,
                  output='ultranest',
                  prior={'mass': (9., 1.6)})
Running nested sampling with UltraNest...
Creating directory for new run ultranest/run1
[ultranest] Sampling 400 live points from prior ...
/Users/tomasstolker/.pyenv/versions/3.8.7/envs/general3.8/lib/python3.8/site-packages/ultranest/store.py:194: DeprecationWarning: `np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  'points', dtype=np.float,
[ultranest] Explored until L=1e+04  1486.70 [11486.2247..11486.2255]*| it/evals=11377/450358 eff=2.5285% N=400
[ultranest] Likelihood function evaluations: 450397
/Users/tomasstolker/.pyenv/versions/3.8.7/envs/general3.8/lib/python3.8/site-packages/ultranest/utils.py:170: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  idx = np.zeros(N, dtype=np.int)
[ultranest] Writing samples and results to disk ...
/Users/tomasstolker/.pyenv/versions/3.8.7/envs/general3.8/lib/python3.8/site-packages/ultranest/utils.py:170: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  idx = np.zeros(N, dtype=np.int)
[ultranest] Writing samples and results to disk ... done
[ultranest]   logZ = 1.146e+04 +- 0.194
[ultranest] Effective samples strategy satisfied (ESS = 2812.3, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.04 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.46, need <0.5)
[ultranest]   logZ error budget: single: 0.23 bs:0.19 tail:0.01 total:0.19 required:<0.50
[ultranest] done iterating.
Log-evidence = 11462.52 +/- 0.27
Best-fit parameters (mean +/- std):
   - teff = 1.72e+03 +/- 4.97e+00
   - logg = 3.98e+00 +/- 4.81e-02
   - feh = 1.33e-01 +/- 1.39e-02
   - radius = 1.50e+00 +/- 1.19e-02
   - corr_len_GPI\_H = -1.49e+00 +/- 1.83e-01
   - corr_amp_GPI\_H = 8.03e-01 +/- 1.05e-01
   - scaling_GPI\_H = 1.12e+00 +/- 3.01e-02
Maximum likelihood sample:
   - Log-likelihood = 11486.73
   - teff = 1720.59
   - logg = 3.98
   - feh = 0.13
   - radius = 1.50
   - corr_len_GPI\_H = -1.49
   - corr_amp_GPI\_H = 0.85
   - scaling_GPI\_H = 1.12
Integrated autocorrelation time = [0.98751055]

Plotting the posterior samples

The samples from the parameter estimation have been stored in the database. We can now run the plot_posterior function to plot the 1D and 2D projections of the posterior distributions by making use of corner.py. Here we specify the database tag with the results from run_ultranest and we also include the derived posterior for the bolometric luminosity and mass.

[12]:
species.plot_posterior(tag='betapic',
                       offset=(-0.3 , -0.3),
                       title_fmt=['.0f', '.2f', '.2f', '.2f', '.2f', '.2f', '.2f', '.3f', '.1f'],
                       inc_luminosity=True,
                       inc_mass=True,
                       output='posterior.png')
Median sample:
   - teff = 1720.86
   - logg = 3.98
   - feh = 0.13
   - radius = 1.50
   - corr_len_GPI\_H = -1.47
   - corr_amp_GPI\_H = 0.82
   - scaling_GPI\_H = 1.13
   - distance = 19.75
Plotting the posterior: posterior.png... [DONE]

Let’s have a look at the plot!

[13]:
from IPython.display import Image
Image('posterior.png')
[13]:
../_images/tutorials_fitting_model_spectra_39_0.png

Extracting boxes with results

We will now extract Box objects with the results. These will be combined in a plot at the end of the tutorial. We start by selecting 30 random samples from the posterior and use these parameters to interpolate the DRIFT-PHOENIX grid. To do so, we use the get_mcmc_spectra method of the Database. We will smooth the spectra to a resolution of \(R=500\) to match with the resolution of the GRAVITY spectrum.

[14]:
samples = database.get_mcmc_spectra(tag='betapic',
                                    random=30,
                                    wavel_range=None,
                                    spec_res=500.)
Getting MCMC spectra: 100%|██████████| 30/30 [00:00<00:00, 146.53it/s]

The function returned a list with 30 ModelBox objects. Let’s have a look at the content of the first box by running the open_box method. This method can be used for all boxes to easily inspect the attributes of a Box object.

[15]:
samples[0].open_box()
Opening ModelBox...
model = drift-phoenix
type = None
wavelength = [ 0.1         0.10002476  0.10004953 ... 49.97524871 49.98762282
 50.        ]
flux = [2.42584190e-69 2.23707323e-69 1.99404456e-69 ... 2.06877619e-19
 2.06600127e-19 2.06398308e-19]
parameters = {'teff': 1722.4622616917563, 'logg': 3.971918644296093, 'feh': 0.1319168469622567, 'radius': 1.5023336843398334, 'distance': 19.75, 'mass': 8.162536214856585, 'luminosity': 0.00018074753024335536}
quantity = flux

Instead of using the plotting functionalities of species, the attributes from a Box can also be directly extracted for further processing. For example, to print the wavelengths of the first ModelBox that was returned by get_mcmc_spectra:

[16]:
print(samples[0].wavelength)
[ 0.1         0.10002476  0.10004953 ... 49.97524871 49.98762282
 50.        ]

Next, we get the median parameter values with the get_median_sample method of the Database and adopt these as best-fit parameters. Similarly, the get_probable_sample method can be used for extracting the sample with the maximum likelihood.

[17]:
best = database.get_median_sample(tag='betapic')

The best-fit model spectrum is now extracted from the DRIFT-PHOENIX grid by using the functionalities of the ReadModel class.

[18]:
read_model = species.ReadModel(model='drift-phoenix',
                               wavel_range=None)

The grid is interpolated at the best-fit parameters with get_model. The argument of model_param contains the parameter dictionary that was returned by get_median_sample. The spectrum will be smoothed to a resolution of \(R=500\). To extract the spectrum without smoothing at the resolution as stored in the database is achieved by setting spec_res=None and smooth=False.

[19]:
modelbox = read_model.get_model(model_param=best,
                                spec_res=500.,
                                smooth=True)
/Users/tomasstolker/applications/species/species/read/read_model.py:595: UserWarning: The 'corr_len_GPI\_H' parameter is not required by 'drift-phoenix' so the parameter will be ignored. The mandatory parameters are ['teff', 'logg', 'feh'].
  warnings.warn(f'The \'{key}\' parameter is not required by \'{self.model}\' so '
/Users/tomasstolker/applications/species/species/read/read_model.py:595: UserWarning: The 'corr_amp_GPI\_H' parameter is not required by 'drift-phoenix' so the parameter will be ignored. The mandatory parameters are ['teff', 'logg', 'feh'].
  warnings.warn(f'The \'{key}\' parameter is not required by \'{self.model}\' so '
/Users/tomasstolker/applications/species/species/read/read_model.py:595: UserWarning: The 'scaling_GPI\_H' parameter is not required by 'drift-phoenix' so the parameter will be ignored. The mandatory parameters are ['teff', 'logg', 'feh'].
  warnings.warn(f'The \'{key}\' parameter is not required by \'{self.model}\' so '

Some warnings are printed because the model_param dictionary contained calibration parameters (that were fitted), but these are not required when reading a model spectrum.

We will also extract all spectra and photometric fluxes of beta Pic b with the get_object method which will return the output in an ObjectBox.

[20]:
objectbox = database.get_object(object_name='beta Pic b',
                                inc_phot=True,
                                inc_spec=True)
Getting object: beta Pic b... [DONE]

We then run the ObjectBox through the update_spectra function with the best-fit parameters. This function applies the flux scaling and could also inflate the errors in case the model_param dictionary contains error inflation parameters.

[21]:
objectbox = species.update_spectra(objectbox=objectbox,
                                   model_param=best)
Scaling the flux of GPI\_H: 1.13... [DONE]

Next, we calculate the residuals of the best-fit model with get_residuals. This function will interpolate the grid at the specified parameters, calculate synthetic photometry for the filters that are provided as argument of inc_phot, and smooth and resample the spectra to the resolution and wavelengths of the spectra that are provided as argument of inc_spec (all spectra of the ObjectBox in this case). The returned residuals will be stored in a ResidualsBox.

[22]:
residuals = species.get_residuals(datatype='model',
                                  spectrum='drift-phoenix',
                                  parameters=best,
                                  objectbox=objectbox,
                                  inc_phot=inc_phot,
                                  inc_spec=True)
Calculating synthetic photometry... [DONE]
Calculating residuals... [DONE]
Residuals (sigma):
   - Paranal/NACO.Lp: -4.69
   - Paranal/NACO.NB374: -1.11
   - Paranal/NACO.NB405: -0.70
   - Paranal/NACO.Mp: -0.78
   - GPI\_H: min: -1.26, max: 2.03
   - GPI\_J: min: -1.94, max: 3.95
   - GPI\_Y: min: -1.37, max: 1.37
   - GRAVITY: min: -5.69, max: 5.14

Finally, we compute synthetic photometry from the best-fit model spectrum with the multi_photometry function. Since we will plot all available photometric data of the ObjectBox, we will calculate synthetic photometry for all the names of the filters attribute of the ObjectBox.

[23]:
synphot = species.multi_photometry(datatype='model',
                                   spectrum='drift-phoenix',
                                   filters=objectbox.filters,
                                   parameters=best)
Calculating synthetic photometry... [DONE]

Let’s have a look inside the SynphotBox.

[24]:
synphot.open_box()
Opening SynphotBox...
name = synphot
flux = {'Gemini/NICI.ED286': 6.324736978089229e-15, 'Magellan/VisAO.Ys': 3.254495050981177e-15, 'Paranal/NACO.H': 6.032863101777227e-15, 'Paranal/NACO.J': 6.3948954907705766e-15, 'Paranal/NACO.Ks': 4.983635923266107e-15, 'Paranal/NACO.Lp': 2.0025302139642507e-15, 'Paranal/NACO.Mp': 8.533908419904405e-16, 'Paranal/NACO.NB374': 2.0928879135235758e-15, 'Paranal/NACO.NB405': 1.6616662765042234e-15}

Plotting the spectral energy distribution

Now that we have gather all the Box objects with the results, we can pass them as list to the boxes parameter of the plot_spectrum function. The ResidualsBox is separately provided as argument of residuals. We also include a list of filter names for which the transmission profiles will be plotted. The arguments of residuals and filters can also be set to None to not include these data in the plot.

The somewhat complex part of plot_spectrum is the optional plot_kwargs parameter. The argument is a list with the same length as the list of boxes. Each item contains a dictionary with keyword arguments that can be used to fine-tune the styling of the plot. The item of the SynphotBox can be set to None because it is automatically adjusted to the styles that are used for the ObjectBox.

[25]:
species.plot_spectrum(boxes=[samples, modelbox, objectbox, synphot],
                      filters=objectbox.filters,
                      residuals=residuals,
                      plot_kwargs=[{'ls': '-', 'lw': 0.2, 'color': 'gray'},
                                   {'ls': '-', 'lw': 1., 'color': 'black'},
                                   {'GRAVITY': {'marker': 'o', 'ms': 5., 'mew': 0., 'color': 'seagreen', 'ls': 'none', 'alpha': 0.5, 'label': 'VLTI/GRAVITY'},
                                    'GPI\\_Y': {'marker': 'o', 'ms': 5., 'mew': 0., 'color': 'darkslateblue', 'ls': 'none', 'alpha': 0.5, 'label': 'Gemini/GPI'},
                                    'GPI\\_J': {'marker': 'o', 'ms': 5., 'mew': 0., 'color': 'darkslateblue', 'ls': 'none', 'alpha': 0.5},
                                    'GPI\\_H': {'marker': 'o', 'ms': 5., 'mew': 0., 'color': 'darkslateblue', 'ls': 'none', 'alpha': 0.5},
                                    'Magellan/VisAO.Ys': {'marker': 's', 'ms': 5., 'color': 'goldenrod', 'ls': 'none', 'label': 'Magellan/VisAO'},
                                    'Gemini/NICI.ED286': {'marker': 's', 'ms': 5., 'color': 'teal', 'ls': 'none', 'label': 'Gemini/NICI'},
                                    'Paranal/NACO.J': {'marker': 's', 'ms': 5., 'color': 'tomato', 'ls': 'none', 'label': 'VLT/NACO'},
                                    'Paranal/NACO.H': {'marker': 's', 'ms': 5., 'color': 'tomato', 'ls': 'none'},
                                    'Paranal/NACO.Ks': {'marker': 's', 'ms': 5., 'color': 'tomato', 'ls': 'none'},
                                    'Paranal/NACO.NB374': {'marker': 's', 'ms': 5., 'color': 'tomato', 'ls': 'none'},
                                    'Paranal/NACO.Lp': {'marker': 's', 'ms': 5., 'color': 'tomato', 'ls': 'none'},
                                    'Paranal/NACO.NB405': {'marker': 's', 'markersize': 5., 'color': 'tomato', 'ls': 'none'},
                                    'Paranal/NACO.Mp': {'marker': 's', 'markersize': 5., 'color': 'tomato', 'ls': 'none'}},
                                    None],
                      xlim=(0.8, 5.2),
                      ylim=(-1.15e-15, 1.15e-14),
                      ylim_res=(-7., 7.),
                      scale=('linear', 'linear'),
                      offset=(-0.4, -0.05),
                      legend=[{'loc': 'lower left', 'frameon': False, 'fontsize': 11.},
                              {'loc': 'upper right', 'frameon': False, 'fontsize': 11.}],
                      figsize=(8., 4.),
                      quantity='flux density',
                      output='spectrum.png')
Plotting spectrum: spectrum.png... [DONE]

Let’s have a look at the result!

[26]:
Image('spectrum.png')
[26]:
../_images/tutorials_fitting_model_spectra_68_0.png