1.1 Atmosphere Structure

#!/usr/bin/env python3
"""
atmosphere_structure.py - Plot standard atmosphere temperature and pressure profiles

Run: python3 atmosphere_structure.py
Requires: pip install numpy matplotlib
"""
import numpy as np
import matplotlib.pyplot as plt

# Constants
g = 9.81      # m/s^2
R = 287.0     # J/(kg·K) for dry air
p0 = 101325   # Pa at sea level
T0 = 288.15   # K at sea level

# Standard atmosphere layers (altitude in km, lapse rate in K/km)
layers = [
    (0, 11, -6.5),      # Troposphere
    (11, 20, 0),        # Tropopause
    (20, 32, 1.0),      # Lower stratosphere
    (32, 47, 2.8),      # Upper stratosphere
    (47, 51, 0),        # Stratopause
    (51, 71, -2.8),     # Mesosphere lower
    (71, 85, -2.0),     # Mesosphere upper
]

def standard_atmosphere(z_km):
    """Calculate T and p at altitude z (km)"""
    z = z_km * 1000  # Convert to meters
    T = T0
    p = p0

    z_prev = 0
    for z_bot, z_top, lapse in layers:
        z_bot_m, z_top_m = z_bot * 1000, z_top * 1000
        if z <= z_bot_m:
            break

        dz = min(z, z_top_m) - z_bot_m
        L = lapse / 1000  # K/m

        if L == 0:
            # Isothermal layer
            p = p * np.exp(-g * dz / (R * T))
        else:
            T_new = T + L * dz
            p = p * (T_new / T) ** (-g / (R * L))
            T = T_new

        if z <= z_top_m:
            break

    return T, p

# Generate profile
altitudes = np.linspace(0, 85, 500)
temps = np.array([standard_atmosphere(z)[0] for z in altitudes])
pressures = np.array([standard_atmosphere(z)[1] for z in altitudes])

# Plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 8))

# Temperature profile
ax1.plot(temps - 273.15, altitudes, 'r-', linewidth=2)
ax1.axhline(y=11, color='gray', linestyle='--', alpha=0.5, label='Tropopause')
ax1.axhline(y=50, color='gray', linestyle='-.', alpha=0.5, label='Stratopause')
ax1.axhline(y=85, color='gray', linestyle=':', alpha=0.5, label='Mesopause')
ax1.set_xlabel('Temperature (°C)', fontsize=12)
ax1.set_ylabel('Altitude (km)', fontsize=12)
ax1.set_title('Temperature Profile', fontsize=14)
ax1.grid(True, alpha=0.3)
ax1.legend()
ax1.set_xlim(-100, 20)

# Pressure profile
ax2.semilogx(pressures / 100, altitudes, 'b-', linewidth=2)
ax2.set_xlabel('Pressure (hPa)', fontsize=12)
ax2.set_ylabel('Altitude (km)', fontsize=12)
ax2.set_title('Pressure Profile', fontsize=14)
ax2.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig('atmosphere_structure.png', dpi=150, bbox_inches='tight')
plt.show()
print("Saved as atmosphere_structure.png")