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()

Speeding Galaxies

Introduction 

Universe is not a static thing. It is dynamic.  We all know galaxies are still moving away as the universe is expanding.  

We took data provided by NASA for educational purposes for 7 galaxies.  It consisted of names, distances,and speeds of the galaxies. 

The data was plotted on XY plot and linear regression analysis were done. Plotting was done using numpy library and matplotlib. 

Coding was done in python language. 

One obvious conclusion from the plots is that more distant galaxies are moving with more speed. 

import numpy as np
import matplotlib.pyplot as plt

# 1. Constants & Data
ONE_MPC_LY = 3.26 
galaxies = {
    "NGC-5357": (0.45, 200), "NGC-3627": (0.9, 650), "NGC-5236": (0.9, 500),
    "NGC-4151": (1.7, 960),  "NGC-4472": (2.0, 850), "NGC-4486": (2.0, 800),
    "NGC-4649": (2.0, 1090),
}

# Convert data for calculation
dist_mpc = np.array([d[0] for d in galaxies.values()])
speed_kms = np.array([d[1] for d in galaxies.values()])

# 2. Linear Regression & R-Squared
a, b = np.polyfit(dist_mpc, speed_kms, 1)
r_squared = np.corrcoef(dist_mpc, speed_kms)[0,1]**2

# 3. NASA Extrapolation (2500 km/s)
speed_nasa = 2500
dist_nasa = (speed_nasa - b) / a

# 4. Unit Conversions (Using min/max distances)
d_min, d_max = dist_mpc.min(), dist_mpc.max()
print(f"Distance Range: {d_min*ONE_MPC_LY:.2f} to {d_max*ONE_MPC_LY:.2f} Million ly")
print(f"Note: There are 3 galaxies at the max distance of {d_max} mpc.")

# 5. Plotting
plt.figure(figsize=(10, 6))

# Plot the 3 overlapping points and the others
plt.scatter(dist_mpc, speed_kms, color='red', label='Galaxy Data', zorder=3)

# Regression Line (Extended)
x_fit = np.linspace(0, dist_nasa + 1, 100)
plt.plot(x_fit, a*x_fit + b, color='blue', alpha=0.8, label=f'Fit: $R^2$={r_squared:.2f}')

# NASA Vertical Line & Annotation
plt.plot([dist_nasa, dist_nasa], [0, speed_nasa], 'g--', alpha=0.5)
plt.plot([0, dist_nasa], [speed_nasa, speed_nasa], 'g--', alpha=0.5)
plt.annotate(f'NASA Target: {dist_nasa:.2f} mpc', (dist_nasa, speed_nasa), 
             textcoords="offset points", xytext=(-15, 10), ha='center', color='green')

# Scale and Labels
plt.xlim(0, dist_nasa + 1)
plt.ylim(0, speed_nasa + 500)
plt.xlabel('Distance (mpc)')
plt.ylabel('Speed (km/s)')
plt.title('Hubble Law Regression with NASA Extrapolation')
plt.legend()
plt.grid(True, linestyle=':', alpha=0.6)

plt.show()

Hope yoy will like the effort.

The Chemistry and Code of Sulfur Recovery

In industrial chemistry, removing sulfur from "sour gas" (natural gas containing H_2S) is a critical safety and environmental task. One of the most elegant ways to do this is the Claus Process. Today, we’ll use Python to model this reaction, balancing the stoichiometry and calculating the physical volumes of our products.

The Reaction

​We are looking at the reaction where sulfur dioxide (SO2)—often from coal burning—is neutralized by hydrogen sulfide (H2S):

SO2 +H2S -> 3S + 2H2O

The Python Script

"""
Project: The 70-Year-Old Coder
Topic: Stoichiometry & Physics of the Claus Process
"""

# 1. Chemical Constants
atwts = {"S": 32.06, "O": 15.999, "H": 1.008}
stoichiometry = {"SO2": 1, "H2S": 2, "S": 3, "H2O": 2}

def mol_weight(chem):
    return sum([v * atwts.get(c.upper()) for c, v in chem.items()])

species_data = {
    "SO2": {"s": 1, "o": 2},
    "H2S": {"h": 2, "s": 1},
    "S":   {"s": 1},
    "H2O": {"h": 2, "o": 1}
}

molecular_wts = {name: mol_weight(comp) for name, comp in species_data.items()}

# 2. Input Parameters
wt_coal_tons = 2
pct_sulfur = 3.5
temp_c = 27
pressure_atm = 1

# 3. Mass & Gas Calculations
sulfur_in_coal_gm = wt_coal_tons * 1_000_000 * (pct_sulfur / 100)
so2_formed = (molecular_wts['SO2'] / atwts['S']) * sulfur_in_coal_gm

# Molar ratio calculation for H2S required
h2s_req = (stoichiometry['H2S'] * molecular_wts['H2S']) / \
          (stoichiometry['SO2'] * molecular_wts['SO2']) * so2_formed

# Volume using Charles's Law (v1/t1 = v2/t2) at 1 atm
R = 0.08206  
temp_k = temp_c + 273.15
moles_h2s = h2s_req / molecular_wts['H2S']
vol_h2s_liters = (moles_h2s * R * temp_k) / pressure_atm

# 4. Physical Recovery (Density Logic)
s_recovered_gm = (stoichiometry['S'] * molecular_wts['S']) / \
                 (stoichiometry['SO2'] * molecular_wts['SO2']) * so2_formed
                 
density_s = 2.07  # g/cm3
vol_s_liters = (s_recovered_gm / density_s) / 1000

# --- OUTPUT ---
print(f"--- Results for {wt_coal_tons} Tons of Coal ({pct_sulfur}% S) ---")
print(f"SO2 Produced:      {so2_formed/1000:.2f} kg")
print(f"H2S Required:      {h2s_req/1000:.2f} kg")
print(f"Volume of H2S:     {vol_h2s_liters:.0f} Liters at {temp_c}°C")
print(f"\n--- Physical Recovery ---")
print(f"Sulfur Recovered:  {s_recovered_gm/1000:.2f} kg")
print(f"Sulfur Volume:     {vol_s_liters:.1f} Liters")

The Cosmic Delivery: Did Icy Invaders Build Our Oceans?

Introduction

​For decades, scientists have looked at the vast, blue expanse of our oceans and asked one simple, profound question: Where did all this water come from? While some believe water was "baked" into the rocks of Earth during its formation, a competing theory suggests a much more violent and spectacular origin. The Comet Impact Hypothesis proposes that Earth was once a dry, rocky wasteland until it was bombarded by a relentless "cosmic delivery service"—trillions of tons of pure water-ice arriving via meteors and comet nuclei.

