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

New function: 'rad_Zfrac' for radial metallicity maps.

parent d8e591fd
......@@ -11,7 +11,7 @@ Utils to compute and plot several useful quantities:
- Radial plots
- half_mass_radius,
- rad_prof, rad_col,
- rad_prof, rad_Zfrac, rad_col,
- check_folder.
- Stars-related
......@@ -1201,6 +1201,87 @@ def rad_prof(genpath, folders, labels, timesteps, var_str='H', overwhat='V', r_m
if savefig:
fig.savefig('mass_cgm.svg')
def rad_Zfrac(genpath, folders, labels, timesteps, r_max=150, factor=2, nbins=100, split_flow=False, xlims=None, ylims=None, saveinfile=False, savefig=False):
"""
Plots the average density of the CGM and it's average mass-weighted temperature through radial chunks. It is also possible to remove satellites from the selection.
Stacks timesteps up to the rvir of the first timestep (the smaller one).
Parameters:
-----------
var_str: 'gas' or 'HI'
The variable to plot: all the gas or neutral gas.
overwhat: 'v' or 'r'
Choose if either sum(var*dV)/sum(dV) or sum(var*dV)/sum(r) (a la Mitchell) is plotted
nbins: int
Number of bins along x axis (number of points in the graph).
Example:
--------
timestep = 17
genpath = '/mnt/lyoccf/scratch/mrey/outputs/2_G8/'
folders = ['1_base/','2_nstar/1_n_star=1/', '2_nstar/2_n_star=50/']
labels = ['1', '2', '3']
plot_cgm_time(genpath, folders, labels, timestep, shape="no_rebloch", xlims=None, ylims=None, saveinfile=False, savefig=False)
"""
from ratatouille import readNsave as ras
import matplotlib.pyplot as plt
import numpy as np
# If there is just one folder
if (not isinstance(folders,list)):
folders=[folders]
# timesteps must be an array
if isinstance(timesteps,float) or isinstance(timesteps,int):
timesteps = [timesteps]
kpc2cm = 3.086e21
r_max_cm = r_max*kpc2cm
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(10,6))
prop_cycle = plt.rcParams['axes.prop_cycle']
colors = prop_cycle.by_key()['color']
for index, folder in enumerate(folders):
RamsesDir = genpath+folder
print('Folder', folder)
var_arr = []
for timestep in timesteps:
info = ras.extract_info(timestep, RamsesDir, saveinfile=saveinfile)
cu2cm = info['unit_l'] * info['boxlen']
center, _ = ras.get_rad(RamsesDir, timestep, 1)
center_cm = [cen*cu2cm for cen in center]
_, cells, cell_pos, _ = ras.extract_cells(RamsesDir, timestep, factor=factor, saveinfile=saveinfile)
_, _, _, _, _, _, _, _, _, _, _, _, mass, Z = cells
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)
# Bin radially sum(mass*Z)/sum(mass)
var_hist, rad = np.histogram(norm_r[norm_r<r_max_cm], range=[0,r_max_cm], weights=Z[norm_r<r_max_cm]*mass[norm_r<r_max_cm], bins=nbins)
nrmlz, _ = np.histogram(norm_r[norm_r<r_max_cm], range=[0,r_max_cm], weights=mass[norm_r<r_max_cm], bins=nbins)
var_arr.append(var_hist/nrmlz)
perc_min, perc_max = 15.9, 84.1
x_axis = (rad[:-1]+rad[1:])/2/kpc2cm
var_arr = np.array(var_arr)
median_vartot = np.median(var_arr, axis=0)
ax.plot(x_axis, median_vartot, label=labels[index], color=colors[index])
ax.fill_between(x_axis,np.percentile(var_arr,perc_min, axis=0), np.percentile(var_arr,perc_max, axis=0), color=colors[index], alpha=0.25)
ax.legend(frameon=False)
ax.set_yscale('log')
ax.set_xlabel('$r$ [kpc]')
ax.set_ylabel('Metal mass fraction')
ax.set_xlim(xlims)
ax.set_ylim(ylims)
plt.tight_layout()
plt.show()
if savefig:
fig.savefig(f'rad_Zprof.svg')
def rad_col(genpath, folders, timesteps, labels, r_max=150, lmax=20, var_str='HI', split_flows=False, nbins = 1000, factor = 2,
plot=True, xlims=None, ylims=None, saveinfile=True, no_err=True, logx=True):
"""
......
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