Galaxy Cluster Simulation

import matplotlib.pyplot as plt
import numpy as np

np.random.seed(42)

# Generate galaxy cluster data
n_galaxies = 300

# Central cluster
r1 = np.random.exponential(1, n_galaxies // 2)
theta1 = np.random.uniform(0, 2*np.pi, n_galaxies // 2)
phi1 = np.random.uniform(0, np.pi, n_galaxies // 2)
x1 = r1 * np.sin(phi1) * np.cos(theta1)
y1 = r1 * np.sin(phi1) * np.sin(theta1)
z1 = r1 * np.cos(phi1)

# Satellite cluster
r2 = np.random.exponential(0.5, n_galaxies // 2)
x2 = r2 * np.sin(np.random.uniform(0, np.pi, n_galaxies // 2)) * np.cos(np.random.uniform(0, 2*np.pi, n_galaxies // 2)) + 3
y2 = r2 * np.sin(np.random.uniform(0, np.pi, n_galaxies // 2)) * np.sin(np.random.uniform(0, 2*np.pi, n_galaxies // 2)) + 2
z2 = r2 * np.cos(np.random.uniform(0, np.pi, n_galaxies // 2)) + 1

x = np.concatenate([x1, x2])
y = np.concatenate([y1, y2])
z = np.concatenate([z1, z2])

# Size based on "mass"
sizes = np.random.exponential(30, n_galaxies)

fig = plt.figure(figsize=(10, 8), facecolor='#0a0a0f')
ax = fig.add_subplot(111, projection='3d', facecolor='#0a0a0f')

scatter = ax.scatter(x, y, z, c=np.sqrt(x**2 + y**2 + z**2), 
                     cmap='cool', s=sizes, alpha=0.6, edgecolors='none')

ax.set_xlabel('X (Mpc)', color='white', fontsize=10)
ax.set_ylabel('Y (Mpc)', color='white', fontsize=10)
ax.set_zlabel('Z (Mpc)', color='white', fontsize=10)
ax.set_title('Galaxy Cluster Simulation', color='white', fontsize=14, fontweight='bold', pad=20)

ax.tick_params(colors='#64748b', labelsize=8)
ax.xaxis.pane.fill = False
ax.yaxis.pane.fill = False
ax.zaxis.pane.fill = False
ax.xaxis.pane.set_edgecolor('#1e293b')
ax.yaxis.pane.set_edgecolor('#1e293b')
ax.zaxis.pane.set_edgecolor('#1e293b')

ax.view_init(elev=20, azim=45)
plt.tight_layout()
plt.show()