%?function?[t?p?rho?a]?=?atmosphere（height?unit）
%?Richard?Rieber
%?rrieber@gmail.com
%?Updated?3/17/2006
%
%?Function?to?determine?temperature?pressure?density?and?speed?of?sound
%?as?a?function?of?height.??Function?based?on?US?Standard?Atmosphere?of
%?1976.??All?calculations?performed?in?metric?units?but?converted?to
%?english?if?the?user?chooses.??Assuming?constant?gravitational
%?acceleration.
%
%?Input??o?height?（feet?or?meters）?can?be?a?1-dimensional?vector
%????????o?unit?-?optional;?default?is?metric;
%?????????????????boolean?True?for?english?units?False?for?metric
%?Output?o?t?-?Temperature?（K?（default）?or?R）
%????????o?p?-?Pressure?（pa?（default）?or?psi）
%????????o?rho?-?Density?（kg/m^3?（default）?or?lbm/ft^3）
%????????o?a?-?speed?of?sound?（m/s?（default）?or?ft/s）

function?[t?p?rho?a]?=?atmosphere（heightunit）

if?nargin?????error（‘Error:??Not?enough?inputs:?see?help?atmosphere.m‘）
elseif?nargin?>?2
????error（‘Error:??Too?few?inputs:?see?help?atmosphere.m‘）
elseif?nargin?==?2
????unit?=?logical（unit）;
elseif?nargin?==?1
???unit?=?0;
end

for?j?=?1:length（height）
	if?height（j）?????????error（‘Height?should?be?greater?than?0‘）
	end
	
	m2ft?=?3.2808;
	K2R?=?1.8;
	pa2psi?=?1.450377377302092e-004;
	kg2lbm?=?2.20462262184878;
	
	if?unit
????????height（j）?=?height（j）/m2ft;
	end
	
	g?=?9.81;???????%Acceleration?of?gravity?（m/s/s）
	gamma?=?1.4;????%Ratio?of?specific?heats
	R?=?287;????????%Gas?constant?for?air?（J/kg-K）
	
	%Altitudes?（m）
	Start?=?0;
	H1?=?11000;
	H2?=?20000;
	H3?=?32000;
	H4?=?47000;
	H5?=?51000;
	H6?=?71000;
	H7?=?84852;
	
	%Lapse?Rates?（K/m）
	L1?=?-0.0065;
	L2?=?0;
	L3?=?.001;
	L4?=?.0028;
	L5?=?0;
	L6?=?-.0028;
	L7?=?-.002;
	
	%Initial?Values
	T0?=?288.16;?????%（k）
	P0?=?1.01325e5;??%（pa）
	Rho0?=?1.225;????%（kg/m^3）
	
	if?height（j）?<=?H1
????????[TNew?PNew?RhoNew]?=?Gradient（Start?height（j）?T0?P0?Rho0?L1）;????
????????
	elseif?height（j）?>?H1?&?height（j）?<=?H2
????????[TNew?PNew?RhoNew]?=?Gradient（Start?H1?T0?P0?Rho0?L1）;
????????[TNew?PNew?RhoNew]?=?IsoThermal（H1?height（j）?TNew?PNew?RhoNew）;
????????
	elseif?height（j）?>?H2?&?height（j）?<=?H3
????????[TNew?PNew?RhoNew]?=?Gradient（Start?H1?T0?P0?Rho0?L1）;
????????[TNew?PNew?RhoNew]?=?IsoThermal（H1?H2?TNew?PNew?RhoNew）;
????????[TNew?PNew?RhoNew]?=?Gradient（H2?height（j）?TNew?PNew?RhoNew?L3）;
????????
	elseif?height（j）?>?H3?&?height（j）?<=?H4
????????[TNew?PNew?RhoNew]?=?Gradient（Start?H1?T0?P0?Rho0?L1）;
????????[TNew?PNew?RhoNew]?=?IsoThermal（H1?H2?TNew?PNew?RhoNew）;
??????