Flux calibration#
Photometric and spectroscopic measurements of directly imaged planets typically provide the flux contrast between the companion and star. To calibrate the contrast of the companion to a flux or apparent magnitude requires an absolute measurement of the stellar flux.
In this tutorial, we will we will fit the Gaia DR3 and 2MASS magnitudes of the G7 type star PZ Tel with the BT-NextGen model spectra. From the posterior samples, we will then compute synthetic photometry for the VLT/ERIS M’ filter and a synthetic spectrum for a given instrument resolution and wavelength binning.
Getting started#
We start by importing the required Python modules.
[29]:
import calistar
import numpy as np
[30]:
from species import SpeciesInit
from species.data.database import Database
from species.fit.fit_model import FitModel
from species.read.read_model import ReadModel
from species.plot.plot_mcmc import plot_posterior, plot_mag_posterior
from species.plot.plot_spectrum import plot_spectrum
from species.util.box_util import update_objectbox
from species.util.fit_util import get_residuals, multi_photometry
Next, we initiate the workflow by calling the SpeciesInit class. This will create both the HDF5 database and the configuration file in the working folder.
[31]:
SpeciesInit()
==============
species v0.7.4
==============
Working folder: /Users/tomasstolker/applications/species/docs/tutorials
Creating species_config.ini... [DONE]
Creating species_database.hdf5... [DONE]
Creating data folder... [DONE]
Configuration settings:
- Database: /Users/tomasstolker/applications/species/docs/tutorials/species_database.hdf5
- Data folder: /Users/tomasstolker/applications/species/docs/tutorials/data
- Magnitude of Vega: 0.03
Multiprocessing: mpi4py installed
Process number 1 out of 1...
[31]:
<species.core.species_init.SpeciesInit at 0x15e59f950>
We then create an instance of Database, which provides read and write access to the HDF5 database.
[32]:
database = Database()
We will use the calistar tool to retrieve the Gaia DR3 parallax of PZ Tel and the Gaia, 2MASS, and WISE magnitudes. Similar to species, the calistar tool uses also the filter names as defined by the SVO Filter Profile Service. We start by creating an instance of the
CaliStar class for which we provide the Gaia DR3 source ID of PZ Tel as input (see Simbad).
[33]:
cal_star = calistar.CaliStar(gaia_source=6655168686921108864, gaia_release='DR3')
===============
calistar v0.0.3
===============
Next, we run the target_star to retrieve the astrometric and photometric properties of PZ Tel, which are returned in a dictionary.
[34]:
target_dict = cal_star.target_star(write_json=False)
-> Querying GAIA DR3...
INFO: Query finished. [astroquery.utils.tap.core]
GAIA DR3 source ID = 6655168686921108864
Reference epoch = 2016.0
Parallax = 21.16 +/- 0.02 mas
RA = 283.274586 deg +/- 0.0162 mas
Dec = -50.180907 deg +/- 0.0174 mas
Coordinates = +18h53m05.90s -50d10m51.27s
Proper motion RA = 16.27 +/- 0.02 mas/yr
Proper motion Dec = -85.52 +/- 0.02 mas/yr
Radial velocity = -3.59 +/- 1.55 km/s
G mag = 8.101660 +/- 0.003166
BP mag = 8.496099 +/- 0.005965
RP mag = 7.526747 +/- 0.005561
GRVS mag = 7.288631 +/- 0.009597
Effective temperature = 5214 K
Surface gravity = 4.25
Metallicity = -0.82
G-band extinction = 0.00
Astrometric excess noise = 0.11
RUWE = 0.95
Non single star = 0
Single star probability from DSC-Combmod = 1.00
XP continuous = True
XP sampled = True
RVS spectrum = False
Reading input file...me...
0/1 [00:00<?, ?spec/s]
Storing Gaia XP plot: gaiaxp_6655168686921108864_0.jpg
Storing Gaia XP data: gaiaxp_6655168686921108864.dat
-> Querying Simbad...
Simbad ID = V* PZ Tel
Spectral type = G9IV
-> Querying VizieR...
2MASS source ID = 18530587-5010499
Separation between Gaia and 2MASS source = 1.2 mas
2MASS J mag = 6.856 +/- 0.021
2MASS H mag = 6.486 +/- 0.049
2MASS Ks mag = 6.366 +/- 0.024
ALLWISE source ID = J185305.88-501050.8
Separation between Gaia and ALLWISE source = 16.2 mas
WISE W1 mag = 6.294 +/- 0.098
WISE W2 mag = 6.247 +/- 0.027
WISE W3 mag = 6.275 +/- 0.015
WISE W4 mag = 6.257 +/- 0.052
We will then assign the parralax to a separate variable and write the magnitudes to a new dictionary.
[35]:
parallax = target_dict['Gaia parallax']
[36]:
magnitudes = {'GAIA/GAIA3.G': target_dict['GAIA/GAIA3.G'],
'GAIA/GAIA3.Gbp': target_dict['GAIA/GAIA3.Gbp'],
'GAIA/GAIA3.Grp': target_dict['GAIA/GAIA3.Grp'],
'GAIA/GAIA3.Grvs': target_dict['GAIA/GAIA3.Grvs'],
'2MASS/2MASS.J': target_dict['2MASS/2MASS.J'],
'2MASS/2MASS.H': target_dict['2MASS/2MASS.H'],
'2MASS/2MASS.Ks': target_dict['2MASS/2MASS.Ks']}
We also create a list of the filter names for use later on.
[37]:
filters = list(magnitudes.keys())
Adding stellar photometry#
We can now store the parallax and magnitudes of PZ Tel in the database by using the add_object method. This will also download a flux-calibrated spectrum of Vega and convert the magnitudes into fluxes.
[38]:
database.add_object(object_name='PZ Tel',
parallax=parallax,
app_mag=magnitudes,
spectrum=None)
----------
Add object
----------
Object name: PZ Tel
Units: None
Deredden: None
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, 173MB/s]
Adding spectrum: VegaReference: Bohlin et al. 2014, PASP, 126
URL: https://ui.adsabs.harvard.edu/abs/2014PASP..126..711B/abstract
Parallax (mas) = 21.16 +/- 0.02
Magnitudes:
- GAIA/GAIA3.G:
- Mean wavelength (um) = 2.1656e+00
- Apparent magnitude = 8.10 +/- 0.00
- Flux (W m-2 um-1) = 1.48e-11 +/- 4.32e-14
- GAIA/GAIA3.Gbp:
- Mean wavelength (um) = 2.1656e+00
- Apparent magnitude = 8.50 +/- 0.01
- Flux (W m-2 um-1) = 1.68e-11 +/- 9.21e-14
- GAIA/GAIA3.Grp:
- Mean wavelength (um) = 2.1656e+00
- Apparent magnitude = 7.53 +/- 0.01
- Flux (W m-2 um-1) = 1.27e-11 +/- 6.52e-14
- GAIA/GAIA3.Grvs:
- Mean wavelength (um) = 2.1656e+00
- Apparent magnitude = 7.29 +/- 0.01
- Flux (W m-2 um-1) = 1.13e-11 +/- 9.97e-14
- 2MASS/2MASS.J:
- Mean wavelength (um) = 2.1656e+00
- Apparent magnitude = 6.86 +/- 0.02
- Flux (W m-2 um-1) = 5.74e-12 +/- 1.11e-13
- 2MASS/2MASS.H:
- Mean wavelength (um) = 2.1656e+00
- Apparent magnitude = 6.49 +/- 0.05
- Flux (W m-2 um-1) = 2.93e-12 +/- 1.32e-13
- 2MASS/2MASS.Ks:
- Mean wavelength (um) = 2.1656e+00
- Apparent magnitude = 6.37 +/- 0.02
- Flux (W m-2 um-1) = 1.23e-12 +/- 2.72e-14
Adding a grid of model spectra#
Next, we will download the BT-NextGen grid and add the spectra of a specified \(T_\mathrm{eff}\) range to the database.
[39]:
database.add_model('bt-nextgen', teff_range=(3000., 8000.))
Downloading data from 'https://home.strw.leidenuniv.nl/~stolker/species/bt-nextgen.tgz' to file '/Users/tomasstolker/applications/species/docs/tutorials/data/bt-nextgen.tgz'.
-------------------------
Add grid of model spectra
-------------------------
Database tag: bt-nextgen
Model name: BT-NextGen
100%|████████████████████████████████████████| 385M/385M [00:00<00:00, 174GB/s]
SHA256 hash of downloaded file: d07093865a3632a45d8463c30fb885bca94d00cb40c524d433165175a2c49c6a
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 414/757 model spectra from BT-NextGen (368 MB)... [DONE]
Wavelength range (um) = 0.1 - 50
Sampling (lambda/d_lambda) = 4000
Teff range (K) = 3000.0 - 8000.0
Adding BT-NextGen model spectra... bt-nextgen_teff_8000_logg_5.0_feh_0.5_spec.dat
Grid points stored in the database:
- Teff = [3000. 3100. 3200. 3300. 3400. 3500. 3600. 3700. 3800. 3900. 4000. 4100.
4200. 4300. 4400. 4500. 4600. 4700. 4800. 4900. 5000. 5100. 5200. 5300.
5400. 5500. 5600. 5700. 5800. 5900. 6000. 6100. 6200. 6300. 6400. 6500.
6600. 6700. 6800. 6900. 7000. 7200. 7400. 7600. 7800. 8000.]
- log(g) = [3. 4. 5.]
- [Fe/H] = [0. 0.3 0.5]
Number of grid points per parameter:
- teff: 46
- logg: 3
- feh: 3
Number of stored grid points: 414
Number of interpolated grid points: 0
Number of missing grid points: 0
Fitting the 2MASS fluxes with the calibration spectrum#
Now that we have prepared the database, we can fit the photometric fluxes with the model grid. To do so, we use the FitModel class, which provides a Bayesian framework for parameter estimation. The argument of bounds contains a dictionary with the uniform priors that are used. In this example, we fix \(T_\mathrm{eff}\), \(\log(g)\), and \([\mathrm{Fe}/\mathrm{H}]\), so we effectively
scale the model spectrum to the data by adjusting the distance and radius. We will also account for a instrument-specific error inflation, relative to the actual uncertainties on the 2MASS fluxes (so allowing them to increase up to a factor of 10). Finally, the parallax is automatically included with a normal prior.
[40]:
fit = FitModel(object_name='PZ Tel',
model='bt-nextgen',
bounds={'teff': None,
'logg': None,
'feh': (0., 0.),
'radius': (1., 20.),
'GAIA/GAIA3_error': (1., 20.),
'2MASS/2MASS_error': (1., 10.)},
inc_phot=True,
inc_spec=False,
fit_corr=None,
apply_weights=False,
normal_prior=None)
-----------------
Fit model spectra
-----------------
Available photometric data for PZ Tel:
- 2MASS/2MASS.H
- 2MASS/2MASS.J
- 2MASS/2MASS.Ks
- GAIA/GAIA3.G
- GAIA/GAIA3.Gbp
- GAIA/GAIA3.Grp
- GAIA/GAIA3.Grvs
Interpolating 2MASS/2MASS.H... [DONE]
Interpolating 2MASS/2MASS.J... [DONE]
Interpolating 2MASS/2MASS.Ks... [DONE]
Interpolating GAIA/GAIA3.G... [DONE]
Interpolating GAIA/GAIA3.Gbp... [DONE]
Interpolating GAIA/GAIA3.Grp... [DONE]
Interpolating GAIA/GAIA3.Grvs... [DONE]
Fixing 1 parameters:
- feh = 0.0
Fitting 6 parameters:
- teff
- logg
- radius
- parallax
- 2MASS/2MASS_error
- GAIA/GAIA3_error
Uniform priors (min, max):
- teff = (3000.0, 8000.0)
- logg = (3.0, 5.0)
- radius = (1.0, 20.0)
- GAIA/GAIA3_error = (1.0, 20.0)
- 2MASS/2MASS_error = (1.0, 10.0)
Normal priors (mean, sigma):
- parallax = (21.162138650630606, 0.022285202518105507)
Weights for the log-likelihood function:
- 2MASS/2MASS.H = 1.00
- 2MASS/2MASS.J = 1.00
- 2MASS/2MASS.Ks = 1.00
- GAIA/GAIA3.G = 1.00
- GAIA/GAIA3.Gbp = 1.00
- GAIA/GAIA3.Grp = 1.00
- GAIA/GAIA3.Grvs = 1.00
We will sample the posterior distribution with the run_dynesty, which uses the nested sampling implementation of Dynesty. The samples will be stored in the database by the tag name. Let’s run the sampler with 500 live points!
[41]:
fit.run_dynesty(tag='pztel',
n_live_points=500,
output='dynesty/',
resume=False,
evidence_tolerance=0.5)
----------------------------
Nested sampling with Dynesty
----------------------------
Database tag: pztel
Number of live points: 500
Resume previous fit: False
Output folder: dynesty/
6852it [01:08, 100.63it/s, +500 | bound: 55 | nc: 1 | ncall: 59110 | eff(%): 12.544 | loglstar: -inf < 198.197 < inf | logz: 185.092 +/- 0.167 | dlogz: 0.001 > 0.500]
Samples shape: (7352, 6)
Number of iterations: 6852
Storing samples: dynesty/retrieval_post_equal_weights.dat
Nested sampling log-evidence: 185.09 +/- 0.25
Sample with the highest probability:
- Log-likelihood = 198.20
- teff = 5351.02
- logg = 3.45
- radius = 11.49
- parallax = 21.16
- 2MASS/2MASS_error = 1.60
- GAIA/GAIA3_error = 3.25
---------------------
Add posterior samples
---------------------
Database tag: pztel
Sampler: dynesty
Samples shape: (7352, 7)
Integrated autocorrelation time:
- teff: 1.03
- logg: 1.06
- radius: 1.00
- parallax: 1.08
- 2MASS/2MASS_error: 0.95
- GAIA/GAIA3_error: 0.96
- feh: fixed
Plotting the posterior distribution#
After running the sampler, we can plot the posterior distribution of the 3 free parameters with the plot_posterior by simply pointing to the database tag that was specified with run_dynesty.
[42]:
fig = plot_posterior(tag='pztel',
offset=(-0.3, -0.3),
title_fmt=['.0f', '.1f', '.1f', '.2f', '.1f', '.1f'],
object_type='star',
output=None)
---------------------
Get posterior samples
---------------------
Database tag: pztel
Random samples: None
Samples shape: (7352, 7)
Parameters:
- teff
- logg
- radius
- parallax
- 2MASS/2MASS_error
- GAIA/GAIA3_error
- feh
Uniform priors (min, max):
- 2MASS/2MASS_error = (1.0, 10.0)
- GAIA/GAIA3_error = (1.0, 20.0)
- logg = (3.0, 5.0)
- radius = (1.0, 20.0)
- teff = (3000.0, 8000.0)
Normal priors (mean, sigma):
- parallax = (21.162138650630606, 0.022285202518105507)
----------------------------
Plot posterior distributions
----------------------------
Database tag: pztel
Object type: star
Manual parameters: None
Median parameters:
- teff = 5.30e+03
- logg = 3.79e+00
- radius = 1.17e+01
- parallax = 2.12e+01
- 2MASS/2MASS_error = 2.08e+00
- GAIA/GAIA3_error = 4.31e+00
- feh = 0.00e+00
Sample with highest probability:
- teff = 5.35e+03
- logg = 3.45e+00
- radius = 1.15e+01
- parallax = 2.12e+01
- 2MASS/2MASS_error = 1.60e+00
- GAIA/GAIA3_error = 3.25e+00
- feh = 0.00e+00
The corner plot shows that the surface gravity, \(\log{(g)}\), is unconstrained. That is to be expected when fitting only photometric fluxes. The plot_posterior 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.
Extracting spectral samples#
Later on, we will create a plot of the data, best-fit spectrum, and random spectral samples of the posterior distribution. We start by drawing 30 random samples the posterior distribution and calculate spectra at \(R = 100\). The get_mcmc_spectra method of Database returns a list with ModelBox objects.
[43]:
samples = database.get_mcmc_spectra(tag='pztel',
random=30,
wavel_range=(0.2, 10.),
spec_res=100.)
---------------------
Get posterior spectra
---------------------
Database tag: pztel
Number of samples: 30
Wavelength range (um): (0.2, 10.0)
Resolution: 100.0
/Users/tomasstolker/applications/species/species/read/read_model.py:799: UserWarning: The '2MASS/2MASS_error' parameter is not required by 'bt-nextgen' so the parameter will be ignored. The mandatory parameters are ['teff', 'logg', 'feh'].
warnings.warn(
/Users/tomasstolker/applications/species/species/read/read_model.py:799: UserWarning: The 'GAIA/GAIA3_error' parameter is not required by 'bt-nextgen' so the parameter will be ignored. The mandatory parameters are ['teff', 'logg', 'feh'].
warnings.warn(
Let’s have a look at the content of the first ModelBox by using the open_box method.
[44]:
samples[0].open_box()
Opening ModelBox...
model = bt-nextgen
type = None
wavelength = [ 0.19996262 0.20001213 0.20006165 ... 9.99710369 9.99957902
10.00205496]
flux = [6.09985870e-15 6.13431738e-15 6.16877050e-15 ... 3.49595033e-15
3.49421933e-15 3.49256521e-15]
parameters = {'teff': 5258.833287074302, 'logg': 3.69725492683077, 'radius': 11.949839039499988, 'parallax': 21.14746182640784, '2MASS/2MASS_error': 2.0656904665000653, 'GAIA/GAIA3_error': 2.9082813777814804, 'feh': 0.0, 'luminosity': 1.0390735478708635, 'mass': 286.92176283542005}
quantity = flux
contribution = None
bol_flux = None
Next, we extract the best-fit parameters from the posterior distribution, for which we adopt the sample with the highest likelihood value. The get_probable_sample function returns a dictionary with the parameters, including the ones that were fixed with FitModel.
[45]:
best_sample = database.get_probable_sample(tag='pztel')
print(best_sample)
-----------------------------------
Get sample with highest probability
-----------------------------------
Database tag: pztel
Parameters:
- teff = 5351.02
- logg = 3.45
- radius = 11.49
- parallax = 21.16
- 2MASS/2MASS_error = 1.60
- GAIA/GAIA3_error = 3.25
- feh = 0.00
{'teff': 5351.017388301505, 'logg': 3.4545974786514826, 'radius': 11.492505237995987, 'parallax': 21.164968841565905, '2MASS/2MASS_error': 1.5960808597977203, 'GAIA/GAIA3_error': 3.2456675816583602, 'feh': 0.0}
Next, we interpolate the the BT-NextGen grid at the best-fit parameters. To do so, we first create an instance of ReadModel and then use the get_model method to interpolate the grid with the best-fit parameters that are provided as argument of model_par. Similar to
get_mcmc_spectra, we also smooth this spectrum to \(R = 100\). The model spectrum is again returned in a ModelBox.
[46]:
read_model = ReadModel(model='bt-nextgen', wavel_range=(0.2, 10.))
model_box = read_model.get_model(model_param=best_sample, spec_res=100.)
Each Box is a Python object so the content is easily accessed as attributes.
[47]:
print(model_box.parameters)
{'teff': 5351.017388301505, 'logg': 3.4545974786514826, 'radius': 11.492505237995987, 'parallax': 21.164968841565905, '2MASS/2MASS_error': 1.5960808597977203, 'GAIA/GAIA3_error': 3.2456675816583602, 'feh': 0.0, 'luminosity': 1.0302424358894176, 'mass': 151.77888286492623}
Extracting the object data#
For the plot, we also require the data of PZ Tel. These can be extracted with the get_object, which returns the data (e.g. photometry, spectra, and parallax) for a given object_name.
[48]:
object_box = database.get_object(object_name='PZ Tel')
----------
Get object
----------
Object name: PZ Tel
Include photometry: True
Include spectra: True
The data are stored in an ObjectBox. Let’s have a look at the content by using open_box.
[49]:
object_box.open_box()
Opening ObjectBox...
name = PZ Tel
filters = ['2MASS/2MASS.H', '2MASS/2MASS.J', '2MASS/2MASS.Ks', 'GAIA/GAIA3.G', 'GAIA/GAIA3.Gbp', 'GAIA/GAIA3.Grp', 'GAIA/GAIA3.Grvs']
mean_wavel = {'2MASS/2MASS.H': 1.6513664598414621, '2MASS/2MASS.J': 1.24105170694321, '2MASS/2MASS.Ks': 2.1656311232670684, 'GAIA/GAIA3.G': 0.6390220344836051, 'GAIA/GAIA3.Gbp': 0.5182576434561539, 'GAIA/GAIA3.Grp': 0.7825078823058362, 'GAIA/GAIA3.Grvs': 0.8579038886678563}
magnitude = {'2MASS/2MASS.H': array([6.48600006, 0.049 ]), '2MASS/2MASS.J': array([6.85599995, 0.021 ]), '2MASS/2MASS.Ks': array([6.36600018, 0.024 ]), 'GAIA/GAIA3.G': array([8.10165977e+00, 3.16641662e-03]), 'GAIA/GAIA3.Gbp': array([8.49609852e+00, 5.96450345e-03]), 'GAIA/GAIA3.Grp': array([7.52674723e+00, 5.56069163e-03]), 'GAIA/GAIA3.Grvs': array([7.28863096, 0.00959738])}
flux = {'2MASS/2MASS.H': array([2.92837677e-12, 1.32204464e-13]), '2MASS/2MASS.J': array([5.74224978e-12, 1.11071881e-13]), '2MASS/2MASS.Ks': array([1.22822262e-12, 2.71518473e-14]), 'GAIA/GAIA3.G': array([1.47960664e-11, 4.31509760e-14]), 'GAIA/GAIA3.Gbp': array([1.67616818e-11, 9.20809414e-14]), 'GAIA/GAIA3.Grp': array([1.27300626e-11, 6.51983988e-14]), 'GAIA/GAIA3.Grvs': array([1.12809246e-11, 9.97191629e-14])}
spectrum = None
parallax = [21.16213865 0.0222852 ]
distance = None
With FitModel it is possible to account for systematic uncertainties, either in the data or the model, for example to scale individual spectra or inflate the uncertainties. Since we inflated the 2MASS uncertainties during the fit, we will use the update_objectbox function to adjust the photometric precision with the best-fit error inflation as extracted with get_probable_sample.
[50]:
object_box = update_objectbox(object_box, best_sample)
Inflating the uncertainty of 2MASS/2MASS.H by a factor 1.60 to 2.11e-13 (W m-2 um-1)... [DONE]
Inflating the uncertainty of 2MASS/2MASS.J by a factor 1.60 to 1.77e-13 (W m-2 um-1)... [DONE]
Inflating the uncertainty of 2MASS/2MASS.Ks by a factor 1.60 to 4.33e-14 (W m-2 um-1)... [DONE]
Inflating the uncertainty of GAIA/GAIA3.G by a factor 3.25 to 1.40e-13 (W m-2 um-1)... [DONE]
Inflating the uncertainty of GAIA/GAIA3.Gbp by a factor 3.25 to 2.99e-13 (W m-2 um-1)... [DONE]
Inflating the uncertainty of GAIA/GAIA3.Grp by a factor 3.25 to 2.12e-13 (W m-2 um-1)... [DONE]
Inflating the uncertainty of GAIA/GAIA3.Grvs by a factor 3.25 to 3.24e-13 (W m-2 um-1)... [DONE]
Residuals and synthetic photometry#
Before creating the plot, there are two more boxes that we will create. First, we use the get_residuals function to calculate the residuals (i.e. data minus model, relative to the data uncertainties), together with the reduced \(\chi^2\). The residuals and mean wavelengths are stored in a ResidualsBox for each of the 2MASS filters.
[51]:
residuals = get_residuals(datatype='model',
spectrum='bt-nextgen',
parameters=best_sample,
objectbox=object_box,
inc_phot=True,
inc_spec=False)
-------------------
Calculate residuals
-------------------
Data type: model
Spectrum name: bt-nextgen
Include photometry: True
Include spectra: False
Parameters:
- teff = 5351.02
- logg = 3.45
- radius = 11.49
- parallax = 21.16
- 2MASS/2MASS_error = 1.60
- GAIA/GAIA3_error = 3.25
- feh = 0.00
- luminosity = 1.03e+00
- mass = 151.78
--------------------------
Calculate multi-photometry
--------------------------
Data type: model
Spectrum name: bt-nextgen
Parameters:
- teff = 5351.02
- logg = 3.45
- radius = 11.49
- parallax = 21.16
- 2MASS/2MASS_error = 1.60
- GAIA/GAIA3_error = 3.25
- feh = 0.00
- luminosity = 1.03e+00
- mass = 151.78
Magnitudes:
- 2MASS/2MASS.H = 6.47
- 2MASS/2MASS.J = 6.90
- 2MASS/2MASS.Ks = 6.42
- GAIA/GAIA3.G = 8.10
- GAIA/GAIA3.Gbp = 8.50
- GAIA/GAIA3.Grp = 7.52
- GAIA/GAIA3.Grvs = 7.33
Fluxes (W m-2 um-1):
- 2MASS/2MASS.H = 2.97e-12
- 2MASS/2MASS.J = 5.51e-12
- 2MASS/2MASS.Ks = 1.17e-12
- GAIA/GAIA3.G = 1.49e-11
- GAIA/GAIA3.Gbp = 1.67e-11
- GAIA/GAIA3.Grp = 1.28e-11
- GAIA/GAIA3.Grvs = 1.08e-11
Residuals (sigma):
- 2MASS/2MASS.H = -0.15
- 2MASS/2MASS.J = 1.10
- 2MASS/2MASS.Ks = 1.10
- GAIA/GAIA3.G = -0.49
- GAIA/GAIA3.Gbp = 0.11
- GAIA/GAIA3.Grp = -0.13
- GAIA/GAIA3.Grvs = 1.35
Reduced chi2 = inf
Number of degrees of freedom = 0
[52]:
residuals.open_box()
Opening ResidualsBox...
name = PZ Tel
photometry = {'2MASS/2MASS.H': array([ 1.65136646, -0.15401567]), '2MASS/2MASS.J': array([1.24105171, 1.10181046]), '2MASS/2MASS.Ks': array([2.16563112, 1.09984114]), 'GAIA/GAIA3.G': array([ 0.63902203, -0.48712528]), 'GAIA/GAIA3.Gbp': array([0.51825764, 0.11020271]), 'GAIA/GAIA3.Grp': array([ 0.78250788, -0.12578623]), 'GAIA/GAIA3.Grvs': array([0.85790389, 1.35339561])}
spectrum = None
chi2_red = inf
Secondly, we will also use the dictionary with best-fit parameters to calculate synthetic photometry for the three filters that were used in the fit. The multi_photometry returns a SynphotBox, which includes a dictionary with the photometric fluxes.
[53]:
syn_phot = multi_photometry(datatype='model',
spectrum='bt-nextgen',
filters=filters,
parameters=best_sample)
--------------------------
Calculate multi-photometry
--------------------------
Data type: model
Spectrum name: bt-nextgen
Parameters:
- teff = 5351.02
- logg = 3.45
- radius = 11.49
- parallax = 21.16
- 2MASS/2MASS_error = 1.60
- GAIA/GAIA3_error = 3.25
- feh = 0.00
- luminosity = 1.03e+00
- mass = 151.78
Magnitudes:
- GAIA/GAIA3.G = 8.10
- GAIA/GAIA3.Gbp = 8.50
- GAIA/GAIA3.Grp = 7.52
- GAIA/GAIA3.Grvs = 7.33
- 2MASS/2MASS.J = 6.90
- 2MASS/2MASS.H = 6.47
- 2MASS/2MASS.Ks = 6.42
Fluxes (W m-2 um-1):
- GAIA/GAIA3.G = 1.49e-11
- GAIA/GAIA3.Gbp = 1.67e-11
- GAIA/GAIA3.Grp = 1.28e-11
- GAIA/GAIA3.Grvs = 1.08e-11
- 2MASS/2MASS.J = 5.51e-12
- 2MASS/2MASS.H = 2.97e-12
- 2MASS/2MASS.Ks = 1.17e-12
[54]:
syn_phot.open_box()
Opening SynphotBox...
name = synphot
wavelength = {'GAIA/GAIA3.G': 0.6390220344836051, 'GAIA/GAIA3.Gbp': 0.5182576434561539, 'GAIA/GAIA3.Grp': 0.7825078823058362, 'GAIA/GAIA3.Grvs': 0.8579038886678563, '2MASS/2MASS.J': 1.24105170694321, '2MASS/2MASS.H': 1.6513664598414621, '2MASS/2MASS.Ks': 2.1656311232670684}
flux = {'GAIA/GAIA3.G': 1.48674548408908e-11, 'GAIA/GAIA3.Gbp': 1.6727218358268868e-11, 'GAIA/GAIA3.Grp': 1.275791525883396e-11, 'GAIA/GAIA3.Grvs': 1.0822571530401296e-11, '2MASS/2MASS.J': 5.511749934471628e-12, '2MASS/2MASS.H': 2.966727237241917e-12, '2MASS/2MASS.Ks': 1.1719769623100563e-12}
app_mag = {'GAIA/GAIA3.G': (8.09643388241963, None), 'GAIA/GAIA3.Gbp': (8.498333182854331, None), 'GAIA/GAIA3.Grp': (7.524374291602889, None), 'GAIA/GAIA3.Grvs': (7.333666537005324, None), '2MASS/2MASS.J': (6.900481381689703, None), '2MASS/2MASS.H': (6.4718733966714295, None), '2MASS/2MASS.Ks': (6.4168952127647385, None)}
abs_mag = {'GAIA/GAIA3.G': (4.724522049695395, None), 'GAIA/GAIA3.Gbp': (5.126421350130096, None), 'GAIA/GAIA3.Grp': (4.152462458878654, None), 'GAIA/GAIA3.Grvs': (3.9617547042810886, None), '2MASS/2MASS.J': (3.528569548965468, None), '2MASS/2MASS.H': (3.0999615639471942, None), '2MASS/2MASS.Ks': (3.0449833800405033, None)}
Plotting the data and model spectra#
We have now prepared all the boxes with data so we are ready to combine them in a plot of the spectral energy distribution (SED) of PZ Tel! The plot_spectrum function requires a list of Box objects as argument of boxes. For each box we can set the plot style, by providing a list with dictionaries as
argument of plot_kwargs, in the same order as the list of boxes. Items in the list can be set to None, in which case some default values are used. The ResidualsBox is passed as argument of residuals and will also plot the filter profiles by providing the list with names as argument to filters. Finally, there is a handful of parameters that can be adjusted for the appearance of the plot
(see the API documentation of plot_spectrum for details). Let’s have a look at the plot!
[26]:
fig = plot_spectrum(boxes=[samples, model_box, object_box],
filters=filters,
residuals=residuals,
plot_kwargs=[{'ls': '-', 'lw': 0.2, 'color': 'gray'},
{'ls': '-', 'lw': 0.7, 'color': 'black'},
{'GAIA/GAIA3.G': {'marker': 's', 'ms': 5., 'color': 'tab:orange', 'ls': 'none', 'label': 'Gaia DR3'},
'GAIA/GAIA3.Gbp': {'marker': 's', 'ms': 5., 'color': 'tab:orange', 'ls': 'none'},
'GAIA/GAIA3.Grp': {'marker': 's', 'ms': 5., 'color': 'tab:orange', 'ls': 'none'},
'GAIA/GAIA3.Grvs': {'marker': 's', 'ms': 5., 'color': 'tab:orange', 'ls': 'none'},
'2MASS/2MASS.J': {'marker': 's', 'ms': 5., 'color': 'tab:pink', 'ls': 'none', 'label': '2MASS'},
'2MASS/2MASS.H': {'marker': 's', 'ms': 5., 'color': 'tab:pink', 'ls': 'none'},
'2MASS/2MASS.Ks': {'marker': 's', 'ms': 5., 'color': 'tab:pink', 'ls': 'none'}}],
xlim=(0.2, 2.5),
ylim=(-1.e-12, 1.3e-11),
ylim_res=(-5., 5.),
scale=('linear', 'linear'),
offset=(-0.5, -0.08),
figsize=(6., 4.),
object_type='star',
quantity='flux',
legend=[{'loc': 'upper right', 'fontsize': 9.},
{'loc': (0.15, 0.2), 'fontsize': 10.}],
output=None)
Plotting spectrum...
[DONE]
Note that the photometric uncertainties are inflated in the plot by the best-fit parameters that were fitted. These parameters had been applied to the data when running the update_objectbox on the ObjectBox.
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.
[27]:
fig.axes
[27]:
[<Axes: ylabel='$\\lambda$$F_\\lambda$ (10$^{-11}$ W m$^{-2}$)'>,
<Axes: ylabel='$T_\\lambda$'>,
<Axes: xlabel='Wavelength (μm)', ylabel='$\\Delta$$F_\\lambda$ ($\\sigma$)'>]
Photometric calibration#
Now that we have the posterior samples of the atmospheric parameters, we can calculated synthetic photometry (including uncertainties) for any other filter. As an example, we will calculate the magnitude and flux for the VLT/ERIS \(M'\) filter. The plot_mag_posterior function plots the distribution of the magnitudes, by propagating the posterior of the atmospheric parameters, and
returns an array with the samples. We simply need to specify the database tag and provide the filter_name name as listed on the SVO Filter Profile Service). The function returns an ndarray with the samples and the Figure of the plot that can be used for further customization.
[28]:
phot_mag, fig = plot_mag_posterior(tag='pztel',
filter_name='Paranal/ERIS.Mp',
xlim=(6.0, 6.8),
output=None)
Adding filter: Paranal/ERIS.Mp... [DONE]
Plotting photometry samples...
[DONE]
There is also the get_mcmc_photometry method of Database, which works in a similar way as plot_mag_posterior, but can also return the posterior of the flux density instead of the magnitude.
[29]:
phot_flux = database.get_mcmc_photometry(tag='pztel',
filter_name='Paranal/ERIS.Mp',
phot_type='flux')
/Users/tomasstolker/applications/species/species/read/read_model.py:799: UserWarning: The '2MASS/2MASS_error' parameter is not required by 'bt-nextgen' so the parameter will be ignored. The mandatory parameters are ['teff', 'logg', 'feh'].
warnings.warn(
/Users/tomasstolker/applications/species/species/read/read_model.py:799: UserWarning: The 'GAIA/GAIA3_error' parameter is not required by 'bt-nextgen' so the parameter will be ignored. The mandatory parameters are ['teff', 'logg', 'feh'].
warnings.warn(
To calculate the flux density in the \(M'\)-band, we simply adopt the mean and the standard deviation of the distribution.
[30]:
print(f'Flux density (W m-2 um-1) = {np.mean(phot_flux):.2e} +/- {np.std(phot_flux):.2e}')
Flux density (W m-2 um-1) = 5.78e-14 +/- 2.25e-15
Spectral calibration#
We can also compute a stellar, synthetic spectrum with uncertainties for a specific instrument. To do so, we use again the get_mcmc_spectra method from earlier, but this time we draw a larger number of spectra. Furthermore, for our hypothetical instrument, we assume a resolving power of \(R = 100\) and set the wavelength points (20 linearly-spaced points between 1 and 2.5
\(\mu\)m) as argument of wavel_resample.
[26]:
samples = database.get_mcmc_spectra(tag='pztel',
random=100,
wavel_range=(0.2, 10.),
spec_res=100.,
wavel_resample=np.linspace(1., 2.5, 20))
Next, we calculate the mean and standard deviation of the 100 samples, which we adopt as the synthetic spectrum of PZ Tel.
[27]:
spec_flux = np.mean([box.flux for box in samples], axis=0)
spec_sigma = np.std([box.flux for box in samples], axis=0)
Let’s write the synthetic spectrum to a text file.
[28]:
calib_spec = np.column_stack([samples[0].wavelength, spec_flux, spec_sigma])
np.savetxt('calib_spec.dat', calib_spec)
Now we can use the add_object method to append the spectrum to the data of PZ Tel in the database. The parallax and magnitudes were already provided previously so we can set the arguments of parallax and app_mag to None.
[29]:
database.add_object(object_name='PZ Tel',
parallax=None,
app_mag=None,
spectrum={'calibration': ('calib_spec.dat', None, 100.)})
Adding object: PZ Tel
- Spectrum:
- Database tag: calibration
- Filename: calib_spec.dat
- Data shape: (20, 3)
- Wavelength range (um): 1.00 - 2.50
- Mean flux (W m-2 um-1): 3.20e-12
- Mean error (W m-2 um-1): 1.00e-13
- Spectral resolution:
- calibration: 100.0
Next, we use the get_object method for extracting all the data of PZ Tel from the database and storing these in an ObjectBox. We need to apply again the update_objectbox function to account for the inflated errors of the 2MASS fluxes.
[30]:
object_box = database.get_object(object_name='PZ Tel')
object_box = update_objectbox(object_box, best_sample)
Getting object: PZ Tel... [DONE]
Inflating the uncertainty of 2MASS/2MASS.H by a factor 1.50 to 1.98e-13 (W m-2 um-1)... [DONE]
Inflating the uncertainty of 2MASS/2MASS.J by a factor 1.50 to 1.66e-13 (W m-2 um-1)... [DONE]
Inflating the uncertainty of 2MASS/2MASS.Ks by a factor 1.50 to 4.07e-14 (W m-2 um-1)... [DONE]
Inflating the uncertainty of GAIA/GAIA3.G by a factor 5.96 to 2.57e-13 (W m-2 um-1)... [DONE]
Inflating the uncertainty of GAIA/GAIA3.Gbp by a factor 5.96 to 5.49e-13 (W m-2 um-1)... [DONE]
Inflating the uncertainty of GAIA/GAIA3.Grp by a factor 5.96 to 3.89e-13 (W m-2 um-1)... [DONE]
Inflating the uncertainty of GAIA/GAIA3.Grvs by a factor 5.96 to 5.94e-13 (W m-2 um-1)... [DONE]
If we now open the ObjectBox with open_box then we see that the synthetic spectrum is indeed included together with the previously added data.
[31]:
object_box.open_box()
Opening ObjectBox...
name = PZ Tel
filters = ['2MASS/2MASS.H', '2MASS/2MASS.J', '2MASS/2MASS.Ks', 'GAIA/GAIA3.G', 'GAIA/GAIA3.Gbp', 'GAIA/GAIA3.Grp', 'GAIA/GAIA3.Grvs']
mean_wavel = {'2MASS/2MASS.H': 1.6513664598414621, '2MASS/2MASS.J': 1.24105170694321, '2MASS/2MASS.Ks': 2.1656311232670684, 'GAIA/GAIA3.G': 0.6390220344836051, 'GAIA/GAIA3.Gbp': 0.5182576434561539, 'GAIA/GAIA3.Grp': 0.7825078823058362, 'GAIA/GAIA3.Grvs': 0.8579038886678563}
magnitude = {'2MASS/2MASS.H': array([6.48600006, 0.049 ]), '2MASS/2MASS.J': array([6.85599995, 0.021 ]), '2MASS/2MASS.Ks': array([6.36600018, 0.024 ]), 'GAIA/GAIA3.G': array([8.10165977e+00, 3.16641662e-03]), 'GAIA/GAIA3.Gbp': array([8.49609852e+00, 5.96450345e-03]), 'GAIA/GAIA3.Grp': array([7.52674723e+00, 5.56069163e-03]), 'GAIA/GAIA3.Grvs': array([7.28863096, 0.00959738])}
flux = {'2MASS/2MASS.H': array([2.92837677e-12, 2.38102242e-13]), '2MASS/2MASS.J': array([5.74224978e-12, 2.00042140e-13]), '2MASS/2MASS.Ks': array([1.22822262e-12, 4.89008883e-14]), 'GAIA/GAIA3.G': array([1.47960664e-11, 2.60789423e-13]), 'GAIA/GAIA3.Gbp': array([1.67616818e-11, 5.56505039e-13]), 'GAIA/GAIA3.Grp': array([1.27300626e-11, 3.94036343e-13]), 'GAIA/GAIA3.Grvs': array([1.12809246e-11, 6.02667781e-13])}
spectrum = {'calibration': (array([[1.00000000e+00, 8.53705861e-12, 1.92172921e-13],
[1.07894737e+00, 7.22693556e-12, 1.72791615e-13],
[1.15789474e+00, 6.37201036e-12, 1.61693227e-13],
[1.23684211e+00, 5.60777597e-12, 1.55672402e-13],
[1.31578947e+00, 4.94908369e-12, 1.46302016e-13],
[1.39473684e+00, 4.47532837e-12, 1.44860258e-13],
[1.47368421e+00, 3.93655436e-12, 1.32863078e-13],
[1.55263158e+00, 3.58643326e-12, 1.34067922e-13],
[1.63157895e+00, 3.17442502e-12, 1.25354728e-13],
[1.71052632e+00, 2.72390996e-12, 1.08215990e-13],
[1.78947368e+00, 2.34909262e-12, 9.43669862e-14],
[1.86842105e+00, 1.98385854e-12, 8.09913122e-14],
[1.94736842e+00, 1.69641995e-12, 6.89867541e-14],
[2.02631579e+00, 1.50677407e-12, 6.08331482e-14],
[2.10526316e+00, 1.31568077e-12, 5.31711059e-14],
[2.18421053e+00, 1.14473785e-12, 4.65322991e-14],
[2.26315789e+00, 1.01282844e-12, 4.11077608e-14],
[2.34210526e+00, 8.75838244e-13, 3.38857642e-14],
[2.42105263e+00, 7.67685204e-13, 2.93872731e-14],
[2.50000000e+00, 6.81197444e-13, 2.61683775e-14]]), None, None, 100.0)}
parallax = [21.16213865 0.0222852 ]
distance = None
Let’s have a look at the content of the returned ObjectBox by simply passing it to the plot_spectrum function!
[37]:
fig = plot_spectrum(boxes=[object_box],
plot_kwargs=[{'calibration': {'marker': 'o', 'ms': 5., 'color': 'black', 'ls': 'none', 'label': 'Synthetic spectrum'},
'2MASS/2MASS.J': {'marker': 'o', 'ms': 5., 'color': 'tab:blue', 'ls': 'none', 'label': '2MASS'},
'2MASS/2MASS.H': {'marker': 'o', 'ms': 5., 'color': 'tab:blue', 'ls': 'none'},
'2MASS/2MASS.Ks': {'marker': 'o', 'ms': 5., 'color': 'tab:blue', 'ls': 'none'}}],
xlim=(0.9, 2.5),
ylim=(0., 1e-11),
scale=('linear', 'linear'),
offset=(-0.1, -0.07),
figsize=(6., 3.),
legend={'loc': 'upper right', 'fontsize': 12.},
output=None)
Plotting spectrum...
[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.
[38]:
fig.axes
[38]:
[<Axes: xlabel='Wavelength (μm)', ylabel='$F_\\lambda$ (10$^{-11}$ W m$^{-2}$ μm$^{-1}$)'>]