Why Galaxies Don’t Fly Apart: Modeling M101 with Python

In this post, I use Python to model the orbital mechanics of stars within a galaxy (using M-101 as a template). By defining a custom velocity function and using calculus to find the peak velocity point, we can visualize how star speeds change with distance and calculate the billions of years required for a single galactic orbit.

Stars inside a galaxy continuously orbit to prevent falling into the center of galaxy.  The galaxy has a radius of about 90000 light years.

We used numpy arrays to hold distances, speeds of stars. The data was also plotted using matplotlib. 

Here is the code

import numpy as np
import matplotlib.pyplot as plt

# Constants
LY_TO_KM = 9.461e12  # More precise Light Year to KM
YR_TO_SEC = 365.25 * 24 * 60 * 60  # Include leap year fraction

def vel_star(x_frac):
    """
    Models orbital speed v(x) = 350x / (1 + x^2)^(3/4)
    x_frac: distance in units of 10,000 ly
    """
    return (350 * x_frac) / np.power((1 + x_frac**2), 0.75)

def orbit_time(r_km, v_kms):
    """Calculates orbital period in seconds."""
    return (2 * np.pi * r_km) / v_kms

def format_time(t_sec):
    """Converts seconds into the most readable astronomical unit."""
    if t_sec >= YR_TO_SEC * 1e9:
        return t_sec / (YR_TO_SEC * 1e9), "By"
    elif t_sec >= YR_TO_SEC * 1e6:
        return t_sec / (YR_TO_SEC * 1e6), "My"
    elif t_sec >= YR_TO_SEC:
        return t_sec / YR_TO_SEC, "yr"
    else:
        return t_sec / (24 * 3600), "days"

# Galaxy Parameters
galaxy = "M-101"
radius_ly = 90000
max_dis_x = radius_ly // 10000

print(f"\nAnalysis of Galaxy: {galaxy}")
print("=" * 45)

# Calculus: Find peak velocity
# The derivative v'(x) = 0 at x = sqrt(2)
max_vel_pt = np.sqrt(2)
max_vel = vel_star(max_vel_pt)

print(f"Velocity Peak at: {max_vel_pt * 10000:,.0f} ly")
print(f"Maximum Velocity: {max_vel:.2f} km/s")
print("-" * 45)

# Calculations for discrete points
dist_x = np.arange(1, max_dis_x + 1)
actual_dist_ly = dist_x * 10000
act_dist_km = actual_dist_ly * LY_TO_KM
velocities = vel_star(dist_x)
times_sec = orbit_time(act_dist_km, velocities)

# Display Table
print(f"{'Dist (ly)':<12 45="" actual_dist_ly="" d:="" d="" el="" f="" for="" in="" km="" period="" print="" rbit="" s="" t="" times_sec="" unit="format_time(t)" v:="" v="" val:="" val="" velocities="" zip="">8.2f} {unit}")
print("=" * 45)

# Visualization
plot_x = np.linspace(0.1, max_dis_x, 200)
plot_y = vel_star(plot_x)

plt.figure(figsize=(10, 6))
plt.plot(plot_x, plot_y, color="indigo", linewidth=2, label="Orbital Velocity Curve")
plt.axvline(max_vel_pt, color="crimson", linestyle="--", alpha=0.7, 
            label=f"Max Velocity Radius (x={max_vel_pt:.2f})")
plt.scatter(max_vel_pt, max_vel, color="crimson") # Mark the peak

plt.title(f'Star Orbital Speeds in {galaxy}', fontsize=14)
plt.xlabel('Distance from Center ($10,000$ ly)', fontsize=12)
plt.ylabel('Orbital Speed (km/s)', fontsize=12)
plt.grid(True, linestyle=':', alpha=0.6)
plt.legend()
plt.show()