"""
Module with reading functionalities for isochrones and cooling curves.
"""
import os
import warnings
from configparser import ConfigParser
from typing import Dict, List, Optional, Tuple, Union
import h5py
import numpy as np
from scipy import interpolate
from typeguard import typechecked
from species.core.box import (
ColorMagBox,
ColorColorBox,
CoolingBox,
IsochroneBox,
ModelBox,
PhotometryBox,
create_box,
)
from species.read.read_model import ReadModel
from species.util.convert_util import apparent_to_absolute
from species.util.plot_util import update_labels
[docs]
class ReadIsochrone:
"""
Class for reading isochrone data from the database.
This class interpolates the evolutionary track or
isochrone data. Please carefully check for interpolation
effects. Setting ``masses=None`` in
:func:`~species.read.read_isochrone.ReadIsochrone.get_isochrone`
extracts the isochrones at the masses of the original
grid, so using that option helps with comparing results
for which the masses have been interpolated. Similarly, by
setting ``ages=None`` with the
:func:`~species.read.read_isochrone.ReadIsochrone.get_isochrone`
method will fix the ages to those of the original grid.
"""
@typechecked
def __init__(
self,
tag: Optional[str] = None,
create_regular_grid: bool = False,
extrapolate: bool = False,
) -> None:
"""
Parameters
----------
tag : str, None
Database tag of the isochrone data (e.g. 'ames-cond',
'sonora+0.0', 'atmo-ceq'). A list with the isochrone
data that are stored in the current ``species``
database is printed if the argument of ``tag`` is
set to ``None``.
create_regular_grid : bool
Evolutionary grids can be irregular in the (age, mass)
space. By setting ``create_regular_grid=True``, the
grid will be interpolated and extrapolate onto
a regular grid that is based on all the unique age
and mass values of the original grid. The
resampling of the evolutionary parameters
(i.e. :math:`T_\\mathrm{eff}` (K),
:math:`\\log{(L/L_\\odot)}`, :math:`\\log{(g)}`, and
:math:`R`) is done as function of mass, separately
for each age. Setting ``create_regular_grid=True``
can be helpful if some values are missing at the edge
of the grid. However, please carefully check results
since there might be inaccuracies in the extrapolated
parts of the parameter space, in particular for the
cooling curves extracted with
:func:`~species.read.read_isochrone.ReadIsochrone.get_cooling_curve`.
extrapolate : str
DEPRECATED: This parameter has been renamed to
``create_regular_grid`` and will be removed in a future
release. Please use the ``create_regular_grid``
parameter instead.
Returns
-------
NoneType
None
"""
self.tag = tag
self.extrapolate = extrapolate
self.create_regular_grid = create_regular_grid
if self.extrapolate:
warnings.warn(
"The 'extrapolate' parameter has been "
"renamed to 'create_regular_grid' and "
"will be removed in a future release.",
DeprecationWarning,
)
if not self.create_regular_grid:
warnings.warn(
"Setting 'create_regular_grid=True' since 'extrapolate=True'."
)
config_file = os.path.join(os.getcwd(), "species_config.ini")
config = ConfigParser()
config.read(config_file)
self.database = config["species"]["database"]
if self.tag is None:
with h5py.File(self.database, "r") as hdf5_file:
tag_list = list(hdf5_file["isochrones"])
self.tag = input(
"Please select one of the following "
"isochrone tags that are stored in "
"the database or use 'add_isochrones' "
"to add another model to the database:"
f"\n{tag_list}:\n"
)
with h5py.File(self.database, "r") as hdf5_file:
if f"isochrones/{self.tag}" not in hdf5_file:
tag_list = list(hdf5_file["isochrones"])
raise ValueError(
f"There is no isochrone data stored with the "
f"selected tag '{tag}'. The following isochrone "
f"tags are found in the database: {tag_list}"
)
self.mag_models = [
"ames",
"atmo",
"baraffe",
"bt-settl",
"linder2019",
"manual",
"nextgen",
]
# Connect isochrone model with atmosphere model
# key = isochrone model, value = atmosphere model
self.match_model = {
"ames-cond": "ames-cond",
"ames-dusty": "ames-dusty",
"atmo-ceq": "atmo-ceq",
"atmo-neq-strong": "atmo-neq-strong",
"atmo-neq-weak": "atmo-neq-weak",
"bt-settl": "bt-settl",
"linder2019-petitCODE-metal_-0.4": "petitcode-linder2019-clear",
"linder2019-petitCODE-metal_0.0": "petitcode-linder2019-clear",
"linder2019-petitCODE-metal_0.4": "petitcode-linder2019-clear",
"linder2019-petitCODE-metal_0.8": "petitcode-linder2019-clear",
"linder2019-petitCODE-metal_1.2": "petitcode-linder2019-clear",
"linder2019-petitCODE-metal_-0.4-fsed_1.0": "petitcode-linder2019-cloudy",
"linder2019-petitCODE-metal_0.0-fsed_1.0": "petitcode-linder2019-cloudy",
"linder2019-petitCODE-metal_0.4-fsed_1.0": "petitcode-linder2019-cloudy",
"linder2019-petitCODE-metal_0.8-fsed_1.0": "petitcode-linder2019-cloudy",
"linder2019-petitCODE-metal_1.2-fsed_1.0": "petitcode-linder2019-cloudy",
"saumon2008-nc_solar": "saumon2008-clear",
"saumon2008-f2_solar": "saumon2008-cloudy",
"sonora+0.0": "sonora-bobcat",
}
# Set attribute with extra parameters
# The attribute will be overwritten by any of
# the methods that have the extra_param parameter
self.extra_param = None
if self.tag.split("-")[0] == "linder2019":
self.extra_param = {}
tag_split = self.tag.split("_")
if len(tag_split) == 2:
self.extra_param["feh"] = float(tag_split[1])
elif len(tag_split) == 3:
self.extra_param["feh"] = float(tag_split[1][:-5])
self.extra_param["fsed"] = float(tag_split[2])
print(f"Setting 'extra_param' attribute: {self.extra_param}")
@typechecked
def _read_data(
self,
) -> Tuple[
str,
np.ndarray,
np.ndarray,
np.ndarray,
np.ndarray,
np.ndarray,
np.ndarray,
Optional[np.ndarray],
]:
"""
Internal function for reading the evolutionary
data from the database.
Returns
-------
str
Model name.
np.ndarray
Array with the age (Myr).
np.ndarray
Array with the mass (:math:`M_\\mathrm{J}`).
np.ndarray
Array with the :math:`T_\\mathrm{eff}` (K).
np.ndarray
Array with the :math:`\\log{(L/L_\\odot)}`.
np.ndarray
Array with the :math:`\\log{(g)}`.
np.ndarray
Array with the radius (:math:`R_\\mathrm{J}`).
np.ndarray, None
Optional array with the absolute magnitudes. The
array has two axes with the length of the second
axis equal to the number of filters for which
there are magnitudes available.
"""
with h5py.File(self.database, "r") as hdf5_file:
model_name = hdf5_file[f"isochrones/{self.tag}/age"].attrs["model"]
iso_age = np.asarray(hdf5_file[f"isochrones/{self.tag}/age"])
iso_mass = np.asarray(hdf5_file[f"isochrones/{self.tag}/mass"])
iso_teff = np.asarray(hdf5_file[f"isochrones/{self.tag}/teff"])
iso_loglum = np.asarray(hdf5_file[f"isochrones/{self.tag}/log_lum"])
iso_logg = np.asarray(hdf5_file[f"isochrones/{self.tag}/log_g"])
iso_radius = np.asarray(hdf5_file[f"isochrones/{self.tag}/radius"])
if f"isochrones/{self.tag}/magnitudes" in hdf5_file:
iso_mag = np.asarray(hdf5_file[f"isochrones/{self.tag}/magnitudes"])
else:
iso_mag = None
if self.create_regular_grid:
age_unique = np.unique(iso_age)
mass_unique = np.unique(iso_mass)
n_ages = age_unique.shape[0]
n_masses = mass_unique.shape[0]
new_age = np.zeros((n_ages * n_masses))
new_mass = np.zeros((n_ages * n_masses))
new_teff = np.zeros((n_ages * n_masses))
new_loglum = np.zeros((n_ages * n_masses))
new_logg = np.zeros((n_ages * n_masses))
new_radius = np.zeros((n_ages * n_masses))
if iso_mag is not None:
new_mag = np.zeros(((n_ages * n_masses, iso_mag.shape[1])))
for j, age_item in enumerate(age_unique):
age_select = iso_age == age_item
ages_tmp = np.full(n_masses, age_item)
new_age[j * n_masses : (j + 1) * n_masses] = ages_tmp
new_mass[j * n_masses : (j + 1) * n_masses] = mass_unique
interp_teff = interpolate.interp1d(
iso_mass[age_select],
iso_teff[age_select],
fill_value="extrapolate",
)
new_teff[j * n_masses : (j + 1) * n_masses] = interp_teff(mass_unique)
interp_loglum = interpolate.interp1d(
iso_mass[age_select],
iso_loglum[age_select],
fill_value="extrapolate",
)
new_loglum[j * n_masses : (j + 1) * n_masses] = interp_loglum(
mass_unique
)
interp_logg = interpolate.interp1d(
iso_mass[age_select],
iso_logg[age_select],
fill_value="extrapolate",
)
new_logg[j * n_masses : (j + 1) * n_masses] = interp_logg(mass_unique)
interp_radius = interpolate.interp1d(
iso_mass[age_select],
iso_radius[age_select],
fill_value="extrapolate",
)
new_radius[j * n_masses : (j + 1) * n_masses] = interp_radius(
mass_unique
)
for k in range(iso_mag.shape[1]):
interp_mag = interpolate.interp1d(
iso_mass[age_select],
iso_mag[age_select, k],
fill_value="extrapolate",
)
new_mag[j * n_masses : (j + 1) * n_masses, k] = interp_mag(
mass_unique
)
iso_age = new_age.copy()
iso_mass = new_mass.copy()
iso_teff = new_teff.copy()
iso_loglum = new_loglum.copy()
iso_logg = new_logg.copy()
iso_radius = new_radius.copy()
if iso_mag is not None:
iso_mag = new_mag.copy()
return (
model_name,
iso_age,
iso_mass,
iso_teff,
iso_loglum,
iso_logg,
iso_radius,
iso_mag,
)
@typechecked
def _check_model(self, atmospheric_model: Optional[str]) -> str:
"""
Internal function for matching the atmospheric model with
the evolutionary model and checking if the expected
atmospheric model is used.
Parameters
----------
atmospheric_model : str, None
Name of the atmospheric model. By setting the argument to
``None``, the atmospheric model associated with the
evolutionary model is automatically selected.
Returns
-------
str
Name of the atmospheric model.
"""
if atmospheric_model is None:
if self.tag in self.match_model:
atmospheric_model = self.match_model[self.tag]
else:
raise ValueError(
"Can not find the atmosphere model "
f"associated with the '{self.tag}' "
"evolutionary model. Please contact "
"the code maintainer."
)
elif self.tag in self.match_model:
if atmospheric_model != self.match_model[self.tag]:
warnings.warn(
"Please note that you have selected "
f"'{atmospheric_model}' as "
f"atmospheric model for '{self.tag}' "
f"while '{self.match_model[self.tag]}'"
" is the atmospheric model that is "
f"self-consistently associated with "
f"'{self.tag}'. It is recommended "
"to set 'atmospheric_model=None' to "
"automatically select the correct "
"grid with model spectra."
)
return atmospheric_model
@typechecked
def _update_param(
self,
atmospheric_model: str,
model_param: Dict[str, float],
param_bounds: Dict[str, Tuple[float, float]],
extra_param: Dict[str, float],
) -> Tuple[Dict[str, float], Dict[str, float]]:
"""
Internal function for updating the dictionary with model
parameters for the atmospheric model. Parameters that are
not part of the evolutionary model but are required for the
atmospheric model will be included by either adopting them
from the ``extra_param`` dictionary or asking for manual
input in case a parameter is missing. In the latter case,
the ``extra_param`` dictionary will also be updated.
Parameters
----------
atmospheric_model : str
Name of the atmospheric model.
model_param : dict
Dictionary with the parameters at which the atmospheric
model spectra will be interpolated.
param_bounds : dict
Dictionary with the parameter boundaries of grid with
atmospheric model spectra.
extra_param : dict
Dictionary with additional parameters that are optionally
required for the atmospheric model but are not part of
the evolutionary model grid. In case additional
parameters are required for the atmospheric model but
they are not provided in ``extra_param`` then a manual
input will be requested when running the
``get_photometry`` method. Typically the ``extra_param``
parameter is not needed so the argument can be set to
``None``. It will only be required if a non-self-consistent
approach will be tested, that is, the calculation of
synthetic photometry from an atmospheric model that is
not associated with the evolutionary model.
Returns
-------
dict
Updated dictionary with parameters for the interpolation
of atmospheric model grid.
dict
Updated dictionary with only the parameters that
are required for the atmospheric model but not for
the evolutionary model.
"""
if len(extra_param) == 0 and self.extra_param is not None:
extra_param = self.extra_param.copy()
for key, value in param_bounds.items():
if key not in model_param:
if key in extra_param:
model_param[key] = extra_param[key]
else:
param_name = update_labels([key])[0]
input_value = input(
f"The '{atmospheric_model}' atmospheric model "
f"requires the '{key}' parameter (i.e. "
f"{param_name}) as input, while it is not "
f"part of the '{self.tag}' evolutionary model. "
"Please provide a value within the available "
f"range from {value[0]} to {value[1]}: "
)
model_param[key] = float(input_value)
extra_param[key] = float(input_value)
return model_param, extra_param
[docs]
@typechecked
def grid_points(self) -> Dict[str, np.ndarray]:
"""
Function for returning a dictionary with the unique grid points
of the parameters of the evolutionary data. The grid may not
be regularly sampled unless ``create_regular_grid=True``.
Returns
-------
dict(str, np.array)
Dictionary with the parameter names and the arrays with
the unique values in the grid of evolutionary data.
"""
(
model,
iso_age,
iso_mass,
iso_teff,
iso_loglum,
iso_logg,
iso_radius,
iso_mag,
) = self._read_data()
grid_points = {
"age": iso_age,
"mass": iso_mass,
"radius": iso_radius,
"log_lum": iso_loglum,
"teff": iso_teff,
"logg": iso_logg,
}
return grid_points
[docs]
@typechecked
def get_isochrone(
self,
age: float,
masses: Optional[np.ndarray] = None,
filter_mag: Optional[str] = None,
filters_color: Optional[Tuple[str, str]] = None,
) -> IsochroneBox:
"""
Function for interpolating an isochrone.
Parameters
----------
age : float
Age (Myr) at which the isochrone data is interpolated.
masses : np.ndarray, None
Masses (:math:`M_\\mathrm{J}`) at which the isochrone
data is interpolated. The masses are not interpolated
if the argument is set to ``None``, in which case the
mass sampling from the evolutionary data is used.
filter_mag : str, None
Filter name for the absolute magnitude as listed in the
file with the isochrone data. Not selected if set to
``None`` or if only evolutionary tracks are available.
filters_color : tuple(str, str), None
Filter names for the color as listed in the file with the
isochrone data. Not selected if set to ``None`` or if only
evolutionary tracks are available.
Returns
-------
species.core.box.IsochroneBox
Box with the isochrone.
"""
color = None
mag_abs = None
# Read isochrone data
(
model,
iso_age,
iso_mass,
iso_teff,
iso_loglum,
iso_logg,
iso_radius,
iso_mag,
) = self._read_data()
if masses is None:
idx_min = (np.abs(iso_age - age)).argmin()
age_select = iso_age == iso_age[idx_min]
masses = np.unique(iso_mass[age_select]) # (Mjup)
age_points = np.full(masses.shape[0], age) # (Myr)
grid_points = np.column_stack([iso_age, iso_mass])
filters = self.get_filters()
if model in self.mag_models:
if filters_color is not None:
if filters_color[0] in filters:
index_color_1 = filters.index(filters_color[0])
else:
raise ValueError(
f"Magnitudes for the selected "
f"'{filters_color[0]}' filter "
f"are not found in the "
f"'{self.tag}' data. Please "
f"select one of the following "
f"filters: {filters}"
)
if filters_color[1] in filters:
index_color_2 = filters.index(filters_color[1])
else:
raise ValueError(
f"Magnitudes for the selected "
f"'{filters_color[1]}' filter "
f"are not found in the "
f"'{self.tag}' data. Please "
f"select one of the following "
f"filters: {filters}"
)
if filter_mag is not None:
if filter_mag in filters:
index_mag = filters.index(filter_mag)
else:
raise ValueError(
f"Magnitudes for the selected "
f"'{filter_mag}' filter are not "
f"found in the '{self.tag}' data. "
f"Please select one of the "
f"following filters: {filters}"
)
if filters_color is not None:
mag_color_1 = interpolate.griddata(
points=grid_points,
values=iso_mag[:, index_color_1],
xi=np.stack((age_points, masses), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
mag_color_2 = interpolate.griddata(
points=grid_points,
values=iso_mag[:, index_color_2],
xi=np.stack((age_points, masses), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
color = mag_color_1 - mag_color_2
if filter_mag is not None:
mag_abs = interpolate.griddata(
points=grid_points,
values=iso_mag[:, index_mag],
xi=np.stack((age_points, masses), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
teff = interpolate.griddata(
points=grid_points,
values=iso_teff,
xi=np.stack((age_points, masses), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
log_lum = interpolate.griddata(
points=grid_points,
values=iso_loglum,
xi=np.stack((age_points, masses), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
logg = interpolate.griddata(
points=grid_points,
values=iso_logg,
xi=np.stack((age_points, masses), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
radius = interpolate.griddata(
points=grid_points,
values=iso_radius,
xi=np.stack((age_points, masses), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
if mag_abs is None and filter_mag is not None:
warnings.warn(
f"The isochrones of {self.tag} do not have "
f"magnitudes for the {filter_mag} filter so "
f"setting the argument of 'filter_mag' to None."
)
filter_mag = None
if color is None and filters_color is not None:
warnings.warn(
f"The isochrones of {self.tag} do not have "
f"magnitudes for the {filters_color} filters so "
f"setting the argument of 'filter_color' to None."
)
filters_color = None
return create_box(
boxtype="isochrone",
model=self.tag,
age=age,
filters_color=filters_color,
filter_mag=filter_mag,
color=color,
magnitude=mag_abs,
log_lum=log_lum,
teff=teff,
logg=logg,
radius=radius,
masses=masses,
)
[docs]
@typechecked
def get_cooling_curve(
self,
mass: float,
ages: Optional[np.ndarray] = None,
filter_mag: Optional[str] = None,
filters_color: Optional[Tuple[str, str]] = None,
) -> CoolingBox:
"""
Function for interpolating a cooling curve.
Parameters
----------
mass : float
Mass (:math:`M_\\mathrm{J}`) for which the cooling
curve will be interpolated.
ages : np.ndarray, None
Ages (Myr) at which the cooling curve will be
interpolated. The ages are not interpolated
if the argument is set to ``None``, in which case the
age sampling from the evolutionary data is used.
filter_mag : str, None
Filter name for the absolute magnitude as listed in the
file with the isochrone data. Not selected if set to
``None`` or if only evolutionary tracks are available.
filters_color : tuple(str, str), None
Filter names for the color as listed in the file with the
isochrone data. Not selected if set to ``None`` or if only
evolutionary tracks are available.
Returns
-------
species.core.box.CoolingBox
Box with the cooling curve.
"""
color = None
mag_abs = None
# Read isochrone data
(
model,
iso_age,
iso_mass,
iso_teff,
iso_loglum,
iso_logg,
iso_radius,
iso_mag,
) = self._read_data()
if ages is None:
idx_min = (np.abs(iso_mass - mass)).argmin()
mass_select = iso_mass == iso_mass[idx_min]
ages = np.unique(iso_age[mass_select]) # (Myr)
mass_points = np.full(ages.shape[0], mass) # (Mjup)
grid_points = np.column_stack([iso_age, iso_mass])
filters = self.get_filters()
if model in self.mag_models:
if filters_color is not None:
index_color_1 = filters.index(filters_color[0])
index_color_2 = filters.index(filters_color[1])
if filter_mag is not None:
index_mag = filters.index(filter_mag)
if filters_color is not None:
mag_color_1 = interpolate.griddata(
points=grid_points,
values=iso_mag[:, index_color_1],
xi=np.stack((ages, mass_points), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
mag_color_2 = interpolate.griddata(
points=grid_points,
values=iso_mag[:, index_color_2],
xi=np.stack((ages, mass_points), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
color = mag_color_1 - mag_color_2
if filter_mag is not None:
mag_abs = interpolate.griddata(
points=grid_points,
values=iso_mag[:, index_mag],
xi=np.stack((ages, mass_points), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
teff = interpolate.griddata(
points=grid_points,
values=iso_teff,
xi=np.stack((ages, mass_points), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
log_lum = interpolate.griddata(
points=grid_points,
values=iso_loglum,
xi=np.stack((ages, mass_points), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
logg = interpolate.griddata(
points=grid_points,
values=iso_logg,
xi=np.stack((ages, mass_points), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
radius = interpolate.griddata(
points=grid_points,
values=iso_radius,
xi=np.stack((ages, mass_points), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
if mag_abs is None and filter_mag is not None:
warnings.warn(
f"The isochrones of {self.tag} do not have "
f"magnitudes for the {filter_mag} filter so "
f"setting the argument of 'filter_mag' to None."
)
filter_mag = None
if mag_abs is None and filters_color is not None:
warnings.warn(
f"The isochrones of {self.tag} do not have "
f"magnitudes for the {filters_color} filters so "
f"setting the argument of 'filter_color' to None."
)
filters_color = None
return create_box(
boxtype="cooling",
model=self.tag,
mass=mass,
filters_color=filters_color,
filter_mag=filter_mag,
color=color,
magnitude=mag_abs,
ages=ages,
log_lum=log_lum,
teff=teff,
logg=logg,
radius=radius,
)
[docs]
@typechecked
def get_color_magnitude(
self,
age: float,
masses: Optional[np.ndarray],
filters_color: Tuple[str, str],
filter_mag: str,
adapt_logg: bool = False,
atmospheric_model: Optional[str] = None,
extra_param: Optional[Dict[str, float]] = None,
) -> ColorMagBox:
"""
Function for calculating color-magnitude pairs
from a selected isochrone. The function selects the
corresponding atmosphere model and computes synthetic
photometry by interpolating and integrating the
spectra for any given filters.
Parameters
----------
age : float
Age (Myr) at which the isochrone data is interpolated.
masses : np.ndarray, None
Masses (:math:`M_\\mathrm{J}`) at which the isochrone
data is interpolated. The masses at the nearest age
in the grid with evolutionary data are selected if
the argument is set to ``None``.
filters_color : tuple(str, str)
Filter names for the color as listed in the file with the
isochrone data. The filter names should be provided in the
format of the SVO Filter Profile Service.
filter_mag : str
Filter name for the absolute magnitude as listed in the
file with the isochrone data. The value should be equal
to one of the ``filters_color`` values.
adapt_logg : bool
Adapt :math:`\\log(g)` to the upper or lower boundary of
the atmospheric model grid whenever the :math:`\\log(g)`
that has been calculated from the isochrone mass and
radius lies outside the available range of the synthetic
spectra. Typically :math:`\\log(g)` has only a minor
impact on the broadband magnitudes and colors.
atmospheric_model : str, None
Atmospheric model used to compute the synthetic photometry.
The argument can be set to ``None`` such that the correct
atmospheric model is automatically selected that is
associated with the evolutionary model. If the user
nonetheless wants to test a non-self-consistent approach
by using a different atmospheric model, then the argument
can be set to any of the models that can be added with
:func:`~species.data.database.Database.add_model`.
extra_param : dict, None
Optional dictionary with additional parameters that are
required for the atmospheric model but are not part of
the evolutionary model grid. In case additional
parameters are required for the atmospheric model but
they are not provided in ``extra_param`` then a manual
input will be requested when running the
``get_photometry`` method. Typically the ``extra_param``
parameter is not needed so the argument can be set to
``None``. It will only be required if a non-self-consistent
approach will be tested, that is, the calculation of
synthetic photometry from an atmospheric model that is
not associated with the evolutionary model.
Returns
-------
species.core.box.ColorMagBox
Box with the color-magnitude data.
"""
if extra_param is None:
extra_param = {}
atmospheric_model = self._check_model(atmospheric_model)
isochrone = self.get_isochrone(
age=age, masses=masses, filters_color=None, filter_mag=None
)
model_1 = ReadModel(model=atmospheric_model, filter_name=filters_color[0])
model_2 = ReadModel(model=atmospheric_model, filter_name=filters_color[1])
param_bounds = model_1.get_bounds()
mag1 = np.zeros(isochrone.mass.shape[0])
mag2 = np.zeros(isochrone.mass.shape[0])
for i in range(isochrone.mass.size):
model_param = {
"teff": isochrone.teff[i],
"logg": isochrone.logg[i],
"radius": isochrone.radius[i],
"distance": 10.0,
}
if atmospheric_model == "sonora-bobcat":
model_param["feh"] = float(self.tag[-4:])
# The get_bounds of model_1 and model_2 are the
# same since the same atmospheric_model is used
model_param, extra_param = self._update_param(
atmospheric_model, model_param, model_1.get_bounds(), extra_param
)
if np.isnan(isochrone.teff[i]):
mag1[i] = np.nan
mag2[i] = np.nan
warnings.warn(
f"The value of Teff is NaN for the following "
f"isochrone sample: {model_param}. Setting "
f"the magnitudes to NaN."
)
else:
for item_bounds in param_bounds:
if model_param[item_bounds] < param_bounds[item_bounds][0]:
if adapt_logg and item_bounds == "logg":
warnings.warn(
f"The log(g) is {model_param[item_bounds]} but the "
f"lower boundary of the model grid is "
f"{param_bounds[item_bounds][0]}. Adapting "
f"log(g) to {param_bounds[item_bounds][0]} since "
f"adapt_logg=True."
)
model_param["logg"] = param_bounds["logg"][0]
else:
mag1[i] = np.nan
mag2[i] = np.nan
warnings.warn(
f"The value of {item_bounds} is "
f"{model_param[item_bounds]}, which is below "
f"the lower bound of the model grid "
f"({param_bounds[item_bounds][0]}). Setting the "
f"magnitudes to NaN for the following isochrone "
f"sample: {model_param}."
)
elif model_param[item_bounds] > param_bounds[item_bounds][1]:
if adapt_logg and item_bounds == "logg":
warnings.warn(
f"The log(g) is {model_param[item_bounds]} but "
f"the upper boundary of the model grid is "
f"{param_bounds[item_bounds][1]}. Adapting "
f"log(g) to {param_bounds[item_bounds][1]} "
f"since adapt_logg=True."
)
model_param["logg"] = param_bounds["logg"][1]
else:
mag1[i] = np.nan
mag2[i] = np.nan
warnings.warn(
f"The value of {item_bounds} is "
f"{model_param[item_bounds]}, which is above "
f"the upper bound of the model grid "
f"({param_bounds[item_bounds][1]}). Setting the "
f"magnitudes to NaN for the following isochrone "
f"sample: {model_param}."
)
if not np.isnan(mag1[i]):
mag1[i], _ = model_1.get_magnitude(model_param)
mag2[i], _ = model_2.get_magnitude(model_param)
if filter_mag == filters_color[0]:
abs_mag = mag1
elif filter_mag == filters_color[1]:
abs_mag = mag2
else:
raise ValueError(
"The argument of filter_mag should be equal to "
"one of the two filter values of filters_color."
)
return create_box(
boxtype="colormag",
library=atmospheric_model,
object_type="model",
filters_color=filters_color,
filter_mag=filter_mag,
color=mag1 - mag2,
magnitude=abs_mag,
mass=isochrone.mass,
radius=isochrone.radius,
iso_tag=self.tag,
)
[docs]
@typechecked
def get_color_color(
self,
age: float,
masses: Optional[np.ndarray],
filters_colors: Tuple[Tuple[str, str], Tuple[str, str]],
atmospheric_model: Optional[str] = None,
extra_param: Optional[Dict[str, float]] = None,
) -> ColorColorBox:
"""
Function for calculating color-color pairs
from a selected isochrone. The function selects the
corresponding atmosphere model and computes synthetic
photometry by interpolating and integrating the spectra
for any given filters.
Parameters
----------
age : float
Age (Myr) at which the isochrone data is interpolated.
masses : np.ndarray, None
Masses (:math:`M_\\mathrm{J}`) at which the isochrone
data is interpolated. The masses at the nearest age
in the grid with evolutionary data are selected if
the argument is set to ``None``.
filters_colors : tuple(tuple(str, str), tuple(str, str))
Filter names for the colors as listed in the file with the
isochrone data. The filter names should be provided in the
format of the SVO Filter Profile Service.
atmospheric_model : str, None
Atmospheric model used to compute the synthetic photometry.
The argument can be set to ``None`` such that the correct
atmospheric model is automatically selected that is
associated with the evolutionary model. If the user
nonetheless wants to test a non-self-consistent approach
by using a different atmospheric model, then the argument
can be set to any of the models that can be added with
:func:`~species.data.database.Database.add_model`.
extra_param : dict, None
Optional dictionary with additional parameters that are
required for the atmospheric model but are not part of
the evolutionary model grid. In case additional
parameters are required for the atmospheric model but
they are not provided in ``extra_param`` then a manual
input will be requested when running the
``get_photometry`` method. Typically the ``extra_param``
parameter is not needed so the argument can be set to
``None``. It will only be required if a non-self-consistent
approach will be tested, that is, the calculation of
synthetic photometry from an atmospheric model that is
not associated with the evolutionary model.
Returns
-------
species.core.box.ColorColorBox
Box with the color-color data.
"""
if extra_param is None:
extra_param = {}
atmospheric_model = self._check_model(atmospheric_model)
isochrone = self.get_isochrone(
age=age, masses=masses, filters_color=None, filter_mag=None
)
model_1 = ReadModel(model=atmospheric_model, filter_name=filters_colors[0][0])
model_2 = ReadModel(model=atmospheric_model, filter_name=filters_colors[0][1])
model_3 = ReadModel(model=atmospheric_model, filter_name=filters_colors[1][0])
model_4 = ReadModel(model=atmospheric_model, filter_name=filters_colors[1][1])
mag1 = np.zeros(isochrone.mass.shape[0])
mag2 = np.zeros(isochrone.mass.shape[0])
mag3 = np.zeros(isochrone.mass.shape[0])
mag4 = np.zeros(isochrone.mass.shape[0])
for i in range(isochrone.mass.size):
model_param = {
"teff": isochrone.teff[i],
"logg": isochrone.logg[i],
"radius": isochrone.radius[i],
"distance": 10.0,
}
if atmospheric_model == "sonora-bobcat":
model_param["feh"] = float(self.tag[-4:])
# The get_bounds of model_1 and model_2 are the
# same since the same atmospheric_model is used
model_param, extra_param = self._update_param(
atmospheric_model, model_param, model_1.get_bounds(), extra_param
)
if np.isnan(isochrone.teff[i]):
mag1[i] = np.nan
mag2[i] = np.nan
mag3[i] = np.nan
mag4[i] = np.nan
warnings.warn(
"The value of Teff is NaN for the following "
f"isochrone sample: {model_param}. Setting "
"the magnitudes to NaN."
)
else:
for item_bounds in model_1.get_bounds():
if model_param[item_bounds] < model_1.get_bounds()[item_bounds][0]:
mag1[i] = np.nan
mag2[i] = np.nan
mag3[i] = np.nan
mag4[i] = np.nan
warnings.warn(
f"The value of {item_bounds} is "
f"{model_param[item_bounds]}, which is "
f"below the lower bound of the model grid "
f" ({model_1.get_bounds()[item_bounds][0]}). "
f"Setting the magnitudes to NaN for the "
f"following isochrone sample: {model_param}."
)
elif (
model_param[item_bounds] > model_1.get_bounds()[item_bounds][1]
):
mag1[i] = np.nan
mag2[i] = np.nan
mag3[i] = np.nan
mag4[i] = np.nan
warnings.warn(
f"The value of {item_bounds} is "
f"{model_param[item_bounds]}, which is above "
f"the upper bound of the model grid "
f"({model_1.get_bounds()[item_bounds][1]}). "
f"Setting the magnitudes to NaN for the "
f"following isochrone sample: {model_param}."
)
if (
not np.isnan(mag1[i])
and not np.isnan(mag2[i])
and not np.isnan(mag3[i])
and not np.isnan(mag4[i])
):
mag1[i], _ = model_1.get_magnitude(model_param)
mag2[i], _ = model_2.get_magnitude(model_param)
mag3[i], _ = model_3.get_magnitude(model_param)
mag4[i], _ = model_4.get_magnitude(model_param)
return create_box(
boxtype="colorcolor",
library=atmospheric_model,
object_type="model",
filters=filters_colors,
color1=mag1 - mag2,
color2=mag3 - mag4,
mass=isochrone.mass,
radius=isochrone.radius,
iso_tag=self.tag,
)
[docs]
@typechecked
def get_mass(
self,
age: float,
log_lum: np.ndarray,
) -> np.ndarray:
"""
Function for interpolating a mass for a given
age and array with bolometric luminosities.
Parameters
----------
age : float
Age (Myr) at which the masses will be interpolated.
log_lum : np.ndarray
Array with the bolometric luminosities,
:math:`\\log{(L/L_\\odot)}`, for which the
masses will be interpolated.
Returns
-------
np.ndarray
Array with masses (:math:`M_\\mathrm{J}`).
"""
# Read isochrone data
(
_,
iso_age,
iso_mass,
_,
iso_loglum,
_,
_,
_,
) = self._read_data()
# Interpolate masses
grid_points = np.stack((iso_age, iso_loglum), axis=1)
age_points = np.full(log_lum.size, age) # (Myr)
mass = interpolate.griddata(
points=grid_points,
values=iso_mass,
xi=np.stack((age_points, log_lum), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
return mass
[docs]
@typechecked
def get_radius(
self,
age: float,
log_lum: np.ndarray,
) -> np.ndarray:
"""
Function for interpolating a radius for a given
age and array with bolometric luminosities.
Parameters
----------
age : float
Age (Myr) at which the masses will be interpolated.
log_lum : np.ndarray
Array with the bolometric luminosities,
:math:`\\log{(L/L_\\odot)}`, for which the
masses will be interpolated.
Returns
-------
np.ndarray
Array with radii (:math:`R_\\mathrm{J}`).
"""
# Read isochrone data
(
_,
iso_age,
_,
_,
iso_loglum,
_,
iso_radius,
_,
) = self._read_data()
# Interpolate radius
grid_points = np.stack((iso_age, iso_loglum), axis=1)
age_points = np.full(log_lum.size, age) # (Myr)
radius = interpolate.griddata(
points=grid_points,
values=iso_radius,
xi=np.stack((age_points, log_lum), axis=1),
method="linear",
fill_value="nan",
rescale=False,
)
return radius
[docs]
@typechecked
def get_filters(self) -> Optional[List[str]]:
"""
Function for get a list with filter names for which there
are are magnitudes stored with the isochrone data.
Returns
-------
list(str), None
List with filter names. A ``None`` is returned if
there are no filters and magnitudes stored with
the isochrone data.
"""
with h5py.File(self.database, "r") as hdf5_file:
if "filters" in hdf5_file[f"isochrones/{self.tag}"]:
filters = list(hdf5_file[f"isochrones/{self.tag}/filters/"])
# Convert from bytes to strings
for i, item in enumerate(filters):
if isinstance(item, bytes):
filters[i] = item.decode("utf-8")
else:
filters = None
return filters
[docs]
@typechecked
def get_photometry(
self,
age: float,
mass: float,
distance: float,
filter_name: str,
atmospheric_model: Optional[str] = None,
extra_param: Optional[Dict[str, float]] = None,
) -> PhotometryBox:
"""
Function for computing synthetic photometry by interpolating
and integrating the associated spectra. Bulk and atmosphere
parameters are interpolated from the evolutionary data for
the requested age and mass. The output from the evolutionary
data is then used as input for the atmospheric model. This
function is useful if the required magnitudes or fluxes
are not part of the available filters of the evolutionary
data (i.e. the filters returned by
:func:`~species.read.read_isochrone.ReadIsochrone.get_filters`).
The atmospheric model that is associated with the evolutionary
model is by default automatically selected and added to the
database if needed.
Parameters
----------
age : float
Age (Myr) at which the bulk parameters will be
interpolated from the grid with evolutionary data.
mass : float
Mass (:math:`M_\\mathrm{J}`) at which the bulk
parameters will be interpolated from the grid with
evolutionary data.
distance : float
Distance (pc) that is used for scaling the fluxes
from the atmosphere to the observer.
filter_name : tuple(str, str), None
Filter name for which the synthetic photometry will be
computed. Any filter name from the `SVO Filter Profile
Service <http://svo2.cab.inta-csic.es/svo/theory/fps/>`_
can be used as argument.
atmospheric_model : str, None
Atmospheric model used to compute the synthetic photometry.
The argument can be set to ``None`` such that the correct
atmospheric model is automatically selected that is
associated with the evolutionary model. If the user
nonetheless wants to test a non-self-consistent approach
by using a different atmospheric model, then the argument
can be set to any of the models that can be added with
:func:`~species.data.database.Database.add_model`.
extra_param : dict, None
Optional dictionary with additional parameters that are
required for the atmospheric model but are not part of
the evolutionary model grid. In case additional
parameters are required for the atmospheric model but
they are not provided in ``extra_param`` then a manual
input will be requested when running the
``get_photometry`` method. Typically the ``extra_param``
parameter is not needed so the argument can be set to
``None``. It will only be required if a non-self-consistent
approach will be tested, that is, the calculation of
synthetic photometry from an atmospheric model that is
not associated with the evolutionary model.
Returns
-------
species.core.box.PhotometryBox
Box with the synthetic photometry (magnitude and flux).
"""
if extra_param is None:
extra_param = {}
atmospheric_model = self._check_model(atmospheric_model)
iso_box = self.get_isochrone(age=age, masses=np.array([mass]))
model_param = {
"teff": iso_box.teff[0],
"logg": iso_box.logg[0],
"radius": iso_box.radius[0],
"distance": distance,
}
model_reader = ReadModel(model=atmospheric_model, filter_name=filter_name)
model_param, _ = self._update_param(
atmospheric_model, model_param, model_reader.get_bounds(), extra_param
)
phot_box = model_reader.get_flux(model_param=model_param, return_box=True)
return phot_box
[docs]
@typechecked
def get_spectrum(
self,
age: float,
mass: float,
distance: float,
wavel_range: Optional[Tuple[float, float]] = None,
spec_res: Optional[float] = None,
atmospheric_model: Optional[str] = None,
extra_param: Optional[Dict[str, float]] = None,
) -> ModelBox:
"""
Function for interpolating the model spectrum at a specified
age and mass. Bulk and atmosphere parameters are interpolated
from the evolutionary data for the requested age and mass.
The output from the evolutionary data is then used as input
for the atmospheric model. The atmospheric model that is
associated with the evolutionary model is by default
automatically selected and added to the database if needed.
Parameters
----------
age : float
Age (Myr) at which the bulk parameters will be
interpolated from the grid with evolutionary data.
mass : float
Mass (:math:`M_\\mathrm{J}`) at which the bulk
parameters will be interpolated from the grid with
evolutionary data.
distance : float
Distance (pc) that is used for scaling the fluxes
from the atmosphere to the observer.
wavel_range : tuple(float, float), None
Wavelength range (um). Full spectrum is selected if
the argument is set to ``None``.
spec_res : float, None
Spectral resolution that is used for smoothing the spectrum
with a Gaussian kernel. No smoothing is applied when the
argument is set to ``None``.
atmospheric_model : str, None
Atmospheric model used to compute the synthetic photometry.
The argument can be set to ``None`` such that the correct
atmospheric model is automatically selected that is
associated with the evolutionary model. If the user
nonetheless wants to test a non-self-consistent approach
by using a different atmospheric model, then the argument
can be set to any of the models that can be added with
:func:`~species.data.database.Database.add_model`.
extra_param : dict, None
Optional dictionary with additional parameters that are
required for the atmospheric model but are not part of
the evolutionary model grid. In case additional
parameters are required for the atmospheric model but
they are not provided in ``extra_param`` then a manual
input will be requested when running the
``get_photometry`` method. Typically the ``extra_param``
parameter is not needed so the argument can be set to
``None``. It will only be required if a non-self-consistent
approach will be tested, that is, the calculation of
synthetic photometry from an atmospheric model that is
not associated with the evolutionary model.
Returns
-------
species.core.box.ModelBox
Box with the model spectrum.
"""
if extra_param is None:
extra_param = {}
atmospheric_model = self._check_model(atmospheric_model)
iso_box = self.get_isochrone(age=age, masses=np.array([mass]))
model_param = {
"teff": iso_box.teff[0],
"logg": iso_box.logg[0],
"radius": iso_box.radius[0],
"distance": distance,
}
model_reader = ReadModel(model=atmospheric_model, wavel_range=wavel_range)
model_param, _ = self._update_param(
atmospheric_model, model_param, model_reader.get_bounds(), extra_param
)
model_box = model_reader.get_model(
model_param=model_param, spec_res=spec_res, smooth=True
)
return model_box
[docs]
@typechecked
def contrast_to_mass(
self,
age: float,
distance: float,
filter_name: str,
star_mag: float,
contrast: Union[List[float], np.ndarray],
use_mag: bool = True,
atmospheric_model: Optional[str] = None,
extra_param: Optional[Dict[str, float]] = None,
calc_phot: bool = False,
) -> np.ndarray:
"""
Function for converting contrast values into masses. This
can be used to convert a list/array with detection
limits of companions into mass limits. Either
one of the available filter names from the isochrone grid
can be selected (i.e. the filters returned by
:func:`~species.read.read_isochrone.ReadIsochrone.get_filters`),
or any of the filters from the `SVO Filter Profile
Service <http://svo2.cab.inta-csic.es/svo/theory/fps/>`_. For
the first case, the magnitudes will be directly interpolated
from the grid of evolution data. For the second case, the
associated model spectra will be used for calculating
synthetic photometry for the isochrone age and selected
filter. These will then be interpolated to the requested
contrast values. The atmospheric model that is associated
with the evolutionary model is by default automatically
selected and added to the database if needed.
Parameters
----------
age : float
Age (Myr) at which the bulk parameters will be
interpolated from the grid with evolutionary data.
distance : float
Distance (pc) that is used for scaling the fluxes
from the atmosphere to the observer.
filter_name : str
Filter name for which the magnitudes will be interpolated,
either directly from the isochrone grid or by calculating
synthetic photometry from the associated model spectra.
The first case only works for the filters that are
returned by the
:func:`~species.read.read_isochrone.ReadIsochrone.get_filters`
method of :class:`~species.read.read_isochrone.ReadIsochrone`
because these will have pre-calculated magnitudes. The
second case will work for any of the filter names from the
`SVO Filter Profile Service
<http://svo2.cab.inta-csic.es/svo/theory/fps/>`_. This will
require more disk space and a bit more computation time.
star_mag : float
Stellar apparent magnitude for the filter that is set
as argument of `filter_name`.
contrast : list(float), np.ndarray
List or array with the contrast values between a companion
and the star. The magnitude of the star should be provided
as argument of ``star_mag``. The contrast values will be
converted into masses, while taking into account the
stellar magnitude. The values should be provided
either as ratio (e.g. ``[1e-2, 1e-3, 1e-4]``) or as
magnitudes (e.g. ``[5.0, 7.5, 10.0]``). For ratios,
it is important to set ``use_mag=False``.
use_mag : bool
Set to ``True`` if the values of ``contrast`` are given as
magnitudes. Set to ``False`` if the values of ``contrast``
are given as ratios. The default is set to ``True``.
atmospheric_model : str, None
Atmospheric model used to compute the synthetic photometry
in case the ``filter_name`` is set to a value from the
SVO Filter Profile Service. The argument can be set to
``None`` such that the correct atmospheric model is
automatically selected that is associated with the
evolutionary model. If the user nonetheless wants to test
a non-self-consistent approach by using a different
atmospheric model, then the argument can be set to any of
the models that can be added with
:func:`~species.data.database.Database.add_model`.
extra_param : dict, None
Optional dictionary with additional parameters that are
required for the atmospheric model but are not part of
the evolutionary model grid, for example because they
were implicitly set by the evolution model (e.g.
solar metallicity). In case additional parameters are
required for the atmospheric model but they are not
provided in ``extra_param`` then a manual input will
be requested when running the ``get_photometry`` method.
Typically the ``extra_param`` parameter is not needed so
the argument can be set to ``None``. It will only be
required if a non-self-consistent approach will be tested,
that is, the calculation of synthetic photometry from an
atmospheric model that is not associated with the
evolutionary model.
calc_phot : bool
Calculate synthetic photometry from the model spectra
regardless if pre-calculated magnitudes for the
``filter_name`` are already available with the isochrone
data. Typically the argument can be set to ``False``,
but to force the calculation of synthetic photometry the
argument can be set to ``True``.
Returns
-------
np.ndarray
Array with the masses (in :math:`M_\\mathrm{J}`) for the
requested contrast values.
"""
if isinstance(contrast, list):
contrast = np.array(contrast)
if use_mag and np.all(contrast < 1.0):
warnings.warn(
"All values in the array of 'contrast' are "
"smaller than 1.0 but the argument of "
"'use_mag' is set to True. Please set the "
"argument of 'magnitude' to False in case "
"the values of 'contrast' are given as "
"ratios instead of magnitudes."
)
if not use_mag:
# Convert contrast from ratio to magnitude
contrast = -2.5 * np.log10(contrast)
if extra_param is None:
extra_param = {}
atmospheric_model = self._check_model(atmospheric_model)
app_mag = star_mag + contrast
abs_mag = apparent_to_absolute((app_mag, None), (distance, None))[0]
filter_list = self.get_filters()
if filter_name in filter_list and not calc_phot:
print(
f"The '{filter_name}' filter is found in the list "
"of available filters from the isochrone data of "
f"'{self.tag}'.\nThe requested contrast values "
"will be directly interpolated from the grid with "
"pre-calculated magnitudes."
)
iso_box = self.get_isochrone(age=age, masses=None, filter_mag=filter_name)
# x (=iso_box.magnitude) must be increasing in np.interp
mass_array = np.interp(
abs_mag,
iso_box.magnitude[::-1],
iso_box.mass[::-1],
left=np.nan,
right=np.nan,
)
else:
if not filter_name in filter_list:
print(
f"The '{filter_name}' filter is not found in the "
"list of available filters from the isochrone "
f"data of '{self.tag}'.\nIt will be tried to "
"download the filter profile (if needed) and to "
"use the associated atmospheric model spectra "
"for calculating synthetic photometry."
)
model_reader = ReadModel(model=atmospheric_model, filter_name=filter_name)
iso_box = self.get_isochrone(age=age, masses=None, filter_mag=None)
model_abs_mag = np.zeros(iso_box.mass.size)
for i in range(iso_box.mass.size):
model_param = {
"teff": iso_box.teff[i],
"logg": iso_box.logg[i],
"radius": iso_box.radius[i],
"distance": distance,
}
model_param, _ = self._update_param(
atmospheric_model,
model_param,
model_reader.get_bounds(),
extra_param,
)
# The get_magnitude method returns the
# apparent magnitude and absolute magnitude
_, model_abs_mag[i] = model_reader.get_magnitude(
model_param=model_param, return_box=False
)
# x (=model_abs_mag) must be increasing in np.interp
mass_array = np.interp(
abs_mag,
model_abs_mag[::-1],
iso_box.mass[::-1],
left=np.nan,
right=np.nan,
)
return mass_array