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.