Commit c538e386 authored by Maxime Rey's avatar Maxime Rey
Browse files

New function 'o_fraction' that plots the ionisation fraction of O in each state.

parent 4c8e66db
......@@ -33,6 +33,9 @@ Utils to compute and plot several useful quantities:
- Things against time
- Ms_Mh_relation, Ms_Mh_t, gas_fraction
- Ions
- o_fraction
- Nicolas
- azimuth_frac
......@@ -2974,6 +2977,102 @@ def gas_fraction(genpath, folders, labels, is_ISM=True, var_str='mass'):
ax.set_xlabel('Time [Gyr]')
ax.set_yscale('log')
###################################################################################################################################
############################################################## Ions ##############################################################
###################################################################################################################################
def o_fraction(genpath, folders, timesteps, r200max=2, factor=2, err=False, savefig=False):
"""Computes the ionisation fraction of several ionisation states for oxygen."""
from ratatouille import readNsave as ras
import matplotlib.pyplot as plt
import numpy as np
import h5py, os
str_ions=['OI','OII','OIII','OIV','OV','OVI','OVII','OVIII','OIX'] #, # Input all ions you have computed +1.
if isinstance(timesteps,float) or isinstance(timesteps,int):
timesteps = [timesteps]
len_fold = len(folders)
tot_ox = np.zeros(len_fold)
bottom, xaxis = np.zeros(len_fold), np.arange(len_fold)
oxygen = np.zeros((len(str_ions),len_fold,len(timesteps))) # [ion, folder, timesstep]
for indfold, fold in enumerate(folders):
RamsesDir = genpath+fold
# Create file to store data.
filepath = f'{RamsesDir}/ratadat/'
if not os.path.exists(filepath):
os.mkdir(filepath)
filename = f'{filepath}fracO_within_xxrvir.h5'
if os.path.exists(filename):
file = h5py.File(filename,'r+')
else:
file = h5py.File(filename,'w')
# Computation
for istep, timestep in enumerate(timesteps):
basename = f'fraco_{timestep:05d}_{r200max}'
try:
oxygen[:,indfold,istep] = file[basename][:]
continue
except KeyError:
pass
print(f'\t Computing... \t {fold}: timestep {timestep}')
# Load cell_dx, total oxygen mass and cell_pos
_, cells, cell_pos, _ = ras.extract_cells(RamsesDir,timestep,factor=factor,saveinfile=False,verbose=False)
cell_dx, _, nH, _, _, _, _, _, _, _, _, _, _, Z = cells
tot_ox_cell = 4.90e-4/0.0134*Z*nH*cell_dx**3 # nO,sun = 4.90e-4/Z_sun * Z*nH = 4.90e-4/0.0134 * nH,sun
# Compute norm to select CGM cells and not have OI overtaking all other fractions
info = ras.extract_info(timestep, RamsesDir, saveinfile=True)
cu2cm = info['unit_l'] * info['boxlen']
center, r200 = ras.get_rad(RamsesDir, timestep, 1)
center_cm, r200_cm = [cen*cu2cm for cen in center], r200*cu2cm
vec_r = np.array([cell_pos[:,0]-center_cm[0], cell_pos[:,1]-center_cm[1], cell_pos[:,2]-center_cm[2]]).transpose()
norm_r = np.linalg.norm(vec_r,2, axis=1)
# Compute ion fractions
sel = (norm_r>0.1*r200_cm) & (norm_r<r200max*r200_cm) # Restrict within 0.1 R200 and r200max*r200.
tot_ox[indfold] = np.sum(tot_ox_cell[sel])
for indion, str_ion in enumerate(str_ions[:-1]):
oxy = ras.extract_ions(RamsesDir, timestep, str_ion, factor)*cell_dx**3 # [Nb of atoms]
oxygen[indion,indfold,istep] = np.sum(oxy[sel])/tot_ox[indfold]
oxygen[-1,indfold,:] = 1 - np.sum(oxygen[:-1,indfold,:], axis=0) # Compute remaining ion(s) fraction.
file.create_dataset(basename,data=oxygen[:,indfold,istep])
# Plot
perc_min, perc_max = 15.9, 84.1
oxy_stack = np.median(oxygen,axis=2)
if err:
yerr = np.array([oxy_stack-np.percentile(oxygen,perc_min, axis=2), np.percentile(oxygen,perc_max, axis=2)-oxy_stack]) # Added or substracted value from main value.
fig, ax = plt.subplots(figsize=(7,7))
for indion2, frac_ox in enumerate(oxy_stack):
if err=='OVI':
if indion2==5:
rects = ax.bar(xaxis, frac_ox, yerr=yerr[:, indion2, :], bottom=bottom, label=str_ions[indion2], error_kw=dict(lw=1, capsize=5, capthick=1))
else:
rects = ax.bar(xaxis, frac_ox, yerr=None, bottom=bottom, label=str_ions[indion2], error_kw=dict(lw=2, capsize=5, capthick=1.5))
elif err:
rects = ax.bar(xaxis, frac_ox, yerr=yerr[:, indion2, :], bottom=bottom, label=str_ions[indion2], error_kw=dict(lw=2, capsize=5, capthick=1.5))
else:
rects = ax.bar(xaxis, frac_ox, bottom=bottom, label=str_ions[indion2])
bottom += np.array(frac_ox)
# Add percentage on bars.
for indfold, rect in enumerate(rects):
height = np.round(float(rect.get_height()*100), 1)
if (height>5) or (indion2==5): # Write fraction for OVI and those higher than 0.05.
xloc = xaxis[indfold]
yloc = rect.get_y() + rect.get_height() / 2
if yloc < 0.127: # Avoid text being behind legend.
yloc=0.1277
if (err=='OVI') and (indion2==5): # If errorbar for OVI, move the text slightly up.
yloc+=0.02
ax.annotate(str(height)+' %',xy=(xloc,yloc),xytext=(xloc,yloc),textcoords="offset points",ha='center',va='center',color='k',clip_on=True)
plt.ylabel('Ion mass fraction')
plt.xticks(ticks=xaxis, labels=[' '.join(fold.split('_')[1:])[:-1] for fold in folders])
plt.legend(loc="lower center", ncol=int(len(str_ions)/2+1), fontsize='small')
plt.tight_layout()
plt.show()
savefig and fig.savefig(f'{str_ions[0][:-1]}_fraction.svg')
###################################################################################################################################
############################################################## M*-Mh ##############################################################
###################################################################################################################################
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment