Skip to content

Commit

Permalink
Some improvements
Browse files Browse the repository at this point in the history
  • Loading branch information
@tomasstolker
tomasstolker committed Mar 14, 2021
1 parent c5524fc commit 1cf96c0
Show file tree
Hide file tree
Showing 4 changed files with 46 additions and 30 deletions.
4 changes: 3 additions & 1 deletion species/analysis/retrieval.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1313,7 +1313,9 @@ def loglike(cube,
elif quenching == 'diffusion':
# Calculate the quenching pressure from timescales
p_quench = retrieval_util.quench_pressure(
cube, cube_index, self.pressure, temp)
self.pressure, temp, cube[cube_index['metallicity']],
cube[cube_index['c_o_ratio']], cube[cube_index['logg']],
cube[cube_index['log_kzz']])

else:
p_quench = None
Expand Down
12 changes: 10 additions & 2 deletions species/data/database.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2008,12 +2008,16 @@ def add_retrieval(self,
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']
dset.attrs['chemistry'] = radtrans['chemistry']
dset.attrs['wavel_min'] = radtrans['wavel_range'][0]
dset.attrs['wavel_max'] = radtrans['wavel_range'][1]

if radtrans['quenching'] is None:
dset.attrs['quenching'] = 'None'
else:
dset.attrs['quenching'] = radtrans['quenching']

if 'pt_smooth' in radtrans:
dset.attrs['pt_smooth'] = radtrans['pt_smooth']

Expand Down Expand Up @@ -2249,7 +2253,11 @@ def get_retrieval_spectra(tag: str,
# Get chemistry attributes

chemistry = dset.attrs['chemistry']
quenching = dset.attrs['quenching']

if dset.attrs['quenching'] == 'None':
quenching = None
else:
quenching = dset.attrs['quenching']

# Get P-T profile attributes

Expand Down
6 changes: 3 additions & 3 deletions species/plot/plot_mcmc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -238,9 +238,9 @@ def plot_posterior(tag: str,
if 'H2O' in samples_box.parameters:
samples_box.parameters.append('c_h_ratio')
samples_box.parameters.append('o_h_ratio')
# samples_box.parameters.append('c_o_ratio')
samples_box.parameters.append('c_o_ratio')

ndim += 2
ndim += 3

abund_index = {}
for i, item in enumerate(samples_box.parameters):
Expand Down Expand Up @@ -403,7 +403,7 @@ def plot_posterior(tag: str,
print(f'Plotting the posterior: {output}...', end='', flush=True)

if 'H2O' in samples_box.parameters:
samples = np.column_stack((samples, c_h_ratio, o_h_ratio))
samples = np.column_stack((samples, c_h_ratio, o_h_ratio, c_o_ratio))

if inc_luminosity:
if 'teff' in samples_box.parameters and 'radius' in samples_box.parameters:
Expand Down
54 changes: 30 additions & 24 deletions species/util/retrieval_util.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -584,7 +584,7 @@ def create_abund_dict(abund_in: dict,
def calc_spectrum_clear(rt_object,
pressure: np.ndarray,
temperature: np.ndarray,
logg: float,
log_g: float,
c_o_ratio: Optional[float],
metallicity: Optional[float],
p_quench: Optional[float],
Expand All @@ -604,7 +604,7 @@ def calc_spectrum_clear(rt_object,
Array with the pressure points (bar).
temperature : np.ndarray
Array with the temperature points (K) corresponding to ``pressure``.
logg : float
log_g : float
Log10 of the surface gravity (cm s-2).
c_o_ratio : float, None
Carbon-to-oxygen ratio.
Expand Down Expand Up @@ -686,7 +686,7 @@ def calc_spectrum_clear(rt_object,
indices=None)

# calculate the emission spectrum
rt_object.calc_flux(temperature, abundances, 10.**logg, mmw, contribution=contribution)
rt_object.calc_flux(temperature, abundances, 10.**log_g, mmw, contribution=contribution)

# convert frequency (Hz) to wavelength (cm)
wavel = constants.LIGHT*1e2/rt_object.freq
Expand All @@ -711,8 +711,8 @@ def calc_spectrum_clouds(rt_object,
log_x_abund: Optional[dict],
log_x_base: dict,
fsed: float,
Kzz: float,
logg: float,
log_kzz: float,
log_g: float,
sigma_lnorm: float,
chemistry: str,
pressure_grid: str = 'smaller',
Expand Down Expand Up @@ -744,9 +744,9 @@ def calc_spectrum_clouds(rt_object,
Dictionary with the log10 of the mass fractions at the cloud base.
fsed : float
Sedimentation parameter.
Kzz : float
Log 10 of the eddy diffusion coefficient (cm2 s-1).
logg : float
log_kzz : float
Log10 of the eddy diffusion coefficient (cm2 s-1).
log_g : float
Log10 of the surface gravity (cm s-2).
sigma_lnorm : float
Geometric standard deviation of the log-normal size distribution.
Expand Down Expand Up @@ -843,7 +843,7 @@ def calc_spectrum_clouds(rt_object,
pressure_grid=pressure_grid,
indices=indices)

Kzz_use = np.full(pressure.shape, 10.**Kzz)
Kzz_use = np.full(pressure.shape, 10.**log_kzz)

if pressure_grid == 'smaller':
temperature = temperature[::3]
Expand Down Expand Up @@ -901,7 +901,7 @@ def calc_spectrum_clouds(rt_object,
plt.ylim(1e3, 1e-6)
plt.xlim(1e-10, 1.)
log_x_base_item = log_x_base[item]
plt.title(f'fsed = {fsed:.2f}, lgK = {Kzz:.2f}, X_b = {log_x_base_item:.2f}')
plt.title(f'fsed = {fsed:.2f}, log(Kzz) = {log_kzz:.2f}, X_b = {log_x_base_item:.2f}')
plt.savefig(f'{item.lower()}_clouds.pdf', bbox_inches='tight')
plt.clf()

Expand All @@ -916,7 +916,7 @@ def calc_spectrum_clouds(rt_object,
# Calculate the emission spectrum
check_scaling = rt_object.calc_flux(temperature,
abundances,
10.**logg,
10.**log_g,
mmw,
sigma_lnorm=sigma_lnorm,
Kzz=Kzz_use,
Expand Down Expand Up @@ -2012,24 +2012,30 @@ def convolve(input_wavel: np.ndarray,


@typechecked
def quench_pressure(cube,
cube_index: Dict[str, int],
pressure: np.ndarray,
temperature: np.ndarray) -> Optional[float]:
def quench_pressure(pressure: np.ndarray,
temperature: np.ndarray,
metallicity: float,
c_o_ratio: float,
log_g: float,
log_kzz: float) -> Optional[float]:
"""
Function to determine the CO/CH4 quenching pressure by intersecting the pressure-dependent
timescales of the vertical mixing and the CO/CH4 reaction rates.
Parameters
----------
cube : LP_c_double
Unit cube.
cube_index : dict
Dictionary with the index of each parameter in the ``cube``.
pressure : np.ndarray
Array with the pressures (bar).
temperature : np.ndarray
Array with the temperatures (K) corresponding to ``pressure``.
metallicity : float
Metallicity [Fe/H].
c_o_ratio : float
Carbon-to-oxygen ratio.
log_g : float
Log10 of the surface gravity (cm s-2).
log_kzz : float
Log10 of the eddy diffusion coefficient (cm2 s-1).
Returns
-------
Expand All @@ -2039,8 +2045,8 @@ def quench_pressure(cube,

# Interpolate the equilibbrium abundances

co_array = np.full(pressure.shape[0], cube[cube_index['c_o_ratio']])
feh_array = np.full(pressure.shape[0], cube[cube_index['metallicity']])
co_array = np.full(pressure.shape[0], c_o_ratio)
feh_array = np.full(pressure.shape[0], metallicity)

from poor_mans_nonequ_chem_FeH.poor_mans_nonequ_chem.poor_mans_nonequ_chem import \
interpol_abundances
Expand All @@ -2049,7 +2055,7 @@ def quench_pressure(cube,
co_array, feh_array, temperature, pressure, Pquench_carbon=None)

# Surface gravity (m s-2)
gravity = 1e-2 * 10.**cube[cube_index['logg']]
gravity = 1e-2 * 10.**log_g

# Mean molecular weight (kg)
mmw = abund_eq['MMW'] * constants.ATOMIC_MASS
Expand All @@ -2058,13 +2064,13 @@ def quench_pressure(cube,
h_scale = constants.BOLTZMANN * temperature / (mmw * gravity)

# Diffusion coefficient (m2 s-1)
kzz = 1e-4 * 10.**cube[cube_index['log_kzz']]
kzz = 1e-4 * 10.**log_kzz

# Mixing timescale (s)
t_mix = h_scale**2 / kzz

# Chemical timescale (see Eq. 12 from Zahnle & Marley 2014)
metal = 10.**cube[cube_index['metallicity']]
metal = 10.**metallicity
t_chem = 1.5e-6 * pressure**-1. * metal**-0.7 * np.exp(42000./temperature)

# Determine pressure at which t_mix = t_chem
Expand Down

0 comments on commit 1cf96c0

Please sign in to comment.