1
$\begingroup$

enter image description here

How can I make the plot shown above? It gives pressure as a function of altitude according to the barometric formula.

$\endgroup$
5

1 Answer 1

5
$\begingroup$

There are two different equations for computing pressure at various height regimes below 86 km.

The first equation is used when the value of standard temperature lapse rate is not equal to zero.

ClearAll["Global`*"]

Pb0 = 101325  (*Pa*);
Pb1 = 22632.1 (*Pa*);
Pb2 = 5474.89 (*Pa*);
Pb3 = 868.02  (*Pa*);
Pb4 = 110.91  (*Pa*);
Pb5 = 66.94   (*Pa*);
Pb6 = 3.96    (*Pa*);
Tb0 = 288.15  (*K*);
Tb1 = 216.65  (*K*);
Tb2 = 216.65  (*K*);
Tb3 = 228.65  (*K*);
Tb4 = 270.65  (*K*);
Tb5 = 270.65  (*K*);
Tb6 = 214.15  (*K*);
g0 = 9.860665 (*m/s^2*);
M = 0.0289644 (*kg/mol*);
hb0 = 0       (*m*);
hb1 = 11000   (*m*);
hb2 = 20000   (*m*);
hb3 = 32000   (*m*);
hb4 = 47000   (*m*);
hb5 = 51000   (*m*);
hb6 = 71000   (*m*);
R = 8.3144598 (* J/mol/K*);
Lb0 = -0.0065 (*k/m*);
Lb1 = 10^-9   (*k/m. I must put a small value not a Zero *);
Lb2 = 0.001   (*k/m*);
Lb3 = 0.0028  (*k/m*);
Lb4 = 10^-9   (*k/m. I must put a small value not a Zero *);
Lb5 = -0.0028 (*k/m*);
Lb6 = -0.002  (*k/m*);

P0[h_] := Pb0*(Tb0/(Tb0 + Lb0*(h - hb0)))^((g0*M)/(R*Lb0));
P1[h_] := Pb1*(Tb1/(Tb1 + Lb1*(h - hb1)))^((g0*M)/(R*Lb1));
P2[h_] := Pb2*(Tb2/(Tb2 + Lb2*(h - hb2)))^((g0*M)/(R*Lb2));
P3[h_] := Pb3*(Tb3/(Tb3 + Lb3*(h - hb3)))^((g0*M)/(R*Lb3));
P4[h_] := Pb4*(Tb4/(Tb4 + Lb4*(h - hb4)))^((g0*M)/(R*Lb4));
P5[h_] := Pb5*(Tb5/(Tb5 + Lb5*(h - hb5)))^((g0*M)/(R*Lb5));
P6[h_] := Pb6*(Tb6/(Tb6 + Lb6*(h - hb6)))^((g0*M)/(R*Lb6));

PP[h_] :=  Piecewise[{{P0[h], h <= 11000}, {P1[h], 11000 <= h <= 20000}, {P2[h],
20000 <= h <= 32000}, {P3[h], 32000 <= h <= 47000}, {P4[h], 
47000 <= h <= 51000}, {P5[h], 51000 <= h <= 71000}, {P6[h], 71000 <= h <= 86000}}]

Plot[Labeled[{PP[h]}, "Pressure as a function of the height above the sea level", 
9000], {h, 0, 71000}, PlotRange -> All, PlotStyle -> Red,AxesLabel -> {"h[m]", "P[Pa]"}]
Plot[Labeled[{PP[h]}, "Pressure as a function of the height above the sea level", 
100], {h, 0, 86000}, ScalingFunctions -> {"Log", None}, 
PlotRange -> All, PlotStyle -> Red, AxesLabel -> {"h[m]", "P[Pa]"}]

Plot[Labeled[{PP[h]/100}, "Pressure as a function of the height above the sea level", 
5000], {h, 0, 21000}, PlotRange -> {Automatic, {0, 1013.25}}, 
PlotStyle -> Red, AxesLabel -> {"h[km]", "P[hPa]"}, 
Ticks -> {Table[{1000 i, i}, {i, 1, 21, 2}], Automatic}]

enter image description here

I'm used the second equation is used when standard temperature lapse rate equals zero.

ClearAll["Global`*"]

Pb0 = 101325  (*Pa*);
Pb1 = 22632.10(*Pa*);
Pb2 = 5474.89 (*Pa*);
Pb3 = 868.02  (*Pa*);
Pb4 = 110.91  (*Pa*);
Pb5 = 66.94   (*Pa*);
Pb6 = 3.96    (*Pa*);
Tb0 = 288.15  (*K*);
Tb1 = 216.65  (*K*);
Tb2 = 216.65  (*K*);
Tb3 = 228.65  (*K*);
Tb4 = 270.65  (*K*);
Tb5 = 270.65  (*K*);
Tb6 = 214.15  (*K*);
g0 = 9.860665 (*m/s^2*);
M = 0.0289644 (*kg/mol*);
hb0 = 0       (*m*);
hb1 = 11000   (*m*);
hb2 = 20000   (*m*);
hb3 = 32000   (*m*);
hb4 = 47000   (*m*);
hb5 = 51000   (*m*);
hb6 = 71000   (*m*);
R = 8.3144598 (*J/mol/K*);

P0[h_] := Pb0*Exp[-g0*M*(h - hb0)/(R*Tb0)];
P1[h_] := Pb1*Exp[-g0*M*(h - hb1)/(R*Tb1)];
P2[h_] := Pb2*Exp[-g0*M*(h - hb2)/(R*Tb2)];
P3[h_] := Pb3*Exp[-g0*M*(h - hb3)/(R*Tb3)];
P4[h_] := Pb4*Exp[-g0*M*(h - hb4)/(R*Tb4)];
P5[h_] := Pb5*Exp[-g0*M*(h - hb5)/(R*Tb5)];
P6[h_] := Pb6*Exp[-g0*M*(h - hb6)/(R*Tb6)];

