Note

This example notebook demonstrates NSphere visualization techniques. The original notebook can be found in the examples directory of the NSphere project.

Copyright 2025 Marc Kamionkowski and Kris Sigurdson. Licensed under Apache-2.0 (http://www.apache.org/licenses/LICENSE-2.0).

This notebook demonstrates how to load, visualize, and create animations from the output data generated by the main NSphere simulation (nsphere). It focuses on plotting the phase-space evolution and the radial density profile.

You must have successfully installed NSphere by following the instructions in the main README. Ensure the nsphere Python virtual environment is activated (source ./activate_nsphere from the project root).

This notebook requires output data from an NSphere simulation run created by running the following command from the project root:

./nsphere --save all

To visualize the tidal stripping scenario shown in the pre-rendered animations when first loading this notebook, run:

./nsphere --save all --ftidal 0.4

and then rerun this notebook.

To save animations, uncomment the ani.save(...) lines in the code cells (they will save to nsphere/results/).

For more information on NSphere, see the full documentation.

(This notebook was last tested with Python 3.12.x)

### THIS CELL MAKES AN ANIMATION OF THE EVOLUTION OF THE RADIUS-RADIAL VELOCITY PHASE SPACE



import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.ndimage import gaussian_filter  # 2D smoothing filter
from IPython.display import HTML

# Define the path to the data files and the time steps
# Adjust the path to your data files as necessary
path_template = "../data/Rank_Mass_Rad_VRad_unsorted_t{timestep:05d}_100000_10001_5.dat"
time_steps = range(0, 100, 1)  # Time steps from 0 to 99

# Define bin edges for radius (R) and radial velocity (Vrad)
# Adjust these ranges as needed for your data.
bin_edges_R = np.arange(0, 300, 5)      # For radius: 0 to 300 in steps of 5
bin_edges_V = np.arange(-300, 310, 10)     # For Vrad: -300 to 300 in steps of 10

# Define the record dtype (as in nsphere.c; See project documentation.)
record_dtype = np.dtype([
    ('rank',  np.int32),
    ('mass',  np.float32),
    ('R',     np.float32),
    ('Vrad',  np.float32),
    ('PsiA',  np.float32),
    ('E',     np.float32),
    ('L',     np.float32)
])

# Storage for 2D histograms
all_histograms = []

# Loop over time steps and process each file
for time_step in time_steps:
    file_path = path_template.format(timestep=time_step)
    try:
        # Read the binary file
        alldata = np.fromfile(file_path, dtype=record_dtype)
        # Extract radius and radial velocity data
        data_R = alldata['R']
        data_Vrad = alldata['Vrad']/1.023e-3  # Convert to km/s

        if np.isnan(data_R).any() or np.isnan(data_Vrad).any():
            print(f"Warning: NaN values found in file at time step {time_step}")
        if not (np.isfinite(data_R).all() and np.isfinite(data_Vrad).all()):
            print(f"Warning: Non-finite values found in file at time step {time_step}")

        # Compute 2D histogram: counts for each (R, Vrad) bin
        hist, _, _ = np.histogram2d(data_R, data_Vrad, bins=[bin_edges_R, bin_edges_V])

        # Apply 2D Gaussian smoothing to reduce noise
        # Sigma value of 0.1 provides light smoothing while preserving structure
        smoothed_hist = gaussian_filter(hist, sigma=(0.1, 0.1))

        all_histograms.append(smoothed_hist)
    except Exception as e:
        print(f"Error reading file {file_path}: {e}")
        # Append a zero histogram if file is missing or cannot be read
        all_histograms.append(np.zeros((len(bin_edges_R)-1, len(bin_edges_V)-1)))

# Convert list to numpy array for easier handling
all_histograms = np.array(all_histograms)
print("Histogram data shape:", all_histograms.shape)  # (time_steps, n_R_bins, n_V_bins)

# Determine common color scale for the animation
vmin = 0
vmax = np.max(all_histograms)  # Use the maximum value across all time steps

# Create the animation
fig, ax = plt.subplots()
# Transpose the histogram so that the horizontal axis corresponds to R and vertical to Vrad
im = ax.imshow(all_histograms[0].T, origin='lower',
               extent=[bin_edges_R[0], bin_edges_R[-1], bin_edges_V[0], bin_edges_V[-1]],
               aspect='auto', cmap='viridis', vmin=vmin, vmax=vmax)

ax.set_xlabel("Radius")
ax.set_ylabel("Radial Velocity")
ax.set_title(f"Time Step: {time_steps[0]}")
fig.colorbar(im, ax=ax, label="Count")

def update(frame):
    im.set_data(all_histograms[frame].T)  # update image data
    ax.set_title(f"Time Step: {time_steps[frame]}")
    return im,

ani = FuncAnimation(fig, update, frames=len(time_steps), repeat=True, interval=200)

# Uncomment this line to show a static plot of the initial conditions
#plt.show()

plt.close(fig)

# Uncomment this line to save the animation as a GIF
# ani.save("../results/nsphere_phasespace_example.gif", writer="pillow", fps=10)

# For displaying in a Jupyter Notebook:
HTML(ani.to_jshtml())

Histogram data shape: (100, 59, 60)

Once Loop Reflect

Visualization of NSphere Phase Space and Density Evolution

The animations above and below show two different views of the same N-body simulation:

1. Phase Space Animation (above): This animation shows the evolution of the radius-radial velocity phase space over time. The horizontal axis represents the particle radius (distance from center), while the vertical axis shows the radial velocity. The color intensity indicates the particle density in each region of phase space.

   A stable self-gravitating system will show coherent patterns in this phase space, typically maintaining its overall structure. In contrast, a disturbed system (like one undergoing tidal stripping) will show more dramatic evolution.

2. Density Profile Animation (below): This animation shows how the radial density distribution \(\frac{dM}{dr} = 4\pi r^2 \rho\) evolves over time. The horizontal axis represents the radius, and the vertical axis shows the particle count at each radius.

   This view makes it easier to see how the overall structure of the halo changes, particularly whether particles are moving inward or outward over time.

# This code creates an animation of the evolution of the radius distribution
# and displays it in a Jupyter Notebook.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.ndimage import gaussian_filter1d
from IPython.display import HTML

path_template = "../data/Rank_Mass_Rad_VRad_unsorted_t{timestep:05d}_100000_10001_5.dat"
time_steps = range(0, 100, 1)  # Time steps from 0 to 99

# Using a larger range for the radius distribution plot compared to the phase space plot
# 500 kpc captures more distant particles that might be present in the outer regions
bin_edges = np.arange(0, 500, 5)      # For radius: 0 to 500 in steps of 5
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2  # Calculate bin centers for plotting

# Storage for smoothed histograms
all_smoothed_counts = []

# Same data structure definition as in the first cell
record_dtype = np.dtype([
    ('rank',  np.int32),
    ('mass',  np.float32),
    ('R',     np.float32),
    ('Vrad',  np.float32),
    ('PsiA',  np.float32),
    ('E',     np.float32),
    ('L',     np.float32)
])

# Process each file
for time_step in time_steps:
    file_path = path_template.format(timestep=time_step)
    try:
        # Read the binary file and extract radius data
        alldata = np.fromfile(file_path, dtype=record_dtype)
        data = alldata['R']

        if np.isnan(data).any():
            print(f"Warning: NaN values found in data at Time Step: {time_step}")
        if not np.isfinite(data).all():
            print(f"Warning: Infinite values found in data at Time Step: {time_step}")

        # Count particles in bins
        counts, _ = np.histogram(data, bins=bin_edges)

        # Apply 1D Gaussian smoothing to reduce noise in the histogram
        # Sigma=1 gives gentle smoothing that preserves important features
        smoothed_counts = gaussian_filter1d(counts, sigma=1)

        # Note: Uncomment the line below to normalize by r² (number density) instead of showing total counts
        # This would convert from dM/dr (mass per radius) to ρ(r) (mass density)
        # smoothed_counts = smoothed_counts / bin_centers**2

        all_smoothed_counts.append(smoothed_counts)
    except Exception as e:
        print(f"Error reading file {file_path}: {e}")
        all_smoothed_counts.append(np.zeros(len(bin_centers)))  # Append zeros if file is missing

# Convert to numpy array for easier handling
all_smoothed_counts = np.array(all_smoothed_counts)

# Create animation
fig, ax = plt.subplots()
line, = ax.plot(bin_centers, all_smoothed_counts[0], color="blue")

ax.set_xlim(0, 500)      # Lower limit is > 0
# Compute the global maximum of all frames
global_max = np.max(all_smoothed_counts)
# Set the y-axis limit to be 5% higher than this maximum for better visualization
ax.set_ylim(0, global_max * 1.05)

ax.set_xlabel("Radius")
ax.set_ylabel("Number of Particles")
ax.set_title("Evolution of Particle Distribution")

def update(frame):
    line.set_ydata(all_smoothed_counts[frame])  # Update the y-data for the line
    ax.set_title(f"Time Step: {time_steps[frame]}")
    return line,

ani = FuncAnimation(fig, update, frames=len(time_steps), repeat=True, interval=200)

# Uncomment this line to show a static plot of the initial conditions
#plt.show()

plt.close(fig)

# Uncomment this line to save the animation as a GIF
# ani.save("../results/nsphere_density_example.gif", writer="pillow", fps=10)

HTML(ani.to_jshtml())

Once Loop Reflect