2022-05-05 10:40:25 +02:00
|
|
|
import matplotlib.pyplot as plt
|
|
|
|
import numpy as np
|
|
|
|
import pandas as pd
|
|
|
|
from matplotlib.axes import Axes
|
|
|
|
from matplotlib.figure import Figure
|
|
|
|
|
2022-05-05 11:10:07 +02:00
|
|
|
from readfiles import ParticlesMeta
|
|
|
|
|
2022-05-05 10:40:25 +02:00
|
|
|
|
|
|
|
def V(r):
|
|
|
|
return 4 * np.pi * r ** 3 / 3
|
|
|
|
|
|
|
|
|
2022-05-05 11:12:37 +02:00
|
|
|
def cumulative_mass_profile(particles_in_halos: pd.DataFrame, halo: pd.Series,
|
|
|
|
particles_meta: ParticlesMeta, plot=False):
|
2022-05-05 10:40:25 +02:00
|
|
|
print(type(particles_in_halos))
|
|
|
|
centre = np.array([halo.X, halo.Y, halo.Z])
|
|
|
|
positions = particles_in_halos[["X", "Y", "Z"]].to_numpy()
|
|
|
|
print(positions)
|
|
|
|
print(positions.shape)
|
|
|
|
distances = np.linalg.norm(positions - centre, axis=1)
|
|
|
|
group_radius = distances.max()
|
|
|
|
normalized_distances = distances / group_radius
|
|
|
|
|
|
|
|
num_bins = 30
|
|
|
|
log_radial_bins = np.geomspace(0.01, 2, num_bins)
|
|
|
|
|
2022-05-05 11:10:07 +02:00
|
|
|
bin_masses = []
|
2022-05-05 10:40:25 +02:00
|
|
|
bin_densities = []
|
|
|
|
for k in range(num_bins - 1):
|
|
|
|
bin_start = log_radial_bins[k]
|
|
|
|
bin_end = log_radial_bins[k + 1]
|
|
|
|
in_bin = np.where((bin_start < normalized_distances) & (normalized_distances < bin_end))[0]
|
|
|
|
count = in_bin.shape[0]
|
2022-05-05 11:10:07 +02:00
|
|
|
mass = count * particles_meta.particle_mass
|
2022-05-05 10:40:25 +02:00
|
|
|
volume = V(bin_end * group_radius) - V(bin_start * group_radius)
|
2022-05-05 11:10:07 +02:00
|
|
|
bin_masses.append(mass)
|
2022-05-05 10:40:25 +02:00
|
|
|
density = count / volume
|
|
|
|
bin_densities.append(density)
|
2022-05-05 11:10:07 +02:00
|
|
|
print(bin_masses)
|
2022-05-05 10:40:25 +02:00
|
|
|
print(bin_densities)
|
|
|
|
|
2022-05-05 11:12:37 +02:00
|
|
|
if plot:
|
|
|
|
fig: Figure = plt.figure()
|
|
|
|
ax: Axes = fig.gca()
|
2022-05-05 10:40:25 +02:00
|
|
|
|
2022-05-05 11:12:37 +02:00
|
|
|
ax2 = ax.twinx()
|
2022-05-05 10:40:25 +02:00
|
|
|
|
2022-05-05 11:12:37 +02:00
|
|
|
ax.loglog(log_radial_bins[:-1], bin_masses, label="counts")
|
|
|
|
ax2.loglog(log_radial_bins[:-1], bin_densities, label="densities", c="C1")
|
|
|
|
ax.set_xlabel(r'R / R$_\mathrm{group}$')
|
|
|
|
ax.set_ylabel(r'M [$10^{10} \mathrm{M}_\odot$]')
|
|
|
|
ax2.set_ylabel("density [$\\frac{10^{10} \mathrm{M}_\odot}{Mpc^3}$]")
|
2022-05-05 10:40:25 +02:00
|
|
|
|
2022-05-05 11:12:37 +02:00
|
|
|
plt.show()
|
2022-05-05 10:40:25 +02:00
|
|
|
|
2022-05-05 11:10:07 +02:00
|
|
|
return bin_masses, bin_densities
|