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

Improvement: added a rmin for the 'o_fraction' plot.

parent c538e386
......@@ -2981,7 +2981,7 @@ def gas_fraction(genpath, folders, labels, is_ISM=True, var_str='mass'):
############################################################## Ions ##############################################################
###################################################################################################################################
def o_fraction(genpath, folders, timesteps, r200max=2, factor=2, err=False, savefig=False):
def o_fraction(genpath, folders, timesteps, r200min=0.3, 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
......@@ -2995,10 +2995,11 @@ def o_fraction(genpath, folders, timesteps, r200max=2, factor=2, err=False, save
len_fold = len(folders)
tot_ox = np.zeros(len_fold)
bottom, xaxis = np.zeros(len_fold), np.arange(len_fold)
# Compute
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):
......@@ -3008,7 +3009,7 @@ def o_fraction(genpath, folders, timesteps, r200max=2, factor=2, err=False, save
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:
......@@ -3016,11 +3017,11 @@ def o_fraction(genpath, folders, timesteps, r200max=2, factor=2, err=False, save
continue
except KeyError:
pass
# Else, compute it
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
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']
......@@ -3029,7 +3030,7 @@ def o_fraction(genpath, folders, timesteps, r200max=2, factor=2, err=False, save
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.
sel = (norm_r>r200min*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]
......
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