大气参数计算MATLAB代码

% 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 < 1
    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） < 0
        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）;