​But does the math actually hold up? If Earth were struck by icy bodies of varying sizes—from small 2 km "scouts" to massive 200 km "titans"—how long would it take to fill the basins of our world? In this post, we use Python to simulate the deposition rates of three distinct types of comet nuclei to see if 400 million years of bombardment is enough to create the world as we know it.

​The Hypothesis: Three Tiers of Impact

​To test this, we categorize our icy visitors into three types based on their diameter and arrival frequency:

Body Type Diameter Arrival Frequency

Type I 2 km Twice per year

Type II 20 km Once every 600 years

Type III 200 km Once every 1,000,000 years

Why 400 Million Years?

​Geologically speaking, the "Late Heavy Bombardment" is a period where impact rates were significantly higher. By calculating the Total Ice Deposition Rate, we can determine if the current volume of the oceans fits within this window of planetary development.

"""
BLOG POST HYPOTHESIS: HOW METEOR IMPACTS CREATED EARTH'S WATER
------------------------------------------------------------
Problem Statement & Parameters:
We explore the hypothesis that Earth’s oceans were filled by icy bodies (comets/meteors).
We assume these bodies consisted of three major spherical types made of pure water-ice:

1. **Type I**:   2 km diameter, arriving once every 6 months.
2. **Type II**: 20 km diameter, arriving once every 600 years.
3. **Type III**: 200 km diameter, arriving every one million years.

Goal: Calculate the volumes and deposition rates of these bodies to determine how long it
would take to create Earth’s current oceans.
"""

import math

# --- 1. Define Core Calculation Function ---
def vol_sphere(radius):
    """
    Calculates the volume of a sphere.
    Formula: (4/3) * pi * r^3
    """
    return (4/3) * math.pi * (radius**3)

# --- 2. Initialize Data and Constants ---

# Dictionary storing raw data for each comet type:
# Structure: "type-key": (diameter_in_km, arrival_interval_in_years)
comet_nuclei_data = {
    "type-I":   (2,   0.5),    # Arrives every 0.5 years (6 months)
    "type-II":  (20,  600),    # Arrives every 600 years
    "type-III": (200, 1e06)    # Arrives every 1,000,000 years
}

# Real-world scientific estimates used as a benchmark:
TOTAL_EARTH_WATER_LIQUID_KM3 = 1.33e09  # Total liquid water volume on Earth in km^3
ICE_TO_WATER_VOL_RATIO = 6              # Assuming 6 units of ice required for 1 unit of liquid water

# --- 3. Run Calculations & Process Data ---

# Create a new dictionary to store expanded results for each type.
# We iterate through the raw data and perform calculations on the fly.
processed_comet_results = {}

for k, (dia, interval) in comet_nuclei_data.items():
    # Calculate radius (needed for volume)
    radius = dia / 2
    
    # PROBLEM 1: What are the volumes of the three types of comet nuclei?
    volume_km3 = vol_sphere(radius)
    
    # Calculate arrival frequency per year
    frequency_per_yr = 1 / interval
    
    # Calculate ice deposition rate per year (Volume * Frequency)
    deposit_rate_km3_yr = volume_km3 * frequency_per_yr
    
    # Store calculated values in the new dictionary
    # Structure: "type-key": (dia, frequency, volume, deposit_rate)
    processed_comet_results[k] = (dia, frequency_per_yr, volume_km3, deposit_rate_km3_yr)


# --- 4. Generate Blog-Worthy Output Report ---

print("\n" + "="*80)
print("             METEOR IMPACT HYPOTHESIS REPORT: EARTH'S WATER ORIGIN")
print("="*80)

print(f"\n[BENCHMARK DATA]")
print(f"Current Estimated Total Liquid Water on Earth: {TOTAL_EARTH_WATER_LIQUID_KM3:,.2f} km³")
print(f"Assumed Ice-to-Liquid Water Expansion Ratio:  {ICE_TO_WATER_VOL_RATIO}:1\n")


print("-" * 73)
print(f"{'Comet Type':<12} | {'Dia (km)':>12} | {'Freq (/yr)':>12} | {'Volume (km³)':>16} | {'Dep. Rate (km³/yr)':>20}")
print("-" * 73)

# Print the processed table rows
for k, (dia, freq, vol, rate) in processed_comet_results.items():
    print(f"{k:<12} | {dia:>12.1f} | {freq:>12.6f} | {vol:>16,.0f} | {rate:>20,.2f}")

print("-" * 73)

# --- 5. Summary and Final Analysis ---

# Calculate Total Ice Deposition Rate by summing the individual rates
total_ice_deposition_rate = sum(item[3] for item in processed_comet_results.values())

# Calculate total ice required to fill current oceans
total_ice_required_for_water_km3 = ICE_TO_WATER_VOL_RATIO * TOTAL_EARTH_WATER_LIQUID_KM3

# Calculate time required (Ice Required / Deposition Rate)
time_required_years = total_ice_required_for_water_km3 / total_ice_deposition_rate

print(f"\n[FINAL ANALYSIS SUMMARY]")
print(f"1. Total Average Ice Deposition Rate:  {total_ice_deposition_rate:,.2f} km³ / year")
print(f"2. Total Ice Required for Current Water: {total_ice_required_for_water_km3:,.2e} km³")
print("-" * 80)
print(f"CONCLUSION: Time required to fill present oceans: {time_required_years/1e06:,.2f} million years")
print("="*80 + "\n")

Physical Properties of Planets Using python

Our solar system is truly amazing. Planets have different composition and different environments. Those near the sun, for example are hotter than Earth and the ones far away are cold. 

In this python code we have calculated many physical parameters using radii, mass, distance from sun, and albedo data. 

The physics used is mostly classical mechanics and thermodynamics.  

Parameters calculated are volumes, densities, orbital speeds, orbital periods, equilibrium temperatures and these properties were compared to Earth. 

We used dictionary with names as keys, and basic parameters are tuples. All calculations were done using functions. Numpy arrays were used to make the work easier. Function argmax was used to tie the names of planet with it's properties. 

Code is given below:

import numpy as np

"""
Solar System Physical & Orbital Dynamics
A comparison of planetary characteristics relative to Earth.
"""

# --- CONSTANTS ---
G = 6.674e-11  # Universal Gravitational Constant (m^3 kg^-1 s^-2)
SOLAR_MASS = 1.989e30  # Mass of the Sun (kg)
METER_IN_KM = 1000.0  # Float conversion to prevent integer overflow errors
SUN_LUMINOSITY = 3.828e26  # Watts
SB_CONST=5.67e-08


# --- PHYSICS FUNCTIONS ---


def force_gravity(mb, r):
    """Calculates gravitational force between the Sun and a body."""
    return G * SOLAR_MASS * mb / r**2


def acc_gravity(mp, r):
    """Calculates acceleration due to gravity at the surface."""
    return G * mp / r**2


def vol_sphere(r):
    """Calculates volume of a sphere."""
    return (4 / 3) * np.pi * r**3


def density(m, v):
    """Calculates density (mass/volume)."""
    return m / v


def esc_vel(gr, r):
    """Calculates escape velocity."""
    return np.sqrt(2 * gr * r)


def orbital_speed(r):
    """Calculates the speed required for a stable circular orbit around the Sun."""
    return np.sqrt(G * SOLAR_MASS / r)


def orbit_period(r):
    """Calculates the orbital period (T) using Kepler's Third Law."""
    return (2 * np.pi * np.power(r, 1.5)) / np.sqrt(G * SOLAR_MASS)

def solar_constants(d):
    return SUN_LUMINOSITY/(4*np.pi*d**2)
    
def eq_tk(sc, alb):
    return np.power(sc*(1-alb)/(4*SB_CONST),1/4)
    
    
# --- DATA INITIALIZATION ---

# solar_system = { Name: (Mass(kg), Radius(km), Distance_from_Sun(mkm)), albedos}
solar_system = {
    "Mercury": (3.301e23, 2440, 57.9, (0.088 + 0.12) / 2),
    "Venus": (4.867e24, 6052, 108.2, 0.76),
    "Earth": (5.972e24, 6371, 149.6, 0.30),
    "Mars": (6.417e23, 3390, 227.9, 0.25),
    "Jupiter": (1.898e27, 69911, 778.6, (0.34 + 0.5) / 2),
    "Saturn": (5.683e26, 58232, 1433.5, 0.34),
    "Uranus": (8.681e25, 25362, 2872.5, 0.30),
    "Neptune": (1.024e26, 24622, 4495.1, 0.29),
}

# Convert dictionary to parallel NumPy arrays
planets = np.array(list(solar_system.keys()))
pl_mass = np.array([data[0] for data in solar_system.values()])
pl_rad = np.array([data[1] * METER_IN_KM for data in solar_system.values()])
pl_dis = np.array([data[2] * 1e6 * METER_IN_KM for data in solar_system.values()])
pl_alb=np.array([data[3] for data in solar_system.values()])
# Identify Earth's index for ratio calculations
earth_idx = np.argmax(planets == "Earth")

# --- CALCULATIONS & OUTPUT ---


def print_table(title, earth_val, unit, labels, ratios):
    """Helper function to print formatted results."""
    print(f"\n{title}: {earth_val:.2e} {unit}")
    print(f"{'Planet':<12 arth="" atio="">15}")
    print("-" * 30)
    for name, ratio in zip(labels, ratios):
        print(f"{name:<12 ratio:="">15.3f}")
    print("-" * 30)


# 1. Surface Gravity
surf_gravities = acc_gravity(pl_mass, pl_rad)
g_ratios = surf_gravities / surf_gravities[earth_idx]
print_table(
    "Earth Surface Gravity", surf_gravities[earth_idx], "m/s²", planets, g_ratios
)

# 2. Volume
pl_vols = vol_sphere(pl_rad)
vol_ratios = pl_vols / pl_vols[earth_idx]
print_table("Earth Volume", pl_vols[earth_idx], "m³", planets, vol_ratios)

# 3. Density
pl_densities = density(pl_mass, pl_vols)
den_ratios = pl_densities / pl_densities[earth_idx]
print_table("Earth Density", pl_densities[earth_idx], "kg/m³", planets, den_ratios)

# 4. Escape Velocity
pl_ev = esc_vel(surf_gravities, pl_rad)
ev_ratios = pl_ev / pl_ev[earth_idx]
print_table("Earth Escape Velocity", pl_ev[earth_idx], "m/s", planets, ev_ratios)

# 5. Sun's Gravitational Pull
pl_sun_f = force_gravity(pl_mass, pl_dis)
f_ratios = pl_sun_f / pl_sun_f[earth_idx]
print_table("Sun's Pull on Earth", pl_sun_f[earth_idx], "N", planets, f_ratios)

# 6. Orbital Speed
pl_os = orbital_speed(pl_dis)
os_ratios = pl_os / pl_os[earth_idx]
print_table("Earth Orbital Speed", pl_os[earth_idx], "m/s", planets, os_ratios)

# 7. Orbital Period
pl_period = orbit_period(pl_dis)
period_ratios = pl_period / pl_period[earth_idx]

earth_days = pl_period[earth_idx] / (60 * 60 * 24)
print(
    f"\nEarth Year Length: {pl_period[earth_idx]:.0f} seconds (~{earth_days:.2f} days)"
)
print(f"{'Planet':<12 atio="" ears="">15}")
print("-" * 30)
for name, ratio in zip(planets, period_ratios):
    print(f"{name:<12 ratio:="">15.3f}")
print("-" * 30)


#solar constants
solar_consts=solar_constants(pl_dis)
#print(solar_consts)
earth_sol_con=solar_consts[earth_idx]
sol_con_ratios=solar_consts/earth_sol_con
print(f"\nEarths Solar Constant:{earth_sol_con:.2f} W/m2")
print(f"{'Planet':<12 const="" olar="" ratio="">15}")
print("-" * 30)
for name, ratio in zip(planets, sol_con_ratios):
    print(f"{name:<12 ratio:="">15.5f}")
print("-" * 30)


#eq temps
eq_temps=eq_tk(solar_consts,pl_alb)

eq_temp_c=eq_temps-273.15
#print(solar_consts)
#earth_eq_c=eq_temp_c[earth_idx]
#eq_tc_ratios=eq_temp_c/earth_eq_c
#print(f"\nEarths Equilibrium Temp:{earth_eq_c:.2f} °K")
print(f"{'Planet':<12 q="" temp="">12} ")
print("-" * 30)
for name,et in zip(planets, eq_temp_c):
    print(f"{name:<12 et:="">10.1f}")
print("-" * 30)
print()
print(f"{'Greenhouse Effect Demo':^45}")
# Actual measured average temperatures in Celsius
actual_temps_c = np.array([167, 464, 15, -65, -110, -140, -195, -200])

# Greenhouse Delta calculation (Actual - Equilibrium)
greenhouse_delta = actual_temps_c - eq_temp_c

print(f"{'Planet':<12 q="" temp="">12} {'Actual (C)':>12} {'GH Delta (C)':>15}")
print("-" * 55)
for name, et, at, delta in zip(planets, eq_temp_c, actual_temps_c, greenhouse_delta):
    print(f"{name:<12 et:="">12.1f} {at:>12.1f} {delta:>15.1f}")
print("-" * 55)

Please feel free to use it if it is of use to you.

The Evaporating Giant: Modeling Atmosphere Loss on HD 189733b

We are looking at HD 189733b, one of the most well-studied "Hot Jupiters". Since it orbits so close to its parent star, its atmosphere is essentially being "boiled off" by intense X-ray and UV radiation. 

HD 189733b looks like a peaceful blue marble from a distance, but don't let the color fool you. Located 63 light-years away, this "Hot Jupiter" is a nightmare world where winds howl at seven times the speed of sound and the rain is made of molten glass.

1. We used the following physics and mathematics 
2. The Shell Method: To find the volume of just the "air" around the planet, we subtract the volume of the solid core (r - t) from the total volume (r).
3. Density Context: The density used (50) gm/cc is a simplified estimate for the upper thermosphere where the escape happens; for comparison, Earth's sea-level density is about 1,225 g/m3. ​
4. The Big Picture: While 6 times 10^8 kg/sec sounds like a catastrophic amount of lost gas, your code reveals the scale of gas giants. Even at that rate, HD 189733b will likely survive for trillions of years—long after its parent star has finished its lifecycle.

import math

def calculate_shell_volume(radius, thickness):
    """
    Calculates the volume of a spherical shell.
    Formula: V = 4/3 * π * (R³ - (R - t)³)
    """
    inner_radius = radius - thickness
    volume = (4/3) * math.pi * (radius**3 - inner_radius**3)
    return volume

# --- Planetary Constants & Data ---
PLANET_NAME = "HD 189733b"
MASS_RATIO_TO_JUPITER = 1.13
ATM_DENSITY = 50           # grams per cubic meter (estimated)
LOSS_RATE_KG_PER_SEC = 6.0e08  # mass loss rate in kg/s

# Reference values for Jupiter
MASS_JUPITER_KG = 1.9e27 
ATM_THICKNESS_KM = 1000  
RADIUS_JUPITER_KM = 70000 

# --- Calculations ---

# 1. Determine the mass of HD 189733b
mass_planet_kg = MASS_RATIO_TO_JUPITER * MASS_JUPITER_KG

# 2. Estimate atmosphere volume (using Jupiter as a structural proxy)
# Result is in km³, which we convert to m³ for density calculations
vol_atm_km3 = calculate_shell_volume(RADIUS_JUPITER_KM, ATM_THICKNESS_KM)
vol_atm_m3 = vol_atm_km3 * 1e09 

# 3. Calculate total mass of the atmospheric layer
# Density (g/m³) * Volume (m³) / 1000 = Total Mass in kg
mass_atm_kg = (vol_atm_m3 * ATM_DENSITY) / 1000

# 4. Time Conversion Constants
SECONDS_IN_YEAR = 60 * 60 * 24 * 365

# 5. Lifespan Estimations
years_to_lose_atm = (mass_atm_kg / LOSS_RATE_KG_PER_SEC) / SECONDS_IN_YEAR
years_to_lose_planet = (mass_planet_kg / LOSS_RATE_KG_PER_SEC) / SECONDS_IN_YEAR

# --- Output Results ---
print(f"--- Planetary Analysis: {PLANET_NAME} ---")
print(f"Total Planetary Mass:     {mass_planet_kg:.2e} kg")
print(f"Atmospheric Mass:        {mass_atm_kg:.2e} kg")
print("-" * 40)
print(f"Time to lose atmosphere:  {years_to_lose_atm:,.0f} years")
print(f"Time to total evaporation: {years_to_lose_planet:,.0f} years")

The Loneliness of the Long-Distance Atom

In my years as a scientist, I’ve often found that scientific notation (10^25) is too efficient—it hides the "physicality" of the universe. To truly feel the scale of the atmosphere, we need to see the zeros.

​In this post, I use a simple Python script to perform a thought experiment: If we take a massive square, 1,000 kilometers on each side, and divide it into tiny 1 cm² tiles, how many molecules would we find in each tile at different altitudes? The transition from the crowded "metropolis" of our breath to the desolate "vacuum" of the Van Allen radiation belt is a humbling reminder of how thin our veil of life really is.

In my years as a scientist, I’ve often found that scientific notation (10^{25}) is too efficient—it hides the "physicality" of the universe. To truly feel the scale of the atmosphere, we need to see the zeros.

​In this post, I use a simple Python script to perform a thought experiment: If we take a massive square, 1,000 kilometers on each side, and divide it into tiny 1 cm² tiles, how many molecules would we find in each tile at different altitudes? The transition from the crowded "metropolis" of our breath to the desolate "vacuum" of the Van Allen radiation belt is a humbling reminder of how thin our veil of life really is.

import math

# ==========================================================
# SCRIPT: The Thinning Universe - From Soil to Space
# Author: Ranjit Singh (The 70-Year-Old Coder)
# Logic: Visualizing density via 1cm2 tiles on a 1000km square
# ==========================================================

# --- 1. Constants & Inputs ---
n_surf = 2.5e25           # molecules/m3 (Surface)
n_vanallen = 900          # atoms/m3 (Van Allen Belt)
ratio_meso = 1e-05        # Density drop at ~80km (Mesosphere)

# --- 2. Geometry Scale ---
side_km = 1000            
tiles_total = (side_km * 100000)**2  # Total 1cm2 tiles in the square

# --- 3. Calculations ---

# Problem 1: Surface Density per tile
mol_tile_surf = n_surf / tiles_total

# Problem 2: Mesosphere Density per tile
mol_tile_meso = mol_tile_surf * ratio_meso

# Problem 3: The Van Allen "Isolation"
# Finding the "Personal Space" of a single atom
atom_tile_va = n_vanallen / tiles_total
tiles_per_atom = 1 / atom_tile_va
side_per_atom_km = math.sqrt(tiles_per_atom) / 100 / 1000

# ==========================================================
# FINAL REPORT: THE AWE OF THE ZEROS
# ==========================================================

print(f"{'ATMOSPHERIC LAYER':<25} | {'MOLECULES PER 1cm2 TILE'}")
print("-" * 60)

# Surface
print(f"{'Earth Surface':<25} | {mol_tile_surf:,.0f}")

# Mesosphere
print(f"{'Mesosphere (80km)':<25} | {mol_tile_meso:,.0f}")

# Van Allen
print(f"{'Van Allen Belt':<25} | {atom_tile_va:.14f}") 

print("\n" + "=" * 60)
print(f"THE SCALE OF EMPTINESS:")
print(f"In the Van Allen Belt, you would have to search a square")
print(f"with a side of {side_per_atom_km:.2f} km just to find ONE atom.")
print("=" * 60)

Closing Thought As the code shows, the surface is a "crowd" of 2.5 billion molecules per centimeter. But once you reach the Van Allen belt, that same centimeter is a ghost town. You would have to travel from Zirakpur to Chandigarh and back (roughly 33 km) just to bump into a single atom. It makes you appreciate every breath a little more, doesn't it?

