From 24c130448938064189098a7c3a90bbd9e4ffc26d Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Mon, 15 Apr 2024 10:24:32 +0200 Subject: [PATCH] Support for fitting multiple disk_teff and disk_radius parameters --- species/data/spec_data/spec_vega.py | 2 +- species/fit/fit_model.py | 183 +++++++++++++++++++++++----- species/plot/plot_mcmc.py | 68 +++++++++++ species/read/read_model.py | 140 +++++++++++++++++---- species/util/plot_util.py | 63 +++++++++- 5 files changed, 399 insertions(+), 57 deletions(-) diff --git a/species/data/spec_data/spec_vega.py b/species/data/spec_data/spec_vega.py index 4b14e5b..2a94bb2 100644 --- a/species/data/spec_data/spec_vega.py +++ b/species/data/spec_data/spec_vega.py @@ -49,7 +49,7 @@ def add_vega(input_path, database): progressbar=True, ) - except requests.exceptions.HTTPError: + except (requests.exceptions.HTTPError, requests.exceptions.ConnectionError): url = ( "https://home.strw.leidenuniv.nl/~stolker/" "species/alpha_lyr_stis_011.fits" diff --git a/species/fit/fit_model.py b/species/fit/fit_model.py index e60fea9..b4272a6 100644 --- a/species/fit/fit_model.py +++ b/species/fit/fit_model.py @@ -188,9 +188,34 @@ def __init__( which adds a constant flux (in W m-2 um-1) to the model spectrum. - Blackbody parameters (with ``model='planck'``): + Blackbody disk emission: + + - Blackbody parameters can be fitted to account for + thermal emission from one or multiple disk + components, in addition to the atmospheric + emission. These parameters should therefore be + combined with an atmospheric model. - - Parameter boundaries have to be provided for 'teff' + - Parameter boundaries have to be provided for + 'disk_teff' and 'disk_radius'. For example, + ``bounds={'teff': (2000., 3000.), 'radius': (1., 5.), + 'logg': (3.5, 4.5), 'disk_teff': (100., 2000.), + 'disk_radius': (1., 100.)}`` for fitting a single + blackbody component, in addition to the atmospheric + parameters. Or, ``bounds={'teff': (2000., 3000.), + 'radius': (1., 5.), 'logg': (3.5, 4.5), + 'disk_teff': [(2000., 500.), (1000., 20.)], + 'disk_radius': [(1., 100.), (50., 1000.)]}`` for + fitting two blackbody components. Any number of + blackbody components can be fitted by including + additional priors in the lists of ``'disk_teff'`` + and ``'disk_radius'``. + + Blackbody parameters (only with ``model='planck'``): + + - This implementation fits both the atmospheric emission + and possible disk emission with blackbody components. + Parameter boundaries have to be provided for 'teff' and 'radius'. - For a single blackbody component, the values are @@ -363,18 +388,6 @@ def __init__( ``powerlaw_ext``, and a log-uniform prior for ``powerlaw_max``. - Blackbody disk emission: - - - Additional blackbody emission can be added to the - atmospheric spectrum to account for thermal emission - from a disk. - - - Parameter boundaries have to be provided for - 'disk_teff' and 'disk_radius'. For example, - ``bounds={'teff': (2000., 3000.), 'radius': (1., 5.), - 'logg': (3.5, 4.5), 'disk_teff': (100., 2000.), - 'disk_radius': (1., 100.)}``. - inc_phot : bool, list(str) Include photometric data in the fit. If a boolean, either all (``True``) or none (``False``) of the data are @@ -649,8 +662,32 @@ def __init__( self.n_planck = 0 if "disk_teff" in self.bounds and "disk_radius" in self.bounds: - self.modelpar.append("disk_teff") - self.modelpar.append("disk_radius") + if isinstance(bounds["disk_teff"], list) and isinstance( + bounds["disk_radius"], list + ): + # Update temperature and radius parameters + # in case of multiple disk components + self.n_disk = len(bounds["disk_teff"]) + + for i, item in enumerate(bounds["disk_teff"]): + self.modelpar.append(f"disk_teff_{i}") + self.modelpar.append(f"disk_radius_{i}") + + self.bounds[f"disk_teff_{i}"] = bounds["disk_teff"][i] + self.bounds[f"disk_radius_{i}"] = bounds["disk_radius"][i] + + del bounds["disk_teff"] + del bounds["disk_radius"] + + else: + # Fitting a single disk component + self.n_disk = 1 + + self.modelpar.append("disk_teff") + self.modelpar.append("disk_radius") + + else: + self.n_disk = 0 if self.binary: # Update list of model parameters @@ -847,7 +884,7 @@ def __init__( self.diskphot = [] self.diskspec = [] - if "disk_teff" in self.bounds and "disk_radius" in self.bounds: + if self.n_disk > 0: for item in inc_phot: print(f"Interpolating {item}...", end="", flush=True) readmodel = ReadModel("blackbody", filter_name=item) @@ -1182,7 +1219,6 @@ def _lnlike_func( corr_len = {} corr_amp = {} dust_param = {} - disk_param = {} veil_param = {} param_dict = {} rad_vel = {} @@ -1223,12 +1259,6 @@ def _lnlike_func( elif self.ext_filter is not None and item == f"phot_ext_{self.ext_filter}": dust_param[item] = params[self.cube_index[item]] - elif item == "disk_teff": - disk_param["teff"] = params[self.cube_index[item]] - - elif item == "disk_radius": - disk_param["radius"] = params[self.cube_index[item]] - elif item == "veil_a": veil_param["veil_a"] = params[self.cube_index[item]] @@ -1244,6 +1274,41 @@ def _lnlike_func( else: param_dict[item] = params[self.cube_index[item]] + # Disk parameters + + disk_param = {} + + if self.n_disk == 1: + if "disk_teff" in self.fix_param: + disk_param["teff"] = self.fix_param["disk_teff"] + else: + disk_param["teff"] = params[self.cube_index["disk_teff"]] + + if "disk_radius" in self.fix_param: + disk_param["radius"] = self.fix_param["disk_radius"] + else: + disk_param["radius"] = params[self.cube_index["disk_radius"]] + + elif self.n_disk > 1: + for disk_idx in range(self.n_disk): + if f"disk_teff_{disk_idx}" in self.fix_param: + disk_param[f"teff_{disk_idx}"] = self.fix_param[ + f"disk_teff_{disk_idx}" + ] + else: + disk_param[f"teff_{disk_idx}"] = params[ + self.cube_index[f"disk_teff_{disk_idx}"] + ] + + if f"disk_radius_{disk_idx}" in self.fix_param: + disk_param[f"radius_{disk_idx}"] = self.fix_param[ + f"disk_radius_{disk_idx}" + ] + else: + disk_param[f"radius_{disk_idx}"] = params[ + self.cube_index[f"disk_radius_{disk_idx}"] + ] + # Add the parallax manually because it should # not be provided in the bounds dictionary @@ -1285,18 +1350,14 @@ def _lnlike_func( elif item[:9] == "phot_ext_": dust_param[item] = self.fix_param[item] - elif item == "disk_teff": - disk_param["teff"] = self.fix_param[item] - - elif item == "disk_radius": - disk_param["radius"] = self.fix_param[item] - elif item == "spec_weight": pass else: param_dict[item] = self.fix_param[item] + # Check if the blackbody temperatures/radii are decreasing/increasing + if self.model == "planck" and self.n_planck > 1: for i in range(self.n_planck - 1): if param_dict[f"teff_{i+1}"] > param_dict[f"teff_{i}"]: @@ -1305,13 +1366,35 @@ def _lnlike_func( if param_dict[f"radius_{i}"] > param_dict[f"radius_{i+1}"]: return -np.inf - if disk_param: + if self.n_disk == 1: if disk_param["teff"] > param_dict["teff"]: return -np.inf if disk_param["radius"] < param_dict["radius"]: return -np.inf + elif self.n_disk > 1: + for disk_idx in range(self.n_disk): + if disk_idx == 0: + if disk_param["teff_0"] > param_dict["teff"]: + return -np.inf + + if disk_param["radius_0"] < param_dict["radius"]: + return -np.inf + + else: + if ( + disk_param[f"teff_{disk_idx}"] + > disk_param[f"teff_{disk_idx-1}"] + ): + return -np.inf + + if ( + disk_param[f"radius_{disk_idx}"] + < disk_param[f"radius_{disk_idx-1}"] + ): + return -np.inf + if self.model != "powerlaw": if "radius_0" in param_dict and "radius_1" in param_dict: flux_scaling_0 = (param_dict["radius_0"] * constants.R_JUP) ** 2 / ( @@ -1503,7 +1586,9 @@ def _lnlike_func( phot_flux *= flux_scaling phot_flux += flux_offset - if disk_param: + # Add blackbody flux from disk components + + if self.n_disk == 1: phot_tmp = self.diskphot[i].spectrum_interp([disk_param["teff"]])[0][0] phot_flux += ( @@ -1512,6 +1597,20 @@ def _lnlike_func( / (1e3 * constants.PARSEC / parallax) ** 2 ) + elif self.n_disk > 1: + for disk_idx in range(self.n_disk): + phot_tmp = self.diskphot[i].spectrum_interp( + [disk_param[f"teff_{disk_idx}"]] + )[0][0] + + phot_flux += ( + phot_tmp + * (disk_param[f"radius_{disk_idx}"] * constants.R_JUP) ** 2 + / (1e3 * constants.PARSEC / parallax) ** 2 + ) + + # Apply extinction + if "lognorm_ext" in dust_param: cross_tmp = self.cross_sections[phot_filter]( (10.0 ** dust_param["lognorm_radius"], dust_param["lognorm_sigma"]) @@ -1766,7 +1865,9 @@ def _lnlike_func( self.spectrum[item][1] * sigma_i * sigma_j ) - if disk_param: + # Add blackbody flux from disk components + + if self.n_disk == 1: model_tmp = self.diskspec[i].spectrum_interp([disk_param["teff"]])[0, :] model_tmp *= (disk_param["radius"] * constants.R_JUP) ** 2 / ( @@ -1775,6 +1876,20 @@ def _lnlike_func( model_flux += model_tmp + elif self.n_disk > 1: + for disk_idx in range(self.n_disk): + model_tmp = self.diskspec[i].spectrum_interp( + [disk_param[f"teff_{disk_idx}"]] + )[0, :] + + model_tmp *= ( + disk_param[f"radius_{disk_idx}"] * constants.R_JUP + ) ** 2 / (1e3 * constants.PARSEC / parallax) ** 2 + + model_flux += model_tmp + + # Apply extinction + if "lognorm_ext" in dust_param: cross_tmp = self.cross_sections["spectrum"]( ( @@ -1821,6 +1936,8 @@ def _lnlike_func( model_flux *= 10.0 ** (-0.4 * ext_spec) + # Calculate the likelihood + if self.spectrum[item][2] is not None: # Use the inverted covariance matrix ln_like += -0.5 * np.dot( diff --git a/species/plot/plot_mcmc.py b/species/plot/plot_mcmc.py index 62c918d..afe8cf3 100644 --- a/species/plot/plot_mcmc.py +++ b/species/plot/plot_mcmc.py @@ -516,7 +516,19 @@ def plot_posterior( / constants.L_SUN ) + n_disk = 0 + if "disk_teff" in box.parameters and "disk_radius" in box.parameters: + n_disk = 1 + + else: + for disk_idx in range(100): + if f"disk_teff_{disk_idx}" in box.parameters and f"disk_radius_{disk_idx}" in box.parameters: + n_disk += 1 + else: + break + + if n_disk == 1: teff_index = np.argwhere(np.array(box.parameters) == "disk_teff")[0] radius_index = np.argwhere(np.array(box.parameters) == "disk_radius")[0] @@ -547,10 +559,50 @@ def plot_posterior( * samples[..., teff_index] ** 4 ) ) + samples = np.append(samples, radius_bb / constants.R_JUP, axis=-1) box.parameters.append("radius_bb") ndim += 1 + elif n_disk > 1: + lum_disk = 0.0 + + for disk_idx in range(n_disk): + teff_index = np.argwhere(np.array(box.parameters) == f"disk_teff_{disk_idx}")[0] + radius_index = np.argwhere(np.array(box.parameters) == f"disk_radius_{disk_idx}")[0] + + lum_disk += ( + 4.0 + * np.pi + * (samples[..., radius_index] * constants.R_JUP) ** 2 + * constants.SIGMA_SB + * samples[..., teff_index] ** 4.0 + / constants.L_SUN + ) + + radius_bb = np.sqrt( + lum_planet + * constants.L_SUN + / ( + 16.0 + * np.pi + * constants.SIGMA_SB + * samples[..., teff_index] ** 4 + ) + ) + + samples = np.append(samples, radius_bb / constants.R_JUP, axis=-1) + box.parameters.append(f"radius_bb_{disk_idx}") + ndim += 1 + + samples = np.append(samples, np.log10(lum_planet + lum_disk), axis=-1) + box.parameters.append("luminosity") + ndim += 1 + + samples = np.append(samples, lum_disk / lum_planet, axis=-1) + box.parameters.append("luminosity_disk_planet") + ndim += 1 + else: samples = np.append(samples, np.log10(lum_planet), axis=-1) box.parameters.append("luminosity") @@ -690,11 +742,27 @@ def plot_posterior( if object_type == "star": samples[:, radius_index] *= constants.R_JUP / constants.AU + for disk_idx in range(100): + if f"disk_radius_{disk_idx}" in box.parameters: + radius_index = np.argwhere(np.array(box.parameters) == f"disk_radius_{disk_idx}")[0] + if object_type == "star": + samples[:, radius_index] *= constants.R_JUP / constants.AU + else: + break + if "radius_bb" in box.parameters: radius_index = np.argwhere(np.array(box.parameters) == "radius_bb")[0] if object_type == "star": samples[:, radius_index] *= constants.R_JUP / constants.AU + for disk_idx in range(100): + if f"radius_bb_{disk_idx}" in box.parameters: + radius_index = np.argwhere(np.array(box.parameters) == f"radius_bb_{disk_idx}")[0] + if object_type == "star": + samples[:, radius_index] *= constants.R_JUP / constants.AU + else: + break + # Include the log-likelihood if inc_loglike: diff --git a/species/read/read_model.py b/species/read/read_model.py index 6b33102..69439f3 100644 --- a/species/read/read_model.py +++ b/species/read/read_model.py @@ -122,6 +122,10 @@ def __init__( "vsini", ] + for i in range(10): + self.extra_param.append(f"disk_teff_{i}") + self.extra_param.append(f"disk_radius_{i}") + # Test if the spectra are present in the database hdf5_file = self.open_database() hdf5_file.close() @@ -884,7 +888,24 @@ def get_model( # Add blackbody disk component to the spectrum + n_disk = 0 + if "disk_teff" in model_param and "disk_radius" in model_param: + n_disk = 1 + + else: + for disk_idx in range(100): + if f"disk_teff_{disk_idx}" in model_param and f"disk_radius_{disk_idx}" in model_param: + n_disk += 1 + else: + break + + if n_disk == 1: + + readplanck = ReadPlanck( + (0.9 * self.wavel_range[0], 1.1 * self.wavel_range[-1]) + ) + disk_param = { "teff": model_param["disk_teff"], "radius": model_param["disk_radius"], @@ -896,10 +917,6 @@ def get_model( 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( @@ -907,9 +924,30 @@ def get_model( ) flux += flux_interp(self.wl_points) - # flux += spectres.spectres( - # self.wl_points, planck_box.wavelength, planck_box.flux - # ) + elif n_disk > 1: + + readplanck = ReadPlanck( + (0.9 * self.wavel_range[0], 1.1 * self.wavel_range[-1]) + ) + + for disk_idx in range(n_disk): + disk_param = { + "teff": model_param[f"disk_teff_{disk_idx}"], + "radius": model_param[f"disk_radius_{disk_idx}"], + } + + if "parallax" in model_param: + disk_param["parallax"] = model_param["parallax"] + + elif "distance" in model_param: + disk_param["distance"] = model_param["distance"] + + 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) # Create ModelBox with the spectrum @@ -1136,12 +1174,10 @@ def get_model( / constants.L_SUN ) # (Lsun) - # Add the blackbody disk component to the luminosity + # Add the blackbody disk components to the luminosity + + if n_disk == 1: - if ( - "disk_teff" in model_box.parameters - and "disk_radius" in model_box.parameters - ): model_box.parameters["luminosity"] += ( 4.0 * np.pi @@ -1151,6 +1187,18 @@ def get_model( / constants.L_SUN ) # (Lsun) + elif n_disk > 1: + + for disk_idx in range(n_disk): + model_box.parameters["luminosity"] += ( + 4.0 + * np.pi + * (model_box.parameters[f"disk_radius_{disk_idx}"] * constants.R_JUP) ** 2 + * constants.SIGMA_SB + * model_box.parameters[f"disk_teff_{disk_idx}"] ** 4.0 + / constants.L_SUN + ) # (Lsun) + # Add the planet mass to the parameter dictionary if "radius" in model_param and "logg" in model_param: @@ -1338,7 +1386,26 @@ def get_data( # Add blackbody disk component to the spectrum + # Add blackbody disk component to the spectrum + + n_disk = 0 + if "disk_teff" in model_param and "disk_radius" in model_param: + n_disk = 1 + + else: + for disk_idx in range(100): + if "disk_teff_{disk_idx}" in model_param: + n_disk += 1 + else: + break + + if n_disk == 1: + + readplanck = ReadPlanck( + (0.9 * self.wavel_range[0], 1.1 * self.wavel_range[-1]) + ) + disk_param = { "teff": model_param["disk_teff"], "radius": model_param["disk_radius"], @@ -1350,10 +1417,6 @@ def get_data( 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( @@ -1361,7 +1424,30 @@ def get_data( ) flux += flux_interp(wl_points) - # flux += spectres.spectres(wl_points, planck_box.wavelength, planck_box.flux) + elif n_disk > 1: + + readplanck = ReadPlanck( + (0.9 * self.wavel_range[0], 1.1 * self.wavel_range[-1]) + ) + + for disk_idx in range(n_disk): + disk_param = { + "teff": model_param[f"disk_teff_{disk_idx}"], + "radius": model_param[f"disk_radius_{disk_idx}"], + } + + if "parallax" in model_param: + disk_param["parallax"] = model_param["parallax"] + + elif "distance" in model_param: + disk_param["distance"] = model_param["distance"] + + 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) # Create ModelBox with the spectrum @@ -1472,12 +1558,10 @@ def get_data( / constants.L_SUN ) # (Lsun) - # Add the blackbody disk component to the luminosity + # Add the blackbody disk components to the luminosity + + if n_disk == 1: - if ( - "disk_teff" in model_box.parameters - and "disk_radius" in model_box.parameters - ): model_box.parameters["luminosity"] += ( 4.0 * np.pi @@ -1487,6 +1571,18 @@ def get_data( / constants.L_SUN ) # (Lsun) + elif n_disk > 1: + + for disk_idx in range(n_disk): + model_box.parameters["luminosity"] += ( + 4.0 + * np.pi + * (model_box.parameters[f"disk_radius_{disk_idx}"] * constants.R_JUP) ** 2 + * constants.SIGMA_SB + * model_box.parameters[f"disk_teff"] ** 4.0 + / constants.L_SUN + ) # (Lsun) + return model_box @typechecked diff --git a/species/util/plot_util.py b/species/util/plot_util.py index 514bcf7..72acab7 100644 --- a/species/util/plot_util.py +++ b/species/util/plot_util.py @@ -690,6 +690,14 @@ def update_labels(param: List[str], object_type: str = "planet") -> List[str]: index = param.index("disk_teff") param[index] = r"$T_\mathrm{disk}$ (K)" + for i in range(100): + if f"disk_teff_{i}" in param: + index = param.index(f"disk_teff_{i}") + param[index] = r"$T_\mathrm{disk,}$" + rf"$_\mathrm{{{i+1}}}$ (K)" + + else: + break + if "disk_radius" in param: index = param.index("disk_radius") if object_type == "planet": @@ -697,6 +705,19 @@ def update_labels(param: List[str], object_type: str = "planet") -> List[str]: elif object_type == "star": param[index] = r"$R_\mathrm{disk}$ (au)" + for i in range(100): + if f"disk_radius_{i}" in param: + index = param.index(f"disk_radius_{i}") + if object_type == "planet": + param[index] = ( + r"$R_\mathrm{disk,}$" + rf"$_{{{i+1}}}$ ($R_\mathrm{{J}}$)" + ) + elif object_type == "star": + param[index] = r"$R_\mathrm{disk,}$" + rf"$_\mathrm{{{i+1}}}$ (au)" + + else: + break + if "radius_bb" in param: index = param.index("radius_bb") if object_type == "planet": @@ -704,6 +725,19 @@ def update_labels(param: List[str], object_type: str = "planet") -> List[str]: elif object_type == "star": param[index] = r"$R_\mathrm{bb}$ (au)" + for i in range(100): + if f"radius_bb_{i}" in param: + index = param.index(f"radius_bb_{i}") + if object_type == "planet": + param[index] = ( + r"$R_\mathrm{bb,}$" + rf"$_{{{i+1}}}$ ($R_\mathrm{{J}}$)" + ) + elif object_type == "star": + param[index] = r"$R_\mathrm{bb,}$" + rf"$_\mathrm{{{i+1}}}$ (au)" + + else: + break + if "log_powerlaw_a" in param: index = param.index("log_powerlaw_a") param[index] = r"$a_\mathrm{powerlaw}$" @@ -1046,14 +1080,41 @@ def quantity_unit( unit.append("K") label.append(r"$T_\mathrm{disk}$") + for i in range(100): + if item == f"disk_teff_{i}": + quantity.append(f"disk_teff_{i}") + unit.append("K") + label.append(r"$T_\mathrm{disk,}$" + rf"$_\mathrm{{{i+1}}}$") + + else: + break + if item == "disk_radius": quantity.append("disk_radius") - label.append(r"$R_\mathrm{disk}$") + if object_type == "planet": unit.append(r"$R_\mathrm{J}$") + elif object_type == "star": unit.append("au") + label.append(r"$R_\mathrm{disk}$") + + for i in range(100): + if item == f"disk_radius_{i}": + quantity.append(f"disk_radius_{i}") + + if object_type == "planet": + unit.append(r"$R_\mathrm{J}$") + + elif object_type == "star": + unit.append("au") + + label.append(r"$R_\mathrm{disk,}$" + rf"$_\mathrm{{{i+1}}}$") + + else: + break + if item == "flux_scaling": quantity.append("flux_scaling") unit.append(None)