Commit 8ad9531c authored by Maxime Rey's avatar Maxime Rey
Browse files

Improvement: 'Ms_Mh_relation' now with DM from DM only simulation + data taken...

Improvement: 'Ms_Mh_relation' now with DM from DM only simulation + data taken from paper in a nicer way.
parent b15ed34f
......@@ -3112,7 +3112,7 @@ def plot_cgm_time(genpath, folders, labels, timesteps, shape="no_rebloch",factor
###################################################################################################################################
############################################################## M*-Mh ##############################################################
###################################################################################################################################
def Ms_Mh_relation(genpath, folders, timesteps, labels, xmin=None, ymin=None, norm=False, plotgas=False):
def Ms_Mh_relation(genpath, folders, timesteps, labels, xmin=None, ymin=None, norm=False, plotgas=False, savefig=False):
"""
Plots the M_star - M_halo relation for specific snapshots a a simulation (with a different namelist, you need to change the timesteps considered ^^').
Also includes several observational relations read from specific files. The path has to be changed if this is not the original simulation. Very specific one-time use function.
......@@ -3125,135 +3125,103 @@ def Ms_Mh_relation(genpath, folders, timesteps, labels, xmin=None, ymin=None, no
import numpy as np
import os
_, ax = plt.subplots(nrows=1, ncols=1, figsize=(12,7.2))
shape_s = ['o', 'X', 'P', 'D']
shape_g = ['.', 'x', '+', '1']
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 7.2))
shape_s = ['o', 'X', 'P', 'D']
shape_g = ['.', 'x', '+', '1']
g2Msun = 1/1.989e33
vmin, vmax = 0, 7
colors = ['k', 'blue', 'cyan', 'lightgreen', 'yellow', 'orange', 'red', 'darkmagenta']
cm = LinearSegmentedColormap.from_list('custom_cmap', colors, N=100)
mark, lab = [], []
for ind_f, folder in enumerate(folders):
print('Folder ', folder)
RamsesDir = genpath+folder
M_dm, redshift = [], []
M_star_ISM, M_gas_ISM = [], []
M_star_CGM, M_gas_CGM = [], []
M_star_ISM, M_gas_ISM = [], []
M_star_CGM, M_gas_CGM = [], []
for timestep in timesteps:
try:
info = ras.extract_info(timestep, RamsesDir, saveinfile=True)
except:
continue
redshift.append(1/info['aexp']-1)
star_mass,star_x,star_y,star_z,_,_,_,_,_,_,_ = ras.extract_stars(RamsesDir, timestep, factor=1, rmsat=True, saveinfile=True)
center, r200 = ras.get_r200_cen(RamsesDir, timestep)
cu2cm = info['unit_l'] * info['boxlen']
center_cm, r200_cm = [cen*cu2cm for cen in center], r200*cu2cm
# Stars
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_star_ISM.append(np.sum(star_mass[rad_stars<0.1*r200_cm])/1.989e33) # Mass of stars in the ISM.
M_star_CGM.append(np.sum(star_mass[(rad_stars>0.1*r200_cm) & (rad_stars<r200_cm)])/1.989e33) # Mass of stars in the CGM.
M_star_ISM.append(np.sum(star_mass[rad_stars<0.1*r200_cm])/1.989e33)
M_star_CGM.append(np.sum(star_mass[(rad_stars>0.1*r200_cm) & (rad_stars < r200_cm)])/1.989e33)
# Gas
_, cells, cell_pos, cell_l = ras.extract_cells(RamsesDir,timestep,factor=1,saveinfile=True,verbose=False)
mass = cells[12]
rad_gas = np.sqrt((cell_pos[:,0]-center_cm[0])**2+(cell_pos[:,1]-center_cm[1])**2+(cell_pos[:,2]-center_cm[2])**2)
M_gas_ISM.append(np.sum(mass[rad_gas<0.1*r200_cm])*g2Msun) # Gas mass in the ISM.
M_gas_CGM.append(np.sum(mass[(rad_gas>0.1*r200_cm) & (rad_gas<r200_cm)])*g2Msun) # Gas mass in the CGM.
M_gas_ISM.append(np.sum(mass[rad_gas<0.1*r200_cm])*g2Msun)
M_gas_CGM.append(np.sum(mass[(rad_gas>0.1*r200_cm) & (rad_gas<r200_cm)])*g2Msun)
# DM
M_dm.append(ras.get_m200(genpath+folder,timestep))
# M_dm.append(ras.get_m200(genpath+folder,timestep)) # DM from simulation. Should use DM from a DM-only simulation.
M_dm.append(ras.get_m200('/scratch/Cral/mrey/ICs/Zoom-56-35796/dmRun2-LSS1Rvir', timestep))
M_star_ISM, M_gas_ISM = np.array(M_star_ISM), np.array(M_gas_ISM)
M_star_CGM, M_gas_CGM = np.array(M_star_CGM), np.array(M_gas_CGM)
M_dm = np.array(M_dm )
M_halo = (M_dm + M_star_ISM + M_gas_ISM + M_star_CGM + M_gas_CGM)
if norm:
starplot = plt.scatter(M_halo, M_star_ISM/M_halo, vmin=vmin, vmax=vmax, c=redshift, cmap=cm, edgecolors='k', marker=shape_s[ind_f], s=100)
if plotgas:
gasplot = plt.scatter(M_halo, M_gas_ISM /M_halo, vmin=vmin, vmax=vmax, c=redshift, cmap=cm, marker=shape_g[ind_f], s=100)
mark.append((starplot, gasplot))
else:
mark.append((starplot))
M_halo = np.array(M_dm) # M_star_ISM + M_gas_ISM + M_star_CGM + M_gas_CGM
# Plot selected simulation points.
starplot = plt.scatter(M_halo, M_star_ISM/M_halo if norm else M_star_ISM, vmin=vmin, vmax=vmax, c=redshift, cmap=cm, edgecolors='k', marker=shape_s[ind_f], s=100)
if plotgas:
gasplot = plt.scatter(M_halo, M_gas_ISM/M_halo if norm else M_gas_ISM, vmin=vmin, vmax=vmax, c=redshift, cmap=cm, marker=shape_g[ind_f], s=100)
mark.append((starplot, gasplot))
else:
starplot = plt.scatter(M_halo, M_star_ISM, vmin=vmin, vmax=vmax, c=redshift, cmap=cm, edgecolors='k', marker=shape_s[ind_f], s=100)
if plotgas:
gasplot = plt.scatter(M_halo, M_gas_ISM , vmin=vmin, vmax=vmax, c=redshift, cmap=cm, marker=shape_g[ind_f], s=100)
mark.append((starplot, gasplot))
else:
mark.append((starplot))
mark.append((starplot))
lab.append(labels[ind_f])
Mh_beh_z0 = np.array([9.781E9, 1.293E10, 1.710E10, 2.330E10, 3.104E10, 4.585E10, 5.390E10, 6.672E10, 8.078E10, 9.087E10, 1.068E11, 1.413E11, 1.673E11, 1.938E11, 2.399E11, 2.927E11, 3.366E11, 3.731E11, 4.419E11, 5.195E11, 5.975E11, 6.871E11, 7.396E11, 8.320E11, 9.638E11, 1.133E12, 1.322E12, 1.637E12, 1.775E12])
Ms_beh_z0 = np.array([0.001723, 0.001896, 0.002122, 0.002401, 0.002686, 0.003074, 0.003270, 0.003557, 0.003958, 0.004305, 0.004790, 0.006235, 0.007630, 0.009028, 0.01168, 0.01430, 0.01673, 0.01830, 0.02094, 0.02265, 0.02409, 0.02506, 0.02535, 0.02578, 0.02593, 0.02578, 0.02507, 0.02385, 0.02306])
Mh_beh_z0_err = np.array([2.347E10, 5.551E10, 1.322E11, 3.127E11, 7.450E11, 1.775E12])
Beh_z0_err = np.array([[0.001705, 0.002457, 0.004355, 0.01123, 0.01852, 0.01694],[0.003380, 0.004841, 0.008344, 0.02105, 0.03469, 0.03121]])
if not norm:
Ms_beh_z0 *= Mh_beh_z0
Beh_z0_err *= Mh_beh_z0_err
plt.plot(Mh_beh_z0, Ms_beh_z0,'k', label='Behroozi et al. 2013')
plt.fill_between(Mh_beh_z0_err, Beh_z0_err[0], Beh_z0_err[1], color='k', alpha=0.1)
# Read data from files
from os import listdir
path2files = 'data/Behroozi/'
for fileuh in sorted(listdir(path2files)):
colors = ['k', 'b', 'c', 'g', 'y', 'orange', 'r', 'purple']
# Pick style and label depending on paper
paper = fileuh[:-9]
if paper == 'Behroozi et al. 2019' or fileuh[:-11] == 'Behroozi et al. 2019':
label, ls = '', 'solid'
paper = fileuh.split(',')[0]
red = int(fileuh.split(',')[1][3])
if paper == 'Behroozi et al. 2019':
ls = 'solid'
elif paper == 'Behroozi et al. 2013':
ls = 'dashdot'
elif paper == 'Behroozi et al. 2010':
ls = 'dashed'
elif paper == 'Moster et al. 2013':
label, ls = paper, 'dashed'
elif paper == 'Moster et al. 2018':
label, ls = paper, (0, (5,1))
elif paper == 'McCracken et al. 2012':
label, ls = paper, 'dotted'
elif paper == 'Coupon et al. 2015':
label, ls = paper, (0, (1,10))
# Select redshift from filename and legend for Behroozi.
try:
red = int(fileuh[-7:-6])
plt.fill_between('', '', '', color=colors[red], alpha=0.25, label=fileuh[:-11])
except:
red = int(fileuh[-5:-4])
ls = 'dotted'
fileuh2 = path2files+fileuh
dat = np.loadtxt(fileuh2, usecols = (0,1))
dat = np.loadtxt(fileuh2,usecols=(0,1))
try:
errmS = np.loadtxt(fileuh2, usecols = 2, dtype=str)
errm = np.array([float(errstr[2:]) for errstr in errmS])
errp = np.loadtxt(fileuh2, usecols = 3)
data = np.concatenate((dat, errm.reshape((len(errm),1)), errp.reshape((len(errp),1))), axis=1) # x, y, y-, y+
errmS = np.loadtxt(fileuh2,usecols=2, dtype=str)
errm = np.array([float(errstr[2:]) for errstr in errmS])
errp = np.loadtxt(fileuh2,usecols=3)
data = np.concatenate((dat,errm.reshape((len(errm), 1)),errp.reshape((len(errp),1))),axis=1) # x, y, y-, y+
except IndexError:
data = dat
data = dat
try: # If there's errorbars
if not norm:
errm, errp = data[:,0]*(data[:,1]-data[:,2]), data[:,0]*(data[:,1]+data[:,3])
plt.plot(data[:,0], data[:,0]*data[:,1], color=colors[red], ls=ls, label=label)
else:
errm, errp = (data[:,1]-data[:,2]), (data[:,1]+data[:,3])
plt.plot(data[:,0], data[:,1], color=colors[red], ls=ls, label=label)
non_zero = errm!=errp
plt.fill_between(data[:,0][non_zero], errm[non_zero], errp[non_zero], color=colors[red], ls=ls, alpha=0.1)
except:
if not norm:
plt.plot(data[:,0], data[:,0]*data[:,1], color=colors[red], ls=ls, label=label)
else:
plt.plot(data[:,0], data[:,1], color=colors[red], ls=ls, label=label)
# Won't work if no errorbars.
if not norm:
errm, errp = data[:,0]*(data[:,1]-data[:,2]), data[:,0]*(data[:,1]+data[:,3])
plt.plot(data[:,0], data[:,0]*data[:,1],color=colors[red], ls=ls, label=paper if red==0 else '')
else:
errm, errp = (data[:,1]-data[:,2]), (data[:,1]+data[:,3])
plt.plot(data[:,0], data[:,1],color=colors[red], ls=ls, label=paper if red==0 else '')
non_zero = errm!=errp
plt.fill_between(data[:,0][non_zero], errm[non_zero], errp[non_zero], color=colors[red], ls=ls, alpha=0.1)
# Legend
leg_plot = plt.legend(mark, lab ,handler_map={tuple: HandlerTuple(ndivide=None)}, loc='lower right', fontsize='small')
leg_plot = plt.legend(mark, lab, handler_map={tuple: HandlerTuple(ndivide=None)}, loc='lower right', fontsize='small')
ax.add_artist(leg_plot)
if plotgas:
legend_elements = [Line2D([0], [0], marker='X', linestyle='None', color='k', markersize=10, label='M$_{*, ISM}$', markerfacecolor='none'),
Line2D([0], [0], marker='x', linestyle='None', color='k', markersize=10, label='M$_{gas, ISM}$')]
leg_symb = plt.legend(handles=legend_elements, bbox_to_anchor=(0,1,1,0), loc="lower left", ncol=3, fontsize='small')
Line2D([0], [0], marker='x', linestyle='None', color='k', markersize=10, label='M$_{gas, ISM}$')]
leg_symb = plt.legend(handles=legend_elements, bbox_to_anchor=(0, 1, 1, 0), loc="lower left", ncol=3, fontsize='small')
ax.add_artist(leg_symb)
ax.legend(loc="upper left")
......@@ -3269,6 +3237,9 @@ def Ms_Mh_relation(genpath, folders, timesteps, labels, xmin=None, ymin=None, no
plt.xlim(xmin)
plt.ylim(ymin)
if savefig:
fig.savefig('Ms_Mh.svg')
def content_t(genpath, folders, labels, is_ISM=True, DM=False, xlims=None, ylims=None):
"""
......
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