mirror of
https://github.com/Richard-Sti/csiborgtools.git
synced 2024-12-22 17:28:02 +00:00
28 KiB
28 KiB
In [1]:
import numpy as np
%matplotlib notebook
import matplotlib.pyplot as plt
try:
import csiborgtools
except ModuleNotFoundError:
import sys
sys.path.append("../")
import csiborgtools
import utils
%load_ext autoreload
%autoreload 2
import joblib
from os.path import join
In [2]:
from astropy.cosmology import FlatLambdaCDM, z_at_value
from astropy import units
In [6]:
Nsim = 9844
Nsnap = 1016
# data, box = utils.load_processed(Nsim, Nsnap)
data, box = utils.load_processed(Nsim, Nsnap)
In [10]:
X = np.vstack([data["peak_{}".format(p)] for p in ("x", "y", "z")]).T
In [12]:
from sklearn.neighbors import NearestNeighbors
In [13]:
neighbors = NearestNeighbors()
In [14]:
neighbors.fit(X)
Out[14]:
In [19]:
p = X[0, :]
neighbors.kneighbors(p.reshape(-1,3))
Out[19]:
In [16]:
Out[16]:
In [8]:
data["peak_y"]
Out[8]:
In [ ]:
In [ ]:
In [ ]:
zcosmo = box.box2cosmoredshift(data["dist"])
z
In [ ]:
zpec = box.box2pecredshift(*[data[p] for p in ["vx", "vy", "vz", "peak_x", "peak_y", "peak_z"]])
zobs = box.box2obsredshift(*[data[p] for p in ["vx", "vy", "vz", "peak_x", "peak_y", "peak_z"]])
In [ ]:
m = zcosmo < 0.05
plt.figure()
# plt.scatter(zcosmo[m], zcosmo[m] - zobs[m], s=0.05)
plt.scatter(zcosmo[m], zpec[m], s=0.05, rasterized=True)
t = np.linspace(0, zcosmo[m].max())
plt.axhline(0, c="red", ls="--")
plt.xlabel(r"$z_{\rm cosmo}$")
plt.ylabel(r"$z_{\rm pec}$")
plt.tight_layout()
plt.savefig("../plots/redshift.png", dpi=450)
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
data[""]
In [ ]:
In [ ]:
In [ ]:
In [ ]:
cosmo = FlatLambdaCDM(H0=70, Om0=0.3, Ob0=0.05)
cosmo.comoving_distance()
x = 10000 * units.Mpc
z_at_value(cosmo.comoving_distance, x)
In [ ]:
from astropy import constants
In [ ]:
constants.c.value * 1e-8
In [ ]:
n_sims = csiborgtools.io.get_csiborg_ids("/mnt/extraspace/hdesmond")
In [ ]:
np.where(n_sims == 7660)
In [ ]:
for n_sim in n_sims:
simpath = csiborgtools.io.get_sim_path(n_sim)
maxsnap = csiborgtools.io.get_maximum_snapshot(simpath)
box = csiborgtools.units.BoxUnits(maxsnap, simpath)
print(maxsnap, box._aexp)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
simpath = csiborgtools.io.get_sim_path(7660)
In [ ]:
simpath
In [ ]:
box = csiborgtools.units.BoxUnits(999, simpath)
In [ ]:
box._aexp
In [ ]:
In [ ]:
In [ ]:
planck = utils.load_planck2015(200)
mcxc = utils.load_mcxc(200)
indxs = csiborgtools.io.match_planck_to_mcxc(planck, mcxc)
In [ ]:
groups = utils.load_2mpp_groups()
Nsim = 9844
Nsnap = 1016
data = utils.load_processed(Nsim, Nsnap)
In [ ]:
In [ ]:
plt.figure()
m = data["m500"] > 1e13
plt.scatter(data["ra"][m], data["dec"][m], label="CSiBORG", s=3)
plt.scatter(groups["RA"], groups["DEC"], label="2M++ galaxy groups", s=3, marker="x")
plt.show()
In [ ]:
In [ ]:
In [ ]:
RAcoma = (12 + 59/60 + 48.7 / 60**2) * 15
DECcoma = 27 + 58 / 60 + 50 / 60**2
RAvirgo = (12 + 27 / 60) * 15
DECvirgo = 12 + 43/60
plt.figure()
plt.scatter(mcxc["RAdeg"], mcxc["DEdeg"], label="MCXC")
plt.scatter(planck["RA"], planck["DEC"], label="Plank",s=8, c="red")
plt.scatter(RAcoma, DECcoma, label="Coma", s=30, marker="x")
plt.scatter(RAvirgo, DECvirgo, label="Virgo", s=30, marker="x")
plt.legend(framealpha=0.5)
plt.xlabel("RA")
plt.ylabel("DEC")
plt.title("Clusters below 200 Mpc")
plt.savefig("../plots/clusters_radec.png", dpi=450)
plt.show()
In [ ]:
In [ ]:
In [ ]:
plt.figure()
plt.scatter(mcxc["COMDIST"], mcxc["M500"], label="MCXC")
# yerr = np.vstack([planck["MSZ"] - planck["MSZ_ERR_LOW"],
# planck["MSZ_ERR_UP"] - planck["MSZ"]])
yerr = np.vstack([planck["MSZ_ERR_LOW"], planck["MSZ_ERR_UP"]])
plt.errorbar(planck["COMDIST"], planck["MSZ"], yerr, label="Plank", fmt=" ", capsize=3, color="red")
plt.yscale("log")
plt.legend()
plt.title("Clusters below 200 Mpc")
plt.xlabel(r"$D_{\rm c} / \mathrm{Mpc}$")
plt.ylabel(r"$M_{\rm 500c} / M_\odot$")
# plt.savefig("../plots/clusters_mass_dist.png", dpi=450)
plt.show()
In [ ]:
d["RAdeg"]
In [ ]:
In [ ]:
In [ ]:
d.dtype.names
In [ ]:
d["Cat"]
In [ ]:
type(d["MCXC"])
In [ ]:
np.asanyarray(["aasdasdaasdasdasdad", "bsdfadadfasddsgasdg"]).dtype
In [ ]:
In [ ]:
d["MCXC"]
In [ ]:
In [ ]:
plt.figure()
# plt.scatter(d["RAdeg"], d["DEdeg"], s=1)
plt.scatter(d["z"], d["M500"], s=1)
plt.yscale("log")
plt.show()
In [ ]:
In [ ]:
In [ ]:
Nsim = 9844
Nsnap = 1016
data = utils.load_processed(Nsim, Nsnap)
In [ ]:
bins = np.arange(11.8, 15.4, 0.2)
plt.figure()
x, mu, std = csiborgtools.match.number_density(data, "m200", bins, 200, True)
plt.errorbar(x, mu, std, capsize=4, label=r"$M_{200c}$")
x, mu, std = csiborgtools.match.number_density(data, "m500", bins, 200, True)
plt.errorbar(x, mu, std, capsize=4, label=r"$M_{500c}$")
x, mu, std = csiborgtools.match.number_density(data, "totpartmass", bins, 200, True)
plt.errorbar(x, mu, std, capsize=4, label=r"$M_{\rm tot}$")
x, mu, std = csiborgtools.match.number_density(data, "mass_mmain", bins, 200, True)
plt.errorbar(x, mu, std, capsize=4, label=r"$M_{\rm mmain}$")
plt.legend()
plt.yscale("log")
plt.xscale("log")
plt.ylabel(r"$\phi / (\mathrm{Mpc}^{-3})~\mathrm{dex}$")
plt.xlabel(r"$M_{\rm x}$")
plt.tight_layout()
plt.savefig("../plots/HMF.png", dpi=450)
plt.show()
In [ ]:
nfw = csiborgtools.fits.NFWProfile()
m200_nfw = nfw.enclosed_mass(data["r200"], data["Rs"], data["rho0"])
plt.figure()
plt.scatter(data["m200"], m200_nfw, s=1)
t = np.linspace(1e11, 1e15)
plt.plot(t, t, c="red", ls="--", lw=1.5)
plt.xscale("log")
plt.yscale("log")
plt.xlabel(r"$M_{200c}$")
plt.ylabel(r"$M_{\mathrm{NFW}, 200c}$")
plt.tight_layout()
plt.savefig("../plots/enclosed_vs_nfw.png", dpi=450)
plt.show()
In [ ]:
In [ ]:
logm200 = np.log10(data["m200"])
conc = data["conc"]
N = 10
bins = np.linspace(logm200.min(), logm200.max(), N)
x = [0.5*(bins[i] + bins[i + 1]) for i in range(N-1)]
y = np.full((N - 1, 3), np.nan)
for i in range(N - 1):
mask = (logm200 >= bins[i]) & (logm200 < bins[i + 1]) & np.isfinite(conc)
y[i, :] = np.percentile(conc[mask], [14, 50, 84])
fig, axs = plt.subplots(nrows=2, sharex=True, figsize=(6.4, 6.4 * 1))
fig.subplots_adjust(hspace=0)
axs[0].plot(x, y[:, 1], c="red", marker="o")
axs[0].fill_between(x, y[:, 0], y[:, 2], color="red", alpha=0.25)
axs[1].hist(logm200, bins="auto", log=True)
for b in bins:
for i in range(2):
axs[i].axvline(b, c="orange", lw=0.5)
axs[0].set_ylim(2, 10)
axs[1].set_xlabel(r"$M_{200c}$")
axs[0].set_ylabel(r"$c_{200c}$")
axs[1].set_ylabel(r"Counts")
plt.tight_layout(h_pad=0)
plt.savefig("../plots/mass_concentration.png", dpi=450)
plt.show()
In [ ]: