Saturday, 2 May 2026

Visualizing Particle Kinematics with Python

Understanding the relationship between position, velocity, and acceleration is fundamental to physics. However, seeing how these variables evolve simultaneously can be challenging with equations alone. In this post, we use Python, NumPy, and Matplotlib to model a particle's motion. By leveraging lambda functions and vectorized operations, we can generate a high-resolution continuous model while calculating specific discrete data points for a clear, side-by-side comparison.

The Python Script



import numpy as np
import matplotlib.pyplot as plt

# 1. Define Kinematic Functions using Lambda and NumPy
# These handle both single numbers and arrays automatically
pos = lambda t: 2 * np.sin(np.pi * t / 4)
vel = lambda t: 2 * np.pi / 4 * np.cos(np.pi * t / 4)
acc = lambda t: -2 * (np.pi / 4) ** 2 * np.sin(np.pi * t / 4)

# 2. Generate Discrete Data Points
time_steps = [t / 10 for t in range(21)]
positions, velocities, accel = [], [], []

print(f"{'Time(s)':<10 100="" 10="" 1="" 2="" 3.="" 4.="" 45="" a:="" a="" a_curve="" acc="" accel.append="" accel="" acceleration="" article="" ax1.grid="" ax1.plot="" ax1.scatter="" ax1.set_title="" ax1.set_ylabel="" ax1="" ax2.grid="" ax2.plot="" ax2.scatter="" ax2.set_ylabel="" ax2="" ax3.grid="" ax3.plot="" ax3.scatter="" ax3.set_xlabel="" ax3.set_ylabel="" ax3="" cc="" cm="" code="" color="green" curves="" el="" f="" fig="" figsize="(8," for="" gen_tpts="" generate="" ime="" in="" kinematics:="" label="Continuous Acc" linestyle="--" os="" p:="" p="" p_curve="" plot="" plotting="" plt.show="" plt.subplots="" plt.tight_layout="" position="" positions.append="" positions="" print="" rue="" s2="" s="" sharex="True)" smooth="" time_steps:="" time_steps="" ts:="" ts="" v:="" v="" v_curve="" vel="" velocities.append="" velocities="" velocity="" visualization="">
Conclusion As shown in the resulting plots, Python provides a powerful environment for physics simulation. By using numpy.linspace, we created a smooth curve that represents the true continuous nature of the motion, while the scatter points highlight our specific calculated intervals. This "stacked" subplot approach is particularly useful because it allows us to see exactly how a peak in position corresponds to a zero-crossing in velocity, making the abstract calculus of motion much more intuitive.

Thursday, 30 April 2026

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.

Star velocities


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

Wednesday, 29 April 2026

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.

Tuesday, 28 April 2026

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

Saturday, 25 April 2026

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.  
Solar system


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.

Thursday, 23 April 2026

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

Visualizing Particle Kinematics with Python

Understanding the relationship between position, velocity, and acceleration is fundamental to physics. However, seeing how these variables e...