Tracking LADEE: The Physics of a Lunar Journey

On September 6, 2013, NASA’s Lunar Atmosphere and Dust Environment Explorer (LADEE) lit up the night sky over Virginia. While the mission was designed to study the moon’s thin atmosphere, the most intense part of its journey began at the Wallops Flight Facility. ​

Using trajectory data from the launch, we can use Python to "reverse engineer" the rocket's performance. By applying linear regression to the altitude, range, and velocity, we can see exactly how much power it takes to break free from Earth's grip.

The Raw Data

​We analyzed a 20-second window starting at 190 seconds after launch. At this stage, the rocket is already roughly 150 \text{ km} high—well into the thermosphere—and moving at hypersonic speeds.

​1. Vertical Ascent: Climbing the Gravity Well

​First, we look at the Altitude vs. Time. By fitting a linear trendline to the data, we can determine the vertical component of the rocket's velocity.

​The Equation: Altitude = 0.885 {time} + 149.636

​Vertical Velocity: 0.885 { km/s}

​Accuracy (R^2): 0.9926

​The high R^2 value tells us that during this specific 20-second window, the rocket’s vertical climb was incredibly steady.

​2. The Gravity Turn: Gaining Horizontal Speed

​Rockets don’t just go "up"—to stay in space, they must go "sideways" very fast. This is known as the Range.

​The Equation: Range = 5.433{time} + 444.667

​Horizontal Velocity: 5.433 { km/s}

​Accuracy (R^2): 0.9985

​Note how the horizontal velocity is nearly six times faster than the vertical velocity. LADEE was working hard to reach orbital velocity.

​The Flight Profile at T+190 Seconds

​By combining these vectors, we can calculate the exact state of the rocket at the start of our data set.

How We Calculated This (The Python Approach)

​For the developers and data nerds, we used numpy.polyfit to calculate the "Line of Best Fit" for the trajectory. This allows us to ignore minor sensor noise and see the true physics of the engine's thrust.

import numpy as np
import matplotlib.pyplot as plt

# Global Constants
g = 9.8  # Earth's gravity in m/s^2

def curvefit(x, y, deg=1):
    """Calculates the line of best fit."""
    return np.polyfit(x, y, deg)
    
def calculate_r2(y, y_pred):
    """Determines the Coefficient of Determination (Accuracy)."""
    ss_res = np.sum((y - y_pred)**2)
    ss_tot = np.sum((y - np.mean(y))**2)
    return 1 - (ss_res / ss_tot)

# Raw Trajectory Data: Time(s), Altitude(km), Range(km), Speed(m/s)
trajectory = [
    (190, 150, 446, 5339), (192, 151, 456, 5344),
    (194, 153, 467, 5535), (196, 154, 478, 5629),
    (198, 156, 489, 5744), (200, 158, 500, 5853),
    (202, 159, 511, 5962), (204, 161, 522, 5991),
    (206, 163, 533, 5993), (208, 164, 546, 5990),
]

# Data Extraction
start_time = trajectory[0][0]
time_steps = np.array([t for (t, alt, rng, v) in trajectory])
time_normalized = time_steps - start_time
altitudes = np.array([alt for (t, alt, rng, v) in trajectory])
ranges = np.array([rng for (t, alt, rng, v) in trajectory])
velocities = np.array([v for (t, alt, rng, v) in trajectory])

print(f"{'Rocket Dynamics Analysis':^40}")
print("-" * 40)

# --- Analysis: Vertical Performance ---
coeffs_alt = curvefit(time_normalized, altitudes)
v_vert = coeffs_alt[0]
alt_pred = np.polyval(coeffs_alt, time_normalized)
r2_alt = calculate_r2(altitudes, alt_pred)

print(f"Vertical Velocity:   {v_vert:.3f} km/s")
print(f"Altitude Fit (R2):   {r2_alt:.4f}")

# --- Analysis: Horizontal Performance ---
coeffs_rng = curvefit(time_normalized, ranges)
v_horiz = coeffs_rng[0]
rng_pred = np.polyval(coeffs_rng, time_normalized)
r2_rng = calculate_r2(ranges, rng_pred)

print(f"Horizontal Velocity: {v_horiz:.3f} km/s")
print(f"Range Fit (R2):      {r2_rng:.4f}")

# --- Vector Summation ---
total_vel = np.sqrt(v_vert**2 + v_horiz**2)
launch_angle = np.degrees(np.arctan2(v_vert, v_horiz))

print("-" * 40)
print(f"Total Velocity:      {total_vel:.3f} km/s")
print(f"Flight Path Angle:   {launch_angle:.2f} degrees")

# --- Acceleration Analysis ---
coeffs_vel = curvefit(time_normalized, velocities)
accel = coeffs_vel[0]
print(f"Acceleration:        {accel:.3f} m/s^2")
print(f"G-Force Load:        {accel/g:.2f} G")

Summary By using simple linear regression, we can visualize the sheer scale of a lunar launch. LADEE’s journey reminds us that getting to the moon isn't just about height—it’s about the massive horizontal velocity required to stay in the sky.

The Cosmic Tug-of-War

Roche Limit is the point of no return. If a moon gets too close to its planet, the side facing the planet is pulled much harder than the side facing away. Eventually, this "stretch" (tidal force) becomes stronger than the gravity holding the moon together.

​The Math: Why the Cube Root?

​Since we live in a 3D universe, gravitational influence spreads out over a volume. When we convert that volume back into a linear distance (radius), we use the cube root.

​

Real-World Examples:

​Saturn's Rings: These are basically "failed moons" or shattered debris that live inside Saturn's Roche limit. They can never clump together to form a moon because Saturn's gravity keeps shredding them.

​The Moon: Briefly mention that the Moon is safely outside our Roche limit and is actually moving away from us at about 3.8 cm per year.

Python Code:

Code is based on NASA problem which calculates tidal force radii for Earth, Saturn, and Black holes. Respective targets are Moon, Pan and a supermassive star. It uses formulas based on densities or masses. We used dictionaries of tuples, and two functions. Code is given below:

import math

def calculate_roche_limit(primary_val, target_val, target_radius, is_density=True, fluid=True):
    """
    Calculates the Roche Limit (Tidal Radius).
    
    Args:
        primary_val: Density or Mass of the primary body.
        target_val:  Density or Mass of the target (satellite).
        target_radius: Physical radius of the target body.
        is_density: Boolean, True if using density, False if using mass.
        fluid: Boolean, True for 2.44 coefficient, False for 1.26 (rigid).
    """
    # The fluid constant (2.44) assumes the satellite deforms before breaking.
    # The rigid constant (1.26) is for solid, rocky bodies.
    multiplier = 2.44 if fluid else 1.26
    
    # Mathematically, the ratio of masses or densities works the same 
    # under the cube root when calculating the disruption distance.
    ratio = primary_val / target_val
    return multiplier * target_radius * math.pow(ratio, 1/3)

