Standard Atmosphere

Computes Standard Atmospheric Properties

Standard Atmosphere

Computes Standard Atmospheric Properties

function [h_out, H, T, a, P, rho, nu, mu] = ISA_Standard_Atmosphere(h)
% ISA_Standard_Atmosphere  Enhanced ISA model with stratospheric lapse layer
%
%   Computes standard‐atmosphere properties at geometric altitude h (m):
%     • h_out – geometric altitude [m]
%     • H     – geopotential altitude [m]
%     • T     – temperature [K]
%     • a     – speed of sound [m/s]
%     • P     – pressure [Pa]
%     • rho   – density [kg/m³]
%     • nu    – kinematic viscosity [m²/s]
%     • mu    – dynamic viscosity [Pa·s]
%
%   Layer structure:
%     – 0–11 km:    linear lapse L1 = 6.5 K/km
%     – 11–20 km:   isothermal at T_trop
%     – 20–32 km:   inverse lapse L2 = –1 K/km
%     – >32 km:     isothermal at T32

    %% Constants
    Re    = 6.356766e6;      % Earth radius [m]
    T0    = 288.15;          % sea‐level temperature [K]
    P0    = 101325;          % sea‐level pressure [Pa]
    mu0   = 1.7894e-5;       % sea‐level viscosity [Pa·s]
    Rgas  = 287.053;         % specific gas constant [J/(kg·K)]
    g     = 9.80665;         % gravity [m/s²]
    gamma = 1.4;             % ratio of specific heats

    %% Layer breakpoints & lapse rates
    h_trop    = 11000;       % tropopause [m]
    L1        = 6.5e-3;      % tropospheric lapse [K/m]
    h_strat2  = 20000;       % base of stratospheric lapse [m]
    h_strat3  = 32000;       % top of stratospheric lapse [m]
    L2        = -1e-3;       % stratospheric lapse [K/m]

    %% Reference values at layer boundaries
    T_trop = T0 - L1*h_trop;
    P_trop = P0 * (T_trop/T0)^(g/(Rgas*L1));

    T20 = T_trop;
    P20 = P_trop * exp(-g*(h_strat2 - h_trop)/(Rgas*T_trop));

    T32 = T20 - L2*(h_strat3 - h_strat2);
    P32 = P20 * (T32/T20)^(g/(Rgas*L2));

    %% Geopotential altitude
    h_out = h;
    H     = Re .* h ./ (Re + h);

    %% Preallocate outputs
    T = zeros(size(H));
    P = zeros(size(H));

    %% Regions
    idx1 = H <= h_trop;
    idx2 = H > h_trop  & H <= h_strat2;
    idx3 = H > h_strat2 & H <= h_strat3;
    idx4 = H > h_strat3;

    %% Troposphere (0–11 km)
    T(idx1) = T0 - L1 .* H(idx1);
    P(idx1) = P0 .* (T(idx1)/T0).^(g/(Rgas*L1));

    %% Lower Stratosphere (11–20 km): isothermal
    T(idx2) = T_trop;
    P(idx2) = P_trop .* exp(-g*(H(idx2)-h_trop)/(Rgas*T_trop));

    %% Upper Stratosphere (20–32 km): inverse lapse
    T(idx3) = T20 - L2 .* (H(idx3)-h_strat2);
    P(idx3) = P20 .* (T(idx3)/T20).^(g/(Rgas*L2));

    %% Above 32 km: isothermal at T32
    T(idx4) = T32;
    P(idx4) = P32 .* exp(-g*(H(idx4)-h_strat3)/(Rgas*T32));

    %% Density (ideal gas law)
    rho = P ./ (Rgas .* T);

    %% Viscosity (Sutherland-like) & kinematic viscosity
    mu = mu0 .* (T./T0).^(3/2);
    nu = mu ./ rho;

    %% Speed of sound
    a = sqrt(gamma * Rgas .* T);
end