From d0a340a7b6812c61af97ff7431479322ee69a87e Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Wed, 18 Oct 2023 12:38:58 +0200 Subject: [PATCH] Added the contrast_to_mass method to ReadIsochrone for converting contrast values in a given filter into masses --- docs/overview.rst | 2 +- docs/tutorials/read_isochrone.ipynb | 145 ++++++++++++++++-- species/data/model_data.json | 4 +- species/plot/plot_comparison.py | 8 +- species/read/read_isochrone.py | 229 ++++++++++++++++++++++++++-- species/read/read_model.py | 2 +- 6 files changed, 359 insertions(+), 31 deletions(-) diff --git a/docs/overview.rst b/docs/overview.rst index 104e0e0c..943fd859 100644 --- a/docs/overview.rst +++ b/docs/overview.rst @@ -28,7 +28,7 @@ The following data and models are currently supported: - `DRIFT-PHOENIX `_ - `Exo-REM `_ - `Morley et al. (2012) T/Y dwarf spectra `_ -- `petitCODE `_ +- `petitCODE `_ - `petitRADTRANS `_ - `Saumon & Marley (2008) `_ - `Sonora Bobcat `_ diff --git a/docs/tutorials/read_isochrone.ipynb b/docs/tutorials/read_isochrone.ipynb index f8eda91b..a6d789e3 100644 --- a/docs/tutorials/read_isochrone.ipynb +++ b/docs/tutorials/read_isochrone.ipynb @@ -63,7 +63,7 @@ "output_type": "stream", "text": [ "==============\n", - "species v0.7.2\n", + "species v0.7.3\n", "==============\n", "Working folder: /Users/tomasstolker/applications/species/docs/tutorials\n", "Creating species_config.ini... [DONE]\n", @@ -79,7 +79,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 2, @@ -455,7 +455,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 10, "id": "2e93e5de", "metadata": {}, "outputs": [], @@ -476,7 +476,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 11, "id": "606c4fde", "metadata": {}, "outputs": [ @@ -590,7 +590,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 12, "id": "c944ed49", "metadata": {}, "outputs": [ @@ -605,7 +605,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/Users/tomasstolker/applications/species/species/read/read_model.py:145: UserWarning: The 'atmo-ceq' 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 yourself.\n", + "/Users/tomasstolker/applications/species/species/read/read_model.py:145: UserWarning: The 'atmo-ceq' 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.\n", " warnings.warn(\n", "Downloading data from 'https://home.strw.leidenuniv.nl/~stolker/species/atmo-ceq.tgz' to file '/Users/tomasstolker/applications/species/docs/tutorials/data/atmo-ceq.tgz'.\n" ] @@ -621,7 +621,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "100%|████████████████████████████████████████| 476M/476M [00:00<00:00, 147GB/s]\n", + "100%|████████████████████████████████████████| 476M/476M [00:00<00:00, 229GB/s]\n", "SHA256 hash of downloaded file: bab2dc6920d9358bd0a554b9d5554b09e5885fec2daf91ee11ccf340c310fcbf\n", "Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.\n" ] @@ -656,7 +656,7 @@ "output_type": "stream", "text": [ "Downloading data from 'https://archive.stsci.edu/hlsps/reference-atlases/cdbs/current_calspec/alpha_lyr_stis_011.fits' to file '/Users/tomasstolker/applications/species/docs/tutorials/data/alpha_lyr_stis_011.fits'.\n", - "100%|████████████████████████████████████████| 288k/288k [00:00<00:00, 161MB/s]" + "100%|████████████████████████████████████████| 288k/288k [00:00<00:00, 105MB/s]" ] }, { @@ -690,7 +690,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 13, "id": "0c32a006", "metadata": {}, "outputs": [ @@ -731,7 +731,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 14, "id": "83a86b99", "metadata": {}, "outputs": [], @@ -752,7 +752,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 15, "id": "be9a4029", "metadata": {}, "outputs": [ @@ -846,7 +846,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 16, "id": "0ca8aa3b", "metadata": {}, "outputs": [ @@ -876,7 +876,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 17, "id": "b7959461", "metadata": {}, "outputs": [ @@ -1080,6 +1080,125 @@ "source": [ "fig" ] + }, + { + "cell_type": "markdown", + "id": "c43b9f35-3b87-4971-b62e-0d4689d2785e", + "metadata": {}, + "source": [ + "## Converting contrast into mass" + ] + }, + { + "cell_type": "markdown", + "id": "5c999f15-785d-43cc-b31c-d24f937d9b2a", + "metadata": {}, + "source": [ + "The [contrast_to_mass](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_isochrone.ReadIsochrone.contrast_to_mass) method of [ReadIsochrone](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_isochrone.ReadIsochrone) can be used to convert a list or array with companion-to-star contrast values into masses of the companion. This feature can be useful for converting a *contrast curve* (i.e. detection limits for point sources) into mass limits as function of separation from the star. Or, to estimate the mass of a companion from a measured brightness contrast with its star." + ] + }, + { + "cell_type": "markdown", + "id": "a0bc0df6-8795-4c9c-ac92-05f3b9379cb2", + "metadata": {}, + "source": [ + "We can either select one of the filters with pre-calculated magnitudes that come with the isochrone grid (i.e., those returned by [get_filters](https://species.readthedocs.io/en/latest/species.read.html#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/)." + ] + }, + { + "cell_type": "markdown", + "id": "ccf4f4f1-37d8-4e45-b799-1143f94fcb46", + "metadata": {}, + "source": [ + "For example, let's select the ``MKO_H`` filter that is include with the ATMO grid and interpolate the masses for a list of contrast values that are provided as magnitudes." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "cf2e2a56-b81b-44bd-803a-d0c0badefe1c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The 'MKO_H' filter is found in the list of available filters from the isochrone data of 'atmo-ceq'.\n", + "The requested contrast values will be directly interpolated from the grid with pre-calculated magnitudes.\n", + "[12.18523169 6.5719161 3.0361701 ]\n" + ] + } + ], + "source": [ + "masses = read_iso.contrast_to_mass(age=20.,\n", + " distance=50.,\n", + " filter_name='MKO_H',\n", + " star_mag=10.,\n", + " contrast=[5.0, 7.5, 10.],\n", + " use_mag=True)\n", + "print(masses)" + ] + }, + { + "cell_type": "markdown", + "id": "3c0b501f-c4c1-4e55-a35f-d8ecd7451400", + "metadata": {}, + "source": [ + "Alternatively, we can select a filter from the SVO Filter Profile Service. In this case, the function will use the associated model spectra for calculating and interpolating the synthetic photometry. In this example we will provide the contrast values as ratios instead of magnitudes, so we set ``use_mag=False``." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "83877efb-0a3d-4cdc-a29f-6b18c29f0872", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The 'JWST/NIRCam.F162M' filter is not found in the list of available filters from the isochrone data of 'atmo-ceq'.\n", + "It will be tried to download the filter profile (if needed) and to use the associated atmospheric model spectra for calculating synthetic photometry.\n", + "Adding filter: JWST/NIRCam.F162M... [DONE]\n", + "[11.68030143 6.2010683 2.81440063]\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/tomasstolker/applications/species/species/read/read_model.py:1606: UserWarning: The set of model parameters {'teff': 3001.357148379523, 'logg': 4.455553724453863, 'radius': 2.526381394429854, 'distance': 50.0} is outside the grid range {'teff': (200.0, 3000.0), 'logg': (2.5, 5.5)} so returning a NaN.\n", + " warnings.warn(\n", + "/Users/tomasstolker/applications/species/species/read/read_model.py:1606: UserWarning: The set of model parameters {'teff': 3006.482388036101, 'logg': 4.4562069114763485, 'radius': 2.5424416720662784, 'distance': 50.0} is outside the grid range {'teff': (200.0, 3000.0), 'logg': (2.5, 5.5)} so returning a NaN.\n", + " warnings.warn(\n", + "/Users/tomasstolker/applications/species/species/read/read_model.py:1606: UserWarning: The set of model parameters {'teff': 3011.408951792053, 'logg': 4.456811598156431, 'radius': 2.558493215446333, 'distance': 50.0} is outside the grid range {'teff': (200.0, 3000.0), 'logg': (2.5, 5.5)} so returning a NaN.\n", + " warnings.warn(\n", + "/Users/tomasstolker/applications/species/species/read/read_model.py:1606: UserWarning: The set of model parameters {'teff': 3016.267586113234, 'logg': 4.457442944709719, 'radius': 2.5743185127516264, 'distance': 50.0} is outside the grid range {'teff': (200.0, 3000.0), 'logg': (2.5, 5.5)} so returning a NaN.\n", + " warnings.warn(\n", + "/Users/tomasstolker/applications/species/species/read/read_model.py:1606: UserWarning: The set of model parameters {'teff': 3021.073033921284, 'logg': 4.458036630627344, 'radius': 2.5901106044291575, 'distance': 50.0} is outside the grid range {'teff': (200.0, 3000.0), 'logg': (2.5, 5.5)} so returning a NaN.\n", + " warnings.warn(\n", + "/Users/tomasstolker/applications/species/species/read/read_model.py:1606: UserWarning: The set of model parameters {'teff': 3025.8044035429903, 'logg': 4.45866077561258, 'radius': 2.6056736019935696, 'distance': 50.0} is outside the grid range {'teff': (200.0, 3000.0), 'logg': (2.5, 5.5)} so returning a NaN.\n", + " warnings.warn(\n" + ] + } + ], + "source": [ + "masses = read_iso.contrast_to_mass(age=20.,\n", + " distance=50.,\n", + " filter_name='JWST/NIRCam.F162M',\n", + " star_mag=10.,\n", + " contrast=[1e-2, 1e-3, 1e-4],\n", + " use_mag=False)\n", + "print(masses)" + ] + }, + { + "cell_type": "markdown", + "id": "8729edb9-9e60-4457-9a83-b1716917490f", + "metadata": {}, + "source": [ + "Several warnings appeared because some of the effective temperatures from the ischrone data are outside the available range of the associated atmospheric model. Those magnitudes are set to NaN such that converting any contrast values into masses for objects in the $\\sim$3000 K regime will not be possible. Any contrast values in that regime will be returned as NaN as well, which is not the case in this example." + ] } ], "metadata": { diff --git a/species/data/model_data.json b/species/data/model_data.json index 1e01afd5..25ca7363 100644 --- a/species/data/model_data.json +++ b/species/data/model_data.json @@ -214,7 +214,7 @@ }, "petitcode-linder2019-clear": { "parameters": ["teff", "logg", "feh"], - "name": "petitCODE (Linder et al. 2019)", + "name": "petitCODE clear (Linder et al. 2019)", "file size": "225 MB", "wavelength range": [0.1, 250], "resolution": 1000, @@ -224,7 +224,7 @@ }, "petitcode-linder2019-cloudy": { "parameters": ["teff", "logg", "feh", "fsed"], - "name": "petitCODE hot cloudy", + "name": "petitCODE cloudy (Linder et al. 2019)", "file size": "1.4 GB", "wavelength range": [0.1, 250], "resolution": 1000, diff --git a/species/plot/plot_comparison.py b/species/plot/plot_comparison.py index d6abb8fb..c38c7b30 100644 --- a/species/plot/plot_comparison.py +++ b/species/plot/plot_comparison.py @@ -530,6 +530,7 @@ def plot_grid_statistic( h5_file = h5py.File(db_path, "r") dset = h5_file[f"results/comparison/{tag}/goodness_of_fit"] + goodness_fit = np.array(dset) n_param = dset.attrs["n_param"] parallax = dset.attrs["parallax"] @@ -552,8 +553,6 @@ def plot_grid_statistic( for item in read_obj.get_spectrum().values(): n_wavel += item[0].shape[0] - goodness_fit = np.array(dset) - model_param = [] coord_points = [] @@ -702,6 +701,7 @@ def plot_grid_statistic( goodness_fit = np.nanmin(goodness_fit, axis=tuple(ax_list)) extra_map = np.zeros(goodness_fit.shape[:2]) + for i in range(goodness_fit.shape[0]): for j in range(goodness_fit.shape[1]): # Get the indices in the goodness_full array @@ -709,6 +709,10 @@ def plot_grid_statistic( # goodness_fit array min_idx = np.argwhere(goodness_fit[i, j] == goodness_full) + if len(min_idx) == 0: + extra_map[i, j] = np.nan + continue + if len(min_idx) > 1: warnings.warn( f"Found {len(min_idx)} positions in the " diff --git a/species/read/read_isochrone.py b/species/read/read_isochrone.py index c5896ea4..84023446 100644 --- a/species/read/read_isochrone.py +++ b/species/read/read_isochrone.py @@ -6,7 +6,7 @@ import os import warnings -from typing import Dict, List, Optional, Tuple +from typing import Dict, List, Optional, Tuple, Union import h5py import numpy as np @@ -16,7 +16,7 @@ from species.core import box from species.read import read_model -from species.util import plot_util +from species.util import phot_util, plot_util class ReadIsochrone: @@ -36,7 +36,10 @@ class ReadIsochrone: @typechecked def __init__( - self, tag: Optional[str] = None, create_regular_grid: bool = False, extrapolate: bool = False + self, + tag: Optional[str] = None, + create_regular_grid: bool = False, + extrapolate: bool = False, ) -> None: """ Parameters @@ -106,11 +109,13 @@ def __init__( with h5py.File(self.database, "r") as h5_file: tag_list = list(h5_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") + 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 h5_file: if f"isochrones/{self.tag}" not in h5_file: @@ -122,7 +127,15 @@ def __init__( f"tags are found in the database: {tag_list}" ) - self.mag_models = ["ames", "atmo", "baraffe", "bt-settl", "linder2019", "manual", "nextgen"] + self.mag_models = [ + "ames", + "atmo", + "baraffe", + "bt-settl", + "linder2019", + "manual", + "nextgen", + ] # Connect isochrone model with atmosphere model # key = isochrone model, value = atmosphere model @@ -159,11 +172,11 @@ def __init__( tag_split = self.tag.split("_") if len(tag_split) == 2: - self.extra_param['feh'] = float(tag_split[1]) + 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]) + 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}") @@ -1547,3 +1560,195 @@ def get_spectrum( ) return model_box + + @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, + ) -> 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 `_. 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 + `_. 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. + + 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 = phot_util.apparent_to_absolute((app_mag, None), (distance, None))[0] + + filter_list = self.get_filters() + + if filter_name in filter_list: + 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) + + mass_interp = interpolate.interp1d( + iso_box.magnitude, iso_box.mass, bounds_error=False, fill_value=np.nan + ) + + mass_array = mass_interp(abs_mag) + + else: + 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 = read_model.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 + ) + + mass_interp = interpolate.interp1d( + model_abs_mag, iso_box.mass, bounds_error=False, fill_value=np.nan + ) + + mass_array = mass_interp(abs_mag) + + return mass_array diff --git a/species/read/read_model.py b/species/read/read_model.py index e4007d00..916eecf3 100644 --- a/species/read/read_model.py +++ b/species/read/read_model.py @@ -1875,7 +1875,7 @@ 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 + 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.