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

Improvement: 'plot_phase_diag' -> included nicer side PDFs.

parent 65a76e3b
......@@ -111,8 +111,7 @@ from ratatouille.ks import *
########################################################## Phase Diagram ##########################################################
###################################################################################################################################
def plot_phase_diag(genpath, folders, timesteps, var_str='H', cgm_gas=False, no_sat=False, bins=300, contour=False, levels=None, sigma=3, \
def plot_phase_diag(genpath, folders, timesteps, var_str='H', cgm_gas=False, no_sat=False, bins=300, plot_PDF=True, contour=False, levels=None, sigma=3, \
xlims=None, ylims=None, vmin=None, vmax=None, factor=1, saveinfile=True, savefig=False):
"""
Plots phase diagram from different output files (for the same timestep) one after another.
......@@ -134,7 +133,6 @@ def plot_phase_diag(genpath, folders, timesteps, var_str='H', cgm_gas=False, no_
# If there is just one folder
if folders==None or folders=='':
folders=['']
labels=['']
if (not isinstance(folders,list)):
folders=[folders]
if (not isinstance(timesteps,list)) and (not isinstance(timesteps, np.ndarray)):
......@@ -228,21 +226,61 @@ def plot_phase_diag(genpath, folders, timesteps, var_str='H', cgm_gas=False, no_
phase_diag_gauss = np.ma.masked_where(phase_diag_gauss <= 0, phase_diag_gauss)
plt.contour(phase_diag_gauss.T,levels=levels,origin='lower',norm=LogNorm(vmin=vmin,vmax=vmax),colors='k', extent=extent)
if plot_PDF:
# Adds T PDF on the left of the plot.
PDF,bin_edges = np.histogram(np.log10(T), bins=bins, weights=weight, density=True)
x_axis = (bin_edges[:-1]+bin_edges[1:])*0.5
if xlims:
shift = min(PDF)-np.log10(xlims[0])
else:
shift = min(PDF)-ax.get_xlim()[0]
ax.plot(PDF-shift, x_axis, '-k')
# Adds nH PDF on the bottom of the plot.
PDF,bin_edges = np.histogram(np.log10(nH), bins=bins, weights=weight, density=True)
x_axis = (bin_edges[:-1]+bin_edges[1:])*0.5
if xlims:
shift = min(PDF)-np.log10(ylims[0])
else:
shift = min(PDF)-ax.get_ylim()[0]
ax.plot(x_axis, PDF-shift, '-k')
if (var_str=='HI') or (var_str=='MgII') or (var_str=='CIV') or (var_str=='OVI'):
# Adds two temperature split chosen for ouflows.
ax.axhline(4.25, color='b', alpha=0.5)
ax.axhline(5.05, color='r', alpha=0.5)
# Adds T PDF on the left of the plot.
PDF,bin_edges = np.histogram(np.log10(T), bins=bins, weights=weight, density=True)
x_axis = (bin_edges[:-1]+bin_edges[1:])*0.5
ax.plot(PDF+ax.get_xlim()[0]+0.5, x_axis, '-k')
# Plot T_vir
if var_str=='OVI':
info = ras.extract_info(timestep, RamsesDir, saveinfile=True)
cu2cm = info['unit_l'] * info['boxlen']
cu2kpc = cu2cm / 3.086e21
center, r200 = ras.get_r200_cen(RamsesDir,timestep)
center_cm, r200_cm = [cen*cu2cm for cen in center], r200*cu2cm
star_mass, star_x, star_y, star_z, _, _, _, _, _, _, _ = ras.extract_stars(RamsesDir, timestep, factor=1, rmsat=True, saveinfile=True)
rad_stars = np.sqrt((star_x-center_cm[0])**2+(star_y-center_cm[1])**2+(star_z-center_cm[2])**2)
M_starhalo = np.sum(star_mass[rad_stars < r200_cm])
M_gashalo = np.sum(rho*(cell_dx)**3) # Reusing loaded data so the result is wrong if factor is not one.
M_DM = ras.get_m200(genpath+folder,timestep)
G = 4.3e-3 # [pc/Msun (km/s)^2]
mp = 1.672e-27 # [kg]
kb = 1.38e-23 # [m²kg/(s²K)]
GM_R = G* ((M_gashalo+M_starhalo)/1.989e33+M_DM)/(r200*cu2kpc*1e3) * 1000**2 # (m/s)^2
T_virh = 2/5 * GM_R * mp/kb
print(T_virh)
ax.axhline(np.log10(T_virh), color='k', alpha=0.5)
plt.xlabel('Density n$_{{H}}$ [atom.cm$^{{-3}}$]')
plt.ylabel('Temperature [K]')
plt.text(0.95, 0.95, ' '.join(folders[index].split('_')[1:])[:-1], horizontalalignment='right',verticalalignment='top', transform=ax.transAxes, size=20)
ax.set_xlim(np.log10(xlims))
ax.set_ylim(np.log10(ylims))
if xlims:
ax.set_xlim(np.log10(xlims))
if ylims:
ax.set_ylim(np.log10(ylims))
ax.xaxis.set_major_formatter(plt.FuncFormatter(axis_format_10pow))
ax.yaxis.set_major_formatter(plt.FuncFormatter(axis_format_10pow))
......@@ -257,7 +295,6 @@ def plot_phase_diag(genpath, folders, timesteps, var_str='H', cgm_gas=False, no_
if savefig:
fig.savefig(f'PhaseDiag_{folder[:-1]}_{timestep}.svg',bbox_inches='tight',pad_inches=0)
def plot_evol_phase_diag(RamsesDir, max_timestep, bins=300, xlims=None, ylims=None, vmin=None, vmax=None, factor=1, savefig=False, saveinfile=True):
"""
Plots phase diagrams for all timesteps until max_timestep (for the same output).
......
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