-
Notifications
You must be signed in to change notification settings - Fork 0
/
EarthAtmosRho.m
52 lines (49 loc) · 1.33 KB
/
EarthAtmosRho.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
function [rho] = EarthAtmosRho(z,R0)
% get geopotential
h = R0*z./(R0+z);
% Pressure
if h <= 11000
P = 101325*(288.15/(288.15-6.5*(h/1000))).^(34.1632/-6.5);
elseif h <= 20000
P = 22632.06*exp(-34.1632*((h/1000)-11)/216.65);
elseif h <= 32000
P = 5474.889*(216.65/(216.65+(h/1000-20))).^(34.1632);
elseif h <= 47000
P = 868.0187*(228.65/(228.65+2.8*(h/1000-32))).^(34.1632/2.8);
elseif h <= 51000
P = 110.9063*exp(-34.1632*((h/1000)-47)/270.65);
elseif h <= 71000
P = 66.93887*(270.65/(270.65-2.8*(h/1000-51))).^(34.1632/-2.8);
elseif z <= 86000
P = 3.956420*(214.65/(214.65-2.0*(h/1000-71))).^(34.1632/-2.0);
else
P = exp(-0.0000000422012*(z/1000).^5+0.0000213489*(z/1000).^4-0.00426388*(z/1000).^3+0.421404*(z/1000).^2-20.8270*(z/1000)+416.225);
end
% Temperature
if h <= 11000
T = 288.15-6.5*(h/1000);
elseif h <= 20000
T = 216.65;
elseif h <= 32000
T = 196.65+(h/1000);
elseif h <= 47000
T = 139.05+2.8*(h/1000);
elseif h <= 51000
T = 270.65;
elseif h <= 71000
T = 413.45-2.8*(h/1000);
elseif z <= 86000
T = 356.65-2.0*(h/1000);
elseif z <= 91000
T = 186.8673;
elseif z <= 110000
T = 263.1905-76.3232*sqrt(1-((z/1000-91)/-19.9429).^2);
else
T = 240+12*(z/1000-110);
end
% Density
rho = P./T/287.053;
if z > 120000
rho = 0;
end
end