% Air.m
% 輸入：幾何高度（海拔Altitude）
% 輸出：大氣溫度、壓力、密度、音速等
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [T,P,D,A] = GetAir(h)

	if(nargin<1)
      disp('請輸入一個高度');
      return;
	end

	ACCELERATION_GRAVITY = 9.80665;		%重力加速度(m/s2)
	RADIUM_EARTH = 6356766.0;				%地球半徑(m)
	CONSTANT_GAS = 287.05287;				%氣體常數(J/K.k)
	
	%設定高度，將幾何高度轉化為位勢高度	
	if (h < 0.0) h = 0.0;end;
	if (h > 51000.0) h = 51000.0;end;
	H = (RADIUM_EARTH * h) / (RADIUM_EARTH + h);

	%計算溫度、壓力、密度
	temp = - (ACCELERATION_GRAVITY / CONSTANT_GAS);
	if (H < 11000.0)
		T = 288.15 - 0.0065 * H;
		P = 101325.0 * (T / 288.15) ^(temp / (-0.0065));
	elseif (H < 20000.0)
		T = 216.65;
		P = 22632.1 * exp(temp * (H - 11000.0) / 216.65);
	elseif (H < 32000.0)
		T = 216.65 + 0.001 * (H - 20000.0);
		P = 5474.88 * (T / 216.65) ^ (temp / 0.001);
	elseif (H < 47000.0)
		T = 228.65 + 0.0028 * (H - 32000.0);
		P = 868.018 * (T / 228.65) ^ (temp / 0.0028);
	else
		T = 270.65;
		P = 110.906 * exp(temp * (H - 47000.0) / 270.65);
	end
    
	D = P / (CONSTANT_GAS * T);
	A = 20.046796 * sqrt(T);