From d1c24befac05a6d94179cd438df848e5d284c1f0 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Wed, 6 May 2020 10:49:39 +0200 Subject: [PATCH] Optical constants and reddening vector in color-mag/color diagrams --- docs/overview.rst | 1 + docs/species.data.rst | 8 ++ requirements.txt | 3 +- species/analysis/photometry.py | 2 +- species/data/companions.py | 6 +- species/data/database.py | 39 ++++++- species/data/dust.py | 92 +++++++++++++++ species/plot/plot_color.py | 199 +++++++++++++++++++++++++++------ species/plot/plot_spectrum.py | 28 ++++- species/read/read_color.py | 8 +- species/util/plot_util.py | 170 +++++++++++++++++++++++++++- 11 files changed, 502 insertions(+), 54 deletions(-) create mode 100644 species/data/dust.py diff --git a/docs/overview.rst b/docs/overview.rst index 657d236c..9268fffa 100644 --- a/docs/overview.rst +++ b/docs/overview.rst @@ -30,6 +30,7 @@ The following data are currently supported by *species* (support for other data - Photometry from `Sandy Leggett `_ - Photometry of directly imaged planets and brown dwarfs (see dictionary in :class:`~species.data.companions`) - Calibration spectrum of `Vega `_ +- Optical constants compiled by `Mollière et al. (2019) `_ (see :func:`~species.data.database.Database.add_dust`) Please give credit to the relevant authors when using these data in a publication. More information is available on the respective websites. diff --git a/docs/species.data.rst b/docs/species.data.rst index de2b9fca..ec891d32 100644 --- a/docs/species.data.rst +++ b/docs/species.data.rst @@ -60,6 +60,14 @@ species.data.drift\_phoenix module :undoc-members: :show-inheritance: +species.data.dust module +------------------------ + +.. automodule:: species.data.dust + :members: + :undoc-members: + :show-inheritance: + species.data.exo\_rem module ---------------------------- diff --git a/requirements.txt b/requirements.txt index a9210f65..b43ff2d2 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,9 +6,10 @@ h5py ~= 2.10 matplotlib ~= 3.2 numpy ~= 1.18 pandas ~= 0.25 +PyMieScatt ~= 1.7 pymultinest ~= 2.9 scipy ~= 1.4 spectres ~= 2.1 -tqdm ~= 4.45 +tqdm ~= 4.46 typeguard ~= 2.7 xlrd ~= 1.2 diff --git a/species/analysis/photometry.py b/species/analysis/photometry.py index 2a6ef319..db97cfc9 100644 --- a/species/analysis/photometry.py +++ b/species/analysis/photometry.py @@ -39,7 +39,7 @@ def __init__(self, self.filter_interp = None self.wavel_range = None - self.vega_mag = 0.03 # [mag] + self.vega_mag = 0.03 # (mag) config_file = os.path.join(os.getcwd(), 'species_config.ini') diff --git a/species/data/companions.py b/species/data/companions.py index 2a07d516..aa0c4f1d 100644 --- a/species/data/companions.py +++ b/species/data/companions.py @@ -104,9 +104,9 @@ def get_data(): 'Paranal/SPHERE.IRDIS_D_H23_3': (17.95, 0.17), # Keppler et al. 2018 'Paranal/SPHERE.IRDIS_D_K12_1': [(16.78, 0.31), (16.72, 0.05)], # Stolker et al. in prep. 'Paranal/SPHERE.IRDIS_D_K12_2': [(16.23, 0.32), (16.38, 0.06)], # Stolker et al. in prep. - # 'Paranal/NACO.Lp': (14.08, 0.33), # Stolker et al. in prep. - # 'Paranal/NACO.NB405': (13.91, 0.34), # Stolker et al. in prep. - # 'Paranal/NACO.Mp': (13.64, 0.22), # Stolker et al. in prep. + 'Paranal/NACO.Lp': (14.08, 0.33), # Stolker et al. in prep. + 'Paranal/NACO.NB405': (13.91, 0.34), # Stolker et al. in prep. + 'Paranal/NACO.Mp': (13.64, 0.22), # Stolker et al. in prep. 'Keck/NIRC2.Lp': (14.64, 0.18)}}, # Wang et al. 2020 'PDS 70 c': {'distance': (113.43, 0.52), diff --git a/species/data/database.py b/species/data/database.py index 186cfd2c..787e3b97 100644 --- a/species/data/database.py +++ b/species/data/database.py @@ -25,7 +25,7 @@ from species.core import box, constants from species.data import drift_phoenix, btnextgen, vega, irtf, spex, vlm_plx, leggett, \ companions, filters, btsettl, ames_dusty, ames_cond, \ - isochrones, petitcode, exo_rem + isochrones, petitcode, exo_rem, dust from species.read import read_model, read_calibration, read_planck from species.util import data_util # from species.util import data_util, retrieval_util @@ -164,6 +164,43 @@ def add_companion(self, distance=data[item]['distance'], app_mag=data[item]['app_mag']) + @typechecked + def add_dust(self) -> None: + """ + Function for adding optical constants of MgSiO3 and Fe to the database. The optical + constants have been compiled by Mollière et al. (2019) for petitRADTRANS from the + following sources: + + - MgSiO3, crystalline + - Scott & Duley (1996), ApJS, 105, 401 + - Jäger et al. (1998), A&A, 339, 904 + + - MgSiO3, amorphous + - Jäger et al. (2003), A&A, 408, 193 + + - Fe, crystalline + - Henning & Stognienko (1996), A&A, 311, 291 + + - Fe, amorphous + - Pollack et al. (1994), ApJ, 421, 615 + + Returns + ------- + NoneType + None + """ + + h5_file = h5py.File(self.database, 'a') + + if 'dust' in h5_file: + del h5_file['dust'] + + h5_file.create_group('dust') + + dust.add_optical_constants(self.input_path, h5_file) + + h5_file.close() + def add_filter(self, filter_name, filename=None): diff --git a/species/data/dust.py b/species/data/dust.py new file mode 100644 index 00000000..a2b0e966 --- /dev/null +++ b/species/data/dust.py @@ -0,0 +1,92 @@ +""" +Module for optical constants of dust grains. +""" + +import os +import zipfile +import urllib.request + +import h5py +import numpy as np + +from typeguard import typechecked + + +@typechecked +def add_optical_constants(input_path: str, + database: h5py._hl.files.File) -> None: + """ + Function for adding the optical constants of crystalline and amorphous MgSiO3 and Fe to the + database. + + Parameters + ---------- + input_path : str + Folder where the data is located. + database : h5py._hl.files.File + Database. + + Returns + ------- + None + NoneType + """ + + if not os.path.exists(input_path): + os.makedirs(input_path) + + url = 'https://people.phys.ethz.ch/~stolkert/species/optical_constants.zip' + + data_file = os.path.join(input_path, 'optical_constants.zip') + + if not os.path.isfile(data_file): + print('Downloading optical constants (87 kB)...', end='', flush=True) + urllib.request.urlretrieve(url, data_file) + print(' [DONE]') + + print('Unpacking optical constants...', end='', flush=True) + + with zipfile.ZipFile(data_file, 'r') as zip_ref: + zip_ref.extractall(input_path) + + print(' [DONE]') + + print('Adding optical constants of MgSiO3...', end='') + + nk_file = os.path.join(input_path, 'optical_constants/mgsio3/crystalline/' + 'mgsio3_jaeger_98_scott_96_axis1.dat') + + data = np.loadtxt(nk_file) + database.create_dataset('dust/mgsio3/crystalline/axis_1/', data=data) + + nk_file = os.path.join(input_path, 'optical_constants/mgsio3/crystalline/' + 'mgsio3_jaeger_98_scott_96_axis2.dat') + + data = np.loadtxt(nk_file) + database.create_dataset('dust/mgsio3/crystalline/axis_2/', data=data) + + nk_file = os.path.join(input_path, 'optical_constants/mgsio3/crystalline/' + 'mgsio3_jaeger_98_scott_96_axis3.dat') + + data = np.loadtxt(nk_file) + database.create_dataset('dust/mgsio3/crystalline/axis_3/', data=data) + + nk_file = os.path.join(input_path, 'optical_constants/mgsio3/amorphous/' + 'mgsio3_jaeger_2003_reformat.dat') + + data = np.loadtxt(nk_file) + database.create_dataset('dust/mgsio3/amorphous', data=data) + + print(' [DONE]') + + print('Adding optical constants of Fe...', end='') + + nk_file = os.path.join(input_path, 'optical_constants/fe/crystalline/fe_henning_1996.dat') + data = np.loadtxt(nk_file) + database.create_dataset('dust/fe/crystalline', data=data) + + nk_file = os.path.join(input_path, 'optical_constants/fe/amorphous/fe_pollack_1994.dat') + data = np.loadtxt(nk_file) + database.create_dataset('dust/fe/amorphous', data=data) + + print(' [DONE]') diff --git a/species/plot/plot_color.py b/species/plot/plot_color.py index d3fc9dbd..e8a84e7c 100644 --- a/species/plot/plot_color.py +++ b/species/plot/plot_color.py @@ -5,10 +5,14 @@ import os import math +from typing import Union, Optional, Tuple, List + +import PyMieScatt import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt +from typeguard import typechecked from scipy.interpolate import interp1d from matplotlib.colorbar import Colorbar @@ -23,19 +27,24 @@ plt.rc('axes', edgecolor='black', linewidth=2.2) -def plot_color_magnitude(boxes, - objects=None, - mass_labels=None, - teff_labels=None, - companion_labels=False, - field_range=None, - label_x='Color (mag)', - label_y='Absolute magnitude (mag)', - xlim=None, - ylim=None, - offset=None, - legend='upper left', - output='color-magnitude.pdf'): +@typechecked +def plot_color_magnitude(boxes: list, + objects: Optional[Union[List[Tuple[str, str, str, str]], + List[Tuple[str, str, str, str, Optional[dict], + Optional[dict]]]]] = None, + mass_labels: Optional[Union[List[float], List[Tuple[float, str]]]] = None, + teff_labels: Optional[Union[List[float], List[Tuple[float, str]]]] = None, + companion_labels: bool = False, + reddening: Optional[List[Tuple[Tuple[str, str], Tuple[str, float], str, + float, Tuple[float, float]]]] = None, + field_range: Optional[Tuple[str, str]] = None, + label_x: str = 'Color', + label_y: str = 'Absolute magnitude', + xlim: Optional[Tuple[float, float]] = None, + ylim: Optional[Tuple[float, float]] = None, + offset: Optional[Tuple[float, float]] = None, + legend: Union[str, dict, Tuple[float, float]] = 'upper left', + output: str = 'color-magnitude.pdf') -> None: """ Function for creating a color-magnitude diagram. @@ -46,8 +55,8 @@ def plot_color_magnitude(boxes, libraries, and/or atmospheric models. The synthetic data have to be created with :func:`~species.read.read_isochrone.ReadIsochrone.get_color_magnitude`. These boxes contain synthetic colors and magnitudes for a given age and a range of masses. - objects : tuple(tuple(str, str, str, str), ), - tuple(tuple(str, str, str, str, str, str, dict, dict), ), None + objects : list(tuple(str, str, str, str), ), + list(tuple(str, str, str, str, dict, dict), ), None Tuple with individual objects. The objects require a tuple with their database tag, the two filter names for the color, and the filter names for the absolute magnitude. Optionally, a dictionary with keyword arguments can be provided for the object's marker and label, @@ -66,6 +75,12 @@ def plot_color_magnitude(boxes, to None. companion_labels : bool Plot labels with the names of the directly imaged companions. + reddening : list(tuple(str, str, str, float, str, float, tuple(float, float)), None + Include reddening arrows by providing a list with tuples. Each tuple contains the filter + names for the color, the filter name for the magnitude, the particle radius (um), and the + start position (color, mag) of the arrow in the plot, so (filter_color_1, filter_color_2, + filter_mag, composition, radius, (x_pos, y_pos)). The composition can be either 'Fe' or + 'MgSiO3' (both with crystalline structure). The parameter is not used if set to ``None``. field_range : tuple(str, str), None Range of the discrete colorbar for the field dwarfs. The tuple should contain the lower and upper value ('early M', 'late M', 'early L', 'late L', 'early T', 'late T', 'early Y). @@ -80,8 +95,8 @@ def plot_color_magnitude(boxes, Limits for the y-axis. Not used if set to None. offset : tuple(float, float), None Offset of the x- and y-axis label. - legend : str, None - Legend position. Not shown if set to None. + legend : str, tuple(float, float), dict, None + Legend position or keyword arguments. No legend is shown if set to ``None``. output : str Output filename. @@ -92,8 +107,6 @@ def plot_color_magnitude(boxes, """ - print(f'Plotting color-magnitude diagram: {output}...', end='', flush=True) - model_color = ('#234398', '#f6a432', 'black') model_linestyle = ('-', '--', ':', '-.') @@ -350,6 +363,53 @@ def plot_color_magnitude(boxes, for item in isochrones: ax1.plot(item.color, item.magnitude, linestyle='-', linewidth=1, color='black') + if reddening is not None: + for item in reddening: + ext_1, ext_2 = plot_util.calc_reddening(item[0], + item[1], + composition=item[2], + structure='crystalline', + radius=item[3]) + + delta_x = ext_1 - ext_2 + delta_y = item[1][1] + + x_pos = item[4][0] + delta_x + y_pos = item[4][1] + delta_y + + ax1.annotate(s='', xy=(x_pos, y_pos), xytext=(item[4][0], item[4][1]), + fontsize=8, arrowprops={'arrowstyle': '->'}, color='black', zorder=3.) + + x_pos_text = item[4][0] + delta_x/2. + y_pos_text = item[4][1] + delta_y/2. + + vector_len = math.sqrt(delta_x**2+delta_y**2) + + if item[2] == 'MgSiO3': + dust_species = r'MgSiO$_{3}$' + elif item[2] == 'Fe': + dust_species = 'Fe' + + if (item[3]).is_integer(): + red_label = f'{dust_species} ({item[3]:.0f} $\mu$m)' + else: + red_label = f'{dust_species} ({item[3]:.1f} $\mu$m)' + + text = ax1.annotate(red_label, xy=(x_pos_text, y_pos_text), + xytext=(7.*delta_y/vector_len, 7.*delta_x/vector_len), + textcoords='offset points', fontsize=8., color='black', + ha='center', va='center') + + line, = ax1.plot([item[4][0], x_pos], [item[4][1], y_pos], '-', color='white') + + sp1 = ax1.transData.transform_point((item[4][0], item[4][1])) + sp2 = ax1.transData.transform_point((x_pos, y_pos)) + + angle = np.degrees(np.arctan2(sp2[1]-sp1[1], sp2[0]-sp1[0])) + text.set_rotation(angle) + + print(f'Plotting color-magnitude diagram: {output}...', end='', flush=True) + if objects is not None: for i, item in enumerate(objects): objdata = read_object.ReadObject(item[0]) @@ -408,19 +468,25 @@ def plot_color_magnitude(boxes, print(' [DONE]') -def plot_color_color(boxes, - objects=None, - mass_labels=None, - teff_labels=None, - companion_labels=False, - field_range=None, - label_x='Color (mag)', - label_y='Color (mag)', - xlim=None, - ylim=None, - offset=None, - legend='upper left', - output='color-color.pdf'): +@typechecked +def plot_color_color(boxes: list, + objects: Optional[Union[List[Tuple[str, Tuple[str, str], Tuple[str, str]]], + List[Tuple[str, Tuple[str, str], Tuple[str, str], Optional[dict], + Optional[dict]]]]] = None, + mass_labels: Optional[Union[List[float], List[Tuple[float, str]]]] = None, + teff_labels: Optional[Union[List[float], List[Tuple[float, str]]]] = None, + companion_labels: bool = False, + reddening: Optional[List[Tuple[Tuple[str, str], Tuple[str, str], + Tuple[str, float], str, float, + Tuple[float, float]]]] = None, + field_range: Optional[Tuple[str, str]] = None, + label_x: str = 'Color', + label_y: str = 'Color', + xlim: Optional[Tuple[float, float]] = None, + ylim: Optional[Tuple[float, float]] = None, + offset: Optional[Tuple[float, float]] = None, + legend: Union[str, dict, Tuple[float, float]] = 'upper left', + output: str = 'color-color.pdf') -> None: """ Function for creating a color-color diagram. @@ -431,8 +497,8 @@ def plot_color_color(boxes, libraries, and/or atmospheric models. The synthetic data have to be created with :func:`~species.read.read_isochrone.ReadIsochrone.get_color_color`. These boxes contain synthetic colors for a given age and a range of masses. - objects : tuple(tuple(str, str, str, str), ), - tuple(tuple(str, str, str, str, str, str, dict, dict), ), None + objects : tuple(tuple(str, tuple(str, str), tuple(str, str)), ), + tuple(tuple(str, tuple(str, str), tuple(str, str), dict, dict), ), None Tuple with individual objects. The objects require a tuple with their database tag, the two filter names for the first color, and the two filter names for the second color. Optionally, a dictionary with keyword arguments can be provided for the object's marker and @@ -451,6 +517,12 @@ def plot_color_color(boxes, to None. companion_labels : bool Plot labels with the names of the directly imaged companions. + reddening : list(tuple(tuple(str, str), tuple(str, str), tuple(str, float), str, float, tuple(float, float)), None + Include reddening arrows by providing a list with tuples. Each tuple contains the filter + names for the color, the filter name for the magnitude, the particle radius (um), and the + start position (color, mag) of the arrow in the plot, so (filter_color_1, filter_color_2, + filter_mag, composition, radius, (x_pos, y_pos)). The composition can be either 'Fe' or + 'MgSiO3' (both with crystalline structure). The parameter is not used if set to ``None``. field_range : tuple(str, str), None Range of the discrete colorbar for the field dwarfs. The tuple should contain the lower and upper value ('early M', 'late M', 'early L', 'late L', 'early T', 'late T', 'early Y). @@ -476,8 +548,6 @@ def plot_color_color(boxes, None """ - print(f'Plotting color-color diagram: {output}...', end='', flush=True) - model_color = ('#234398', '#f6a432', 'black') model_linestyle = ('-', '--', ':', '-.') @@ -588,8 +658,8 @@ def plot_color_color(boxes, mass_val = mass_item mass_pos = 'right' - if j == 0: # if j == 0 or (j > 0 and mass_val < 20.): + if j == 0: pos_color1 = interp_color1(mass_val) pos_color2 = interp_color2(mass_val) @@ -722,6 +792,61 @@ def plot_color_color(boxes, for item in isochrones: ax1.plot(item.colors[0], item.colors[1], linestyle='-', linewidth=1, color='black') + if reddening is not None: + for item in reddening: + ext_1, ext_2 = plot_util.calc_reddening(item[0], + item[2], + composition=item[3], + structure='crystalline', + radius=item[4]) + + ext_3, ext_4 = plot_util.calc_reddening(item[1], + item[2], + composition=item[3], + structure='crystalline', + radius=item[4]) + + delta_x = ext_1 - ext_2 + delta_y = ext_3 - ext_4 + + x_pos = item[5][0] + delta_x + y_pos = item[5][1] + delta_y + + ax1.annotate(s='', xy=(x_pos, y_pos), xytext=(item[5][0], item[5][1]), + fontsize=8, arrowprops={'arrowstyle': '->'}, color='black', zorder=3.) + + x_pos_text = item[5][0] + delta_x/2. + y_pos_text = item[5][1] + delta_y/2. + + vector_len = math.sqrt(delta_x**2+delta_y**2) + + if item[3] == 'MgSiO3': + dust_species = r'MgSiO$_{3}$' + + elif item[3] == 'Fe': + dust_species = 'Fe' + + if item[4].is_integer(): + red_label = f'{dust_species} ({item[4]:.0f} $\mu$m)' + + else: + red_label = f'{dust_species} ({item[4]:.1f} $\mu$m)' + + text = ax1.annotate(red_label, xy=(x_pos_text, y_pos_text), + xytext=(-7.*delta_y/vector_len, 7.*delta_x/vector_len), + textcoords='offset points', fontsize=8., color='black', + ha='center', va='center') + + line, = ax1.plot([item[5][0], x_pos], [item[5][1], y_pos], '-', color='white') + + sp1 = ax1.transData.transform_point((item[5][0], item[5][1])) + sp2 = ax1.transData.transform_point((x_pos, y_pos)) + + angle = np.degrees(np.arctan2(sp2[1]-sp1[1], sp2[0]-sp1[0])) + text.set_rotation(angle) + + print(f'Plotting color-color diagram: {output}...', end='', flush=True) + if objects is not None: for i, item in enumerate(objects): objdata = read_object.ReadObject(item[0]) diff --git a/species/plot/plot_spectrum.py b/species/plot/plot_spectrum.py index ee11b551..689443d6 100644 --- a/species/plot/plot_spectrum.py +++ b/species/plot/plot_spectrum.py @@ -46,10 +46,28 @@ def plot_spectrum(boxes, Filter IDs for which the transmission profile is plotted. Not plotted if set to None. residuals : species.core.box.ResidualsBox, None Box with residuals of a fit. Not plotted if set to None. - plot_kwargs : list(str, ), None - TODO Colors to be used for the different boxes. Note that a box with residuals requires a tuple - with two colors (i.e., for the photometry and spectrum). Automatic colors are used if set - to None. + plot_kwargs : list(dict, ), None + List with dictionaries of keyword arguments for each box. For example, if the ``boxes`` + are a ``ModelBox`` and ``ObjectBox``: + + .. code-block:: python + + plot_kwargs=[{'ls': '-', 'lw': 1., 'color': 'black'}, + {'spectrum_1': {'marker': 'o', 'ms': 3., 'color': 'tab:brown', 'ls': 'none'}, + 'spectrum_2': {'marker': 'o', 'ms': 3., 'color': 'tab:blue', 'ls': 'none'}, + 'Paranal/SPHERE.IRDIS_D_H23_3': {'marker': 's', 'ms': 4., 'color': 'tab:cyan', 'ls': 'none'}, + 'Paranal/SPHERE.IRDIS_D_K12_1': [{'marker': 's', 'ms': 4., 'color': 'tab:orange', 'ls': 'none'}, + {'marker': 's', 'ms': 4., 'color': 'tab:red', 'ls': 'none'}], + 'Paranal/NACO.Lp': {'marker': 's', 'ms': 4., 'color': 'tab:green', 'ls': 'none'}, + 'Paranal/NACO.Mp': {'marker': 's', 'ms': 4., 'color': 'tab:green', 'ls': 'none'}}] + + For an ``ObjectBox``, the dictionary contains items for the different spectrum and filter + names stored with :func:`~species.data.database.Database.add_object`. In case both + and ``ObjectBox`` and a ``SynphotBox`` are provided, then the latter can be set to ``None`` + in order to use the same (but open) symbols as the data from the ``ObjectBox``. Note that + if a filter name is duplicated in an ``ObjectBox`` (Paranal/SPHERE.IRDIS_D_K12_1 in the + example) then a list with two dictionaries should be provided. Colors are automatically + chosen if ``plot_kwargs`` is set to ``None``. xlim : tuple(float, float) Limits of the wavelength axis. ylim : tuple(float, float) @@ -82,8 +100,6 @@ def plot_spectrum(boxes, None """ - marker = itertools.cycle(('o', 's', '*', 'p', '<', '>', 'P', 'v', '^')) - if plot_kwargs is None: plot_kwargs = [] diff --git a/species/read/read_color.py b/species/read/read_color.py index 5f6cecfa..2e009c38 100644 --- a/species/read/read_color.py +++ b/species/read/read_color.py @@ -219,9 +219,9 @@ def get_color_color(self, Box with the colors. """ - h5_file = h5py.File(self.database, 'r') - if self.lib_type == 'phot_lib': + h5_file = h5py.File(self.database, 'r') + sptype = np.asarray(h5_file[f'photometry/{self.library}/sptype']) flag = np.asarray(h5_file[f'photometry/{self.library}/flag']) @@ -259,6 +259,8 @@ def get_color_color(self, color2=color2[indices], sptype=sptype[indices]) + h5_file.close() + elif self.lib_type == 'spec_lib': read_spec_0 = read_spectrum.ReadSpectrum(spec_library=self.library, filter_name=self.filters_colors[0][0]) @@ -285,6 +287,4 @@ def get_color_color(self, color2=phot_box_2.app_mag[:, 0]-phot_box_3.app_mag[:, 0], sptype=phot_box_0.sptype) - h5_file.close() - return colorbox diff --git a/species/util/plot_util.py b/species/util/plot_util.py index f8a2b82d..91c052ec 100644 --- a/species/util/plot_util.py +++ b/species/util/plot_util.py @@ -2,8 +2,21 @@ Utility functions for plotting data. """ +import os +import configparser + +from typing import Tuple + +import h5py +import PyMieScatt import numpy as np +from typeguard import typechecked +from scipy.interpolate import interp1d + +from species.data import database +from species.read import read_filter + def sptype_substellar(sptype, shape): @@ -214,8 +227,23 @@ def model_name(key): elif key == 'bt-nextgen': name = 'BT-NextGen' + elif key == 'petitcode-cool-clear': + name = 'petitCODE' + + elif key == 'petitcode-cool-cloudy': + name = 'petitCODE' + + elif key == 'petitcode-hot-clear': + name = 'petitCODE' + + elif key == 'petitcode-hot-cloudy': + name = 'petitCODE' + + elif key == 'exo-rem': + name = 'Exo-REM' + elif key == 'planck': - name = 'Planck radiation' + name = 'Blackbody radiation' elif key == 'zhu2015': name = 'Zhu (2015)' @@ -360,3 +388,143 @@ def field_bounds_ticks(field_range): labels = spectral_ranges[index_start:index_end] return bounds, ticks, labels + + +@typechecked +def dust_cross_section(wavelengths: np.ndarray, + n_index: np.ndarray, + k_index: np.ndarray, + radius: float) -> np.ndarray: + """ + Function for calculating the extinction cross section of a dust grain. + + Parameters + ---------- + wavelengths : np.ndarray + Wavelengths (um). + n_index : np.ndarray + Real part of the refractive index. + k_index : np.ndarray + Imaginary part of the refractive index. + radius : np.ndarray + Radius of the dust grain (um). + + Returns + ------- + np.ndarray + Extinction cross section (um2) + """ + + sigma = np.zeros(wavelengths.shape) + + for i, item in enumerate(wavelengths): + # PyMieScatt units are in nm and the grain size is provided as the diameter + mie = PyMieScatt.MieQ(complex(n_index[i], k_index[i]), item*1e3, 2.*radius*1e3, + asDict=True, asCrossSection=True) + + if 'Cext' in mie: + sigma[i] = mie['Cext'] # (nm2) + + return sigma*1e-6 # (um2) + + +@typechecked +def calc_reddening(filters_color: Tuple[str, str], + extinction: Tuple[str, float], + composition: str = 'MgSiO3', + structure: str = 'crystalline', + radius: float = 1.) -> Tuple[float, float]: + """ + Function for calculating the reddening of a color given the extinction for a given filter. + + Parameters + ---------- + filters_color : tuple(str, str) + Filter names for which the extinction is calculated. + extinction : str + Filter name and extinction (mag). + composition : str + Dust composition ('MgSiO3' or 'Fe'). + structure : str + Grain structure ('crystalline' or 'amorphous'). + radius : float + Radius of the dust grain (um). + + Returns + ------- + float + Extinction (mag) for ``filters_color[0]``. + float + Extinction (mag) for ``filters_color[1]``. + """ + + config_file = os.path.join(os.getcwd(), 'species_config.ini') + + config = configparser.ConfigParser() + config.read_file(open(config_file)) + + database_path = config['species']['database'] + + h5_file = h5py.File(database_path, 'r') + + try: + h5_file['dust'] + + except KeyError: + h5_file.close() + species_db = database.Database() + species_db.add_dust() + h5_file = h5py.File(database_path, 'r') + + if composition == 'MgSiO3' and structure == 'crystalline': + for i in range(3): + data = h5_file[f'dust/mgsio3/crystalline/axis_{i+1}'] + + # Average cross section of the three axes + if i == 0: + sigma = dust_cross_section(data[:, 0], data[:, 1], data[:, 2], radius) / 3. + + else: + sigma += dust_cross_section(data[:, 0], data[:, 1], data[:, 2], radius) / 3. + + else: + if composition == 'MgSiO3' and structure == 'amorphous': + data = h5_file[f'dust/mgsio3/amorphous/'] + + elif composition == 'Fe' and structure == 'crystalline': + data = h5_file[f'dust/fe/crystalline/'] + + elif composition == 'Fe' and structure == 'amorphous': + data = h5_file[f'dust/fe/amorphous/'] + + sigma = dust_cross_section(data[:, 0], data[:, 1], data[:, 2], radius) + + interp_sigma = interp1d(data[:, 0], sigma, kind='linear') + + h5_file.close() + + read_filt = read_filter.ReadFilter(extinction[0]) + transmission = read_filt.get_filter() + + # Weighted average of the cross section for extinction[0] + sigma_mag = np.trapz(interp_sigma(transmission[0, ])*transmission[1, ], + transmission[0, ]) / np.trapz(transmission[1, ], transmission[0, ]) + + read_filt = read_filter.ReadFilter(filters_color[0]) + transmission = read_filt.get_filter() + + # Weighted average of the cross section for filters_color[0] + sigma_color_0 = np.trapz(interp_sigma(transmission[0, ])*transmission[1, ], + transmission[0, ]) / np.trapz(transmission[1, ], transmission[0, ]) + + read_filt = read_filter.ReadFilter(filters_color[1]) + transmission = read_filt.get_filter() + + # Weighted average of the cross section for filters_color[1] + sigma_color_1 = np.trapz(interp_sigma(transmission[0, ])*transmission[1, ], + transmission[0, ]) / np.trapz(transmission[1, ], transmission[0, ]) + + density = extinction[1]/sigma_mag/2.5/np.log10(np.exp(1.)) + + return 2.5 * np.log10(np.exp(1.)) * sigma_color_0 * density, \ + 2.5 * np.log10(np.exp(1.)) * sigma_color_1 * density