# Plotting pressure as a function of altitude

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

• Meet us halfway. Can you give us the equation representing the barometric formula? Commented Oct 14, 2017 at 16:13
• @J.M. Here: en.wikipedia.org/wiki/Barometric_formula Commented Oct 14, 2017 at 17:18
• @Mariusz, I know perfectly well what the barometric formula looks like, being a chemist. :) I was giving the OP a chance to improve his question by including the formula. In its current form I was sorely tempted to close it, but I decided to be nice. Commented Oct 14, 2017 at 17:23
• Welcome to Mathematica StackExchange. In order to learn how to use this site take the tour. Users will appreciate you typing in your equations using formatted inline code. See the Markdown help on how to do this. Commented Oct 14, 2017 at 18:26
• reference.wolfram.com/language/ref/StandardAtmosphereData.html Commented Oct 14, 2017 at 21:04

## 1 Answer

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}]


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}]


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}]


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]"}]


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