Source code for species.read.read_model

"""
Module with reading functionalities for atmospheric model spectra.
"""

import os
import warnings

from configparser import ConfigParser
from typing import Dict, List, Optional, Tuple, Union

import h5py
import numpy as np

from PyAstronomy.pyasl import rotBroad, fastRotBroad
from typeguard import typechecked
from scipy.integrate import simps
from scipy.interpolate import interp1d, RegularGridInterpolator

from species.core import constants
from species.core.box import (
    ColorColorBox,
    ColorMagBox,
    ModelBox,
    PhotometryBox,
    create_box,
)
from species.data.model_data.model_spectra import add_model_grid
from species.data.spec_data.spec_vega import add_vega
from species.phot.syn_phot import SyntheticPhotometry
from species.read.read_filter import ReadFilter
from species.read.read_planck import ReadPlanck
from species.util.convert_util import logg_to_mass
from species.util.dust_util import check_dust_database, ism_extinction, convert_to_av
from species.util.model_util import binary_to_single
from species.util.spec_util import smooth_spectrum


[docs] class ReadModel: """ Class for reading a model spectrum from the database. """ @typechecked def __init__( self, model: str, wavel_range: Optional[Tuple[float, float]] = None, filter_name: Optional[str] = None, ): """ Parameters ---------- model : str Name of the atmospheric model. wavel_range : tuple(float, float), None Wavelength range (um). Full spectrum is selected if set to ``None``. Not used if ``filter_name`` is not ``None``. filter_name : str, None Filter name that is used for the wavelength range. The ``wavel_range`` is used if set to ``None``. Returns ------- NoneType None """ self.model = model if self.model == "bt-settl": warnings.warn( "It is recommended to use the CIFIST " "grid of the BT-Settl, because it is " "a newer version. In that case, set " "model='bt-settl-cifist' when using " "add_model of Database." ) self.spectrum_interp = None self.wl_points = None self.wl_index = None self.filter_name = filter_name self.wavel_range = wavel_range if self.filter_name is not None: transmission = ReadFilter(self.filter_name) self.wavel_range = transmission.wavelength_range() self.mean_wavelength = transmission.mean_wavelength() else: self.mean_wavelength = None config_file = os.path.join(os.getcwd(), "species_config.ini") config = ConfigParser() config.read(config_file) self.database = config["species"]["database"] self.data_folder = config["species"]["data_folder"] self.extra_param = [ "radius", "distance", "parallax", "mass", "luminosity", "lognorm_radius", "lognorm_sigma", "lognorm_ext", "ism_ext", "ism_red", "powerlaw_max", "powerlaw_exp", "powerlaw_ext", "disk_teff", "disk_radius", "veil_a", "veil_b", "veil_ref", "vsini", ] # Test if the spectra are present in the database hdf5_file = self.open_database() hdf5_file.close()
[docs] @typechecked def open_database(self) -> h5py._hl.files.File: """ Internal function for opening the HDF5 database. Returns ------- h5py._hl.files.File The HDF5 database. """ with h5py.File(self.database, "r") as hdf5_file: if f"models/{self.model}" in hdf5_file: model_found = True else: model_found = False warnings.warn( f"The '{self.model}' model spectra are not present " "in the database. Will try to add the model grid. " "If this does not work (e.g. currently without an " "internet connection), then please use the " "'add_model' method of 'Database' to add the " "grid of spectra at a later moment." ) if not model_found: # This will not work when using multiprocessing. # Model spectra should be added to the database # before running FitModel with MPI with h5py.File(self.database, "a") as hdf5_file: add_model_grid(self.model, self.data_folder, hdf5_file) return h5py.File(self.database, "r")
[docs] @typechecked def wavelength_points(self) -> Tuple[np.ndarray, np.ndarray]: """ Internal function for extracting the wavelength points and indices that are used. Returns ------- np.ndarray Wavelength points (um). np.ndarray Array with the size of the original wavelength grid. The booleans indicate if a wavelength point was used. """ with self.open_database() as hdf5_file: wl_points = np.array(hdf5_file[f"models/{self.model}/wavelength"]) if self.wavel_range is None: wl_index = np.ones(wl_points.shape[0], dtype=bool) else: wl_index = (wl_points > self.wavel_range[0]) & ( wl_points < self.wavel_range[1] ) index = np.where(wl_index)[0] if index[0] - 1 >= 0: wl_index[index[0] - 1] = True if index[-1] + 1 < wl_index.size: wl_index[index[-1] + 1] = True return wl_points[wl_index], wl_index
[docs] @typechecked def interpolate_model(self) -> None: """ Internal function for linearly interpolating the full grid of model spectra. Returns ------- NoneType None """ grid_points = list(self.get_points().values()) if self.wl_points is None: self.wl_points, self.wl_index = self.wavelength_points() with self.open_database() as hdf5_file: grid_flux = np.array(hdf5_file[f"models/{self.model}/flux"]) grid_flux = grid_flux[..., self.wl_index] self.spectrum_interp = RegularGridInterpolator( grid_points, grid_flux, method="linear", bounds_error=False, fill_value=np.nan, )
[docs] @typechecked def interpolate_grid( self, wavel_resample: Optional[np.ndarray] = None, spec_res: Optional[float] = None, **kwargs, ) -> None: """ Internal function for linearly interpolating the grid of model spectra for a given filter or wavelength sampling. wavel_resample : np.ndarray, None Wavelength points for the resampling of the spectrum. The ``filter_name`` is used if set to ``None``. spec_res : float, None Spectral resolution that is used for the Gaussian filter. No smoothing is applied if the argument is set to ``None``. Returns ------- NoneType None """ self.interpolate_model() if "smooth" in kwargs: warnings.warn( "The 'smooth' parameter has been " "deprecated. Please set only the " "'spec_res' argument, which can be set " "to None for not applying a smoothing.", DeprecationWarning, ) if not kwargs["smooth"] and spec_res is not None: spec_res = None points = [] for item in self.get_points().values(): points.append(list(item)) param_list = self.get_parameters() n_param = len(param_list) dim_size = [] for item in points: dim_size.append(len(item)) if self.filter_name is not None: dim_size.append(1) else: dim_size.append(wavel_resample.size) flux_new = np.zeros(dim_size) if n_param == 1: model_param = {} for i, item_0 in enumerate(points[0]): model_param[param_list[0]] = item_0 if self.filter_name is not None: flux_new[i] = self.get_flux(model_param)[0] else: flux_new[i, :] = self.get_model( model_param, spec_res=spec_res, wavel_resample=wavel_resample, ).flux elif n_param == 2: model_param = {} for i, item_0 in enumerate(points[0]): for j, item_1 in enumerate(points[1]): model_param[param_list[0]] = item_0 model_param[param_list[1]] = item_1 if self.filter_name is not None: flux_new[i, j] = self.get_flux(model_param)[0] else: flux_new[i, j, :] = self.get_model( model_param, spec_res=spec_res, wavel_resample=wavel_resample, ).flux elif n_param == 3: model_param = {} for i, item_0 in enumerate(points[0]): for j, item_1 in enumerate(points[1]): for k, item_2 in enumerate(points[2]): model_param[param_list[0]] = item_0 model_param[param_list[1]] = item_1 model_param[param_list[2]] = item_2 if self.filter_name is not None: flux_new[i, j, k] = self.get_flux(model_param)[0] else: flux_new[i, j, k, :] = self.get_model( model_param, spec_res=spec_res, wavel_resample=wavel_resample, ).flux elif n_param == 4: model_param = {} for i, item_0 in enumerate(points[0]): for j, item_1 in enumerate(points[1]): for k, item_2 in enumerate(points[2]): for m, item_3 in enumerate(points[3]): model_param[param_list[0]] = item_0 model_param[param_list[1]] = item_1 model_param[param_list[2]] = item_2 model_param[param_list[3]] = item_3 if self.filter_name is not None: flux_new[i, j, k, m] = self.get_flux(model_param)[0] else: flux_new[i, j, k, m, :] = self.get_model( model_param, spec_res=spec_res, wavel_resample=wavel_resample, ).flux elif n_param == 5: model_param = {} for i, item_0 in enumerate(points[0]): for j, item_1 in enumerate(points[1]): for k, item_2 in enumerate(points[2]): for m, item_3 in enumerate(points[3]): for n, item_4 in enumerate(points[4]): model_param[param_list[0]] = item_0 model_param[param_list[1]] = item_1 model_param[param_list[2]] = item_2 model_param[param_list[3]] = item_3 model_param[param_list[4]] = item_4 if self.filter_name is not None: flux_new[i, j, k, m, n] = self.get_flux( model_param )[0] else: flux_new[i, j, k, m, n, :] = self.get_model( model_param, spec_res=spec_res, wavel_resample=wavel_resample, ).flux if self.filter_name is not None: transmission = ReadFilter(self.filter_name) self.wl_points = [transmission.mean_wavelength()] else: self.wl_points = wavel_resample self.spectrum_interp = RegularGridInterpolator( points, flux_new, method="linear", bounds_error=False, fill_value=np.nan, )
[docs] @typechecked def apply_lognorm_ext( self, wavelength: np.ndarray, flux: np.ndarray, radius_interp: float, sigma_interp: float, v_band_ext: float, ) -> np.ndarray: """ Internal function for applying extinction by dust to a spectrum. wavelength : np.ndarray Wavelengths (um) of the spectrum. flux : np.ndarray Fluxes (W m-2 um-1) of the spectrum. radius_interp : float Logarithm of the mean geometric radius (um) of the log-normal size distribution. sigma_interp : float Geometric standard deviation (dimensionless) of the log-normal size distribution. v_band_ext : float The extinction (mag) in the V band. Returns ------- np.ndarray Fluxes (W m-2 um-1) with the extinction applied. """ check_dust_database() with h5py.File(self.database, "r") as hdf5_file: # Read array with cross sections dust_cross = np.array( hdf5_file["dust/lognorm/mgsio3/crystalline/cross_section"] ) # Read array with wavelengths dust_wavel = np.array( hdf5_file["dust/lognorm/mgsio3/crystalline/wavelength"] ) # Read array with geometric radius dust_radius = np.array( hdf5_file["dust/lognorm/mgsio3/crystalline/radius_g"] ) # Read array with geometric standard deviation dust_sigma = np.array(hdf5_file["dust/lognorm/mgsio3/crystalline/sigma_g"]) if wavelength[0] < dust_wavel[0]: raise ValueError( f"The shortest wavelength ({wavelength[0]:.2e} um) for which the " f"spectrum will be calculated is smaller than the shortest " f"wavelength ({dust_wavel[0]:.2e} um) of the grid with dust cross " f"sections." ) if wavelength[-1] > dust_wavel[-1]: raise ValueError( f"The longest wavelength ({wavelength[-1]:.2e} um) for which the " f"spectrum will be calculated is larger than the longest wavelength " f"({dust_wavel[-1]:.2e} um) of the grid with dust cross sections." ) # Interpolate cross sections as function of wavelength, # geometric radius, and geometric standard deviation dust_interp = RegularGridInterpolator( (dust_wavel, dust_radius, dust_sigma), dust_cross, method="linear", bounds_error=True, ) # Read Bessell V-band filter profile read_filt = ReadFilter("Generic/Bessell.V") filt_trans = read_filt.get_filter() # Create empty array for V-band extinction cross_phot = np.zeros((dust_radius.shape[0], dust_sigma.shape[0])) for i in range(dust_radius.shape[0]): for j in range(dust_sigma.shape[0]): # Calculate the filter-weighted average of the extinction cross sections cross_interp = interp1d( dust_wavel, dust_cross[:, i, j], kind="linear", bounds_error=True ) cross_tmp = cross_interp(filt_trans[:, 0]) integral1 = np.trapz(filt_trans[:, 1] * cross_tmp, x=filt_trans[:, 0]) integral2 = np.trapz(filt_trans[:, 1], x=filt_trans[:, 0]) cross_phot[i, j] = integral1 / integral2 # Interpolate the grid with the V-band extinction as # function of geometric radius and standard deviation cross_interp = RegularGridInterpolator( (dust_radius, dust_sigma), cross_phot, method="linear", bounds_error=True ) cross_v_band = cross_interp((10.0**radius_interp, sigma_interp)) # Create arrays with fixed particle properties radius_full = np.full(wavelength.shape[0], 10.0**radius_interp) sigma_full = np.full(wavelength.shape[0], sigma_interp) # Interpolate the cross-sections for the input wavelengths # and the geometric radius and standard deviation cross_new = dust_interp(np.column_stack((wavelength, radius_full, sigma_full))) # Calculate number of grains for the requested V-band extinction n_grains = v_band_ext / cross_v_band / 2.5 / np.log10(np.exp(1.0)) return flux * np.exp(-cross_new * n_grains)
[docs] @typechecked def apply_powerlaw_ext( self, wavelength: np.ndarray, flux: np.ndarray, r_max_interp: float, exp_interp: float, v_band_ext: float, ) -> np.ndarray: """ Internal function for applying extinction by dust to a spectrum. wavelength : np.ndarray Wavelengths (um) of the spectrum. flux : np.ndarray Fluxes (W m-2 um-1) of the spectrum. r_max_interp : float Maximum radius (um) of the power-law size distribution. exp_interp : float Exponent of the power-law size distribution. v_band_ext : float The extinction (mag) in the V band. Returns ------- np.ndarray Fluxes (W m-2 um-1) with the extinction applied. """ check_dust_database() with h5py.File(self.database, "r") as hdf5_file: # Read array with cross sections dust_cross = np.array( hdf5_file["dust/powerlaw/mgsio3/crystalline/cross_section"] ) # Read array with wavelengths dust_wavel = np.array( hdf5_file["dust/powerlaw/mgsio3/crystalline/wavelength"] ) # Read array with maximum particle radii dust_r_max = np.array( hdf5_file["dust/powerlaw/mgsio3/crystalline/radius_max"] ) # Read array with power-law exponents dust_exp = np.array(hdf5_file["dust/powerlaw/mgsio3/crystalline/exponent"]) # Interpolate cross sections as function of wavelength, # geometric radius, and geometric standard deviation dust_interp = RegularGridInterpolator( (dust_wavel, dust_r_max, dust_exp), dust_cross, method="linear", bounds_error=True, ) if wavelength[0] < dust_wavel[0]: raise ValueError( f"The shortest wavelength ({wavelength[0]:.2e} um) for which the " f"spectrum will be calculated is smaller than the shortest " f"wavelength ({dust_wavel[0]:.2e} um) of the grid with dust cross " f"sections." ) if wavelength[-1] > dust_wavel[-1]: raise ValueError( f"The longest wavelength ({wavelength[-1]:.2e} um) for which the " f"spectrum will be calculated is larger than the longest wavelength " f"({dust_wavel[-1]:.2e} um) of the grid with dust cross sections." ) # Read Bessell V-band filter profile read_filt = ReadFilter("Generic/Bessell.V") filt_trans = read_filt.get_filter() # Create empty array for V-band extinction cross_phot = np.zeros((dust_r_max.shape[0], dust_exp.shape[0])) for i in range(dust_r_max.shape[0]): for j in range(dust_exp.shape[0]): # Calculate the filter-weighted average of the extinction cross sections cross_interp = interp1d( dust_wavel, dust_cross[:, i, j], kind="linear", bounds_error=True ) cross_tmp = cross_interp(filt_trans[:, 0]) integral1 = np.trapz(filt_trans[:, 1] * cross_tmp, x=filt_trans[:, 0]) integral2 = np.trapz(filt_trans[:, 1], x=filt_trans[:, 0]) cross_phot[i, j] = integral1 / integral2 # Interpolate the grid with the V-band extinction as # function of geometric radius and standard deviation cross_interp = RegularGridInterpolator( (dust_r_max, dust_exp), cross_phot, method="linear", bounds_error=True ) cross_v_band = cross_interp((10.0**r_max_interp, exp_interp)) # Create arrays with fixed particle properties r_max_full = np.full(wavelength.shape[0], 10.0**r_max_interp) exp_full = np.full(wavelength.shape[0], exp_interp) # Interpolate the cross-sections for the input wavelengths # and the geometric radius and standard deviation cross_new = dust_interp(np.column_stack((wavelength, r_max_full, exp_full))) # Calculate number of grains for the requested V-band extinction n_grains = v_band_ext / cross_v_band / 2.5 / np.log10(np.exp(1.0)) return flux * np.exp(-cross_new * n_grains)
[docs] @staticmethod @typechecked def apply_ext_ism( wavelengths: np.ndarray, flux: np.ndarray, v_band_ext: float, v_band_red: float ) -> Tuple[np.ndarray, np.ndarray]: """ Internal function for applying ISM extinction to a spectrum. wavelengths : np.ndarray Wavelengths (um) of the spectrum. flux : np.ndarray Fluxes (W m-2 um-1) of the spectrum. v_band_ext : float Extinction (mag) in the V band. v_band_red : float Reddening in the V band. Returns ------- np.ndarray Fluxes (W m-2 um-1) with the extinction applied. np.ndarray Extinction (mag) as function of wavelength. """ ext_mag = ism_extinction(v_band_ext, v_band_red, wavelengths) return flux * 10.0 ** (-0.4 * ext_mag), ext_mag
[docs] @typechecked def get_model( self, model_param: Dict[str, float], spec_res: Optional[float] = None, wavel_resample: Optional[np.ndarray] = None, magnitude: bool = False, fast_rot_broad: bool = True, ext_filter: Optional[str] = None, **kwargs, ) -> ModelBox: """ Function for extracting a model spectrum by linearly interpolating the model grid. Parameters ---------- model_param : dict Dictionary with the model parameters and values. The values should be within the boundaries of the grid. The grid boundaries of the spectra in the database can be obtained with :func:`~species.read.read_model.ReadModel.get_bounds()`. spec_res : float, None Spectral resolution that is used for smoothing the spectrum with a Gaussian kernel. No smoothing is applied if the argument is set to ``None``. wavel_resample : np.ndarray, None Wavelength points (um) to which the spectrum is resampled. Optional smoothin with ``spec_res`` is applied for resampling with ``wavel_resample``. The wavelength points as stored in the database are used if the argument is set to ``None``. magnitude : bool Normalize the spectrum with a flux calibrated spectrum of Vega and return the magnitude instead of flux density. fast_rot_broad : bool Apply fast algorithm for the rotational broadening if set to ``True``, otherwise a slow but more accurate broadening is applied if set to ``False``. The fast algorithm will only provide an accurate broadening if the wavelength range of the spectrum is somewhat narrow (e.g. only the :math:`K` band). The argument is only used if the ``vsini`` parameter is included in the ``model_param`` dictionary. ext_filter : str, None Filter that is associated with the (optional) extinction parameter, ``ism_ext``. When the argument of ``ext_filter`` is set to ``None``, the extinction is defined in the visual as usual (i.e. :math:`A_V`). By providing a filter name from the `SVO Filter Profile Service <http://svo2.cab. inta-csic.es/svo/theory/fps/>`_ as argument then the extinction ``ism_ext`` is defined in that filter instead of the $V$ band. Returns ------- species.core.box.ModelBox Box with the model spectrum. """ if "smooth" in kwargs: warnings.warn( "The 'smooth' parameter has been " "deprecated. Please set only the " "'spec_res' argument, which can be set " "to None for not applying a smoothing.", DeprecationWarning, ) if not kwargs["smooth"] and spec_res is not None: spec_res = None # Get grid boundaries grid_bounds = self.get_bounds() # Check if all parameters are present and within the grid boundaries param_list = self.get_parameters() for key in param_list: if key not in model_param.keys(): raise ValueError( f"The '{key}' parameter is required by '{self.model}'. " f"The mandatory parameters are {param_list}." ) if model_param[key] < grid_bounds[key][0]: raise ValueError( f"The input value of '{key}' is smaller than the lower " f"boundary of the model grid ({model_param[key]} < " f"{grid_bounds[key][0]})." ) if model_param[key] > grid_bounds[key][1]: raise ValueError( f"The input value of '{key}' is larger than the upper " f"boundary of the model grid ({model_param[key]} > " f"{grid_bounds[key][1]})." ) # Print a warning if redundant parameters are included in the dictionary ignore_param = [] param_list = self.get_parameters() for key in model_param.keys(): if ( key not in param_list and key not in self.extra_param and not key.startswith("phot_ext_") ): warnings.warn( f"The '{key}' parameter is not required by " f"'{self.model}' so the parameter will be " f"ignored. The mandatory parameters are " f"{param_list}." ) ignore_param.append(key) # Interpolate the model grid if self.spectrum_interp is None: self.interpolate_model() # Set the wavelength range if self.wavel_range is None: wl_points = self.get_wavelengths() self.wavel_range = (wl_points[0], wl_points[-1]) # Create a list with the parameter values check_param = [ "teff", "logg", "feh", "c_o_ratio", "fsed", "log_kzz", "ad_index", ] parameters = [] for item in check_param: if item in model_param and item not in ignore_param: parameters.append(model_param[item]) # Check if the ext_filter should be adjusted # to the name that is extracted from the # phot_ext_{ext_filter} parameter if ext_filter is None: for param_item in model_param: if param_item.startswith("phot_ext_"): ext_filter = param_item[9:] # Interpolate the spectrum from the grid flux = self.spectrum_interp(parameters)[0] # Add the radius to the parameter dictionary if the mass if given if "mass" in model_param and "radius" not in model_param: mass = 1e3 * model_param["mass"] * constants.M_JUP # (g) radius = np.sqrt( 1e3 * constants.GRAVITY * mass / (10.0 ** model_param["logg"]) ) # (cm) model_param["radius"] = 1e-2 * radius / constants.R_JUP # (Rjup) # Apply (radius/distance)^2 scaling if "radius" in model_param and "parallax" in model_param: scaling = (model_param["radius"] * constants.R_JUP) ** 2 / ( 1e3 * constants.PARSEC / model_param["parallax"] ) ** 2 flux *= scaling elif "radius" in model_param and "distance" in model_param: scaling = (model_param["radius"] * constants.R_JUP) ** 2 / ( model_param["distance"] * constants.PARSEC ) ** 2 flux *= scaling elif "flux_scaling" in model_param: flux *= model_param["flux_scaling"] elif "log_flux_scaling" in model_param: flux *= 10.0 ** model_param["log_flux_scaling"] # Add optional offset to the flux if "flux_offset" in model_param: flux += model_param["flux_offset"] # Add blackbody disk component to the spectrum if "disk_teff" in model_param and "disk_radius" in model_param: disk_param = { "teff": model_param["disk_teff"], "radius": model_param["disk_radius"], } if "parallax" in model_param: disk_param["parallax"] = model_param["parallax"] elif "distance" in model_param: disk_param["distance"] = model_param["distance"] readplanck = ReadPlanck( (0.9 * self.wavel_range[0], 1.1 * self.wavel_range[-1]) ) planck_box = readplanck.get_spectrum(disk_param, spec_res=spec_res) flux_interp = interp1d( planck_box.wavelength, planck_box.flux, bounds_error=False ) flux += flux_interp(self.wl_points) # flux += spectres.spectres( # self.wl_points, planck_box.wavelength, planck_box.flux # ) # Create ModelBox with the spectrum model_box = create_box( boxtype="model", model=self.model, wavelength=self.wl_points, flux=flux, parameters=model_param, quantity="flux", spec_res=spec_res, ) # Apply rotational broadening vsin(i) in km/s if "vsini" in model_param: spec_interp = interp1d( model_box.wavelength, model_box.flux, bounds_error=False ) wavel_new = np.linspace( model_box.wavelength[0], model_box.wavelength[-1], 2 * model_box.wavelength.size, ) if fast_rot_broad: flux_broad = fastRotBroad( wvl=wavel_new, flux=spec_interp(wavel_new), epsilon=0.0, vsini=model_param["vsini"], effWvl=None, ) else: flux_broad = rotBroad( wvl=wavel_new, flux=spec_interp(wavel_new), epsilon=0.0, vsini=model_param["vsini"], edgeHandling="firstlast", ) spec_interp = interp1d(wavel_new, flux_broad, bounds_error=False) model_box.flux = spec_interp(model_box.wavelength) # Apply veiling if ( "veil_a" in model_param and "veil_b" in model_param and "veil_ref" in model_param ): lambda_ref = 0.5 # (um) veil_flux = model_param["veil_ref"] + model_param["veil_b"] * ( model_box.wavelength - lambda_ref ) model_box.flux = model_param["veil_a"] * model_box.flux + veil_flux # Apply extinction if ( "lognorm_radius" in model_param and "lognorm_sigma" in model_param and "lognorm_ext" in model_param ): model_box.flux = self.apply_lognorm_ext( model_box.wavelength, model_box.flux, model_param["lognorm_radius"], model_param["lognorm_sigma"], model_param["lognorm_ext"], ) if ( "powerlaw_max" in model_param and "powerlaw_exp" in model_param and "powerlaw_ext" in model_param ): model_box.flux = self.apply_powerlaw_ext( model_box.wavelength, model_box.flux, model_param["powerlaw_max"], model_param["powerlaw_exp"], model_param["powerlaw_ext"], ) if "ism_ext" in model_param or ext_filter is not None: ism_reddening = model_param.get("ism_red", 3.1) if ext_filter is not None: ism_ext_av = convert_to_av( filter_name=ext_filter, filter_ext=model_param[f"phot_ext_{ext_filter}"], v_band_red=ism_reddening, ) else: ism_ext_av = model_param["ism_ext"] model_box.flux, ext_mag = self.apply_ext_ism( model_box.wavelength, model_box.flux, ism_ext_av, ism_reddening, ) idx_select = ext_mag >= 0.0 model_box.wavelength = model_box.wavelength[idx_select] model_box.flux = model_box.flux[idx_select] # Smooth the spectrum if spec_res is not None: model_box.flux = smooth_spectrum( model_box.wavelength, model_box.flux, spec_res ) # Resample the spectrum if wavel_resample is not None: flux_interp = interp1d( model_box.wavelength, model_box.flux, bounds_error=False ) model_box.flux = flux_interp(wavel_resample) # model_box.flux = spectres.spectres( # wavel_resample, # model_box.wavelength, # model_box.flux, # spec_errs=None, # fill=np.nan, # verbose=True, # ) model_box.wavelength = wavel_resample # elif spec_res is not None and not smooth: # index = np.where(np.isnan(model_box.flux))[0] # # if index.size > 0: # raise ValueError( # "Flux values should not contains NaNs. Please make sure that " # "the parameter values and the wavelength range are within " # "the grid boundaries as stored in the database." # ) # # wavel_resample = create_wavelengths( # (self.wl_points[0], self.wl_points[-1]), spec_res # ) # # indices = np.where( # (wavel_resample > self.wl_points[0]) # & (wavel_resample < self.wl_points[-2]) # )[0] # # wavel_resample = wavel_resample[indices] # # flux_interp = interp1d( # model_box.wavelength, model_box.flux, bounds_error=False # ) # model_box.flux = flux_interp(wavel_resample) # # # model_box.flux = spectres.spectres( # # wavel_resample, # # model_box.wavelength, # # model_box.flux, # # spec_errs=None, # # fill=np.nan, # # verbose=True, # # ) # # model_box.wavelength = wavel_resample # Convert flux to magnitude if magnitude: with h5py.File(self.database, "a") as hdf5_file: if "spectra/calibration/vega" not in hdf5_file: add_vega(self.data_folder, hdf5_file) vega_spec = np.array(hdf5_file["spectra/calibration/vega"]) flux_interp = interp1d(vega_spec[0,], vega_spec[1,]) flux_vega = flux_interp(model_box.wavelength) # flux_vega, _ = spectres.spectres( # model_box.wavelength, # calib_box.wavelength, # calib_box.flux, # spec_errs=calib_box.error, # fill=np.nan, # verbose=True, # ) model_box.flux = -2.5 * np.log10(model_box.flux / flux_vega) model_box.quantity = "magnitude" # Check if the contains NaNs if np.isnan(np.sum(model_box.flux)): warnings.warn( "The resampled spectrum contains " f"{np.sum(np.isnan(model_box.flux))} " "NaNs, probably because the original " "wavelength range does not fully encompass " "the new wavelength range. This happened with " f"the following parameters: {model_param}." ) # Add the luminosity to the parameter dictionary if "radius" in model_box.parameters: model_box.parameters["luminosity"] = ( 4.0 * np.pi * (model_box.parameters["radius"] * constants.R_JUP) ** 2 * constants.SIGMA_SB * model_box.parameters["teff"] ** 4.0 / constants.L_SUN ) # (Lsun) # Add the blackbody disk component to the luminosity if ( "disk_teff" in model_box.parameters and "disk_radius" in model_box.parameters ): model_box.parameters["luminosity"] += ( 4.0 * np.pi * (model_box.parameters["disk_radius"] * constants.R_JUP) ** 2 * constants.SIGMA_SB * model_box.parameters["disk_teff"] ** 4.0 / constants.L_SUN ) # (Lsun) # Add the planet mass to the parameter dictionary if "radius" in model_param and "logg" in model_param: model_param["mass"] = logg_to_mass( model_param["logg"], model_param["radius"] ) return model_box
[docs] @typechecked def get_data( self, model_param: Dict[str, float], spec_res: Optional[float] = None, wavel_resample: Optional[np.ndarray] = None, ext_filter: Optional[str] = None, ) -> ModelBox: """ Function for selecting a model spectrum (without interpolation) for a set of parameter values that coincide with the grid points. The stored grid points can be inspected with :func:`~species.read.read_model.ReadModel.get_points`. Parameters ---------- model_param : dict Model parameters and values. Only discrete values from the original grid are possible. Else, the nearest grid values are selected. spec_res : float, None Spectral resolution that is used for smoothing the spectrum with a Gaussian kernel. No smoothing is applied to the spectrum if the argument is set to ``None``. wavel_resample : np.ndarray, None Wavelength points (um) to which the spectrum will be resampled. In that case, ``spec_res`` can still be used for smoothing the spectrum with a Gaussian kernel. The original wavelength points are used if the argument is set to ``None``. ext_filter : str, None Filter that is associated with the (optional) extinction parameter, ``ism_ext``. When the argument of ``ext_filter`` is set to ``None``, the extinction is defined in the visual as usual (i.e. :math:`A_V`). By providing a filter name from the `SVO Filter Profile Service <http://svo2.cab. inta-csic.es/svo/theory/fps/>`_ as argument then the extinction ``ism_ext`` is defined in that filter instead of the $V$ band. Returns ------- species.core.box.ModelBox Box with the model spectrum. """ # Check if all parameters are present param_list = self.get_parameters() for key in param_list: if key not in model_param.keys(): raise ValueError( f"The '{key}' parameter is required by '{self.model}'. " f"The mandatory parameters are {param_list}." ) # Print a warning if redundant parameters are included in the dictionary ignore_param = [] for key in model_param.keys(): if key not in param_list and key not in self.extra_param: warnings.warn( f"The '{key}' parameter is not required by " f"'{self.model}' so the parameter will be " f"ignored. The mandatory parameters are " f"{param_list}." ) ignore_param.append(key) # Set the wavelength range if self.wavel_range is None: wl_points = self.get_wavelengths() self.wavel_range = (wl_points[0], wl_points[-1]) # Model parameters to check check_param = [ "teff", "logg", "feh", "c_o_ratio", "fsed", "log_kzz", "ad_index", ] # Create lists with the parameter names and values param_key = [] param_val = [] for param_item in check_param: if param_item in model_param and param_item not in ignore_param: param_key.append(param_item) param_val.append(model_param[param_item]) # Check if the ext_filter should be adjusted # to the name that is extracted from the # phot_ext_{ext_filter} parameter if ext_filter is None: for param_item in model_param: if param_item.startswith("phot_ext_"): ext_filter = param_item[9:] # Read the wavelength grid and the indices that will # be used for the wavelength range wl_points, wl_index = self.wavelength_points() # Open de HDF5 database with self.open_database() as hdf5_file: # Read the grid of fluxes from the database flux = np.array(hdf5_file[f"models/{self.model}/flux"]) # Find the indices of the grid points for which the spectrum will be extracted indices = [] for i, item in enumerate(param_key): data = np.array(hdf5_file[f"models/{self.model}/{item}"]) data_index = np.argwhere( np.round(data, 4) == np.round(model_param[item], 4) ) if len(data_index) == 0: raise ValueError( f"The parameter {item}={param_val[i]} is not found." ) data_index = data_index[0] indices.append(data_index[0]) # Append the wavelength indices to the list of grid indices indices.append(wl_index) # Extract the spectrum at the requested grid point and # wavelength range flux = flux[tuple(indices)] # Apply (radius/distance)^2 scaling if "radius" in model_param and "parallax" in model_param: scaling = (model_param["radius"] * constants.R_JUP) ** 2 / ( 1e3 * constants.PARSEC / model_param["parallax"] ) ** 2 flux *= scaling elif "radius" in model_param and "distance" in model_param: scaling = (model_param["radius"] * constants.R_JUP) ** 2 / ( model_param["distance"] * constants.PARSEC ) ** 2 flux *= scaling elif "flux_scaling" in model_param: flux *= model_param["flux_scaling"] elif "log_flux_scaling" in model_param: flux *= 10.0 ** model_param["log_flux_scaling"] # Add optional offset to the flux if "flux_offset" in model_param: flux += model_param["flux_offset"] # Add blackbody disk component to the spectrum if "disk_teff" in model_param and "disk_radius" in model_param: disk_param = { "teff": model_param["disk_teff"], "radius": model_param["disk_radius"], } if "parallax" in model_param: disk_param["parallax"] = model_param["parallax"] elif "distance" in model_param: disk_param["distance"] = model_param["distance"] readplanck = ReadPlanck( (0.9 * self.wavel_range[0], 1.1 * self.wavel_range[-1]) ) planck_box = readplanck.get_spectrum(disk_param, spec_res=spec_res) flux_interp = interp1d( planck_box.wavelength, planck_box.flux, bounds_error=False ) flux += flux_interp(wl_points) # flux += spectres.spectres(wl_points, planck_box.wavelength, planck_box.flux) # Create ModelBox with the spectrum model_box = create_box( boxtype="model", model=self.model, wavelength=wl_points, flux=flux, parameters=model_param, quantity="flux", spec_res=spec_res, ) # Apply extinction if ( "lognorm_radius" in model_param and "lognorm_sigma" in model_param and "lognorm_ext" in model_param ): model_box.flux = self.apply_lognorm_ext( model_box.wavelength, model_box.flux, model_param["lognorm_radius"], model_param["lognorm_sigma"], model_param["lognorm_ext"], ) if ( "powerlaw_max" in model_param and "powerlaw_exp" in model_param and "powerlaw_ext" in model_param ): model_box.flux = self.apply_powerlaw_ext( model_box.wavelength, model_box.flux, model_param["powerlaw_max"], model_param["powerlaw_exp"], model_param["powerlaw_ext"], ) if "ism_ext" in model_param or ext_filter is not None: ism_reddening = model_param.get("ism_red", 3.1) if ext_filter is not None: ism_ext_av = convert_to_av( filter_name=ext_filter, filter_ext=model_param[f"phot_ext_{ext_filter}"], v_band_red=ism_reddening, ) else: ism_ext_av = model_param["ism_ext"] model_box.flux, ext_mag = self.apply_ext_ism( model_box.wavelength, model_box.flux, ism_ext_av, ism_reddening, ) idx_select = ext_mag >= 0.0 model_box.wavelength = model_box.wavelength[idx_select] model_box.flux = model_box.flux[idx_select] # Smooth the spectrum if spec_res is not None: model_box.flux = smooth_spectrum( model_box.wavelength, model_box.flux, spec_res ) # Resample the spectrum if wavel_resample is not None: flux_interp = interp1d( model_box.wavelength, model_box.flux, bounds_error=False ) model_box.flux = flux_interp(wavel_resample) # model_box.flux = spectres.spectres( # wavel_resample, # model_box.wavelength, # model_box.flux, # spec_errs=None, # fill=np.nan, # verbose=True, # ) model_box.wavelength = wavel_resample # Add the planet mass to the parameter dictionary if "radius" in model_param and "logg" in model_param: model_param["mass"] = logg_to_mass( model_param["logg"], model_param["radius"] ) # Add the luminosity to the parameter dictionary if "radius" in model_box.parameters: model_box.parameters["luminosity"] = ( 4.0 * np.pi * (model_box.parameters["radius"] * constants.R_JUP) ** 2 * constants.SIGMA_SB * model_box.parameters["teff"] ** 4.0 / constants.L_SUN ) # (Lsun) # Add the blackbody disk component to the luminosity if ( "disk_teff" in model_box.parameters and "disk_radius" in model_box.parameters ): model_box.parameters["luminosity"] += ( 4.0 * np.pi * (model_box.parameters["disk_radius"] * constants.R_JUP) ** 2 * constants.SIGMA_SB * model_box.parameters["disk_teff"] ** 4.0 / constants.L_SUN ) # (Lsun) return model_box
[docs] @typechecked def get_flux( self, model_param: Dict[str, float], synphot=None, return_box: bool = False ) -> Union[Tuple[Optional[float], Optional[float]], PhotometryBox]: """ Function for calculating the average flux density for the ``filter_name``. Parameters ---------- model_param : dict Model parameters and values. synphot : SyntheticPhotometry, None Synthetic photometry object. The object is created if set to ``None``. return_box : bool Return a :class:`~species.core.box.PhotometryBox` if set to ``True`` or return the two values that are specified below if set to ``False``. By default, the argument is set to ``False``. The advantage of returning the output in a :class:`~species.core.box.PhotometryBox` is that it can directly be provided as input to :func:`~species.plot.plot_spectrum.plot_spectrum`. Returns ------- float Average flux (W m-2 um-1). float, None Uncertainty (W m-2 um-1), which is set to ``None``. """ param_list = self.get_parameters() for key in param_list: if key not in model_param.keys(): raise ValueError( f"The '{key}' parameter is required by '{self.model}'. " f"The mandatory parameters are {param_list}." ) if self.spectrum_interp is None: self.interpolate_model() spectrum = self.get_model(model_param) if synphot is None: synphot = SyntheticPhotometry(self.filter_name) model_flux = synphot.spectrum_to_flux(spectrum.wavelength, spectrum.flux) if return_box: model_mag = self.get_magnitude(model_param) phot_box = create_box( boxtype="photometry", name=self.model, wavelength=[self.mean_wavelength], flux=[model_flux], app_mag=[(model_mag[0], None)], abs_mag=[(model_mag[1], None)], filter_name=[self.filter_name], ) return phot_box return model_flux
[docs] @typechecked def get_magnitude( self, model_param: Dict[str, float], return_box: bool = False, ) -> Union[Tuple[Optional[float], Optional[float]], PhotometryBox]: """ Function for calculating the apparent and absolute magnitudes for the ``filter_name``. Parameters ---------- model_param : dict Dictionary with the model parameters. A ``radius`` (Rjup), and ``parallax`` (mas) or ``distance`` (pc) are required for the apparent magnitude (i.e. to scale the flux from the planet to the observer). Only a ``radius`` is required for the absolute magnitude. return_box : bool Return a :class:`~species.core.box.PhotometryBox` if set to ``True`` or return the two values that are specified below if set to ``False``. By default, the argument is set to ``False``. The advantage of returning the output in a :class:`~species.core.box.PhotometryBox` is that it can directly be provided as input to :func:`~species.plot.plot_spectrum.plot_spectrum`. Returns ------- float Apparent magnitude. A ``None`` is returned if the dictionary of ``model_param`` does not contain a ``radius``, and ``parallax`` or ``distance``. float, None Absolute magnitude. A ``None`` is returned if the dictionary of ``model_param`` does not contain a ``radius``. """ param_list = self.get_parameters() for key in param_list: if key not in model_param.keys(): raise ValueError( f"The '{key}' parameter is required by '{self.model}'. " f"The mandatory parameters are {param_list}." ) if self.spectrum_interp is None: self.interpolate_model() try: spectrum = self.get_model(model_param) except ValueError: warnings.warn( f"The set of model parameters {model_param} is outside " f"the grid range {self.get_bounds()} so returning a NaN." ) return np.nan, np.nan if spectrum.wavelength.size == 0: app_mag = (np.nan, None) abs_mag = (np.nan, None) else: synphot = SyntheticPhotometry(self.filter_name) if "radius" in model_param and "parallax" in model_param: app_mag, abs_mag = synphot.spectrum_to_magnitude( spectrum.wavelength, spectrum.flux, parallax=(model_param["parallax"], None), ) elif "radius" in model_param and "distance" in model_param: app_mag, abs_mag = synphot.spectrum_to_magnitude( spectrum.wavelength, spectrum.flux, distance=(model_param["distance"], None), ) else: app_mag = (None, None) abs_mag = (None, None) if "radius" in model_param: distance = 10.0 # (pc) spectrum.flux *= (model_param["radius"] * constants.R_JUP) ** 2 spectrum.flux /= (distance * constants.PARSEC) ** 2 _, abs_mag = synphot.spectrum_to_magnitude( spectrum.wavelength, spectrum.flux, distance=(distance, None) ) if return_box: model_flux = self.get_flux(model_param) phot_box = create_box( boxtype="photometry", name=self.model, wavelength=[self.mean_wavelength], flux=[model_flux], app_mag=[(app_mag[0], None)], abs_mag=[(abs_mag[0], None)], filter_name=[self.filter_name], ) return phot_box return app_mag[0], abs_mag[0]
[docs] @typechecked def get_bounds(self) -> Dict[str, Tuple[float, float]]: """ Function for extracting the grid boundaries. Returns ------- dict Boundaries of parameter grid. """ param_list = self.get_parameters() with self.open_database() as hdf5_file: bounds = {} for param_item in param_list: data = hdf5_file[f"models/{self.model}/{param_item}"] bounds[param_item] = (data[0], data[-1]) return bounds
[docs] @typechecked def get_wavelengths(self) -> np.ndarray: """ Function for extracting the wavelength points. Returns ------- np.ndarray Wavelength points (um). """ with self.open_database() as hdf5_file: wavelength = np.array(hdf5_file[f"models/{self.model}/wavelength"]) return wavelength
[docs] @typechecked def get_points(self) -> Dict[str, np.ndarray]: """ Function for extracting the grid points. Returns ------- dict Parameter points of the model grid. """ param_list = self.get_parameters() points = {} with self.open_database() as hdf5_file: points = {} for param_item in param_list: data = hdf5_file[f"models/{self.model}/{param_item}"] points[param_item] = np.array(data) return points
[docs] @typechecked def get_parameters(self) -> List[str]: """ Function for extracting the parameter names. Returns ------- list(str) Model parameters. """ with self.open_database() as hdf5_file: dset = hdf5_file[f"models/{self.model}"] if "n_param" in dset.attrs: n_param = dset.attrs["n_param"] elif "nparam" in dset.attrs: n_param = dset.attrs["nparam"] param = [] for i in range(n_param): param.append(dset.attrs[f"parameter{i}"]) return param
[docs] @typechecked def get_sampling(self) -> float: """ Function for returning the wavelength sampling, :math:`\\lambda/\\Delta\\lambda`, of the model spectra as stored in the database. Returns ------- float Wavelength sampling, :math:`\\lambda/\\Delta\\lambda`. """ wavel_points = self.get_wavelengths() wavel_mean = (wavel_points[1:] + wavel_points[:-1]) / 2.0 wavel_sampling = wavel_mean / np.diff(wavel_points) return np.mean(wavel_sampling)
[docs] @typechecked def binary_spectrum( self, model_param: Dict[str, float], spec_res: Optional[float] = None, wavel_resample: Optional[np.ndarray] = None, **kwargs, ) -> ModelBox: """ Function for extracting a model spectrum of a binary system. A weighted combination of two spectra will be returned. The ``model_param`` dictionary should contain the parameters for both components (e.g. ``teff_0`` and ``teff_1``, instead of ``teff``). Apart from that, the same parameters are used as with :meth:`~species.read.read_model.ReadModel.get_model`. Parameters ---------- model_param : dict Dictionary with the model parameters and values. The values should be within the boundaries of the grid. The grid boundaries of the spectra in the database can be obtained with :func:`~species.read.read_model.ReadModel.get_bounds()`. spec_res : float, None Spectral resolution that is used for smoothing the spectrum with a Gaussian kernel. No smoothing is applied if the argument is set to ``None``. wavel_resample : np.ndarray, None Wavelength points (um) to which the spectrum is resampled. In that case, ``spec_res`` can still be used for smoothing the spectrum with a Gaussian kernel. The original wavelength points are used if the argument is set to ``None``. Returns ------- species.core.box.ModelBox Box with the model spectrum. """ if "smooth" in kwargs: warnings.warn( "The 'smooth' parameter has been " "deprecated. Please set only the " "'spec_res' argument, which can be set " "to None for not applying a smoothing.", DeprecationWarning, ) if not kwargs["smooth"] and spec_res is not None: spec_res = None # Get grid boundaries param_0 = binary_to_single(model_param, 0) model_box_0 = self.get_model( param_0, spec_res=spec_res, wavel_resample=wavel_resample, ) param_1 = binary_to_single(model_param, 1) model_box_1 = self.get_model( param_1, spec_res=spec_res, wavel_resample=wavel_resample, ) flux_comb = ( model_param["spec_weight"] * model_box_0.flux + (1.0 - model_param["spec_weight"]) * model_box_1.flux ) model_box = create_box( boxtype="model", model=self.model, wavelength=model_box_0.wavelength, flux=flux_comb, parameters=model_param, quantity="flux", spec_res=spec_res, ) return model_box
[docs] @typechecked def integrate_spectrum(self, model_param: Dict[str, float]) -> float: """ Function for calculating the bolometric flux by integrating a model spectrum at the requested parameters. In principle, the calculated luminosity should be approximately the same as the value that can be calculated directly from the :math:`T_\\mathrm{eff}` and radius parameters, unless the atmospheric model had not properly converged. Parameters ---------- model_param : dict Dictionary with the model parameters and values. The values should be within the boundaries of the grid. The grid boundaries of the spectra in the database can be obtained with :func:`~species.read.read_model.ReadModel.get_bounds()`. Returns ------- float Bolometric luminosity (:math:`\\log{(L/L_\\odot)}`). """ if "radius" not in model_param: raise ValueError( "Please include the 'radius' parameter " "in the 'model_param' dictionary, " "which is require for calculating the " "bolometric luminosity." ) wavel_points = self.get_wavelengths() if self.wavel_range is None: self.wavel_range = (wavel_points[0], wavel_points[-1]) if ( self.wavel_range[0] != wavel_points[0] or self.wavel_range[1] != wavel_points[-1] ): warnings.warn( "The 'wavel_range' is not set to the " "maximum available range. To maximize the " "accuracy when calculating the bolometric " "luminosity, it is recommended to set " "'wavel_range=None'." ) if "parallax" in model_param: del model_param["parallax"] if "distance" in model_param: del model_param["distance"] model_box = self.get_model(model_param) bol_lum = ( 4.0 * np.pi * (model_param["radius"] * constants.R_JUP) ** 2 * simps(model_box.flux, model_box.wavelength) ) return np.log10(bol_lum / constants.L_SUN)
[docs] @typechecked def create_color_magnitude( self, model_param: Dict[str, float], filters_color: Tuple[str, str], filter_mag: str, ) -> ColorMagBox: """ Function for creating a :class:`~species.core.box. ColorMagBox` for a given set of filter names and model parameters. The effective temperature, :math:`T_\\mathrm{eff}`, is varied such that the returned :class:`~species.core.box. ColorMagBox` contains the colors as function of :math:`T_\\mathrm{eff}` and can be provide as input to the :func:`~species.plot.plot_color.plot_color_magnitude` function. Parameters ---------- model_param : dict Dictionary with the model parameters and values. The values should be within the boundaries of the grid. The boundaries of the model grid can be inspected by using the :func:`~species.read.read_model.ReadModel.get_bounds()` method. The effective temperature, :math:`T_\\mathrm{eff}`, does not need to be included in the dictionary since it is varied. The values of :math:`T_\\mathrm{eff}` are set to the grid points. The grid points can be inspected with the :func:`~species.read.read_model.ReadModel.get_points()` method. filters_color : tuple(str, str) Filter names that are used for the color. Any of the filter names from the `SVO Filter Profile Service <http://svo2.cab.inta-csic.es/svo/theory/fps/>`_ are compatible. filter_mag : str Filter name that is used for the magnitude. Any of the filter names from the `SVO Filter Profile Service <http://svo2.cab.inta-csic.es/svo/theory/fps/>`_ are compatible. Returns ------- species.core.box.ColorMagBox Box with the colors and magnitudes. """ if "distance" not in model_param: model_param["distance"] = 10.0 if "radius" not in model_param: model_param["radius"] = 1.0 if "parallax" in model_param: del model_param["parallax"] read_model_1 = ReadModel(self.model, filter_name=filters_color[0]) read_model_2 = ReadModel(self.model, filter_name=filters_color[0]) read_model_3 = ReadModel(self.model, filter_name=filter_mag) model_points = self.get_points() param_list = [] color_list = [] mag_list = [] for param_item in model_points["teff"]: model_param["teff"] = param_item mag_1 = read_model_1.get_magnitude(model_param) mag_2 = read_model_2.get_magnitude(model_param) mag_3 = read_model_3.get_magnitude(model_param) param_list.append(param_item) color_list.append(mag_1[0] - mag_2[0]) mag_list.append(mag_3[0]) return create_box( "colormag", library=self.model, object_type="spectra", filters_color=filters_color, filter_mag=filter_mag, color=color_list, magnitude=mag_list, sptype=param_list, )
[docs] @typechecked def create_color_color( self, model_param: Dict[str, float], filters_colors: Tuple[Tuple[str, str], Tuple[str, str]], ) -> ColorColorBox: """ Function for creating a :class:`~species.core.box. ColorColorBox` for a given set of filter names and model parameters. The effective temperature, :math:`T_\\mathrm{eff}`, is varied such that the returned :class:`~species.core.box. ColorColorBox` contains the colors as function of :math:`T_\\mathrm{eff}` and can be provide as input to the :func:`~species.plot.plot_color.plot_color_color` function. Parameters ---------- model_param : dict Dictionary with the model parameters and values. The values should be within the boundaries of the grid. The boundaries of the model grid can be inspected by using the :func:`~species.read.read_model.ReadModel.get_bounds()` method. The effective temperature, :math:`T_\\mathrm{eff}`, does not need to be included in the dictionary since it is varied. The values of :math:`T_\\mathrm{eff}` are set to the grid points. The grid points can be inspected with the :func:`~species.read.read_model.ReadModel.get_points()` method. filters_colors : tuple(tuple(str, str), tuple(str, str)) Filter names that are used for the two colors. Any of the filter names from the `SVO Filter Profile Service <http://svo2.cab.inta-csic.es/svo/theory/fps/>`_ are compatible. Returns ------- species.core.box.ColorColorBox Box with the colors. """ if "distance" not in model_param: model_param["distance"] = 10.0 if "radius" not in model_param: model_param["radius"] = 1.0 if "parallax" in model_param: del model_param["parallax"] read_model_1 = ReadModel(self.model, filter_name=filters_colors[0][0]) read_model_2 = ReadModel(self.model, filter_name=filters_colors[0][1]) read_model_3 = ReadModel(self.model, filter_name=filters_colors[1][0]) read_model_4 = ReadModel(self.model, filter_name=filters_colors[1][1]) model_points = self.get_points() param_list = [] color_1_list = [] color_2_list = [] for param_item in model_points["teff"]: model_param["teff"] = param_item mag_1 = read_model_1.get_magnitude(model_param) mag_2 = read_model_2.get_magnitude(model_param) mag_3 = read_model_3.get_magnitude(model_param) mag_4 = read_model_4.get_magnitude(model_param) param_list.append(param_item) color_1_list.append(mag_1[0] - mag_2[0]) color_2_list.append(mag_3[0] - mag_4[0]) return create_box( "colorcolor", library=self.model, object_type="spectra", filters=filters_colors, color1=color_1_list, color2=color_2_list, sptype=param_list, )