function [rho, Temp, a, mu] = atm(h)

T_o = 518.7;
rho_o = 0.002377;
R = 1716;
gamma = 1.4;

if h < 36089
   theta = 1- 6.875E-6*h;
   Temp = T_o*theta;
   rho = rho_o*theta^4.2561;
else
   Temp = -69.7 + 459.67;
   rho = rho_o*0.2971*exp(-((h-36089)/20806.7));
end

a = sqrt(gamma*R*Temp);	
mu = (1E-10*0.3170)*Temp^1.5*(734.7/(Temp+216));
