Modelo de Atmósfera Estándar Internacional (ISA)

Asignatura: Vehículos aeroespaciales

Departamento: Ciencia y Tecnología de Materiales y Fluidos

Centro: Escuela Universitaria Politécnica de Teruel

Profesor: Adrián Navas Montilla

import numpy as np
import matplotlib.pyplot as plt

# Constants
g = 9.81         # gravity [m/s²]
R = 287.05          # specific gas constant for dry air [J/(kg·K)]
gamma = 1.4

# Data from Table 5.3 (Book: S. Corda)
z_n = np.array([0, 11, 20, 32, 47, 51, 71])                      # [km]
lapse_n = np.array([-6.5, 0.0, 1.0, 2.8, 0.0, -2.8, -2.0])      # [K/km]
T_n = np.array([288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65])  # [K]
p_n = np.array([101325, 22632.1, 5474.8, 868.019, 110.906, 66.9389, 3.95642])  # [Pa]
rho_n = np.array([1.2250, 0.36392, 8.8035e-2, 1.3550e-2, 1.4275e-3, 8.160e-4, 6.4211e-5])  # [kg/m³]

# Arrays to hold the complete profile
z_total = []
T_total = []
p_total = []
rho_total = []
c_total = []

# Loop through each layer
for n in range(len(z_n)-1):  # up to index 6 (mesosphere)
    z0 = z_n[n]*1000 # [m]
    z1 = z_n[n+1]*1000# [m]
    lapse = lapse_n[n]/1000 # [K/m]
    T0 = T_n[n]
    p0 = p_n[n]

    pts = 20

    z = np.linspace(z0, z1, pts)             # [km]

    if lapse == 0:
        # Isothermal layer
        T = T0 * np.ones_like(z)             # Temperature remains constant
        p = p0 * np.exp(-g * (z - z0) / (R * T0))
    else:
        # Gradient layer
        T = T0 + lapse * (z - z0)            # Temperature profile
        p = p0 * (T / T0) ** (-g / (lapse * R))

    rho = p / (R * T)                           # Density from ideal gas law
    c_sound = np.sqrt(gamma * R * T)                    # Sound speed

    #print(n)
    #print(z[0],p[0],rho[0])

    # Append to total profile
    z_total.extend(z)
    T_total.extend(T)
    p_total.extend(p)
    rho_total.extend(rho)
    c_total.extend(c_sound)

# Convert to numpy arrays
z_total = np.array(z_total)
T_total = np.array(T_total)
p_total = np.array(p_total)
rho_total = np.array(rho_total)

# --- Optional: Plotting ---
plt.figure(figsize=(15, 7))

plt.subplot(1, 4, 1)
plt.plot(T_total, z_total)
plt.xlabel("Temperature [K]")
plt.ylabel("Altitude [m]")
plt.title("Temperature Profile")

plt.subplot(1, 4, 2)
plt.plot(p_total*0.001, z_total)
plt.xlabel("Pressure [kPa]")
plt.title("Pressure Profile")

plt.subplot(1, 4, 3)
plt.plot(rho_total, z_total)
plt.xlabel("Density [kg/m³]")
plt.title("Density Profile")

plt.subplot(1, 4, 4)
plt.plot(c_total, z_total)
plt.xlabel("Speed of sound [m/s]")
plt.title("Speed of sound")


plt.tight_layout()
plt.show()

# ----------------------
# Combined plot with shared y-axis
# ----------------------
fig, ax1 = plt.subplots(figsize=(5, 7))

ax1.set_ylabel("Altitude [m]", fontsize=12)
ax1.set_xlabel("Temperature [K]", color='tab:red')
ax1.plot(T_total, z_total, color='tab:red', label='Temperature')
ax1.set_xlim(150, 300)
ax1.tick_params(axis='x', labelcolor='tab:red')

# Add more axes
ax2 = ax1.twiny()
ax2.plot(p_total * 0.001, z_total, color='tab:blue', label='Pressure')
ax2.set_xlabel("Pressure [kPa]", color='tab:blue')
ax2.tick_params(axis='x', labelcolor='tab:blue')

ax3 = ax1.twiny()
ax3.spines['top'].set_position(("axes", 1.1))
ax3.plot(rho_total, z_total, color='tab:green', label='Density')
ax3.set_xlabel("Density [kg/m³]", color='tab:green')
ax3.tick_params(axis='x', labelcolor='tab:green')

ax4 = ax1.twiny()
ax4.spines['top'].set_position(("axes", 1.2))
ax4.plot(c_total, z_total, color='tab:purple', label='Speed of Sound')
ax4.set_xlim(200, 350)
ax4.set_xlabel("Speed of sound [m/s]", color='tab:purple')
ax4.tick_params(axis='x', labelcolor='tab:purple')


# Atmospheric layer names (from Table 5.3, n = 0 to 6)
layer_names = [
    "Troposphere",
    "Tropopause",
    "Stratosphere (1)",
    "Stratosphere (2)",
    "Stratopause",
    "Mesosphere ",
    "Mesosphere (2)"
]

# Altitudes at boundaries z_n
z_bounds = z_n * 1000

# Plot dashed lines and text labels
for i in range(len(z_bounds) - 1):
    z_start = z_bounds[i]
    z_end = z_bounds[i + 1]
    z_mid = (z_start + z_end) / 2

    # Horizontal dashed line at boundary (not at 0 km)
    if i > 0:
        ax1.axhline(z_start, color='gray', linestyle='--', linewidth=0.8)

    # Add layer name at midpoint
    if i == 1:
        ax1.text(180, z_mid, layer_names[i],
                 fontsize=9, color='black', va='center')
    else:
        ax1.text(165, z_mid, layer_names[i],
             fontsize=9, color='black', va='center')

plt.title("International Standard Atmosphere (ISA)", fontsize=14)
plt.grid(True)
plt.tight_layout()
plt.show()