The standard atmosphere¶
The pressure of a world's standard atmosphere at an altitude, the hydrostatic equation over its layers, lusine.laws.layer_pressure through each of the world's rows of lusine.atmosphere_layers the altitude reaches, from its atmosphere's sea-level pressure (lusine.atmospheres) up.
The altitude's layer is the highest of the world's layers (world.atmosphere_layers) whose base is at or below it, the lowest also below its base; an ordered fold over the layers below it, in their order, carries the pressure across each to its top, the next layer's base (next()): the altitude's layer's base pressure, never stored. That layer's law gives the pressure at the altitude. Each layer's scale height is lusine.properties.scale_height at its base temperature, of the air's molar mass, under the gravity the world's standard atmosphere is defined under: a convention, passed in (the Earth's ISO 2533 is computed with the standard gravity, lusine.conventions.g_n; Mars's under its surface gravity). Its temperature is lusine.laws.standard_temperature, and its density lusine.properties.gas_density of the two. The altitude is the one its layers are given in: ISO 2533's geopotential for the Earth, to its row's top. THE GAS CONSTANT: R is lusine.constants.R, CODATA's 8.31446 J/(mol K), where ISO 2533 and the U.S. Standard Atmosphere tabulate with their own R* = 8.31432; the pressures differ by 3e-4 relative at most (the reference's arms: 0.03% at every layer's base, 2e-6 with R* itself). Nothing restates R* here: a law reads one R. Source: ISO 2533:1975, Standard Atmosphere, section 2.4 and Table 4; U.S. Standard Atmosphere, 1976 (NOAA-S/T 76-1562), eqs. 33a and 33b. Checked against a reference in numpy, which has its own tests.
- Source
- ISO 2533:1975, Standard Atmosphere; U.S. Standard Atmosphere, 1976 (NOAA-S/T 76-1562), ch. 1 and Table 4.
- Proof
- reference, curve, invariant, limit.
standard_atmosphere.trait
namespace lusine.laws;
float<Pa> standard_atmosphere(relation<lusine.bodies> world, float<m> altitude, float<m/s2> gravity)
{
let air = world.atmospheres;
let layers = world.atmosphere_layers;
let layer = max(layers.layer) where layers.layer == 0 || layers.base <= altitude;
let base_pressure = fold(p = air.sea_level_pressure)
lusine.laws.layer_pressure(altitude = next(layers.base), base = layers.base, base_pressure = p,
base_temperature = layers.base_temperature, gradient = layers.gradient,
scale_height = lusine.properties.scale_height(temperature = layers.base_temperature,
molar_mass = air.air_molar_mass,
gravity = gravity))
where layers.layer < layer;
return fold(p = base_pressure)
lusine.laws.layer_pressure(altitude, base = layers.base, base_pressure = p,
base_temperature = layers.base_temperature, gradient = layers.gradient,
scale_height = lusine.properties.scale_height(temperature = layers.base_temperature,
molar_mass = air.air_molar_mass,
gravity = gravity))
where layers.layer == layer;
}
The ambient pressure at an altitude: the hydrostatic equation over a world's layers, ISO 2533's seven for the Earth to 86 km. It is the power law where the temperature changes with height and the exponential where it holds, and each layer's base pressure is the layer below's carried across it, derived and never stored.
Every number it runs on is a row's:
- the item's body,
lusine.bodies, whoseatmospherenames its row oflusine.atmospheres; - the gravity its standard atmosphere is defined under: the Earth's ISO 2533 the standard gravity, the convention
lusine.conventions.g_n(not the Earth's own GM/R², 9.798 m/s²,lusine.sky.surface_gravity), Mars's the surface gravity its row keeps; - R,
lusine.constants.R.
The kernel declares no constant.
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import library
from lusine.library_sources import atmospheres, bodies
from lusine.kernels.laws import standard_atmosphere as sa
from lusine.kernels.laws.standard_atmosphere import standard_atmosphere
The reference¶
The kernel's reference, in numpy, which every section below runs. The kernel's Trait form, above, is one kernel over the world's rows of lusine.atmosphere_layers, the gravity passed in; it compiles, and this page will run it once a program can run a fold that reads the world's air through its relation; until then the sections below run the reference. Its temperature is lusine.laws.standard_temperature, proved on its own page.
library.reference(standard_atmosphere)
the reference, standard_atmosphere.py
def standard_atmosphere(altitude, body="earth", unwritten=np.nan):
"""THE STANDARD ATMOSPHERE: the ambient pressure at an altitude, over a
world's layers.
Requires: `Altitude` h (m) above the world's standard datum, from -2 km
to its row's top (geopotential: the body's surface gravity throughout);
the item's `Body` (a row of lusine.bodies, per item), whose
`atmosphere` names a row of lusine.atmospheres or none.
Provides: `AmbientPressure` p (Pa); an item on an airless body is
skipped (here, `unwritten`).
Law: the hydrostatic equation dp/dh = -p g M / (R T) over the layers'
piecewise-linear temperature T = Tb + a (h - hb). In a layer with a
gradient, p = pb (1 + a (h - hb) / Tb)^(-g M / (R a)); in an isothermal
one, p = pb exp(-g M (h - hb) / (R Tb)). Each layer's base pressure pb
is the layer below's carried across it, from the row's sea-level
pressure. M and the layers are the atmosphere's row, g the body's,
R lusine.constants' N_A k. Every power and exponential is lusine.dmath's.
Source: ISO 2533:1975, Standard Atmosphere; U.S. Standard Atmosphere,
1976 (NOAA-S/T 76-1562), ch. 1 and Table 4.
Proof: reference, curve, invariant, limit."""
h = np.asarray(altitude, dtype=np.float64)
names = np.broadcast_to(np.asarray(body, dtype=object), h.shape)
out = np.full(h.shape, unwritten, dtype=np.float64)
for name in set(names.reshape(-1).tolist()):
air = _air(name)
if air is None:
continue # an airless body's items: skipped, the kernel's guard
here = names == name
out[here] = pressure(h[here], air[0], air[1])
return out if out.shape else float(out)
Reference: ISO 2533 at every layer's base¶
The standard's pressure at each of the seven layers' bases and the top, against the kernel's.
The tables were computed with the gas constant of their day, R* = 8.31432 J/(mol K). The kernel uses lusine's R, CODATA's N_A k = 8.314462618, as the law should, so they part by up to about 2e-4 at 85 km.
With R* in its place, the same layers and the same law give the tables to about 2e-7. So the layers and the law are the standard's, and R is the only difference.
hs = np.array(sorted(sa.BOUNDARIES), dtype=float)
want = np.array([sa.BOUNDARIES[int(h)] for h in hs])
got = standard_atmosphere(hs)
g = sa._gravity('earth') # ISO 2533's g0, the standard gravity lusine.conventions.g_n
star = np.array(sa.base_pressures(g, atmospheres.rows()['earth'], sa.R_STAR) + [float(sa.pressure(84852.0, g, atmospheres.rows()['earth'], sa.R_STAR))])
for h, w, p, s in zip(hs, want, got, star):
print(f'{h:8.0f} m ISO {w:12.6g} Pa kernel {p:12.6g} ({p / w - 1:+.1e}) with R* {s:12.6g} ({s / w - 1:+.1e})')
assert np.max(np.abs(got / want - 1)) < 3e-4
assert np.max(np.abs(star / want - 1)) < 2e-6
0 m ISO 101325 Pa kernel 101325 (+0.0e+00) with R* 101325 (+0.0e+00) 11000 m ISO 22632.1 Pa kernel 22632.6 (+2.6e-05) with R* 22632.1 (+1.8e-07) 20000 m ISO 5474.89 Pa kernel 5475.16 (+5.0e-05) with R* 5474.89 (-6.0e-08) 32000 m ISO 868.019 Pa kernel 868.09 (+8.2e-05) with R* 868.019 (-1.8e-08) 47000 m ISO 110.906 Pa kernel 110.919 (+1.2e-04) with R* 110.906 (+5.0e-08) 51000 m ISO 66.9389 Pa kernel 66.9473 (+1.3e-04) with R* 66.9389 (+4.7e-08) 71000 m ISO 3.95642 Pa kernel 3.95711 (+1.7e-04) with R* 3.95642 (+1.1e-07) 84852 m ISO 0.373384 Pa kernel 0.373464 (+2.1e-04) with R* 0.373384 (-2.7e-08)
Curve: inside every layer, the hydrostatic equation itself¶
The closed forms against the hydrostatic equation integrated numerically up the profile, with no closed form used, within 1e-8 inside every layer. Also ISO 2533's tabulated points inside the lower layers, within 0.02%.
earth = atmospheres.rows()['earth']
inside = list(np.linspace(500.0, 84000.0, 85))
ref = sa.integrated(inside, g, earth)
err = max(abs(float(standard_atmosphere(h)) / ref[h] - 1) for h in inside)
print(f'worst against the integration: {err:.1e}')
assert err < 1e-8
tab = np.array(sorted(sa.TABLE), dtype=float)
assert np.max(np.abs(standard_atmosphere(tab) / np.array([sa.TABLE[int(h)] for h in tab]) - 1)) < 2e-4
line = np.linspace(0, 84852, 1200)
fig, ax = plt.subplots(figsize=(7, 4))
ax.semilogx(standard_atmosphere(line), line / 1000, label='kernel')
ax.semilogx(want, hs / 1000, 'o', label='ISO 2533, each layer\'s base')
for b in sa.BASES[1:]:
ax.axhline(b / 1000, lw=0.4, ls='--', color='grey')
ax.set_xlabel('pressure (Pa)')
ax.set_ylabel('geopotential altitude (km)')
ax.legend()
plt.show()
worst against the integration: 2.2e-09
Invariant: the pressure falls with height, and the layers meet¶
Every metre up lowers the pressure, and at every layer's base the two layers' laws agree.
hs = np.linspace(-1000, 84852, 8001)
assert np.all(np.diff(standard_atmosphere(hs)) < 0)
for b in sa.BASES[1:]:
lo, hi = float(standard_atmosphere(b - 1e-6)), float(standard_atmosphere(b + 1e-6))
assert abs(lo / hi - 1) < 1e-9
print(f'{b:7.0f} m: {lo:.6g} Pa below, {hi:.6g} Pa above')
11000 m: 22632.6 Pa below, 22632.6 Pa above 20000 m: 5475.16 Pa below, 5475.16 Pa above 32000 m: 868.09 Pa below, 868.09 Pa above 47000 m: 110.919 Pa below, 110.919 Pa above 51000 m: 66.9473 Pa below, 66.9473 Pa above 71000 m: 3.95711 Pa below, 3.95711 Pa above
Limit: sea level is standard pressure¶
At zero altitude the kernel gives the row's sea-level pressure, exactly.
assert float(standard_atmosphere(0.0)) == sa.P0
print(f'{float(standard_atmosphere(0.0))} Pa')
101325.0 Pa
Worlds: a body per item¶
A game on several worlds supplies each item's body.
- The Earth takes ISO 2533.
- Mars takes its own row: NASA Glenn's two-layer fit and Mars's gravity.
- An item on the Moon is skipped: its body names no atmosphere, and the kernel's guard leaves its pressure unwritten.
h = np.array([0.0, 5000.0, 5000.0, 5000.0, 20000.0])
world = np.array(['earth', 'earth', 'mars', 'moon', 'mars'], dtype=object)
p = standard_atmosphere(h, body=world)
for hh, w, pp in zip(h, world, p):
print(f'{w:6s} {hh:7.0f} m ' + ('skipped: no atmosphere' if np.isnan(pp) else f'{pp:10.4g} Pa'))
assert np.isnan(p[3]) and p[2] < p[1] and p[0] == sa.P0
mars = np.linspace(0, 40000, 400)
fig, ax = plt.subplots(figsize=(7, 3.2))
ax.semilogx(standard_atmosphere(mars, body='earth'), mars / 1000, label='Earth')
ax.semilogx(standard_atmosphere(mars, body='mars'), mars / 1000, label='Mars')
ax.set_xlabel('pressure (Pa)')
ax.set_ylabel('altitude (km)')
ax.legend()
plt.show()
earth 0 m 1.013e+05 Pa earth 5000 m 5.402e+04 Pa mars 5000 m 466.9 Pa moon 5000 m skipped: no atmosphere mars 20000 m 126.1 Pa
The planted fault¶
The reference's planted fault takes every layer as isothermal, dropping the gradients. It misses ISO 2533 by percents at the first layer's top already, which is the page's check refusing it.
bad = sa.WRONG(11000.0)
print(f'at 11 km: {bad:.6g} Pa against ISO 2533\'s {sa.BOUNDARIES[11000]} ({bad / sa.BOUNDARIES[11000] - 1:+.1%})')
assert abs(bad / sa.BOUNDARIES[11000] - 1) > 0.01
at 11 km: 27500.1 Pa against ISO 2533's 22632.06 (+21.5%)