# --- Dataset for the Blog ---
# Format: "Primary Name": (Target Name, Primary Density/Mass, Target Density/Mass, Target Rad, Actual Dist)
planets_density = {
    "Earth": ("Moon", 5514, 3340, 1737, 384400),
    "Saturn": ("Pan", 687, 400, 14, 133584),
}

stars_mass = {
    "Blackhole A": ("Sun-like Star", 10, 1, 695700, 50e6), # Mass in Solar Masses
    "Blackhole B": ("Sun-like Star", 3, 1, 695700, 50e6),
}

def run_simulation():
    print("-" * 50)
    print("ROCHE LIMIT SIMULATION: WILL IT DISRUPT?")
    print("-" * 50)

    # 1. Density Analysis (Planetary)
    print("\n[Section 1: Planetary Systems (Density Based)]")
    for primary, (target, p_den, t_den, t_rad, dist) in planets_density.items():
        limit = calculate_roche_limit(p_den, t_den, t_rad, is_density=True)
        status = "DISRUPTED" if dist < limit else "STABLE"
        print(f"{primary} vs {target}:")
        print(f"  Distance: {dist:,.0f} km | Roche Limit: {limit:,.0f} km")
        print(f"  Result: {status}")

    # 2. Mass Analysis (Astrophysical)
    print("\n[Section 2: High Gravity Systems (Mass Based)]")
    for primary, (target, p_mass, t_mass, t_rad, dist) in stars_mass.items():
        limit = calculate_roche_limit(p_mass, t_mass, t_rad, is_density=False)
        status = "DISRUPTED" if dist < limit else "STABLE"
        print(f"{primary} vs {target}:")
        print(f"  Distance: {dist:,.0f} km | Roche Limit: {limit:,.0f} km")
        print(f"  Result: {status}")

if __name__ == "__main__":
    run_simulation()
    
    

Tracking the Thaw: Modeling the Antarctic Ice Cap with Python

The Antarctic ice cap is one of the most critical barometers for our planet’s climate health. But when we look at 45 years of satellite data, how do we differentiate between yearly "noise" and a definitive trend? Today, we’re using Linear and Quadratic Regression to analyze ice cap area data from 1979 to 2024 and project what the horizon looks like for 2030.

​Starting in 1979, the September ice cap area sat at roughly 7.2 Million sq km. Fast forward to the 2020s, and we are frequently seeing measurements dip below 4.5 Million sq km.

​To make sense of this, we’ve written a script to calculate the "Goodness of Fit" (R^2) and project future values.

import numpy as np
import matplotlib.pyplot as plt
plt.cla()
def r_squared(actual,predicted):
    ss_res=np.sum((actual-predicted)**2)
    ss_tot=np.sum((actual-np.mean(actual))**2)
    return 1-(ss_res/ss_tot)

def curve_fit(x,y,deg=1):
    return np.polyfit(x,y, deg)
    
#antarctic ice cap in septembe4 
ice_cap_area=[
    #Survey Yr,Ice area(M sq km) 
   (1979, 7.2),
   (1980, 7.9), 
   (1981, 7.3), 
   (1982, 7.5), 
   (1983, 7.5),
   (1984, 7.2), 
   (1985, 6.9), 
   (1986, 7.5), 
   (1987, 7.5), 
   (1988, 7.5)
,
  (1989, 7.0),
   (1990,6.2), 
   (1991,6.6), 
   (1992, 7.6), 
  (1993, 6.5), 
  (1994,7.2 ),
  (1995, 6.1), 
  (1996, 7.9), 
  (1997,6.7),
  (1998, 6.6), 
  (1999, 6.2), 
  (2000, 6.3),
  (2001, 6.8), 
  (2002, 5.9),
   (2003, 6.2),
(2004, 6.1), 
(2005 ,5.6), 
  (2006, 5.9),
  (2007,4.3 ),
  (2008, 4.7), 
   (2009, 5.4),
(2010, 4.9), 
  (2011, 4.6),
(2012, 3.6), #record 
(2013, 5.2),
(2014, 5.3),
(2015, 4.6),
(2016, 4.4),
(2017, 4.8),
(2018, 4.7),
(2019, 4.3),
(2020, 4.2),
(2021, 4.9),
(2022, 4.9),
(2023, 4.4),
(2024, 4.4),
  ]
predict_year=2030
#take 1979 as 0
time_1979=np.array([x-1979 for (x,y) in ice_cap_area])
#print(time_1979)

cap_area=np.array([y for (x,y) in ice_cap_area])
#print(cap_area)

#try linear fit
print("Try a linear fit between data points")
deg_lin=1
lin_fit_coeffs=curve_fit(time_1979,cap_area,deg_lin)
#print(lin_fit_coeffs)
slope=lin_fit_coeffs[0]
intercept=lin_fit_coeffs[1]
print(f"Eqn: ice_cap_area={slope:.4f}*time+{intercept:.4f}")
pred_cap_area_lin=slope*time_1979+intercept
r2_fit_lin=r_squared(cap_area,pred_cap_area_lin)
print(f"R2:{r2_fit_lin:.5f}")
ice_area_2030_lin=slope*(predict_year-1979)+intercept
print(f"Predicted Ice cap area in {predict_year}:{ice_area_2030_lin:.3f} Msqkm")



print()
#nasa asks for quadratic regression also
print("Try a quadratic fit between data")
deg_quad=2
quad_fit_coeffs=curve_fit(time_1979,cap_area,deg_quad)
#print(quad_fit_coeffs)
coeff_x2=quad_fit_coeffs[0]
coeff_x=quad_fit_coeffs[1]
coeff_k=quad_fit_coeffs[2]
print(f"ice_cap_area={coeff_x2:.4f}t2{coeff_x2:.4f}t+{coeff_k:.4f}")
pred_cap_area_quad=coeff_x2*time_1979**2+coeff_x*time_1979+coeff_k
r2_fit_quad=r_squared(cap_area,pred_cap_area_quad)
print(f"R2:{r2_fit_quad:.5f}")
ice_area_2030_quad=coeff_x2*(predict_year-1979)**2+coeff_x*(predict_year-1979)+coeff_k

print(f"Predicted Ice cap area in {predict_year}:{ice_area_2030_quad:.3f} Msqmt")



# Create a smooth range of x-values for the curves
x_range = np.linspace(0, time_1979[-1], 100)

# Calculate y-values for both models
y_lin = slope * x_range + intercept
y_quad = coeff_x2 * x_range**2 + coeff_x * x_range + coeff_k

