Commit 7c4806a1 authored by Maxime Rey's avatar Maxime Rey
Browse files

Improvement: Phase diagrams now contain stacked PDFs.

parent 37c38931
......@@ -106,7 +106,7 @@ from ratatouille.ks import *
###################################################################################################################################
###################################################### PDFs & Phase Diagrams ######################################################
###################################################################################################################################
def PDF(what2plot, genpath, folders, timesteps, labels, weight=None, density=True, bins=300, \
def PDF(what2plot, genpath, folders, timesteps, labels, weight=None, density=False, bins=300, \
x_log=True, y_log=True, xlims=None, ylims=None, factor=1, savefig=False, saveinfile=True, rmsat=False):
"""
Plots PDF for 'T', 'Tmu', 'nH', 'rho', 'pres', 'mass', 'star_mass', 'star_age', 'star_mets' or 'star_minit'.
......@@ -278,14 +278,15 @@ def plot_phase_diag(genpath, folders, timesteps, var_str='H', cgm_gas=False, no_
g2Msun = 1./1.9891e33
for index, folder in enumerate(folders):
RamsesDir = genpath + folder
PD_all = np.empty((len(timesteps), bins, bins))
T_all = []
nH_all = []
weight_all = []
T_all, nH_all, weight_all = [], [], [] # Array for method 1
PD_all = np.empty((len(timesteps), bins, bins)) # Array for method 2
PDF_Tall = np.empty((len(timesteps), bins)) # Array for PDFs
PDF_nHall = np.empty((len(timesteps), bins))
for istep, timestep in enumerate(timesteps):
_, cells, cell_pos, _ = ras.extract_cells(RamsesDir,timestep,factor=factor,saveinfile=saveinfile)
cell_dx, rho, nH, _, vx, vy, vz, xHII, _, _, _, T, _, _ = cells
# Select variables to plot
if var_str=='H':
weight = rho*cell_dx**3*g2Msun
if var_str=='HI':
......@@ -319,7 +320,7 @@ def plot_phase_diag(genpath, folders, timesteps, var_str='H', cgm_gas=False, no_
norm_r = np.linalg.norm(vec_r,2, axis=1)
vec_v = np.array([vx,vy,vz]).transpose()
weight = np.array([np.dot(vec_v[ind],vec_r[ind])/norm_r[ind] for ind in range(len(norm_r))])/1e5 # Proj. vec_v vector on vec_r (going out or in) [km/s]
# Compute restrictions on data
restr = True
if no_sat:
sat = is_sat(RamsesDir,timestep, factor=1)
......@@ -334,15 +335,14 @@ def plot_phase_diag(genpath, folders, timesteps, var_str='H', cgm_gas=False, no_
if vmin==None and (np.max(weight)/np.min(weight)>1e16):
vmin = np.max(weight)/1e16
print('Used a default maximum limit for vmin.')
# Phase diagram of stacked data or average of phase diagrams.
# Phase diagrams
if method1:
T_all = np.append(T_all, T)
nH_all = np.append(nH_all, nH)
weight_all = np.append(weight_all, weight)
else:
if xlims and ylims:
sel = tuple([(xlims[0]<nH) & (xlims[-1]>nH) & (ylims[0]<T) & (ylims[-1]>T)]) # If lims inputed, use them to restrict data.
sel = tuple([(xlims[0]<nH) & (xlims[-1]>nH) & (ylims[0]<T) & (ylims[-1]>T)]) # restrict data range by the first snapshot if none given
PD, xedges, yedges = np.histogram2d(np.log10(nH[sel]),np.log10(T[sel]),bins=(bins,bins),weights=weight[sel])
else:
PD, xedges, yedges = np.histogram2d(np.log10(nH),np.log10(T),bins=(bins,bins),weights=weight)
......@@ -352,6 +352,12 @@ def plot_phase_diag(genpath, folders, timesteps, var_str='H', cgm_gas=False, no_
if not ylims:
ylims = 10**yedges[0], 10**yedges[-1]
PD_all[istep,:,:] = PD
# Compute PDFs
if plot_PDF:
PDF_nHall[istep,:], nHbin_edges = np.histogram(np.log10(nH), bins=bins, range=np.log10(xlims), weights=weight, density=False) # True ??
PDF_Tall[istep,:] , Tbin_edges = np.histogram(np.log10(T) , bins=bins, range=np.log10(ylims), weights=weight, density=False) # True ??
# Phase diagram with method 1 or 2.
if method1:
phase_diag,xedges,yedges = np.histogram2d(np.log10(nH_all),np.log10(T_all),bins=(bins,bins),weights=weight_all)
phase_diag /= len(timesteps) # Normalise total mass by the number of timesteps considered.
......@@ -375,52 +381,49 @@ 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)
# Have to replace by stacked.
if plot_PDF:
perc_min, perc_max = 15.9, 84.1
nrm = 2e8 # Handpicked normalisation factor to get a nicer plot. Otherwise density has to be used...
# 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')
PDF_T = np.median(PDF_Tall, axis=0) / nrm
shift_T = -min(PDF_T)+np.log10(xlims[0])
xaxis_T = (Tbin_edges[:-1]+Tbin_edges[1:])*0.5
ax.plot(PDF_T+shift_T, xaxis_T, '-k')
ax.fill_betweenx(xaxis_T,np.percentile(PDF_Tall,perc_min, axis=0)/nrm+shift_T, np.percentile(PDF_Tall,perc_max, axis=0)/nrm+shift_T, color='k', alpha=0.25)
# 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')
PDF_nH = np.median(PDF_nHall, axis=0) / nrm
shift_nH = -min(PDF_nH)+np.log10(ylims[0])
xaxis_nH = (nHbin_edges[:-1]+nHbin_edges[1:])*0.5
ax.plot(xaxis_nH, PDF_nH+shift_nH, '-k')
ax.fill_between(xaxis_nH,np.percentile(PDF_nHall,perc_min, axis=0)/nrm+shift_nH, np.percentile(PDF_nHall,perc_max, axis=0)/nrm+shift_nH, color='k', alpha=0.25)
if (var_str=='HI') or (var_str=='MgII') or (var_str=='CIV') or (var_str=='OVI'):
# 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)
# Plot T_vir
if var_str=='OVI':
info = ras.extract_info(timestep, RamsesDir, saveinfile=saveinfile)
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
# Plot T_vir
if var_str=='OVI':
print("Only the T_vir of the last ouput is included. You might want to 'stack' it.")
info = ras.extract_info(timestep, RamsesDir, saveinfile=saveinfile)
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=saveinfile)
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)
star_mass, star_x, star_y, star_z, _, _, _, _, _, _, _ = ras.extract_stars(RamsesDir, timestep, factor=1, rmsat=True, saveinfile=saveinfile)
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}}$]')
......
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