-
Notifications
You must be signed in to change notification settings - Fork 1
/
EarthAtmos.m
90 lines (84 loc) · 3.06 KB
/
EarthAtmos.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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
function value = EarthAtmos(z,R0,Type)
% EarthAtmos: Coarse model of Earth's atmosphere
% Input:
% h = geometric altitude (m)
% Type = type of data requested, as string ('rho','T','P','a')
% Output:
% rho = density (kg/m^3)
% T = temperature (k)
% P = pressure (Pa)
% a = speed of sound (m/s)
h = R0*z./(R0+z);
% Generate Values:
if h <= 11*10^3
T = 288.15 - 6.5*h/1000;
P = 101325*(288.15/T)^(34.1632/-6.5);
rho = P/T/287.053;
a = sqrt(1.4*287.053*T);
elseif (h > 11*10^3) && (h <= 20*10^3)
T = 216.65;
P = 22632.06*exp(-34.1632*(h/1000-11)/T);
rho = P/T/287.053;
a = sqrt(1.4*287.053*T);
elseif (h > 20*10^3) && (h <= 32*10^3)
T = 196.65 + h/1000;
P = 5474.889*(216.65/(216.65+(h/1000-20)))^34.1632;
rho = P/T/287.053;
a = sqrt(1.4*287.053*T);
elseif (h > 32*10^3) && (h <= 47*10^3)
T = 139.05 + 2.8*h/1000;
P = 868.0187*(228.65/(228.65+2.8*(h/1000-32)))^(34.1632/2.8);
rho = P/T/287.053;
a = sqrt(1.4*287.053*T);
elseif (h > 47*10^3) && (h <= 51*10^3)
T = 270.65;
P = 110.9063*exp(-34.1632*(h/1000-47)/270.65);
rho = P/T/287.053;
a = sqrt(1.4*287.053*T);
elseif (h > 51*10^3) && (h <= 71*10^3)
T = 413.45 - 2.8*h/1000;
P = 66.93887*(270.65/(270.65-2.8*(h/1000-51)))^(34.1632/-2.8);
rho = P/T/287.053;
a = sqrt(1.4*287.053*T);
elseif (h > 71*10^3) && (h <= 86*10^3)
T = 356.65 - 2.0*h/1000;
P = 3.956420*(214.65/(214.65-2*(h/1000-71)))^(34.1632/-2);
rho = P/T/287.053;
a = sqrt(1.4*287.053*T);
elseif (h > 86*10^3) && (h <= 91*10^3)
T = 186.8673;
P = exp(2.159582*10^-6*(h/1000)^3 - 4.836957*10^-4*(h/1000)^2 - 0.1425192*(h/1000) + 13.47530);
rho = exp(-3.322622*10^-6*(h/1000)^3 + 9.111460*10^-4*(h/1000)^2 - 0.2609971*(h/1000) + 5.944694);
a = sqrt(1.4*287.053*T);
elseif (h > 91*10^3) && (h <= 100*10^3)
T = 263.1905-76.3232*sqrt(1-((h/1000-91)/-19.9429)^2);
P = exp(3.304895*10^-5*(h/1000)^3 - 0.009062730*(h/1000)^2 + 0.6516698*(h/1000) - 11.03037);
rho = exp(2.873405*10^-5*(h/1000)^3 - 0.008492037*(h/1000)^2 + 0.6541179*(h/1000) - 23.62010);
a = sqrt(1.4*287.053*T);
elseif (h > 100*10^3) && (h <= 110*10^3)
T = 263.1905-76.3232*sqrt(1-((h/1000-91)/-19.9429)^2);
P = exp(6.693926*10^-5*(h/1000)^3 - 0.01945388*(h/1000)^2 + 1.719080*(h/1000) - 47.75030);
rho = exp(-1.240774*10^-5*(h/1000)^4 + 0.005162063*(h/1000)^3 - 0.8048342*(h/1000)^2 + 55.55996*(h/1000) - 1443.338);
a = sqrt(1.4*287.053*T);
elseif (h > 110*10^3) && (h <= 120*10^3)
T = 240 + 12*(h/1000-110);
P = exp(-6.539316*10^-5*(h/1000)^3 + 0.02485568*(h/1000)^2 - 3.223620*(h/1000) + 135.9355);
rho = exp(-8.854164*10^-5*(h/1000)^3 + 0.03373254*(h/1000)^2 - 4.390837*(h/1000) + 176.5294);
a = sqrt(1.4*287.053*T);
else % height above 120 km
T = 360;
P = 0;
rho = 0;
a = 0;
end
% Return requested type:
if strcmp(Type,'rho')
value = rho;
elseif strcmp(Type,'T')
value = T;
elseif strcmp(Type,'P')
value = P;
else % Type == 'a'
value = a;
end
end