Commit 37c38931 authored by Maxime Rey's avatar Maxime Rey
Browse files

Bugfix: Phase Diagrams can now stack data.

parent 440b6528
......@@ -245,9 +245,10 @@ def PDF(what2plot, genpath, folders, timesteps, labels, weight=None, density=Tru
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):
xlims=None, ylims=None, vmin=None, vmax=None, factor=1, saveinfile=True, savefig=False, method1=True):
"""
Plots phase diagram from different output files (for the same timestep) one after another.
"method1" appends all timesteps together and make a phase diagram out of it if True. If False it averages the phase diagrams of each timesteps.
Example:
--------
......@@ -277,8 +278,12 @@ 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
for timestep in timesteps:
_, cells, cell_pos, _ = ras.extract_cells(RamsesDir,timestep,factor=factor,saveinfile=True)
PD_all = np.empty((len(timesteps), bins, bins))
T_all = []
nH_all = []
weight_all = []
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
if var_str=='H':
......@@ -303,7 +308,7 @@ def plot_phase_diag(genpath, folders, timesteps, var_str='H', cgm_gas=False, no_
nMgII = ras.extract_ions(RamsesDir, timestep, 'MgII', factor=factor) # MgII / cm^3
weight = nMgII/(nMgI+nMgII) # MgII fraction
elif var_str=='v_rad':
info = ras.extract_info(timestep, RamsesDir, saveinfile=True)
info = ras.extract_info(timestep, RamsesDir, saveinfile=saveinfile)
cu2cm = info['unit_l'] * info['boxlen']
if info['is_cosmo']:
center, _ = ras.get_r200_cen(RamsesDir,timestep)
......@@ -329,25 +334,35 @@ 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.')
if timestep==timesteps[0]:
T_all = T
nH_all = nH
weight_all = weight
else:
# Phase diagram of stacked data or average of 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)
phase_diag,_,_ = np.histogram2d(np.log10(nH),np.log10(T),bins=(bins,bins),weights=weight)
extent=(np.log10(nH.min()), np.log10(nH.max()), # To have good values on imshow
np.log10(T.min()), np.log10(T.max())) # Otherwise he counts bins...
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.
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)
if timestep==timesteps[0]:
if not xlims:
xlims = 10**xedges[0], 10**xedges[-1]
if not ylims:
ylims = 10**yedges[0], 10**yedges[-1]
PD_all[istep,:,:] = PD
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.
else:
phase_diag = np.average(PD_all, axis=0)
extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]] # Extent used to have good values instead of the number of bins...
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(10,10))
if var_str=='v_rad':
im = plt.imshow(phase_diag.T,interpolation='nearest',origin='lower',vmin=vmin,vmax=vmax,cmap='seismic', extent=extent)
else:
im = plt.imshow(phase_diag.T,interpolation='nearest',origin='lower',norm=LogNorm(vmin=vmin,vmax=vmax),cmap='BuPu_r', extent=extent)
im = plt.imshow(phase_diag.T,interpolation='nearest',origin='lower',norm=LogNorm(vmin=vmin,vmax=vmax),cmap='BuPu_r', extent=extent) # Enlèves le "_r"
if contour:
if folder==folders[0]:
extent_0 = extent
......@@ -360,6 +375,7 @@ 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:
# Adds T PDF on the left of the plot.
PDF,bin_edges = np.histogram(np.log10(T), bins=bins, weights=weight, density=True)
......@@ -381,18 +397,18 @@ def plot_phase_diag(genpath, folders, timesteps, var_str='H', cgm_gas=False, no_
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)
# 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=True)
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=True)
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.
......
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