From 101dd1147e4ddb9d19736a459a50e2064cd67488 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Thu, 12 Mar 2020 16:47:50 +0100 Subject: [PATCH] Updated unit tests --- species/__init__.py | 6 +- species/analysis/retrieval.py | 758 --------------------------------- species/data/database.py | 539 +++++++++++------------ species/plot/plot_retrieval.py | 163 ------- species/read/read_model.py | 7 +- species/read/read_radtrans.py | 185 -------- species/util/phot_util.py | 30 +- species/util/retrieval_util.py | 685 ----------------------------- test/test_read/test_model.py | 16 +- test/test_read/test_planck.py | 6 +- 10 files changed, 307 insertions(+), 2088 deletions(-) delete mode 100644 species/analysis/retrieval.py delete mode 100644 species/plot/plot_retrieval.py delete mode 100644 species/read/read_radtrans.py delete mode 100644 species/util/retrieval_util.py diff --git a/species/__init__.py b/species/__init__.py index 70787bcb..30ca852b 100644 --- a/species/__init__.py +++ b/species/__init__.py @@ -6,7 +6,7 @@ from species.analysis.photometry import SyntheticPhotometry -from species.analysis.retrieval import AtmosphericRetrieval +# from species.analysis.retrieval import AtmosphericRetrieval from species.read.read_calibration import ReadCalibration @@ -18,7 +18,7 @@ from species.read.read_planck import ReadPlanck -from species.read.read_radtrans import ReadRadtrans +# from species.read.read_radtrans import ReadRadtrans from species.read.read_spectrum import ReadSpectrum @@ -44,7 +44,7 @@ plot_walkers, \ plot_photometry -from species.plot.plot_retrieval import plot_pt_profile +# from species.plot.plot_retrieval import plot_pt_profile from species.plot.plot_spectrum import plot_spectrum diff --git a/species/analysis/retrieval.py b/species/analysis/retrieval.py deleted file mode 100644 index 5b2c7e18..00000000 --- a/species/analysis/retrieval.py +++ /dev/null @@ -1,758 +0,0 @@ -""" -Module with functionalities for atmospheric retrieval with petitRADTRANS (Mollière et al. 2019). -More details on the retrieval code are available at https://petitradtrans.readthedocs.io. -""" - -import os -import json -import warnings - -import pymultinest -import numpy as np -import matplotlib.pyplot as plt - -from scipy.stats import invgamma -from rebin_give_width import rebin_give_width - -from petitRADTRANS import Radtrans -from petitRADTRANS_ck_test_speed import nat_cst as nc -from petitRADTRANS_ck_test_speed import Radtrans as RadtransScatter - -from species.analysis import photometry -from species.data import database -from species.core import constants -from species.read import read_object -from species.util import retrieval_util - - -os.environ['OMP_NUM_THREADS'] = '1' - - -class AtmosphericRetrieval: - """ - Class for atmospheric retrieval with petitRADTRANS. - """ - - def __init__(self, - object_name, - line_species, - cloud_species, - scattering, - output_name): - """ - Parameters - ---------- - object_name : str - Object name in the database. - line_species : list - List with the line species. - cloud_species : list - List with the cloud species. No clouds are used if an empty list is provided. - scattering : bool - Include scattering in the radiative transfer. - output_name : str - Output name that is used for the output files from MultiNest. - - Returns - ------- - NoneType - None - """ - - # input parameters - - self.object_name = object_name - self.line_species = line_species - self.cloud_species = cloud_species - self.scattering = scattering - self.output_name = output_name - - # get object data - - self.object = read_object.ReadObject(self.object_name) - self.distance = self.object.get_distance()[0] # [pc] - - print(f'Object: {self.object_name}') - print(f'Distance: {self.distance}') - - if len(self.line_species) == 0: - print(f'Line species: None') - - else: - print(f'Line species:') - for item in self.line_species: - print(f' - {item}') - - if len(self.cloud_species) == 0: - print(f'Cloud species: None') - - else: - print(f'Cloud species:') - for item in self.cloud_species: - print(f' - {item}') - - print(f'Scattering: {self.scattering}') - - species_db = database.Database() - - objectbox = species_db.get_object(object_name, - filters=None, - inc_phot=True, - inc_spec=True) - - # get photometric data - - self.objphot = [] - self.synphot = [] - - if len(objectbox.filters) != 0: - warnings.warn('Support for photometric data is not yet implemented.') - print('Photometric data:') - - for item in objectbox.filters: - obj_phot = self.object.get_photometry(item) - self.objphot.append((obj_phot[2], obj_phot[3])) - - print(f' - {item} (W m-2 um-1) = {obj_phot[2]:.2e} +/- {obj_phot[3]}') - - sphot = photometry.SyntheticPhotometry(item) - self.synphot.append(sphot) - - if not self.objphot: - self.objphot = None - - if not self.synphot: - self.synphot = None - - # get spectroscopic data - - self.spectrum = self.object.get_spectrum() - - if self.spectrum is None: - raise ValueError('A spectrum is required for atmospheric retrieval.') - - # set wavelength bins and add to spectrum dictionary - - self.wavel_min = [] - self.wavel_max = [] - - print('Spectroscopic data:') - - for key, value in self.spectrum.items(): - dict_val = list(value) - wavel_data = dict_val[0][:, 0] - - wavel_bins = np.zeros_like(wavel_data) - wavel_bins[:-1] = np.diff(wavel_data) - wavel_bins[-1] = wavel_bins[-2] - - dict_val.append(wavel_bins) - self.spectrum[key] = dict_val - - # min and max wavelength for Radtrans object - - self.wavel_min.append(wavel_data[0]) - self.wavel_max.append(wavel_data[-1]) - - print(f' - {key}') - print(f' Wavelength range (um) = {wavel_data[0]:.2f} - {wavel_data[-1]:.2f}') - print(f' Spectral resolution = {self.spectrum[key][3]:.2f}') - - # mock p-t profile for Radtrans object - - temp_params = {} - temp_params['log_delta'] = -6. - temp_params['log_gamma'] = 1. - temp_params['t_int'] = 750. - temp_params['t_equ'] = 0. - temp_params['log_p_trans'] = -3. - temp_params['alpha'] = 0. - - self.pressure, _ = nc.make_press_temp(temp_params) - - print(f'Initiating {self.pressure.size} pressure levels (bar): ' - f'{self.pressure[0]:.2e} - {self.pressure[-1]:.2e}') - - # initiate parameter list and counters - - self.parameters = [] - - def set_parameters(self, - bounds, - quenching, - pt_profile): - """ - Function to set the list with parameters. - - Parameters - ---------- - bounds : dict - Dictionary with the parameter boundaries. - quenching : bool - Fitting a quenching pressure. - pt_profile : str - The parametrization for the pressure-temperature profile ('molliere' or 'line'). - - Returns - ------- - NoneType - None - """ - - # generic parameters - - self.parameters.append('logg') - self.parameters.append('radius') - - # p-t profile parameters - - if pt_profile == 'molliere': - self.parameters.append('tint') - self.parameters.append('t1') - self.parameters.append('t2') - self.parameters.append('t3') - self.parameters.append('alpha') - self.parameters.append('log_delta') - - elif pt_profile == 'line': - for i in range(15): - self.parameters.append(f't{i}') - - self.parameters.append('gamma_r') - - # abundance parameters - - self.parameters.append('feh') - self.parameters.append('co') - - if quenching: - self.parameters.append('log_p_quench') - - # cloud parameters - - if len(self.cloud_species) > 0: - self.parameters.append('fe_fraction') - self.parameters.append('mgsio3_fraction') - self.parameters.append('fsed') - self.parameters.append('kzz') - self.parameters.append('sigma_lnorm') - - # add the flux scaling parameters - - for item in self.spectrum: - if item in bounds: - if bounds[item][0] is not None: - self.parameters.append(f'scaling_{item}') - - # add the error offset parameters - - for item in self.spectrum: - if item in bounds: - if bounds[item][1] is not None: - self.parameters.append(f'error_{item}') - - print(f'Fitting {len(self.parameters)} parameters:') - - for item in self.parameters: - print(f' - {item}') - - def run_multinest(self, - bounds, - quenching=True, - pt_profile='molliere', - live_points=2000, - efficiency=0.05, - resume=False, - plotting=False): - """ - Function to sample the posterior distribution with MultiNest. See also - https://github.com/farhanferoz/MultiNest. - - Parameters - ---------- - bounds : dict - Dictionary with the prior boundaries. - quenching : bool - Fitting a quenching pressure. - pt_profile : str - The parametrization for the pressure-temperature profile ('molliere' or 'line'). - live_points : int - Number of live points. - efficiency : float - Sampling efficiency. - resume : bool - Resume from a previous run. - plotting : bool - Plot sample results for testing. - - Returns - ------- - NoneType - None - """ - - # set initial number of parameters (not including the flux scaling and error offeset) - - if len(self.cloud_species) == 0: - n_param = 11 - else: - n_param = 16 - - # create list with parameters for MultiNest - - self.set_parameters(bounds, quenching, pt_profile) - - # create a dictionary with the cube indices of the parameters - - cube_index = {} - for i, item in enumerate(self.parameters): - cube_index[item] = i - - # delete the cloud parameters from the boundaries dictionary in case of no cloud species - - if len(self.cloud_species) == 0: - if 'fe_fraction' in bounds: - del bounds['fe_fraction'] - - if 'mgsio3_fraction' in bounds: - del bounds['mgsio3_fraction'] - - if 'fsed' in bounds: - del bounds['fsed'] - - if 'kzz' in bounds: - del bounds['kzz'] - - if 'sigma_lnorm' in bounds: - del bounds['sigma_lnorm'] - - # create Ratrans object - - print('Setting up petitRADTRANS...') - - if self.scattering: - rt_object = RadtransScatter(line_species=self.line_species, - rayleigh_species=['H2', 'He'], - cloud_species=self.cloud_species, - continuum_opacities=['H2-H2', 'H2-He'], - wlen_bords_micron=(0.95*min(self.wavel_min), - 1.05*max(self.wavel_max)), - mode='c-k', - test_ck_shuffle_comp=self.scattering, - do_scat_emis=self.scattering) - - else: - rt_object = Radtrans(line_species=self.line_species, - rayleigh_species=['H2', 'He'], - cloud_species=self.cloud_species, - continuum_opacities=['H2-H2', 'H2-He'], - wlen_bords_micron=(0.95*min(self.wavel_min), - 1.05*max(self.wavel_max)), - mode='c-k') - - # create RT arrays of appropriate lengths by using every three pressure points - - print(f'Decreasing the number of pressure levels: {self.pressure.size} -> ' - f'{self.pressure[::3].size}.') - - rt_object.setup_opa_structure(self.pressure[::3]) - - if pt_profile == 'line': - knot_press = np.logspace(np.log10(self.pressure[0]), np.log10(self.pressure[-1]), 15) - - def prior(cube, n_dim, n_param): - """ - Function to transform the unit cube into the parameter cube. - - Parameters - ---------- - cube : pymultinest.run.LP_c_double - Unit cube. - - Returns - ------- - NoneType - None - """ - - # surface gravity (dex) - if 'logg' in bounds: - logg = bounds['logg'][0] + (bounds['logg'][1]-bounds['logg'][0])*cube[cube_index['logg']] - else: - # default: 2 - 5.5 dex - logg = 2. + 3.5*cube[cube_index['logg']] - - cube[cube_index['logg']] = logg - - # planet radius (Rjup) - if 'radius' in bounds: - radius = bounds['radius'][0] + (bounds['radius'][1]-bounds['radius'][0])*cube[cube_index['radius']] - else: - # defaul: 0.8-2 Rjup - radius = 0.8 + 1.2*cube[cube_index['radius']] - - cube[cube_index['radius']] = radius - - if pt_profile == 'molliere': - - # internal temperature (K) of the Eddington model - # see Eq. 2 in Mollière et al. in prep. - if 'tint' in bounds: - tint = bounds['tint'][0] + (bounds['tint'][1]-bounds['tint'][0])*cube[cube_index['tint']] - else: - # default: 500 - 3000 K - tint = 500. + 2500.*cube[cube_index['tint']] - - cube[cube_index['tint']] = tint - - # connection temperature (K) - t_connect = (3./4.*tint**4.*(0.1+2./3.))**0.25 - - # the temperature (K) at temp_3 is scaled down from t_connect - temp_3 = t_connect*(1-cube[cube_index['t3']]) - cube[cube_index['t3']] = temp_3 - - # the temperature (K) at temp_2 is scaled down from temp_3 - temp_2 = temp_3*(1-cube[cube_index['t2']]) - cube[cube_index['t2']] = temp_2 - - # the temperature (K) at temp_1 is scaled down from temp_2 - temp_1 = temp_2*(1-cube[cube_index['t1']]) - cube[cube_index['t1']] = temp_1 - - # alpha: power law index in tau = delta * press_cgs**alpha - # see Eq. 1 in Mollière et al. in prep. - if 'alpha' in bounds: - alpha = bounds['alpha'][0] + (bounds['alpha'][1]-bounds['alpha'][0])*cube[cube_index['alpha']] - else: - # default: 1 - 2 - alpha = 1. + cube[cube_index['alpha']] - - cube[cube_index['alpha']] = alpha - - # photospheric pressure (bar) - # default: 1e-3 - 1e2 bar - p_phot = 10.**(-3. + 5.*cube[cube_index['log_delta']]) - - # delta: proportionality factor in tau = delta * press_cgs**alpha - # see Eq. 1 in Mollière et al. in prep. - delta = (p_phot*1e6)**(-alpha) - log_delta = np.log10(delta) - - cube[cube_index['log_delta']] = log_delta - - elif pt_profile == 'line': - # 15 temperature (K) knots - for i in range(15): - # default: 0 - 4000 K - cube[cube_index[f't{i}']] = 4000.*cube[cube_index[f't{i}']] - - # penalization of wiggles in the P-T profile - # inverse Gamma: a=1, b=5e-5 - gamma_r = invgamma.ppf(cube[cube_index['gamma_r']], a=1., scale=5e-5) - cube[cube_index['gamma_r']] = gamma_r - - # metallicity (dex) for the nabla_ad interpolation - if 'feh' in bounds: - feh = bounds['feh'][0] + (bounds['feh'][1]-bounds['feh'][0])*cube[cube_index['feh']] - else: - # default: -1.5 - 1.5 dex - feh = -1.5 + 3.*cube[cube_index['feh']] - - cube[cube_index['feh']] = feh - - # carbon-to-oxygen ratio for the nabla_ad interpolation - if 'co' in bounds: - co_ratio = bounds['co'][0] + (bounds['co'][1]-bounds['co'][0])*cube[cube_index['co']] - else: - # default: 0.1 - 1.6 - co_ratio = 0.1 + 1.5*cube[cube_index['co']] - - cube[cube_index['co']] = co_ratio - - # quench pressure (bar) - # default: 1e-6 - 1e3 bar - if quenching: - if 'log_p_quench' in bounds: - log_p_quench = bounds['log_p_quench'][0] + (bounds['log_p_quench'][1]-bounds['log_p_quench'][0])*cube[cube_index['log_p_quench']] - else: - # default: -6 - 3. - log_p_quench = -6. + 9.*cube[cube_index['log_p_quench']] - - cube[cube_index['log_p_quench']] = log_p_quench - - if len(self.cloud_species) > 0: - # cloud base mass fractions of Fe (iron) - # relative to the maximum values allowed from elemental abundances - # see Eq. 3 in Mollière et al. in prep. - # default: 0.05 - 1. - fe_fraction = np.log10(0.05)+(np.log10(1.)-np.log10(0.05))*cube[cube_index['fe_fraction']] - cube[cube_index['fe_fraction']] = fe_fraction - - # cloud base mass fractions of MgSiO3 (enstatite) - # relative to the maximum values allowed from elemental abundances - # see Eq. 3 in Mollière et al. in prep. - # default: 0.05 - 1. - mgsio3_fraction = np.log10(0.05)+(np.log10(1.)-np.log10(0.05))*cube[cube_index['mgsio3_fraction']] - cube[cube_index['mgsio3_fraction']] = mgsio3_fraction - - # sedimentation parameter - # ratio of the settling and mixing velocities of the cloud particles - # see Eq. 3 in Mollière et al. in prep. - if 'fsed' in bounds: - fsed = bounds['fsed'][0] + (bounds['fsed'][1]-bounds['fsed'][0])*cube[cube_index['fsed']] - else: - # default: 0 - 10 - fsed = 10.*cube[cube_index['fsed']] - - cube[cube_index['fsed']] = fsed - - # eddy diffusion coefficient, log(Kzz) - if 'kzz' in bounds: - kzz = bounds['kzz'][0] + (bounds['kzz'][1]-bounds['kzz'][0])*cube[cube_index['kzz']] - else: - # default: 5 - 13 - kzz = 5. + 8.*cube[cube_index['kzz']] - - cube[cube_index['kzz']] = kzz - - # width of the log-normal particle size distribution TODO (um?) - if 'sigma_lnorm' in bounds: - sigma_lnorm = bounds['sigma_lnorm'][0] + (bounds['sigma_lnorm'][1] - - bounds['sigma_lnorm'][0])*cube[cube_index['sigma_lnorm']] - else: - # default: 1.05 - 3. TODO um (?) - sigma_lnorm = 1.05 + 1.95*cube[cube_index['sigma_lnorm']] - - cube[cube_index['sigma_lnorm']] = sigma_lnorm - - # add flux scaling parameter if the boundaries are provided - - for item in self.spectrum: - if item in bounds: - if bounds[item][0] is not None: - cube[cube_index[f'scaling_{item}']] = bounds[item][0][0] + \ - (bounds[item][0][1]-bounds[item][0][0])*cube[cube_index[f'scaling_{item}']] - - # add error inflation parameter if the boundaries are provided - - for item in self.spectrum: - if item in bounds: - if bounds[item][1] is not None: - cube[cube_index[f'error_{item}']] = bounds[item][1][0] + \ - (bounds[item][1][1]-bounds[item][1][0]) * \ - cube[cube_index[f'error_{item}']] - - def loglike(cube, n_dim, n_param): - """ - Function for the logarithm of the likelihood, computed from the parameter cube. - - Parameters - ---------- - cube : pymultinest.run.LP_c_double - Unit cube. - - Returns - ------- - float - Sum of the logarithm of the prior and likelihood. - """ - - # initiate the logarithm of the prior and likelihood - - log_prior = 0. - log_likelihood = 0. - - # create dictionary with flux scaling parameters - - scaling = {} - - for item in self.spectrum: - if item in bounds and bounds[item][0] is not None: - scaling[item] = cube[cube_index[f'scaling_{item}']] - else: - scaling[item] = 1. - - # create dictionary with error offset parameters - - err_offset = {} - - for item in self.spectrum: - if item in bounds and bounds[item][1] is not None: - err_offset[item] = cube[cube_index[f'error_{item}']] - else: - err_offset[item] = -100. - - # create a p-t profile - - if pt_profile == 'molliere': - temp, _, _ = retrieval_util.pt_ret_model(np.array([cube[cube_index['t1']], - cube[cube_index['t2']], - cube[cube_index['t3']]]), - 10.**cube[cube_index['log_delta']], - cube[cube_index['alpha']], - cube[cube_index['tint']], - self.pressure, - cube[cube_index['feh']], - cube[cube_index['co']]) - - elif pt_profile == 'line': - knot_temp = [] - for i in range(15): - knot_temp.append(cube[cube_index[f't{i}']]) - - temp = retrieval_util.pt_spline_interp(knot_press, knot_temp, self.pressure) - - knot_temp = np.asarray(knot_temp) - - temp_sum = np.sum((knot_temp[2:] + knot_temp[:-2] - 2.*knot_temp[1:-1])**2.) - - # temp_sum = np.sum((temp[::3][2:] + temp[::3][:-2] - 2.*temp[::3][1:-1])**2.) - - log_prior += -1.*temp_sum/(2.*cube[cube_index['gamma_r']]) - \ - 0.5*np.log(2.*np.pi*cube[cube_index['gamma_r']]) - - # return zero probability if the minimum temperature is negative - - if np.min(temp) < 0.: - return -np.inf - - # set the quenching pressure - if quenching: - log_p_quench = cube[cube_index['log_p_quench']] - else: - log_p_quench = -10. - - # calculate the emission spectrum - - if len(self.cloud_species) > 0: - # cloudy atmosphere - - # mass fraction of Fe - x_fe = retrieval_util.return_XFe(cube[cube_index['feh']], cube[cube_index['co']]) - - # logarithm of the cloud base mass fraction of Fe - log_x_base_fe = np.log10(1e1**cube[cube_index['fe_fraction']]*x_fe) - - # mass fraction of MgSiO3 - x_mgsio3 = retrieval_util.return_XMgSiO3(cube[cube_index['feh']], cube[cube_index['co']]) - - # logarithm of the cloud base mass fraction of MgSiO3 - log_x_base_mgsio3 = np.log10(1e1**cube[cube_index['mgsio3_fraction']]*x_mgsio3) - - # wlen_micron, flux_lambda, Pphot_esti, tau_pow, tau_cloud = \ - wlen_micron, flux_lambda = retrieval_util.calc_spectrum_clouds( - rt_object, self.pressure, temp, cube[cube_index['co']], cube[cube_index['feh']], log_p_quench, - log_x_base_fe, log_x_base_mgsio3, cube[cube_index['fsed']], cube[cube_index['feh']], cube[cube_index['kzz']], - cube[cube_index['logg']], cube[cube_index['sigma_lnorm']], half=True, plotting=plotting) - - else: - # clear atmosphere - wlen_micron, flux_lambda = retrieval_util.calc_spectrum_clear( - rt_object, self.pressure, temp, cube[cube_index['logg']], - cube[cube_index['co']], cube[cube_index['feh']], log_p_quench, half=True) - - # return zero probability if the spectrum contains NaN values - - if np.sum(np.isnan(flux_lambda)) > 0: - if len(flux_lambda) > 1: - warnings.warn('Spectrum with NaN values encountered.') - - return -np.inf - - # scale the emitted spectrum to the observation - flux_lambda *= (cube[cube_index['radius']]*constants.R_JUP / (self.distance*constants.PARSEC))**2. - - for key, value in self.spectrum.items(): - # get spectrum - data_wavel = value[0][:, 0] - data_flux = value[0][:, 1] - data_error = value[0][:, 2] - - # get inverted covariance matrix - data_cov_inv = value[2] - - # get spectral resolution - spec_res = value[3] - - # get wavelength bins - data_wavel_bins = value[4] - - # fitted error component - err_fit = 10.**err_offset[key] - - # convolve with Gaussian LSF - flux_smooth = retrieval_util.convolve(wlen_micron, - flux_lambda, - spec_res) - - # resample to the observation - flux_rebinned = rebin_give_width(wlen_micron, - flux_smooth, - data_wavel, - data_wavel_bins) - - # difference between the observed and modeled spectrum - diff = flux_rebinned - scaling[key]*data_flux - - if data_cov_inv is not None: - # calculate the log-likelihood with the covariance matrix - # TODO include err_fit in the covariance matrix - log_likelihood += -np.dot(diff, data_cov_inv.dot(diff))/2. - - else: - # calculate the log-likelihood without the covariance matrix - var_infl = data_error**2.+err_fit**2 - log_likelihood += -0.5*np.sum(diff**2/var_infl + np.log(2.*np.pi*var_infl)) - - if plotting: - plt.errorbar(data_wavel, scaling[key]*data_flux, yerr=data_error+err_fit, - marker='o', ms=3, color='tab:blue', markerfacecolor='tab:blue') - - plt.plot(data_wavel, flux_rebinned, marker='o', ms=3, color='tab:orange') - - if plotting: - plt.plot(wlen_micron, flux_smooth, color='black', zorder=-20) - plt.xlabel(r'Wavelength ($\mu$m)') - plt.ylabel(r'Flux (W m$^{-2}$ $\mu$m$^{-1}$)') - plt.savefig('spectrum.pdf', bbox_inches='tight') - plt.clf() - - return log_prior + log_likelihood - - # store the model parameters in a JSON file - - print(f'Storing the model parameters: {self.output_name}_params.json') - - with open(f'{self.output_name}_params.json', 'w') as json_file: - json.dump(self.parameters, json_file) - - # store the Radtrans arguments in a JSON file - - print(f'Storing the Radtrans arguments: {self.output_name}_radtrans.json') - - radtrans_dict = {} - radtrans_dict['line_species'] = self.line_species - radtrans_dict['cloud_species'] = self.cloud_species - radtrans_dict['distance'] = self.distance - radtrans_dict['scattering'] = self.scattering - radtrans_dict['quenching'] = quenching - radtrans_dict['pt_profile'] = pt_profile - - with open(f'{self.output_name}_radtrans.json', 'w', encoding='utf-8') as json_file: - json.dump(radtrans_dict, json_file, ensure_ascii=False, indent=4) - - # run the nested sampling with MultiNest - - print('Sampling the posterior distribution with MultiNest...') - - pymultinest.run(loglike, - prior, - len(self.parameters), - outputfiles_basename=f'{self.output_name}_', - resume=resume, - verbose=True, - const_efficiency_mode=True, - sampling_efficiency=efficiency, - n_live_points=live_points, - evidence_tolerance=0.5) diff --git a/species/data/database.py b/species/data/database.py index b886b40a..bf6bae22 100644 --- a/species/data/database.py +++ b/species/data/database.py @@ -14,9 +14,9 @@ from astropy.io import fits -from petitRADTRANS import Radtrans -from petitRADTRANS_ck_test_speed import nat_cst as nc -from petitRADTRANS_ck_test_speed import Radtrans as RadtransScatter +# from petitRADTRANS import Radtrans +# from petitRADTRANS_ck_test_speed import nat_cst as nc +# from petitRADTRANS_ck_test_speed import Radtrans as RadtransScatter from species.analysis import photometry from species.core import box, constants @@ -24,7 +24,8 @@ companions, filters, btsettl, ames_dusty, ames_cond, \ isochrones, petitcode, exo_rem from species.read import read_model, read_calibration, read_planck -from species.util import data_util, retrieval_util +from species.util import data_util +# from species.util import data_util, retrieval_util class Database: @@ -1296,268 +1297,268 @@ def get_samples(self, prob_sample=prob_sample, median_sample=median_sample) - def add_retrieval(self, - tag, - output_name): - """ - Parameters - ---------- - tag : str - Database tag. - output_name : str - Output name that was used for the output files by MultiNest. - - Returns - ------- - NoneType - None - """ - - print('Storing samples in the database...', end='', flush=True) - - with open(f'{output_name}_params.json') as json_file: - parameters = json.load(json_file) - - with open(f'{output_name}_radtrans.json') as json_file: - radtrans = json.load(json_file) - - samples = np.loadtxt(f'{output_name}_post_equal_weights.dat') - - with h5py.File(self.database, 'a') as h5_file: - - if 'results' not in h5_file: - h5_file.create_group('results') - - if 'results/mcmc' not in h5_file: - h5_file.create_group('results/mcmc') - - if f'results/mcmc/{tag}' in h5_file: - del h5_file[f'results/mcmc/{tag}'] - - # remove the column with the log-likelihood value - samples = samples[:, :-1] - - if samples.shape[1] != len(parameters): - raise ValueError('The number of parameters is not equal to the parameter size ' - 'of the samples array.') - - dset = h5_file.create_dataset(f'results/mcmc/{tag}/samples', data=samples) - - dset.attrs['type'] = 'model' - dset.attrs['spectrum'] = 'petitradtrans' - dset.attrs['n_param'] = len(parameters) - dset.attrs['distance'] = radtrans['distance'] - - count_scale = 0 - count_error = 0 - - for i, item in enumerate(parameters): - dset.attrs[f'parameter{i}'] = item - - for i, item in enumerate(parameters): - if item[0:6] == 'scaling_': - dset.attrs[f'scaling{count_scale}'] = item - count_scale += 1 - - for i, item in enumerate(parameters): - if item[0:6] == 'error_': - dset.attrs[f'error{count_error}'] = item - count_error += 1 - - dset.attrs['n_scaling'] = count_scale - dset.attrs['n_error'] = count_error - - for i, item in enumerate(radtrans['line_species']): - dset.attrs[f'line_species{i}'] = item - - for i, item in enumerate(radtrans['cloud_species']): - dset.attrs[f'cloud_species{i}'] = item - - dset.attrs['n_line_species'] = len(radtrans['line_species']) - dset.attrs['n_cloud_species'] = len(radtrans['cloud_species']) - - dset.attrs['scattering'] = radtrans['scattering'] - dset.attrs['quenching'] = radtrans['quenching'] - dset.attrs['pt_profile'] = radtrans['pt_profile'] - - print(' [DONE]') - - def get_retrieval_spectra(self, - tag, - random, - wavel_range, - spec_res=None): - """ - Parameters - ---------- - tag : str - Database tag with the MCMC samples. - random : int - Number of randomly selected samples. - wavel_range : tuple(float, float) or str - Wavelength range (um) or filter name. - spec_res : float - Spectral resolution that is used for the smoothing with a Gaussian kernel. No smoothing - is applied if set to None. - - Returns - ------- - list(species.core.box.ModelBox, ) - Boxes with the randomly sampled spectra. - """ - - 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') - dset = h5_file[f'results/mcmc/{tag}/samples'] - - spectrum_type = dset.attrs['type'] - spectrum_name = dset.attrs['spectrum'] - - if 'n_param' in dset.attrs: - n_param = dset.attrs['n_param'] - elif 'nparam' in dset.attrs: - n_param = dset.attrs['nparam'] - - n_line_species = dset.attrs['n_line_species'] - n_cloud_species = dset.attrs['n_cloud_species'] - - scattering = dset.attrs['scattering'] - quenching = dset.attrs['quenching'] - pt_profile = dset.attrs['pt_profile'] - - if dset.attrs.__contains__('distance'): - distance = dset.attrs['distance'] - else: - distance = None - - samples = np.asarray(dset) - - random_indices = np.random.randint(samples.shape[0], size=random) - samples = samples[random_indices, :] - - parameters = [] - for i in range(n_param): - parameters.append(dset.attrs[f'parameter{i}']) - - parameters = np.asarray(parameters) - - line_species = [] - for i in range(n_line_species): - line_species.append(dset.attrs[f'line_species{i}']) - - line_species = np.asarray(line_species) - - cloud_species = [] - for i in range(n_cloud_species): - cloud_species.append(dset.attrs[f'cloud_species{i}']) - - cloud_species = np.asarray(cloud_species) - - # create mock p-t profile - - temp_params = {} - temp_params['log_delta'] = -6. - temp_params['log_gamma'] = 1. - temp_params['t_int'] = 750. - temp_params['t_equ'] = 0. - temp_params['log_p_trans'] = -3. - temp_params['alpha'] = 0. - - pressure, _ = nc.make_press_temp(temp_params) - - logg_index = np.argwhere(parameters == 'logg')[0] - radius_index = np.argwhere(parameters == 'radius')[0] - feh_index = np.argwhere(parameters == 'feh')[0] - co_index = np.argwhere(parameters == 'co')[0] - - if quenching: - log_p_quench_index = np.argwhere(parameters == 'log_p_quench')[0] - - if pt_profile == 'molliere': - tint_index = np.argwhere(parameters == 'tint')[0] - t1_index = np.argwhere(parameters == 't1')[0] - t2_index = np.argwhere(parameters == 't2')[0] - t3_index = np.argwhere(parameters == 't3')[0] - alpha_index = np.argwhere(parameters == 'alpha')[0] - log_delta_index = np.argwhere(parameters == 'log_delta')[0] - - elif pt_profile == 'line': - temp_index = [] - for i in range(15): - temp_index.append(np.argwhere(parameters == f't{i}')[0]) - - knot_press = np.logspace(np.log10(pressure[0]), np.log10(pressure[-1]), 15) - - if scattering: - rt_object = RadtransScatter(line_species=line_species, - rayleigh_species=['H2', 'He'], - cloud_species=cloud_species, - continuum_opacities=['H2-H2', 'H2-He'], - wlen_bords_micron=wavel_range, - mode='c-k', - test_ck_shuffle_comp=scattering, - do_scat_emis=scattering) - - else: - rt_object = Radtrans(line_species=line_species, - rayleigh_species=['H2', 'He'], - cloud_species=cloud_species, - continuum_opacities=['H2-H2', 'H2-He'], - wlen_bords_micron=wavel_range, - mode='c-k') - - # create RT arrays of appropriate lengths by using every three pressure points - rt_object.setup_opa_structure(pressure[::3]) - - boxes = [] - - for i, item in tqdm.tqdm(enumerate(samples), desc='Getting MCMC spectra'): - - if pt_profile == 'molliere': - temp, _, _ = retrieval_util.pt_ret_model( - np.array([item[t1_index][0], item[t2_index][0], item[t3_index][0]]), - 10.**item[log_delta_index][0], item[alpha_index][0], item[tint_index][0], pressure, - item[feh_index][0], item[co_index][0]) - - elif pt_profile == 'line': - knot_temp = [] - for i in range(15): - knot_temp.append(item[temp_index[i]][0]) - - temp = retrieval_util.pt_spline_interp(knot_press, knot_temp, pressure) - - if quenching: - log_p_quench = item[log_p_quench_index][0] - else: - log_p_quench = -10. - - wavelength, flux = retrieval_util.calc_spectrum_clear( - rt_object, pressure, temp, item[logg_index][0], item[co_index][0], - item[feh_index][0], log_p_quench, half=True) - - flux *= (item[radius_index]*constants.R_JUP/(distance*constants.PARSEC))**2. - - if spec_res is not None: - # convolve with a Gaussian line spread function - flux = retrieval_util.convolve(wavelength, flux, spec_res) - - model_box = box.create_box(boxtype='model', - model='petitradtrans', - wavelength=wavelength, - flux=flux, - parameters=None, - quantity='flux') - - model_box.type = 'mcmc' - - boxes.append(model_box) - - h5_file.close() - - return boxes + # def add_retrieval(self, + # tag, + # output_name): + # """ + # Parameters + # ---------- + # tag : str + # Database tag. + # output_name : str + # Output name that was used for the output files by MultiNest. + # + # Returns + # ------- + # NoneType + # None + # """ + # + # print('Storing samples in the database...', end='', flush=True) + # + # with open(f'{output_name}_params.json') as json_file: + # parameters = json.load(json_file) + # + # with open(f'{output_name}_radtrans.json') as json_file: + # radtrans = json.load(json_file) + # + # samples = np.loadtxt(f'{output_name}_post_equal_weights.dat') + # + # with h5py.File(self.database, 'a') as h5_file: + # + # if 'results' not in h5_file: + # h5_file.create_group('results') + # + # if 'results/mcmc' not in h5_file: + # h5_file.create_group('results/mcmc') + # + # if f'results/mcmc/{tag}' in h5_file: + # del h5_file[f'results/mcmc/{tag}'] + # + # # remove the column with the log-likelihood value + # samples = samples[:, :-1] + # + # if samples.shape[1] != len(parameters): + # raise ValueError('The number of parameters is not equal to the parameter size ' + # 'of the samples array.') + # + # dset = h5_file.create_dataset(f'results/mcmc/{tag}/samples', data=samples) + # + # dset.attrs['type'] = 'model' + # dset.attrs['spectrum'] = 'petitradtrans' + # dset.attrs['n_param'] = len(parameters) + # dset.attrs['distance'] = radtrans['distance'] + # + # count_scale = 0 + # count_error = 0 + # + # for i, item in enumerate(parameters): + # dset.attrs[f'parameter{i}'] = item + # + # for i, item in enumerate(parameters): + # if item[0:6] == 'scaling_': + # dset.attrs[f'scaling{count_scale}'] = item + # count_scale += 1 + # + # for i, item in enumerate(parameters): + # if item[0:6] == 'error_': + # dset.attrs[f'error{count_error}'] = item + # count_error += 1 + # + # dset.attrs['n_scaling'] = count_scale + # dset.attrs['n_error'] = count_error + # + # for i, item in enumerate(radtrans['line_species']): + # dset.attrs[f'line_species{i}'] = item + # + # for i, item in enumerate(radtrans['cloud_species']): + # dset.attrs[f'cloud_species{i}'] = item + # + # dset.attrs['n_line_species'] = len(radtrans['line_species']) + # dset.attrs['n_cloud_species'] = len(radtrans['cloud_species']) + # + # dset.attrs['scattering'] = radtrans['scattering'] + # dset.attrs['quenching'] = radtrans['quenching'] + # dset.attrs['pt_profile'] = radtrans['pt_profile'] + # + # print(' [DONE]') + # + # def get_retrieval_spectra(self, + # tag, + # random, + # wavel_range, + # spec_res=None): + # """ + # Parameters + # ---------- + # tag : str + # Database tag with the MCMC samples. + # random : int + # Number of randomly selected samples. + # wavel_range : tuple(float, float) or str + # Wavelength range (um) or filter name. + # spec_res : float + # Spectral resolution that is used for the smoothing with a Gaussian kernel. No smoothing + # is applied if set to None. + # + # Returns + # ------- + # list(species.core.box.ModelBox, ) + # Boxes with the randomly sampled spectra. + # """ + # + # 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') + # dset = h5_file[f'results/mcmc/{tag}/samples'] + # + # spectrum_type = dset.attrs['type'] + # spectrum_name = dset.attrs['spectrum'] + # + # if 'n_param' in dset.attrs: + # n_param = dset.attrs['n_param'] + # elif 'nparam' in dset.attrs: + # n_param = dset.attrs['nparam'] + # + # n_line_species = dset.attrs['n_line_species'] + # n_cloud_species = dset.attrs['n_cloud_species'] + # + # scattering = dset.attrs['scattering'] + # quenching = dset.attrs['quenching'] + # pt_profile = dset.attrs['pt_profile'] + # + # if dset.attrs.__contains__('distance'): + # distance = dset.attrs['distance'] + # else: + # distance = None + # + # samples = np.asarray(dset) + # + # random_indices = np.random.randint(samples.shape[0], size=random) + # samples = samples[random_indices, :] + # + # parameters = [] + # for i in range(n_param): + # parameters.append(dset.attrs[f'parameter{i}']) + # + # parameters = np.asarray(parameters) + # + # line_species = [] + # for i in range(n_line_species): + # line_species.append(dset.attrs[f'line_species{i}']) + # + # line_species = np.asarray(line_species) + # + # cloud_species = [] + # for i in range(n_cloud_species): + # cloud_species.append(dset.attrs[f'cloud_species{i}']) + # + # cloud_species = np.asarray(cloud_species) + # + # # create mock p-t profile + # + # temp_params = {} + # temp_params['log_delta'] = -6. + # temp_params['log_gamma'] = 1. + # temp_params['t_int'] = 750. + # temp_params['t_equ'] = 0. + # temp_params['log_p_trans'] = -3. + # temp_params['alpha'] = 0. + # + # pressure, _ = nc.make_press_temp(temp_params) + # + # logg_index = np.argwhere(parameters == 'logg')[0] + # radius_index = np.argwhere(parameters == 'radius')[0] + # feh_index = np.argwhere(parameters == 'feh')[0] + # co_index = np.argwhere(parameters == 'co')[0] + # + # if quenching: + # log_p_quench_index = np.argwhere(parameters == 'log_p_quench')[0] + # + # if pt_profile == 'molliere': + # tint_index = np.argwhere(parameters == 'tint')[0] + # t1_index = np.argwhere(parameters == 't1')[0] + # t2_index = np.argwhere(parameters == 't2')[0] + # t3_index = np.argwhere(parameters == 't3')[0] + # alpha_index = np.argwhere(parameters == 'alpha')[0] + # log_delta_index = np.argwhere(parameters == 'log_delta')[0] + # + # elif pt_profile == 'line': + # temp_index = [] + # for i in range(15): + # temp_index.append(np.argwhere(parameters == f't{i}')[0]) + # + # knot_press = np.logspace(np.log10(pressure[0]), np.log10(pressure[-1]), 15) + # + # if scattering: + # rt_object = RadtransScatter(line_species=line_species, + # rayleigh_species=['H2', 'He'], + # cloud_species=cloud_species, + # continuum_opacities=['H2-H2', 'H2-He'], + # wlen_bords_micron=wavel_range, + # mode='c-k', + # test_ck_shuffle_comp=scattering, + # do_scat_emis=scattering) + # + # else: + # rt_object = Radtrans(line_species=line_species, + # rayleigh_species=['H2', 'He'], + # cloud_species=cloud_species, + # continuum_opacities=['H2-H2', 'H2-He'], + # wlen_bords_micron=wavel_range, + # mode='c-k') + # + # # create RT arrays of appropriate lengths by using every three pressure points + # rt_object.setup_opa_structure(pressure[::3]) + # + # boxes = [] + # + # for i, item in tqdm.tqdm(enumerate(samples), desc='Getting MCMC spectra'): + # + # if pt_profile == 'molliere': + # temp, _, _ = retrieval_util.pt_ret_model( + # np.array([item[t1_index][0], item[t2_index][0], item[t3_index][0]]), + # 10.**item[log_delta_index][0], item[alpha_index][0], item[tint_index][0], pressure, + # item[feh_index][0], item[co_index][0]) + # + # elif pt_profile == 'line': + # knot_temp = [] + # for i in range(15): + # knot_temp.append(item[temp_index[i]][0]) + # + # temp = retrieval_util.pt_spline_interp(knot_press, knot_temp, pressure) + # + # if quenching: + # log_p_quench = item[log_p_quench_index][0] + # else: + # log_p_quench = -10. + # + # wavelength, flux = retrieval_util.calc_spectrum_clear( + # rt_object, pressure, temp, item[logg_index][0], item[co_index][0], + # item[feh_index][0], log_p_quench, half=True) + # + # flux *= (item[radius_index]*constants.R_JUP/(distance*constants.PARSEC))**2. + # + # if spec_res is not None: + # # convolve with a Gaussian line spread function + # flux = retrieval_util.convolve(wavelength, flux, spec_res) + # + # model_box = box.create_box(boxtype='model', + # model='petitradtrans', + # wavelength=wavelength, + # flux=flux, + # parameters=None, + # quantity='flux') + # + # model_box.type = 'mcmc' + # + # boxes.append(model_box) + # + # h5_file.close() + # + # return boxes diff --git a/species/plot/plot_retrieval.py b/species/plot/plot_retrieval.py deleted file mode 100644 index cb0ef096..00000000 --- a/species/plot/plot_retrieval.py +++ /dev/null @@ -1,163 +0,0 @@ -import numpy as np -import matplotlib as mpl -import matplotlib.pyplot as plt - -from petitRADTRANS_ck_test_speed import nat_cst as nc - -from species.data import database -from species.util import retrieval_util - - -def plot_pt_profile(tag, - random=100, - xlim=None, - ylim=None, - offset=None, - output='pt_profile.pdf'): - """ - Function to plot the posterior distribution. - - Parameters - ---------- - tag : str - Database tag with the MCMC samples. - random : int - Number of randomly selected samples from the posterior. - xlim : tuple(float, float) - Limits of the wavelength axis. - ylim : tuple(float, float) - Limits of the flux axis. - offset : tuple(float, float), None - Offset of the x- and y-axis label. - output : str - Output filename. - - Returns - ------- - NoneType - None - """ - - print(f'Plotting the P-T profiles: {output}...', end='', flush=True) - - species_db = database.Database() - box = species_db.get_samples(tag, burnin=0) - - parameters = np.asarray(box.parameters) - samples = box.samples - median = box.median_sample - - indices = np.random.randint(samples.shape[0], size=random) - samples = samples[indices, ] - - mpl.rcParams['font.serif'] = ['Bitstream Vera Serif'] - mpl.rcParams['font.family'] = 'serif' - - plt.rc('axes', edgecolor='black', linewidth=2.5) - - plt.figure(1, figsize=(4., 5.)) - gridsp = mpl.gridspec.GridSpec(1, 1) - gridsp.update(wspace=0, hspace=0, left=0, right=1, bottom=0, top=1) - - ax = plt.subplot(gridsp[0, 0]) - - ax.tick_params(axis='both', which='major', colors='black', labelcolor='black', - direction='in', width=1, length=5, labelsize=12, top=True, - bottom=True, left=True, right=True) - - ax.tick_params(axis='both', which='minor', colors='black', labelcolor='black', - direction='in', width=1, length=3, labelsize=12, top=True, - bottom=True, left=True, right=True) - - ax.set_xlabel('Temperature (K)', fontsize=13) - ax.set_ylabel('Pressure (bar)', fontsize=13) - - if xlim: - ax.set_xlim(xlim[0], xlim[1]) - else: - ax.set_xlim(1000., 5000.) - - if ylim: - ax.set_ylim(ylim[0], ylim[1]) - else: - ax.set_ylim(1e3, 1e-6) - - ax.set_yscale('log') - - if offset is not None: - ax.get_xaxis().set_label_coords(0.5, offset[0]) - ax.get_yaxis().set_label_coords(offset[1], 0.5) - - else: - ax.get_xaxis().set_label_coords(0.5, -0.06) - ax.get_yaxis().set_label_coords(-0.14, 0.5) - - # create pressure levels - - temp_params = {} - temp_params['log_delta'] = -6. - temp_params['log_gamma'] = 1. - temp_params['t_int'] = 750. - temp_params['t_equ'] = 0. - temp_params['log_p_trans'] = -3. - temp_params['alpha'] = 0. - - pressure, _ = nc.make_press_temp(temp_params) - - if 'tint' in parameters: - pt_profile = 'molliere' - - tint_index = np.argwhere(parameters == 'tint')[0] - t1_index = np.argwhere(parameters == 't1')[0] - t2_index = np.argwhere(parameters == 't2')[0] - t3_index = np.argwhere(parameters == 't3')[0] - alpha_index = np.argwhere(parameters == 'alpha')[0] - log_delta_index = np.argwhere(parameters == 'log_delta')[0] - - else: - pt_profile = 'line' - - temp_index = [] - for i in range(15): - temp_index.append(np.argwhere(parameters == f't{i}')[0]) - - knot_press = np.logspace(np.log10(pressure[0]), np.log10(pressure[-1]), 15) - - feh_index = np.argwhere(parameters == 'feh')[0] - co_index = np.argwhere(parameters == 'co')[0] - - for item in samples: - if pt_profile == 'molliere': - temp, _, _ = retrieval_util.pt_ret_model( - np.array([item[t1_index][0], item[t2_index][0], item[t3_index][0]]), - 10.**item[log_delta_index][0], item[alpha_index][0], item[tint_index][0], pressure, - item[feh_index][0], item[co_index][0]) - - elif pt_profile == 'line': - knot_temp = [] - for i in range(15): - knot_temp.append(item[temp_index[i]][0]) - - temp = retrieval_util.pt_spline_interp(knot_press, knot_temp, pressure) - - ax.plot(temp, pressure, '-', lw=0.3, color='gray', alpha=0.5, zorder=1) - - if pt_profile == 'molliere': - temp, _, _ = retrieval_util.pt_ret_model( - np.array([median['t1'], median['t2'], median['t3']]), 10.**median['log_delta'], - median['alpha'], median['tint'], pressure, median['feh'], median['co']) - - elif pt_profile == 'line': - knot_temp = [] - for i in range(15): - knot_temp.append(median[f't{i}']) - - temp = retrieval_util.pt_spline_interp(knot_press, knot_temp, pressure) - - ax.plot(temp, pressure, '-', lw=1, color='black', zorder=2) - - plt.savefig(output, bbox_inches='tight') - plt.clf() - plt.close() - - print(' [DONE]') diff --git a/species/read/read_model.py b/species/read/read_model.py index 30fafd62..0337b527 100644 --- a/species/read/read_model.py +++ b/species/read/read_model.py @@ -470,7 +470,12 @@ def get_model(self, readcalib = read_calibration.ReadCalibration('vega', filter_name=None) calibbox = readcalib.get_spectrum() - flux_vega, _ = spectres.spectres(new_spec_wavs=self.wl_points, + if wavel_resample is not None: + new_spec_wavs = wavel_resample + else: + new_spec_wavs = self.wl_points + + flux_vega, _ = spectres.spectres(new_spec_wavs=new_spec_wavs, old_spec_wavs=calibbox.wavelength, spec_fluxes=calibbox.flux, spec_errs=calibbox.error) diff --git a/species/read/read_radtrans.py b/species/read/read_radtrans.py deleted file mode 100644 index c6df001e..00000000 --- a/species/read/read_radtrans.py +++ /dev/null @@ -1,185 +0,0 @@ -""" -Module with reading functionalities for atmospheric models from petitRADTRANS. See -Mollière et al. 2019 for details about the retrieval code. -""" - -import os -import configparser - -import numpy as np - -from petitRADTRANS import Radtrans -from petitRADTRANS_ck_test_speed import nat_cst as nc -from petitRADTRANS_ck_test_speed import Radtrans as RadtransScatter - -from species.core import box, constants -from species.read import read_filter -from species.util import read_util, retrieval_util - - -class ReadRadtrans: - """ - Class for reading a model spectrum from the database. - """ - - def __init__(self, - line_species=['H2O', 'CO', 'CH4'], - cloud_species=[], - scattering=False, - wavel_range=None, - filter_name=None): - """ - Parameters - ---------- - line_species : list - List with the line species. - cloud_species : list - List with the cloud species. No clouds are used if an empty list is provided. - scattering : bool - Include scattering in the radiative transfer. - wavel_range : tuple(float, float), None - Wavelength range (um). The wavelength range is set to 0.8-10 um if set to None or - 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 - ''filter_name`` is set to None. - - Returns - ------- - NoneType - None - """ - - self.filter_name = filter_name - self.wavel_range = wavel_range - self.scattering = scattering - - if self.filter_name is not None: - transmission = read_filter.ReadFilter(self.filter_name) - self.wavel_range = transmission.wavelength_range() - - elif self.wavel_range is None: - self.wavel_range = (0.8, 10.) - - config_file = os.path.join(os.getcwd(), 'species_config.ini') - - config = configparser.ConfigParser() - config.read_file(open(config_file)) - - self.database = config['species']['database'] - - # create mock p-t profile - - temp_params = {} - temp_params['log_delta'] = -6. - temp_params['log_gamma'] = 1. - temp_params['t_int'] = 750. - temp_params['t_equ'] = 0. - temp_params['log_p_trans'] = -3. - temp_params['alpha'] = 0. - - self.pressure, _ = nc.make_press_temp(temp_params) - - # create Radtrans object - - if self.scattering: - self.rt_object = RadtransScatter(line_species=line_species, - rayleigh_species=['H2', 'He'], - cloud_species=cloud_species, - continuum_opacities=['H2-H2', 'H2-He'], - wlen_bords_micron=wavel_range, - mode='c-k', - test_ck_shuffle_comp=self.scattering, - do_scat_emis=self.scattering) - - else: - self.rt_object = Radtrans(line_species=line_species, - rayleigh_species=['H2', 'He'], - cloud_species=cloud_species, - continuum_opacities=['H2-H2', 'H2-He'], - wlen_bords_micron=wavel_range, - mode='c-k') - - # create RT arrays of appropriate lengths by using every three pressure points - self.rt_object.setup_opa_structure(self.pressure[::3]) - - def get_model(self, - model_param, - spec_res=None, - wavel_resample=None): - """ - Function for extracting a model spectrum by linearly interpolating the model grid. The - parameters values should lie within the boundaries of the grid points that are stored - in the database. The stored grid points can be inspected with - :func:`~species.read.read_model.ReadModel.get_points`. - - Parameters - ---------- - model_param : dict - Model parameters and values. The values should be within the boundaries of the grid. - The grid boundaries of the available spectra in the database can be obtained with - :func:`~species.read.read_model.ReadModel.get_bounds()`. - spec_res : float, None - Spectral resolution, achieved by smoothing with a Gaussian kernel. The original - wavelength points are used if both ``spec_res`` and ``wavel_resample`` are set to None. - wavel_resample : numpy.ndarray - Wavelength points (um) to which the spectrum is resampled. Only used if - ``spec_res`` is set to None. - - Returns - ------- - species.core.box.ModelBox - Box with the model spectrum. - """ - - if spec_res is not None and wavel_resample is not None: - raise ValueError('The \'spec_res\' and \'wavel_resample\' parameters can not be used ' - 'simultaneously. Please set one of them to None.') - - if 'tint' in model_param: - temp, _, _ = retrieval_util.pt_ret_model( - np.array([model_param['t1'], model_param['t2'], model_param['t3']]), - 10.**model_param['log_delta'], model_param['alpha'], model_param['tint'], - self.pressure, model_param['feh'], model_param['co']) - - else: - knot_press = np.logspace(np.log10(self.pressure[0]), np.log10(self.pressure[-1]), 15) - - knot_temp = [] - for i in range(15): - knot_temp.append(model_param[f't{i}']) - - temp = retrieval_util.pt_spline_interp(knot_press, knot_temp, self.pressure) - - if 'log_p_quench' in model_param: - log_p_quench = model_param['log_p_quench'] - else: - log_p_quench = -10. - - if self.scattering: - pass - - else: - wavelength, flux = retrieval_util.calc_spectrum_clear( - self.rt_object, self.pressure, temp, model_param['logg'], model_param['co'], - model_param['feh'], log_p_quench, half=True) - - if 'radius' in model_param: - model_param['mass'] = read_util.get_mass(model_param) - - if 'distance' in model_param: - scaling = (model_param['radius']*constants.R_JUP)**2 / \ - (model_param['distance']*constants.PARSEC)**2 - - flux *= scaling - - if spec_res is not None: - # convolve with Gaussian LSF - flux = retrieval_util.convolve(wavelength, flux, spec_res) - - return box.create_box(boxtype='model', - model='petitradtrans', - wavelength=wavelength, - flux=flux, - parameters=model_param, - quantity='flux') diff --git a/species/util/phot_util.py b/species/util/phot_util.py index 65a94bd3..90cd9101 100644 --- a/species/util/phot_util.py +++ b/species/util/phot_util.py @@ -11,7 +11,8 @@ import numpy as np from species.core import box -from species.read import read_model, read_calibration, read_filter, read_planck, read_radtrans +from species.read import read_model, read_calibration, read_filter, read_planck +# from species.read import read_model, read_calibration, read_filter, read_planck, read_radtrans def multi_photometry(datatype, @@ -206,18 +207,21 @@ def get_residuals(datatype, else: if spectrum == 'petitradtrans': - radtrans = read_radtrans.ReadRadtrans(line_species=kwargs_radtrans['line_species'], - cloud_species=kwargs_radtrans['cloud_species'], - scattering=kwargs_radtrans['scattering'], - wavel_range=wavel_range) - - model = radtrans.get_model(parameters, spec_res=None) - - # separate resampling to the new wavelength points - - flux_new = spectres.spectres(new_spec_wavs=wl_new, - old_spec_wavs=model.wavelength, - spec_fluxes=model.flux) + # TODO change back + pass + + # radtrans = read_radtrans.ReadRadtrans(line_species=kwargs_radtrans['line_species'], + # cloud_species=kwargs_radtrans['cloud_species'], + # scattering=kwargs_radtrans['scattering'], + # wavel_range=wavel_range) + # + # model = radtrans.get_model(parameters, spec_res=None) + # + # # separate resampling to the new wavelength points + # + # flux_new = spectres.spectres(new_spec_wavs=wl_new, + # old_spec_wavs=model.wavelength, + # spec_fluxes=model.flux) else: readmodel = read_model.ReadModel(spectrum, wavel_range=wavel_range) diff --git a/species/util/retrieval_util.py b/species/util/retrieval_util.py deleted file mode 100644 index 97c80630..00000000 --- a/species/util/retrieval_util.py +++ /dev/null @@ -1,685 +0,0 @@ -import copy - -import numpy as np -import matplotlib.pyplot as plt - -from scipy.interpolate import interp1d, CubicSpline -from scipy.ndimage.filters import gaussian_filter - -from petitRADTRANS_ck_test_speed import nat_cst as nc -from poor_mans_nonequ_chem_FeH.poor_mans_nonequ_chem.poor_mans_nonequ_chem import interpol_abundances - - -def pt_ret_model(T3, delta, alpha, tint, press, FeH, CO, conv=True): - """ - Self-luminous retrieval P-T model. - - It has 7 free parameters: - - T3 = np.array([t1, t2, t3]): temperature points to be added on top - of the radiative Eddington structure (above tau = 0.1). - Use spline interpolation, t1 < t2 < t3 < tconnect as prior. - - delta: proportionality factor in tau = delta * press_cgs**alpha - - alpha: power law index in tau = delta * press_cgs**alpha - For the tau model: use proximity to kappa_rosseland photosphere - as prior. - - tint: internal temperature of the Eddington model - - press: input pressure profile in bar - - FeH: metallicity for the nabla_ad interpolation - - CO: C/O for the nabla_ad interpolation - - conv: enforce convective adiabat yes/no - """ - - # Go grom bar to cgs - press_cgs = press*1e6 - - # Calculate the optical depth - tau = delta*press_cgs**alpha - - # This is the eddington temperature - tedd = (3./4.*tint**4.*(2./3.+tau))**0.25 - - ab = interpol_abundances(CO*np.ones_like(tedd), FeH*np.ones_like(tedd), tedd, press) - - nabla_ad = ab['nabla_ad'] - - # Enforce convective adiabat - if conv: - # Calculate the current, radiative temperature gradient - nab_rad = np.diff(np.log(tedd))/np.diff(np.log(press_cgs)) - - # Extend to array of same length as pressure structure - nabla_rad = np.ones_like(tedd) - nabla_rad[0] = nab_rad[0] - nabla_rad[-1] = nab_rad[-1] - nabla_rad[1:-1] = (nab_rad[1:]+nab_rad[:-1])/2. - - # Where is the atmosphere convectively unstable? - conv_index = nabla_rad > nabla_ad - - for i in range(10): - - if i == 0: - t_take = copy.copy(tedd) - else: - t_take = copy.copy(tfinal) - - ab = interpol_abundances(CO*np.ones_like(t_take), FeH*np.ones_like(t_take), t_take, press) - - nabla_ad = ab['nabla_ad'] - - # Calculate the average nabla_ad between the layers - nabla_ad_mean = nabla_ad - nabla_ad_mean[1:] = (nabla_ad[1:]+nabla_ad[:-1])/2. - - # What are the increments in temperature due to convection - tnew = nabla_ad_mean[conv_index]*np.mean(np.diff(np.log(press_cgs))) - - # What is the last radiative temperature? - tstart = np.log(t_take[~conv_index][-1]) - - # Integrate and translate to temperature from log(temperature) - tnew = np.exp(np.cumsum(tnew)+tstart) - - # Add upper radiative and - # lower conective part into one single array - tfinal = copy.copy(t_take) - tfinal[conv_index] = tnew - - if np.max(np.abs(t_take-tfinal)/t_take) < 0.01: - # print('n_ad', 1./(1.-nabla_ad[conv_index])) - break - - else: - tfinal = tedd - - # Add the three temperature-point P-T description above tau = 0.1 - def press_tau(tau): - # Returns the pressure at a given tau, in cgs - return (tau/delta)**(1./alpha) - - # Where is the uppermost pressure of the Eddington radiative structure? - p_bot_spline = press_tau(0.1) - - for i_intp in range(2): - - if i_intp == 0: - - # Create the pressure coordinates for the spline support nodes at low pressure - support_points_low = np.logspace(np.log10(press_cgs[0]), np.log10(p_bot_spline), 4) - - # Create the pressure coordinates for the spline support nodes at high pressure, - # the corresponding temperatures for these nodes will be taken from the radiative+convective solution - support_points_high = 1e1**np.arange(np.log10(p_bot_spline), np.log10(press_cgs[-1]), np.diff(np.log10(support_points_low))[0]) - - # Combine into one support node array, don't add the p_bot_spline point twice. - support_points = np.zeros(len(support_points_low)+len(support_points_high)-1) - support_points[:4] = support_points_low - support_points[4:] = support_points_high[1:] - - else: - - # Create the pressure coordinates for the spline support nodes at low pressure - support_points_low = np.logspace(np.log10(press_cgs[0]), np.log10(p_bot_spline), 7) - - # Create the pressure coordinates for the spline support nodes at high pressure, - # the corresponding temperatures for these nodes will be taken from the radiative+convective solution - support_points_high = np.logspace(np.log10(p_bot_spline), np.log10(press_cgs[-1]), 7) - - # Combine into one support node array, don't add the p_bot_spline point twice. - support_points = np.zeros(len(support_points_low)+len(support_points_high)-1) - support_points[:7] = support_points_low - support_points[7:] = support_points_high[1:] - - # Define the temperature values at the node points. - t_support = np.zeros_like(support_points) - - if i_intp == 0: - tfintp = interp1d(press_cgs, tfinal) - - # The temperature at p_bot_spline (from the radiative-convectice solution) - t_support[int(len(support_points_low))-1] = tfintp(p_bot_spline) - - # The temperature at pressures below p_bot_spline (free parameters) - t_support[:(int(len(support_points_low))-1)] = T3 - # t_support[:3] = tfintp(support_points_low) - - # The temperature at pressures above p_bot_spline (from the radiative-convectice solution) - t_support[int(len(support_points_low)):] = tfintp(support_points[(int(len(support_points_low))):]) - - else: - tfintp1 = interp1d(press_cgs, tret) - - t_support[:(int(len(support_points_low))-1)] = tfintp1(support_points[:(int(len(support_points_low))-1)]) - - tfintp = interp1d(press_cgs, tfinal) - - # The temperature at p_bot_spline (from the radiative-convectice solution) - t_support[int(len(support_points_low))-1] = tfintp(p_bot_spline) - - # print('diff', t_connect_calc - tfintp(p_bot_spline)) - t_support[int(len(support_points_low)):] = tfintp(support_points[(int(len(support_points_low))):]) - - # Make the temperature spline interpolation to be returned to the user - # tret = spline(np.log10(support_points), t_support, np.log10(press_cgs), order = 3) - - cs = CubicSpline(np.log10(support_points), t_support) - tret = cs(np.log10(press_cgs)) - - # Return the temperature, the pressure at tau = 1, and the temperature at the connection point. - # The last two are needed for the priors on the P-T profile. - return tret, press_tau(1.)/1e6, tfintp(p_bot_spline) - - -def pt_spline_interp(knot_press, - knot_temp, - pressure): - - pt_interp = CubicSpline(np.log10(knot_press), knot_temp) - - return pt_interp(np.log10(pressure)) - - -def calc_spectrum_clear(rt_object, - press, - temp, - logg, - co, - feh, - log_p_quench, - half=False): - - # create arrays for constant values of C/O and Fe/H - - co_list = np.full(press.shape, co) - feh_list = np.full(press.shape, feh) - - # interpolate the abundances, following chemical equilibrium - - abundances_interp = interpol_abundances(co_list, feh_list, temp, press, Pquench_carbon=10.**log_p_quench) - - # extract the mean molecular weight - - mmw = abundances_interp['MMW'] - - # extract every three levels if half=True - - if half: - temp = temp[::3] - press = press[::3] - mmw = mmw[::3] - - # create a dictionary with the abundances by replacing species ending with _all_iso - - abundances = {} - - if half: - for species in rt_object.line_species: - abundances[species] = abundances_interp[species.replace('_all_iso', '')][::3] - - abundances['H2'] = abundances_interp['H2'][::3] - abundances['He'] = abundances_interp['He'][::3] - - else: - for species in rt_object.line_species: - abundances[species] = abundances_interp[species.replace('_all_iso', '')] - - abundances['H2'] = abundances_interp['H2'] - abundances['He'] = abundances_interp['He'] - - # Corretion for the nuclear spin degeneracy that was not included in the partition function - # See Charnay et al. (2018) - - if 'FeH' in abundances: - abundances['FeH'] = abundances['FeH']/2. - - # calculate the emission spectrum - - rt_object.calc_flux(temp, abundances, 10.**logg, mmw) - - # convert frequency (Hz) to wavelength (cm) - - wlen = nc.c/rt_object.freq - - # return wavelength (micron) and flux (W m-2 um-1) - - return 1e4*wlen, 1e-7*rt_object.flux*nc.c/wlen**2. - - -def calc_spectrum_clouds(rt_object, - press, - temp, - CO, - FeH, - log_p_quench, - log_X_cloud_base_Fe, - log_X_cloud_base_MgSiO3, - fsed_Fe, - fsed_MgSiO3, - Kzz, - logg, - sigma_lnorm, - half=False, - plotting=False): - - COs = CO * np.ones_like(press) - FeHs = FeH * np.ones_like(press) - - abundances_interp = interpol_abundances(COs, FeHs, temp, press, Pquench_carbon=1e1**log_p_quench) - - MMW = abundances_interp['MMW'] - - P_base_Fe = simple_cdf_Fe(press, temp, FeH, CO, np.mean(MMW), plotting) - P_base_MgSiO3 = simple_cdf_MgSiO3(press, temp, FeH, CO, np.mean(MMW), plotting=plotting) - - abundances = {} - - abundances['Fe(c)'] = np.zeros_like(temp) - - abundances['Fe(c)'][press < P_base_Fe] = \ - 1e1**log_X_cloud_base_Fe * (press[press <= P_base_Fe]/P_base_Fe)**fsed_Fe - - abundances['MgSiO3(c)'] = np.zeros_like(temp) - - abundances['MgSiO3(c)'][press < P_base_MgSiO3] = \ - 1e1**log_X_cloud_base_MgSiO3 * (press[press <= P_base_MgSiO3]/P_base_MgSiO3)**fsed_MgSiO3 - - if half: - abundances['Fe(c)'] = abundances['Fe(c)'][::3] - abundances['MgSiO3(c)'] = abundances['MgSiO3(c)'][::3] - - if half: - for species in rt_object.line_species: - abundances[species] = abundances_interp[species.replace('_all_iso', '')][::3] - - abundances['H2'] = abundances_interp['H2'][::3] - abundances['He'] = abundances_interp['He'][::3] - - else: - for species in rt_object.line_species: - abundances[species] = abundances_interp[species.replace('_all_iso', '')] - - abundances['H2'] = abundances_interp['H2'] - abundances['He'] = abundances_interp['He'] - - # Corretion for the nuclear spin degeneracy that was not included in the partition function - # See Charnay et al. (2018) - - if 'FeH' in abundances: - abundances['FeH'] = abundances['FeH']/2. - - Kzz_use = (1e1**Kzz) * np.ones_like(press) - - if half: - temp = temp[::3] - press = press[::3] - MMW = MMW[::3] - Kzz_use = Kzz_use[::3] - - fseds = {} - fseds['Fe(c)'] = fsed_Fe - fseds['MgSiO3(c)'] = fsed_MgSiO3 - - if plotting: - plt.plot(abundances['CO_all_iso'], press, label='CO') - plt.plot(abundances['CH4'], press, label='CH4') - plt.plot(abundances['H2O'], press, label='H2O') - plt.xlim([1e-10, 1.]) - plt.ylim([press[-1], press[0]]) - plt.yscale('log') - plt.xscale('log') - plt.xlabel('Mass fraction') - plt.ylabel('Pressure (bar)') - plt.axhline(1e1**log_p_quench) - plt.legend(loc='best') - plt.savefig('abundances.pdf', bbox_inches='tight') - plt.clf() - - plt.plot(temp, press) - plt.axhline(P_base_Fe, label='Cloud deck Fe') - plt.axhline(P_base_MgSiO3, label='Cloud deck MgSiO3') - plt.yscale('log') - plt.ylim([1e3, 1e-6]) - plt.xlim([0., 4000.]) - plt.savefig('pt_cloud_deck.pdf', bbox_inches='tight') - plt.clf() - - plt.plot(abundances['Fe(c)'], press) - plt.axhline(P_base_Fe) - plt.yscale('log') - if np.count_nonzero(abundances['Fe(c)']) > 0: - plt.xscale('log') - plt.ylim([1e3, 1e-6]) - plt.xlim([1e-10, 1.]) - plt.title('fsed_Fe = '+str(fsed_Fe)+' lgK='+str(Kzz)+' X_b = '+str(log_X_cloud_base_Fe)) - plt.savefig('fe_clouds.pdf', bbox_inches='tight') - plt.clf() - - plt.plot(abundances['MgSiO3(c)'], press) - plt.axhline(P_base_MgSiO3) - plt.yscale('log') - if np.count_nonzero(abundances['MgSiO3(c)']) > 0: - plt.xscale('log') - plt.ylim([1e3, 1e-6]) - plt.xlim([1e-10, 1.]) - plt.title('fsed_MgSiO3 = '+str(fsed_MgSiO3)+' lgK='+str(Kzz)+' X_b = '+str(log_X_cloud_base_MgSiO3)) - plt.savefig('mgsio3_clouds.pdf', bbox_inches='tight') - plt.clf() - - # Turn off clouds - # abundances['MgSiO3(c)'] = np.zeros_like(press) - # abundances['Fe(c)'] = np.zeros_like(press) - - rt_object.calc_flux(temp, abundances, 1e1**logg, MMW, Kzz=Kzz_use, fsed=fseds, sigma_lnorm=sigma_lnorm) - - wlen_micron = nc.c/rt_object.freq/1e-4 - wlen = nc.c/rt_object.freq - flux = rt_object.flux - - # convert flux f_nu to f_lambda - f_lambda = flux*nc.c/wlen**2. - # convert to flux per m^2 (from flux per cm^2) cancels with step below - # f_lambda = f_lambda * 1e4 - # convert to flux per micron (from flux per cm) cancels with step above - # f_lambda = f_lambda * 1e-4 - # convert from ergs to Joule - f_lambda = f_lambda * 1e-7 - - # plt.yscale('log') - # plt.xscale('log') - # plt.ylim([1e2,1e-6]) - # plt.ylabel('P (bar)') - # plt.xlabel('Average particle size of MgSiO3 particles (microns)') - # plt.plot(rt_object.r_g[:,rt_object.cloud_species.index('MgSiO3(c)')]/1e-4, press) - # plt.savefig('mgsio3_size.png') - # plt.show() - # plt.clf() - - # plt.yscale('log') - # plt.xscale('log') - # plt.ylim([1e2,1e-6]) - # plt.ylabel('P (bar)') - # plt.xlabel('Average particle size of Fe particles (microns)') - # plt.plot(rt_object.r_g[:,rt_object.cloud_species.index('Fe(c)')]/1e-4, press) - # plt.savefig('fe_size.png') - # plt.show()rt_object - # plt.clf() - - # return wlen_micron, f_lambda, rt_object.pphot, rt_object.tau_pow, np.mean(rt_object.tau_cloud) - return wlen_micron, f_lambda - - -############################################################# -# To calculate X_Fe from [Fe/H], C/O -############################################################# - -# metal species -metals = ['C', 'N', 'O', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'K', 'Ca', 'Ti', 'V', 'Fe', 'Ni'] - -# solar abundances, [Fe/H] = 0, from Asplund+ 2009 -nfracs = {} -nfracs['H'] = 0.9207539305 -nfracs['He'] = 0.0783688694 -nfracs['C'] = 0.0002478241 -nfracs['N'] = 6.22506056949881e-05 -nfracs['O'] = 0.0004509658 -nfracs['Na'] = 1.60008694353205e-06 -nfracs['Mg'] = 3.66558742055362e-05 -nfracs['Al'] = 2.595e-06 -nfracs['Si'] = 2.9795e-05 -nfracs['P'] = 2.36670201997668e-07 -nfracs['S'] = 1.2137900734604e-05 -nfracs['Cl'] = 2.91167958499589e-07 -nfracs['K'] = 9.86605611925677e-08 -nfracs['Ca'] = 2.01439011429255e-06 -nfracs['Ti'] = 8.20622804366359e-08 -nfracs['V'] = 7.83688694089992e-09 -nfracs['Fe'] = 2.91167958499589e-05 -nfracs['Ni'] = 1.52807116806281e-06 - -# atomic masses -masses = {} -masses['H'] = 1. -masses['He'] = 4. -masses['C'] = 12. -masses['N'] = 14. -masses['O'] = 16. -masses['Na'] = 23. -masses['Mg'] = 24.3 -masses['Al'] = 27. -masses['Si'] = 28. -masses['P'] = 31. -masses['S'] = 32. -masses['Cl'] = 35.45 -masses['K'] = 39.1 -masses['Ca'] = 40. -masses['Ti'] = 47.9 -masses['V'] = 51. -masses['Fe'] = 55.8 -masses['Ni'] = 58.7 - - -def return_XFe(FeH, CO): - - nfracs_use = copy.copy(nfracs) - - for spec in nfracs.keys(): - - if (spec != 'H') and (spec != 'He'): - nfracs_use[spec] = nfracs[spec]*1e1**FeH - - nfracs_use['O'] = nfracs_use['C']/CO - - XFe = masses['Fe']*nfracs_use['Fe'] - - add = 0. - for spec in nfracs_use.keys(): - add += masses[spec]*nfracs_use[spec] - - XFe = XFe / add - - return XFe - - -def return_XMgSiO3(FeH, CO): - - nfracs_use = copy.copy(nfracs) - - for spec in nfracs.keys(): - - if (spec != 'H') and (spec != 'He'): - nfracs_use[spec] = nfracs[spec]*1e1**FeH - - nfracs_use['O'] = nfracs_use['C']/CO - - nfracs_mgsio3 = np.min([nfracs_use['Mg'], nfracs_use['Si'], nfracs_use['O']/3.]) - - masses_mgsio3 = masses['Mg'] + masses['Si'] + 3. * masses['O'] - - Xmgsio3 = masses_mgsio3*nfracs_mgsio3 - - add = 0. - for spec in nfracs_use.keys(): - add += masses[spec]*nfracs_use[spec] - - Xmgsio3 = Xmgsio3 / add - - return Xmgsio3 - - -############################################################# -# Fe saturation pressure, from Ackerman & Marley (2001), including erratum (P_vap is in bar, not cgs!) -############################################################# - -def return_T_cond_Fe(FeH, CO, MMW = 2.33): - - T = np.linspace(100.,10000.,1000) - P_vap = lambda x: np.exp(15.71 - 47664./x) - - XFe = return_XFe(FeH, CO) - - return P_vap(T)/(XFe*MMW/masses['Fe']), T - - -def return_T_cond_Fe_l(FeH, CO, MMW = 2.33): - - T = np.linspace(100.,10000.,1000) - P_vap = lambda x: np.exp(9.86 - 37120./x) - - XFe = return_XFe(FeH, CO) - - return P_vap(T)/(XFe*MMW/masses['Fe']), T - - -def return_T_cond_Fe_comb(FeH, CO, MMW = 2.33): - - P1, T1 = return_T_cond_Fe(FeH, CO, MMW) - P2, T2 = return_T_cond_Fe_l(FeH, CO, MMW) - - retP = np.zeros_like(P1) - - index = P1 1e-8) & (Pc < 1e5) - Pc, Tc = Pc[index], Tc[index] - tcond_p = interp1d(Pc, Tc) - #print(Pc, press) - Tcond_on_input_grid = tcond_p(press) - - Tdiff = Tcond_on_input_grid - temp - diff_vec = Tdiff[1:]*Tdiff[:-1] - ind_cdf = (diff_vec < 0.) - if len(diff_vec[ind_cdf]) > 0: - P_clouds = (press[1:]+press[:-1])[ind_cdf]/2. - P_cloud = P_clouds[-1] - else: - P_cloud = 1e-8 - - if plotting: - plt.plot(temp, press) - plt.plot(Tcond_on_input_grid, press) - plt.axhline(P_cloud, color = 'red', linestyle = '--') - plt.yscale('log') - plt.xlim([0., 3000.]) - plt.ylim([1e2,1e-6]) - plt.savefig('condensation_fe.pdf', bbox_inches='tight') - plt.clf() - - return P_cloud - - -def simple_cdf_MgSiO3(press, temp, FeH, CO, MMW=2.33, plotting=False): - - Pc, Tc = return_T_cond_MgSiO3(FeH, CO, MMW) - index = (Pc > 1e-8) & (Pc < 1e5) - Pc, Tc = Pc[index], Tc[index] - tcond_p = interp1d(Pc, Tc) - # print(Pc, press) - Tcond_on_input_grid = tcond_p(press) - - Tdiff = Tcond_on_input_grid - temp - diff_vec = Tdiff[1:]*Tdiff[:-1] - ind_cdf = (diff_vec < 0.) - if len(diff_vec[ind_cdf]) > 0: - P_clouds = (press[1:]+press[:-1])[ind_cdf]/2. - P_cloud = P_clouds[-1] - else: - P_cloud = 1e-8 - - if plotting: - plt.plot(temp, press) - plt.plot(Tcond_on_input_grid, press) - plt.axhline(P_cloud, color='red', linestyle='--') - plt.yscale('log') - plt.xlim([0., 3000.]) - plt.ylim([1e2, 1e-6]) - plt.savefig('condensation_mgsio3.pdf', bbox_inches='tight') - plt.clf() - - return P_cloud - -# if plotting: -# kappa_IR = 0.01 -# gamma = 0.4 -# T_int = 200. -# T_equ = 1550. -# gravity = 1e1**2.45 -# -# pressures = np.logspace(-6, 2, 100) -# -# temperature = nc.guillot_global(pressures, kappa_IR, gamma, gravity, T_int, T_equ) -# -# simple_cdf_Fe(pressures, temperature, 0., 0.55) -# simple_cdf_MgSiO3(pressures, temperature, 0., 0.55) - -def convolve(input_wavelength, input_flux, instrument_res): - # From talking to Ignas: delta lambda of resolution element - # is FWHM of the LSF's standard deviation, hence: - sigma_lsf = 1./instrument_res/(2.*np.sqrt(2.*np.log(2.))) - - # The input spacing of petitRADTRANS is 1e3, but just compute - # it to be sure, or more versatile in the future. - # Also, we have a log-spaced grid, so the spacing is constant - # as a function of wavelength - spacing = np.mean(2.*np.diff(input_wavelength)/(input_wavelength[1:]+input_wavelength[:-1])) - - # Calculate the sigma to be used in the gauss filter in units - # of input wavelength bins - sigma_lsf_gauss_filter = sigma_lsf/spacing - - return gaussian_filter(input_flux, sigma=sigma_lsf_gauss_filter, mode='nearest') diff --git a/test/test_read/test_model.py b/test/test_read/test_model.py index e8c58c9a..be53cf03 100644 --- a/test/test_read/test_model.py +++ b/test/test_read/test_model.py @@ -43,36 +43,36 @@ def test_get_model(self): magnitude=False, smooth=True) - assert np.sum(model_box.wavelength) == pytest.approx( 45.969308588797304, rel=self.limit, abs=0.) - assert np.sum(model_box.flux) == pytest.approx(8.523746906754028e-13, rel=self.limit, abs=0.) + assert np.sum(model_box.wavelength) == pytest.approx(39.36850660634572, rel=self.limit, abs=0.) + assert np.sum(model_box.flux) == pytest.approx(7.549674442955264e-13, rel=self.limit, abs=0.) model_box = read_model.get_model(self.model_param, spec_res=100., magnitude=True, smooth=True) - assert np.sum(model_box.wavelength) == pytest.approx(45.969308588797304, rel=self.limit, abs=0.) - assert np.sum(model_box.flux) == pytest.approx(324.13066887417983, rel=self.limit, abs=0.) + assert np.sum(model_box.wavelength) == pytest.approx(39.36850660634572, rel=self.limit, abs=0.) + assert np.sum(model_box.flux) == pytest.approx(275.6061491339208, rel=self.limit, abs=0.) def test_get_data(self): read_model = species.ReadModel('ames-cond', filter_name='Paranal/NACO.H') model_box = read_model.get_data(self.model_param) assert np.sum(model_box.wavelength) == pytest.approx(47.859770458276856, rel=self.limit, abs=0.) - assert np.sum(model_box.flux) == pytest.approx(8.65388169807849e-13, rel=self.limit, abs=0.) + assert np.sum(model_box.flux) == pytest.approx(8.275363683007199e-13, rel=self.limit, abs=0.) def test_get_flux(self): read_model = species.ReadModel('ames-cond', filter_name='Paranal/NACO.H') flux = read_model.get_flux(self.model_param) - assert flux[0] == pytest.approx(3.4895271009175577e-14, rel=self.limit, abs=0.) + assert flux[0] == pytest.approx(3.3368963026400554e-14, rel=self.limit, abs=0.) def test_get_magnitude(self): read_model = species.ReadModel('ames-cond', filter_name='Paranal/NACO.H') magnitude = read_model.get_magnitude(self.model_param) - assert magnitude[0] == pytest.approx(11.30856476397497, rel=self.limit, abs=0.) - assert magnitude[1] == pytest.approx(11.30856476397497, rel=self.limit, abs=0.) + assert magnitude[0] == pytest.approx(11.357124426046317, rel=self.limit, abs=0.) + assert magnitude[1] == pytest.approx(11.357124426046317, rel=self.limit, abs=0.) def test_get_bounds(self): read_model = species.ReadModel('ames-cond', filter_name='Paranal/NACO.H') diff --git a/test/test_read/test_planck.py b/test/test_read/test_planck.py index c2dee7b6..c58a37a1 100644 --- a/test/test_read/test_planck.py +++ b/test/test_read/test_planck.py @@ -39,14 +39,14 @@ def test_get_spectrum(self): assert modelbox.flux.shape == (22, ) assert np.sum(modelbox.wavelength) == pytest.approx(27.672469420634147, rel=self.limit, abs=0.) - assert np.sum(modelbox.flux) == pytest.approx(4.568408954870407e-13, rel=self.limit, abs=0.) + assert np.sum(modelbox.flux) == pytest.approx(4.368588209687883e-13, rel=self.limit, abs=0.) def test_get_flux(self): read_planck = species.ReadPlanck(filter_name='MKO/NSFCam.J') flux = read_planck.get_flux({'teff': 2000., 'radius': 1., 'distance': 10.}) - assert flux[0] == pytest.approx(2.0806115243951512e-14, rel=self.limit, abs=0.) + assert flux[0] == pytest.approx(1.9896062423928137e-14, rel=self.limit, abs=0.) synphot = species.SyntheticPhotometry(filter_name='MKO/NSFCam.J') flux = read_planck.get_flux({'teff': 2000., 'radius': 1., 'distance': 10.}, synphot=synphot) - assert flux[0] == pytest.approx(2.0806115243951512e-14, rel=self.limit, abs=0.) + assert flux[0] == pytest.approx(1.9896062423928137e-14, rel=self.limit, abs=0.)