{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Copyright 2025 Marc Kamionkowski and Kris Sigurdson.\n",
"Licensed under Apache-2.0 (http://www.apache.org/licenses/LICENSE-2.0).\n",
"\n",
"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.\n",
"\n",
"You must have successfully installed NSphere by following the instructions in the main [README](https://github.com/kris-sigurdson/NSphere-Development). Ensure the `nsphere` Python virtual environment is activated (`source ./activate_nsphere` from the project root).\n",
"\n",
"This notebook requires output data from an NSphere simulation run created by running the following command from the project root:\n",
"```bash\n",
"./nsphere --save all\n",
"```\n",
"\n",
"To visualize the tidal stripping scenario shown in the pre-rendered animations when first loading this notebook, run:\n",
"```bash\n",
"./nsphere --save all --ftidal 0.4\n",
"```\n",
"and then rerun this notebook.\n",
"\n",
"To save animations, uncomment the `ani.save(...)` lines in the code cells (they will save to `nsphere/results/`).\n",
"\n",
"For more information on NSphere, see the [full documentation](https://kris-sigurdson.github.io/NSphere-Development/).\n",
"\n",
"*(This notebook was last tested with Python 3.12.x)*"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Histogram data shape: (100, 59, 60)\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"### THIS CELL MAKES AN ANIMATION OF THE EVOLUTION OF THE RADIUS-RADIAL VELOCITY PHASE SPACE\n",
"\n",
"\n",
"\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.animation import FuncAnimation\n",
"from scipy.ndimage import gaussian_filter # 2D smoothing filter\n",
"from IPython.display import HTML\n",
"\n",
"# Define the path to the data files and the time steps\n",
"# Adjust the path to your data files as necessary\n",
"path_template = \"../data/Rank_Mass_Rad_VRad_unsorted_t{timestep:05d}_100000_10001_5.dat\"\n",
"time_steps = range(0, 100, 1) # Time steps from 0 to 99\n",
"\n",
"# Define bin edges for radius (R) and radial velocity (Vrad)\n",
"# Adjust these ranges as needed for your data.\n",
"bin_edges_R = np.arange(0, 300, 5) # For radius: 0 to 300 in steps of 5\n",
"bin_edges_V = np.arange(-300, 310, 10) # For Vrad: -300 to 300 in steps of 10\n",
"\n",
"# Define the record dtype (as in nsphere.c; See project documentation.)\n",
"record_dtype = np.dtype([\n",
" ('rank', np.int32),\n",
" ('mass', np.float32),\n",
" ('R', np.float32),\n",
" ('Vrad', np.float32),\n",
" ('PsiA', np.float32),\n",
" ('E', np.float32),\n",
" ('L', np.float32)\n",
"])\n",
"\n",
"# Storage for 2D histograms\n",
"all_histograms = []\n",
"\n",
"# Loop over time steps and process each file\n",
"for time_step in time_steps:\n",
" file_path = path_template.format(timestep=time_step)\n",
" try:\n",
" # Read the binary file\n",
" alldata = np.fromfile(file_path, dtype=record_dtype)\n",
" # Extract radius and radial velocity data\n",
" data_R = alldata['R']\n",
" data_Vrad = alldata['Vrad']/1.023e-3 # Convert to km/s\n",
" \n",
" if np.isnan(data_R).any() or np.isnan(data_Vrad).any():\n",
" print(f\"Warning: NaN values found in file at time step {time_step}\")\n",
" if not (np.isfinite(data_R).all() and np.isfinite(data_Vrad).all()):\n",
" print(f\"Warning: Non-finite values found in file at time step {time_step}\")\n",
" \n",
" # Compute 2D histogram: counts for each (R, Vrad) bin\n",
" hist, _, _ = np.histogram2d(data_R, data_Vrad, bins=[bin_edges_R, bin_edges_V])\n",
" \n",
" # Apply 2D Gaussian smoothing to reduce noise\n",
" # Sigma value of 0.1 provides light smoothing while preserving structure\n",
" smoothed_hist = gaussian_filter(hist, sigma=(0.1, 0.1))\n",
" \n",
" all_histograms.append(smoothed_hist)\n",
" except Exception as e:\n",
" print(f\"Error reading file {file_path}: {e}\")\n",
" # Append a zero histogram if file is missing or cannot be read\n",
" all_histograms.append(np.zeros((len(bin_edges_R)-1, len(bin_edges_V)-1)))\n",
"\n",
"# Convert list to numpy array for easier handling\n",
"all_histograms = np.array(all_histograms)\n",
"print(\"Histogram data shape:\", all_histograms.shape) # (time_steps, n_R_bins, n_V_bins)\n",
"\n",
"# Determine common color scale for the animation\n",
"vmin = 0\n",
"vmax = np.max(all_histograms) # Use the maximum value across all time steps\n",
"\n",
"# Create the animation\n",
"fig, ax = plt.subplots()\n",
"# Transpose the histogram so that the horizontal axis corresponds to R and vertical to Vrad\n",
"im = ax.imshow(all_histograms[0].T, origin='lower',\n",
" extent=[bin_edges_R[0], bin_edges_R[-1], bin_edges_V[0], bin_edges_V[-1]],\n",
" aspect='auto', cmap='viridis', vmin=vmin, vmax=vmax)\n",
"\n",
"ax.set_xlabel(\"Radius\")\n",
"ax.set_ylabel(\"Radial Velocity\")\n",
"ax.set_title(f\"Time Step: {time_steps[0]}\")\n",
"fig.colorbar(im, ax=ax, label=\"Count\")\n",
"\n",
"def update(frame):\n",
" im.set_data(all_histograms[frame].T) # update image data\n",
" ax.set_title(f\"Time Step: {time_steps[frame]}\")\n",
" return im,\n",
"\n",
"ani = FuncAnimation(fig, update, frames=len(time_steps), repeat=True, interval=200)\n",
"\n",
"# Uncomment this line to show a static plot of the initial conditions\n",
"#plt.show()\n",
"\n",
"plt.close(fig)\n",
"\n",
"# Uncomment this line to save the animation as a GIF\n",
"# ani.save(\"../results/nsphere_phasespace_example.gif\", writer=\"pillow\", fps=10)\n",
"\n",
"# For displaying in a Jupyter Notebook:\n",
"HTML(ani.to_jshtml())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Visualization of NSphere Phase Space and Density Evolution\n",
"\n",
"The animations above and below show two different views of the same N-body simulation:\n",
"\n",
"1. **Phase Space Animation (above)**: \n",
" 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.\n",
" \n",
" 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.\n",
"\n",
"2. **Density Profile Animation (below)**:\n",
" 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.\n",
" \n",
" 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."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# This code creates an animation of the evolution of the radius distribution\n",
"# and displays it in a Jupyter Notebook.\n",
"\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.animation import FuncAnimation\n",
"from scipy.ndimage import gaussian_filter1d\n",
"from IPython.display import HTML\n",
"\n",
"path_template = \"../data/Rank_Mass_Rad_VRad_unsorted_t{timestep:05d}_100000_10001_5.dat\"\n",
"time_steps = range(0, 100, 1) # Time steps from 0 to 99\n",
"\n",
"# Using a larger range for the radius distribution plot compared to the phase space plot\n",
"# 500 kpc captures more distant particles that might be present in the outer regions\n",
"bin_edges = np.arange(0, 500, 5) # For radius: 0 to 500 in steps of 5\n",
"bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2 # Calculate bin centers for plotting\n",
"\n",
"# Storage for smoothed histograms\n",
"all_smoothed_counts = []\n",
"\n",
"# Same data structure definition as in the first cell\n",
"record_dtype = np.dtype([\n",
" ('rank', np.int32),\n",
" ('mass', np.float32),\n",
" ('R', np.float32),\n",
" ('Vrad', np.float32),\n",
" ('PsiA', np.float32),\n",
" ('E', np.float32),\n",
" ('L', np.float32)\n",
"])\n",
"\n",
"# Process each file\n",
"for time_step in time_steps:\n",
" file_path = path_template.format(timestep=time_step)\n",
" try:\n",
" # Read the binary file and extract radius data\n",
" alldata = np.fromfile(file_path, dtype=record_dtype)\n",
" data = alldata['R']\n",
"\n",
" if np.isnan(data).any():\n",
" print(f\"Warning: NaN values found in data at Time Step: {time_step}\")\n",
" if not np.isfinite(data).all():\n",
" print(f\"Warning: Infinite values found in data at Time Step: {time_step}\")\n",
" \n",
" # Count particles in bins\n",
" counts, _ = np.histogram(data, bins=bin_edges)\n",
" \n",
" # Apply 1D Gaussian smoothing to reduce noise in the histogram\n",
" # Sigma=1 gives gentle smoothing that preserves important features\n",
" smoothed_counts = gaussian_filter1d(counts, sigma=1)\n",
" \n",
" # Note: Uncomment the line below to normalize by r² (number density) instead of showing total counts\n",
" # This would convert from dM/dr (mass per radius) to ρ(r) (mass density)\n",
" # smoothed_counts = smoothed_counts / bin_centers**2\n",
" \n",
" all_smoothed_counts.append(smoothed_counts)\n",
" except Exception as e:\n",
" print(f\"Error reading file {file_path}: {e}\")\n",
" all_smoothed_counts.append(np.zeros(len(bin_centers))) # Append zeros if file is missing\n",
"\n",
"# Convert to numpy array for easier handling\n",
"all_smoothed_counts = np.array(all_smoothed_counts)\n",
"\n",
"# Create animation\n",
"fig, ax = plt.subplots()\n",
"line, = ax.plot(bin_centers, all_smoothed_counts[0], color=\"blue\")\n",
"\n",
"ax.set_xlim(0, 500) # Lower limit is > 0\n",
"# Compute the global maximum of all frames\n",
"global_max = np.max(all_smoothed_counts)\n",
"# Set the y-axis limit to be 5% higher than this maximum for better visualization\n",
"ax.set_ylim(0, global_max * 1.05)\n",
"\n",
"ax.set_xlabel(\"Radius\")\n",
"ax.set_ylabel(\"Number of Particles\")\n",
"ax.set_title(\"Evolution of Particle Distribution\")\n",
"\n",
"def update(frame):\n",
" line.set_ydata(all_smoothed_counts[frame]) # Update the y-data for the line\n",
" ax.set_title(f\"Time Step: {time_steps[frame]}\")\n",
" return line,\n",
"\n",
"ani = FuncAnimation(fig, update, frames=len(time_steps), repeat=True, interval=200)\n",
"\n",
"# Uncomment this line to show a static plot of the initial conditions\n",
"#plt.show()\n",
"\n",
"plt.close(fig)\n",
"\n",
"# Uncomment this line to save the animation as a GIF\n",
"# ani.save(\"../results/nsphere_density_example.gif\", writer=\"pillow\", fps=10)\n",
"\n",
"HTML(ani.to_jshtml())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"apache_2_0_license_text": "Copyright 2025 Marc Kamionkowski and Kris Sigurdson.\n\nLicensed under the Apache License, Version 2.0 (the \"License\");\nyou may not use this file except in compliance with the License.\nYou may obtain a copy of the License at\n\n http://www.apache.org/licenses/LICENSE-2.0\n\nUnless required by applicable law or agreed to in writing, software\ndistributed under the License is distributed on an \"AS IS\" BASIS,\nWITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\nSee the License for the specific language governing permissions and\nlimitations under the License.",
"copyright": "Copyright 2025 Marc Kamionkowski and Kris Sigurdson. Licensed under Apache-2.0",
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.3"
},
"license": "Apache-2.0"
},
"nbformat": 4,
"nbformat_minor": 4
}