PP[h_] := Piecewise[{{P0[h], h <= 11000}, {P1[h], 11000 <= h <= 20000}, {P2[h],
 20000 <= h <= 32000}, {P3[h], 32000 <= h <= 47000}, {P4[h], 
 47000 <= h <= 51000}, {P5[h], 51000 <= h <= 71000}, {P6[h], 71000 <= h <= 86000}}]

 Plot[Labeled[{PP[h]}, "Pressure as a function of the height above the sea level", 
 5000], {h, 0, 51000}, PlotRange -> All, PlotStyle -> Red, AxesLabel -> {"h[m]", "P[Pa]"}]

 Plot[Labeled[{PP[h]/100}, "Pressure as a function of the height above the sea level", 
 4000], {h, 0, 21000}, PlotRange -> {Automatic, {0, 1013.25}}, 
 PlotStyle -> Red, AxesLabel -> {"h[km]", "P[hPa]"}, 
 Ticks -> {Table[{1000 i, i}, {i, 1, 21, 2}], Automatic}]

enter image description here

We can use equations from Practical Meteorology

p1 = 1013.258*(288.15/(288.15 - 6.5*h))^(-5.255877);(* hPa *)
p2 = 226.32*Exp[-0.1577*(h - 11)];
p3 = 54.749*(216.65/(216.65 + 1*(h - 20)))^34.16319;
p4 = 8.868*(228.65/(228.65 + 2.8*(h - 32)))^12.2011;
p5 = 1.1109*Exp[-0.1262*(h - 47)];
P[h_] := Piecewise[{{p1, h <= 11}, {p2, 11 <= h <= 20}, {p3, 
20 <= h <= 32}, {p4, 32 <= h <= 47}, {p5, 47 <= h <= 51}}]
Plot[Labeled[P[h], "Pressure as a function of the height above the sea level", 8], {h, 
0, 51}, PlotRange -> All, PlotStyle -> Red, AxesLabel -> {h, P}] 

enter image description here

EDITED: 15.10.2017.

U.S. Standard Atmosphere-1976 year from 0 to 1000 km.

p1 = 1013.25*(288.15/(288.15 - 6.5*h))^(-5.255877);
p2 = 226.32*Exp[-0.1577*(h - 11)];
p3 = 54.749*(216.65/(216.65 + 1*(h - 20)))^34.16319;
p4 = 8.868*(228.65/(228.65 + 2.8*(h - 32)))^12.2011;
p5 = 1.1109*Exp[-0.1262*(h - 47)];
p6 = 0.6693887*(270.65/(270.65 - 2.8*(h - 51)))^-12.2011;
p7 = 0.03956420*(214.65/(214.65 - 2*(h - 71)))^-17.0816;
p8 = Exp[2.159582*10^-6*h^3 - 4.836957*10^-4*h^2 - 0.1425192*h + 13.47530];
p9 = Exp[3.304895*10^-5*h^3 - 0.00906273*h^2 + 0.6516698*h - 11.03037];
p10 = Exp[6.693926*10^-5*h^3 - 0.01945388*h^2 + 1.71908*h - 47.75030];
p11 = Exp[-6.539316*10^-5*h^3 + 0.02485568*h^2 - 3.223620*h + 
135.9355];
p12 = Exp[2.283506*10^-7*h^4 - 1.343221*10^-4*h^3 + 0.02999016*h^2 - 
3.055446*h + 113.5764];
p13 = Exp[1.20943*10^-8*h^4 - 9.692458*10^-6*h^3 + 0.003002041*h^2 - 
0.4523015*h + 19.19151];
p14 = Exp[8.113942*10^-10*h^4 - 9.822568*10^-7*h^3 + 4.687616*10^-4*h^2 - 
0.1231710*h + 3.067409];
p15 = Exp[9.814674*10^-11*h^4 - 1.654439*10^-7*h^3 + 1.148115*10^-4*h^2 - 
0.05431334*h - 2.011365];
p16 = Exp[-7.835161*10^-11*h^4 + 1.96489*10^-7*h^3 - 
1.657213*10^-4*h^2 + 0.04305869*h - 14.77132];
p17 = Exp[2.813255*10^-11*h^4 - 1.120689*10^-7*h^3 + 1.695568*10^-4*h^2 - 
0.1188941*h + 14.56718];

P2[h_] := Piecewise[{{p1, h <= 11}, {p2, 11 <= h <= 20}, {p3, 
20 <= h <= 32}, {p4, 32 <= h <= 47}, {p5, 47 <= h <= 51}, {p6, 
51 <= h <= 71}, {p7, 71 <= h <= 86}, {p8, 86 <= h <= 91}, {p9, 
91 <= h <= 100}, {p10, 100 <= h <= 110}, {p11, 
110 <= h <= 120}, {p12, 120 <= h <= 150}, {p13, 
150 <= h <= 200}, {p14, 200 <= h <= 300}, {p15, 
300 <= h <= 500}, {p16, 500 <= h <= 750}, {p17, 750 <= h <= 1000}}]

Plot[Labeled[P2[h],"Pressure as a function of the height above the sea level", 
1.5], {h, 0, 1000}, ScalingFunctions -> {"Log", "Log"}, 
PlotRange -> All, PlotStyle -> Red, AxesLabel -> {"h[km]", "P[hPa]"}]

enter image description here

We can see discontinuity it is caused by two different models.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.