# A graph of various atmospheric values vs. height.

# Conventions:
#   "0"     - At height 0 meters.
#   "rel"   - Relative. Scaled so 1 at height 0 meters.
#   "s"     - Substitute specific values.
#   "dX_dY" - Derivative of X with respect to Y.

# Graph the pressure and density formulas from
#   https://en.wikipedia.org/wiki/Density_of_air

# See above URL for these variables.
var("p0 h T0 g L R M")

# The pressure (p) with respect to height (h).
p = p0*(1 - (L*h)/T0)^((g*M)/(R*L))
p_rel = p/p0 # Relative. 1 at a height of 0.

# Substitute specific values for variables from the above URL.
subs = {p0: 101.325, T0:288.15, g:9.80665, L:0.0065, R:8.31447, M:0.0289644}
p_rel_s = p_rel.subs(subs)

# The derivative of p with respect to h.
dp_dh = diff(p, h)
dp_dh_0 = dp_dh.subs({h:0}) # Relative. 1 at a height of 0.
dp_dh_rel = dp_dh/dp_dh_0
dp_dh_rel_s = dp_dh_rel.subs(subs)

# Formula for density given the above and the URL.
T = T0 - L*h
d = (M/(R*T))*p
d_0 = d.subs({h:0})
d_rel = d/d_0
d_rel_s = d_rel.subs(subs)

# The derivative of d with respect to h.
dd_dh = diff(d, h)
dd_dh_0 = dd_dh.subs({h:0}) # Relative. 1 at a height of 0.
dd_dh_rel = dd_dh/dd_dh_0
dd_dh_rel_s = dd_dh_rel.subs(subs)

# A straight line to compare to the density slope. Assume that the density
# slope goes to 0 at 18 km up. This should be accurate += 0.05 for most
# flights.
d_slope_zero_height = 18000
dsl = 1 - h/d_slope_zero_height

r = (h, 0, 15000) # The range to plot.
plots = \
  plot(    d_rel_s, r, legend_label="density"              , color="red") + \
  plot(dd_dh_rel_s, r, legend_label="density slope"        , color="orange") + \
  plot(        dsl, r, legend_label="density slope linear" , color="magenta") + \
  plot(    p_rel_s, r, legend_label="pressure"             , color="green") + \
  plot(dp_dh_rel_s, r, legend_label="pressure slope"       , color="red")

# Save to the /tmp directory. Adjust for your operating system, or use
# ".show()", but delete the path.
plots.save("/tmp/atmosphere-vs-height.png", axes_labels=["Height (m)",
  "Relative to zero height"], figsize=16, frame=True, gridlines=True,
  title="Atmosphere vs. height")