# Plotting
plt.scatter(time_1979, cap_area, color='gray', alpha=0.5, label='Actual Data')
plt.plot(x_range, y_lin, color='blue', label=f'Linear (R²={r2_fit_lin:.3f})')
plt.plot(x_range, y_quad, color='red', linestyle='--', label=f'Quadratic (R²={r2_fit_quad:.3f})')

plt.xlabel('Years since 1979')
plt.ylabel('Ice Cap Area (M sq km)')
plt.title('Antarctic Ice Cap: Linear vs Quadratic Trend')
plt.legend()
plt.grid(True)
plt.show()

Carbonyl Sulphide COS

The important source of aerosols in the stratosphere is the formation of carbonyl sulfide (COS) droplets. These droplets are more dangerous than carbon dioxide.  

Although volcanism injects millions of tons of SO2 into the atmosphere every few years, scientists have found that during other times, COS is by far the biggest source of sulfur compounds in the stratosphere, leading to the production of sulphuric acid aerosols. 

NASA has tabulated the sources and sinks of COS. At present,  the sinks outnumber the sources. 

From this table, we created a python dictionary with sources names as keys and tuples with rate and type as values. We then created lists of sources, sinks, types using comprehension. We also calculated the number of molecules of COS using avagadro number and volume of atmosphetre. For formatting we used the length of longest name. 

 

# Constants for Atmospheric Calculations
AVOGADRO = 6.023e23          # molecules/mole
MOL_WT_COS = 60              # g/mole
VOL_ATM_KM3 = 4.23e09        # Atmospheric volume in km^3
MT_TO_GRAMS = 1e12           # 1 Million Ton = 10^12 grams
KM3_TO_M3 = 1e09             # 1 km^3 = 10^9 m^3

# Data structure: {Source Name: (Rate in Mt/yr, Type)}
carbonyl_sulfide = {
    "Open ocean": (0.10, "Natural"),
    "Coastal ocean, salt marshes": (0.20, "Natural"),
    "Anoxic soils": (0.02, "Natural"),
    "Wetlands": (0.03, "Natural"),
    "Volcanism": (0.05, "Natural"),
    "Natural CS2 oxidation": (0.21, "Natural"),
    "Oxic soils": (-0.92, "Natural"),
    "Vegetation": (-0.56, "Natural"),
    "Anthropogenic CS2 oxidation": (0.21, "Human"),
    "Biomass burning": (0.07, "Human"),
    "Anthropogenic production": (0.12, "Human"),
    "Precipitation": (0.13, "Other"),
    "DMS oxidation": (0.17, "Other"),
    "Reactions with OH": (-0.24, "Other"),
    "Reactions with oxygen": (-0.02, "Other"),
    "Photodissociation": (-0.05, "Other"),
}

# 1. Calculate dynamic padding for clean table formatting
source_names = list(carbonyl_sulfide.keys())
max_str_len = len(max(source_names, key=len)) + 2 
 
# 2. Categorise into Sources and Sinks
sources = {k: v for k, (v, t) in carbonyl_sulfide.items() if v > 0}
sinks = {k: v for k, (v, t) in carbonyl_sulfide.items() if v < 0}

# 3. Print Sources Table
print(f"{'**** Sources of COS ****':^{max_str_len + 10}}")
print(f"{'Source Name':<{max_str_len}} {'Mt/yr':>8}")
print("-" * (max_str_len + 10))
for name, val in sources.items():
    print(f"{name:<{max_str_len}} {val:>8.2f}")

# 4. Calculate Net Flux (Mass Balance)
sum_sources = sum(sources.values())
sum_sinks = sum(sinks.values())
net_diff = sum_sources + sum_sinks  # Result is in Mt/yr

# 5. Molecule Calculations
# Convert Net Difference (Mt) to Grams, then to Molecules
diff_grams = abs(net_diff) * MT_TO_GRAMS
total_molecules = (diff_grams / MOL_WT_COS) * AVOGADRO

# Convert Atmospheric Volume to m^3
vol_atm_m3 = VOL_ATM_KM3 * KM3_TO_M3
molecules_per_m3 = total_molecules / vol_atm_m3

# 6. Human Impact Analysis
human_sources = {k: v for k, (v, t) in carbonyl_sulfide.items() if t == "Human" and v > 0}
sum_human = sum(human_sources.values())
pct_human_contribution = (sum_human / sum_sources) * 100

# 7. Final Output Results
print("\n" + "="*40)
print(f"{'ATMOSPHERIC ANALYSIS':^40}")
print("="*40)
print(f"Total Annual Sources:    {sum_sources:>6.2f} Mt/yr")
print(f"Total Annual Sinks:      {abs(sum_sinks):>6.2f} Mt/yr")
print(f"Net Flux:                {net_diff:>6.2f} Mt/yr")

if net_diff > 0:
    print("\nRESULT: COS production exceeds reduction.")
else:
    print("\nRESULT: COS reduction exceeds production.")

print(f"Concentration Change:    {molecules_per_m3:.2e} molecules/m3 per year")
print(f"Human Contribution:      {pct_human_contribution:.1f}% of total sources")
print("="*40)

Dwarf Planets

The code is based on finding volume, density and orbit radius of some dwarf planets. Further we tried to match these densities with one's calculated for different mixtures of granite and ice. Based on the results following conclusions can be drawn. Eris is a "rock star"—likely over 70% of granite rock. Pluto is the "balanced" one—roughly a 50/50 split. Ceres is surprisingly light—suggesting either more ice or a lot of trapped salts and clays (brines).

Here is the code

import math

# 1. Helper Functions
def vol_sphere(r):
    """Calculates volume of a sphere given radius."""
    return (4/3) * math.pi * r**3

def calc_density(rock_frac):
    """Calculates theoretical density based on rock/ice mix."""
    rho_rock = 3000  # kg/m3 (Granite-like)
    rho_ice = 900    # kg/m3 (Water ice)
    return rock_frac * rho_rock + (1 - rock_frac) * rho_ice

def orbit(period_yrs):
    """Calculates semi-major axis using Kepler's Third Law."""
    period_secs = period_yrs * secs_in_year
    kepler_const = 4 * math.pi**2 / (G * m_sun)
    return math.pow(period_secs**2 / kepler_const, 1/3)

# 2. Constants & Data
G = 6.674e-11
secs_in_year = 365.25 * 24 * 60 * 60
m_sun = 1.989e30
m_moon = 7.342e22
dia_moon = 3.474e06

dwarfs = {
    # Name: (diameter_%_moon, mass_%_moon, orbit_period_yrs)
    "Ceres": (27, 1.3, 4.6),
    "Pluto": (66, 17.8, 248),
    "Haumea": (36, 5.5, 283),
    "Makemake": (46, 5.4, 310),
    "Eris": (67, 22.7, 557)
}

