function [acousticSpeed,rho] = getRho(state)
% Compute atmospheric properties from mean sea level altitude.
g = 9.8066; % m/sec^2
lapserate = -0.0065; % K/m
R = 287.05; % J/kg-K
gamma = 1.4; % adiabatic constant for air
z0 = 0; % mean sea level reference altitude (m)
T0 = 288.15; % sea level temperature (K)
P0 = 1.01325e5; % sea level pressure (Pa)

z = state(3); % altitude above mean sea level (m)

if z <= 11000
    T = T0 + lapserate*(z-z0);
    P = P0*(T/T0)^(-g/(lapserate*R));
else
    T = T0 + lapserate*(11000-z0);
    P11k = P0*(T/T0)^(-g/(lapserate*R));
    P = P11k*exp(-g*(z-11000)/(R*T));
end

rho = P/(R*T);
acousticSpeed = sqrt(gamma*R*T);
end