# 3. Conversions and Calculations
# Convert relative moon data to absolute SI units
dwarfs_actual = {
    name: (
        (dia_moon / 2) * (d_pct / 100), # Radius (m)
        m_moon * (m_pct / 100),         # Mass (kg)
        p                               # Period (yrs)
    ) for name, (d_pct, m_pct, p) in dwarfs.items()
}

# Calculate Volume, Mass, Actual Density, and Orbital Distance
dwarfs_calc = {
    name: (
        vol_sphere(rad), 
        m, 
        m / vol_sphere(rad), 
        orbit(tp)
    ) for name, (rad, m, tp) in dwarfs_actual.items()
}

# 4. Output Physical Data
print(f"{'Name':<10} | {'Density':>8} | {'Distance (AU)':>12}")
print("-" * 40)
for dwarf, (vol, m, d, a) in dwarfs_calc.items():
    # Convert meters to Astronomical Units (AU) for readability
    dist_au = a / 1.496e11
    print(f"{dwarf:<10} | {d:>6.0f} kg/m³ | {dist_au:>10.2f} AU")

# 5. Composition Simulation
simulation = []  
print("\n--- Running Composition Simulation ---")

for rf in range(0, 101, 10):
    rock_p = rf / 100
    theo_den = calc_density(rock_p)
    
    for nm, (vol, m, dens, orb) in dwarfs_calc.items():
        diff = abs(dens - theo_den)
        simulation.append((nm, rf, diff))

# 6. Analyze Results (Find Best Fit)
print("\n--- Most Likely Composition (Best Fit) ---")
for name in dwarfs.keys():
    # Find the rock percentage with the smallest density difference
    planet_data = [item for item in simulation if item[0] == name]
    best_fit = min(planet_data, key=lambda x: x[2])
    
    print(f"{name:<10}: Approximately {best_fit[1]}% Rock / {100-best_fit[1]}% Ice")

The Universal Geyser: Mapping Energy Across the Solar System

Welcome back to the workbench! One of the most beautiful aspects of physics is its universality—the same math we use to calculate the efficiency of a hydroelectric turbine at Nevada Falls applies to the towering plumes of Enceladus. 

In this session, we’re looking at the "Many Faces of Energy." We start by calculating the massive power output of Earth's falling water, then take that same energy budget to the stars.

By keeping our mass and energy fixed, we can see exactly how the varying gravitational pulls of different celestial bodies dictate the height and duration of a geyser’s arc. Whether you're calculating E = mgh for a local power grid or a moon-hopping adventure, the logic remains as steady as the stars. 

Summary of the Script's Logic The script operates on a few core physical principles: 

• Newton’s Law of Gravitation: We calculate local gravity (g) using the mass and radius of the target world:  Conservation of Energy: By setting Kinetic Energy equal to Potential Energy (1/2 mv^2 = mgh), we find that regardless of the world, a fixed amount of energy results in a fixed initial velocity. 
• The Height Equation: The final height is inversely proportional to the local gravity: On a low-gravity moon like Enceladus, that 1.0 MJ of energy doesn't just push water—it launches it into a slow-motion, five-minute journey into the heavens. Let’s look at the numbers.

"""
THE 70-YEAR-OLD CODER: Many Faces of Energy
This script compares the gravitational potential energy and 
geyser mechanics across different worlds.

We use the conservation of energy (E = mgh = 1/2 mv^2) to see 
how high a 2,000 kg mass would rise if given 1.0 Megajoule.
"""
import math

# --- Constants & Conversions ---
G = 6.6743e-11          # Universal Gravitational Constant
MT_IN_FOOT = 0.3048     # Meters in a foot
M3_IN_CUFT = MT_IN_FOOT**3
SECS_IN_DAY = 24 * 60 * 60

def ht(pe, m, g):
    """Calculates height (h) given Potential Energy (pe), Mass (m), and Gravity (g)."""
    return pe / (m * g)

def calculate_g(m, r):
    """Calculates surface gravity (g) using Mass (M) and Radius (R)."""
    return G * m / r**2

# --- Example 1: Earth-bound Power (Nevada Falls) ---
print("--- EXAMPLE 1: NEVADA FALLS HYDRO-POWER ---")
height_nf = 180         # meters
water_flow_cuft = 500   # cubic feet per second
density = 1000          # kg/m3
eff_pct = 70            # Turbine efficiency

water_rate_m3 = water_flow_cuft * M3_IN_CUFT
water_mass_rate = water_rate_m3 * density
energy_per_sec = water_mass_rate * 9.81 * height_nf  # E = mgh

energy_per_day_mj = (energy_per_sec * SECS_IN_DAY) / 1e06
useful_energy_mj = (eff_pct / 100) * energy_per_day_mj

print(f"Water Mass Flow: {water_mass_rate:.2f} kg/s")
print(f"Total Energy/Day: {energy_per_day_mj:,.0f} MegaJoules")
print(f"Useful Energy (70%): {useful_energy_mj:,.0f} MegaJoules\n")

# --- Example 2: The Universal Geyser ---
print("--- EXAMPLE 2: THE UNIVERSAL GEYSER COMPARISON ---")
energy_geyser = 1e06    # 1.0 MegaJoule
mass_lifted = 2000      # 2000 kg (approx. 2 small cars)

# Initial Velocity is the same for all (Kinetic Energy = 1/2 mv^2)
# v = sqrt(2E/m)
init_vel = math.sqrt(2 * energy_geyser / mass_lifted)

print(f"Fixed Energy: {energy_geyser/1e6:.1f} MJ")
print(f"Fixed Mass: {mass_lifted} kg")
print(f"Shared Launch Velocity: {init_vel:.2f} m/s (~114 km/h)\n")

# Dictionary of celestial bodies: (Mass in kg, Radius in meters)
bodies = {
    "Earth":     (5.972e24, 6.371e06),
    "Io":        (8.932e22, 1.8216e06),
    "Europa":    (4.80e22,  1.5608e06),
    "Triton":    (2.14e22,  1.3534e06),
    "Enceladus": (1.08e20,  0.2521e06)
}

# Header for the table
print(f"{'World':<12 g="" m="" s2="">10} | {'Height (m)':>12} | {'Time to Top (s)':>15}")
print("-" * 58)

for name, (m, r) in bodies.items():
    g_local = calculate_g(m, r)
    height = ht(energy_geyser, mass_lifted, g_local)
    time_to_top = init_vel / g_local
    
    print(f"{name:<12 g_local:="">10.3f} | {height:>12.2f} | {time_to_top:>15.2f}")

print("\nConclusion: On Enceladus, you can watch the geyser rise for nearly 5 minutes!")