The standard temperature¶

The temperature of a world's standard atmosphere at an altitude, its layers' piecewise-linear profile: lusine.laws.layer_temperature over the layer the altitude is in, the highest of the world's rows of lusine.atmosphere_layers whose base is at or below it, the lowest also below its base.

One ordered fold over the world's layers in their order (world.atmosphere_layers), its value the last layer the where passes. A world with no layers (an airless body's) gives the fold's start, its atmosphere's sea-level temperature: none. The altitude is the one its layers are given in: ISO 2533's geopotential for the Earth, to its row's top. Source: ISO 2533:1975, Standard Atmosphere, section 2.3 and Table 4; U.S. Standard Atmosphere, 1976 (NOAA-S/T 76-1562), eq. 23 and Table 4. Checked against a reference in numpy, which has its own tests.

Source

ISO 2533:1975, Standard Atmosphere, section 2.3 and Table 4; U.S. Standard Atmosphere, 1976 (NOAA-S/T 76-1562), eq. 23 and Table 4.

Proof

reference, invariant, limit.

standard_temperature.trait

namespace lusine.laws;

float<K> standard_temperature(relation<lusine.bodies> world, float<m> altitude)
{
  let layers = world.atmosphere_layers;
  return fold(t = world.atmospheres.sea_level_temperature)
           lusine.laws.layer_temperature(altitude, layers.base, layers.base_temperature, layers.gradient)
         where layers.layer == 0 || layers.base <= altitude;
}

The temperature of a world's standard atmosphere at an altitude, its layers' piecewise-linear profile, T = Tb + a (h - hb) in the layer the altitude is in. The kernel takes the world, a row of lusine.bodies, and folds over that world's rows of lusine.atmosphere_layers (world.atmosphere_layers) in one ordered fold whose value is the last layer the altitude reaches: one call, no host and no layer passed in. Every number is a row's; the kernel declares none.

In [1]:

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 trait_run
from lusine import rows
F32 = np.float32

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 trait_run from lusine import rows F32 = np.float32

The kernel¶

The law as the library states it, and as this page runs it.

In [2]:

print(trait_run.source('laws/standard_temperature.trait'))

print(trait_run.source('laws/standard_temperature.trait'))

// lusine.laws.standard_temperature: the temperature of a world's standard atmosphere at an altitude, its layers'
// piecewise-linear profile: lusine.laws.layer_temperature over the layer the altitude is in, the highest of the
// world's rows of lusine.atmosphere_layers whose base is at or below it, the lowest also below its base. One ordered
// fold over the world's layers in their order (world.atmosphere_layers), its value the last layer the where passes.
// A world with no layers (an airless body's) gives the fold's start, its atmosphere's sea-level temperature: none.
// The altitude is the one its layers are given in: ISO 2533's geopotential for the Earth, to its row's top.
// Source: ISO 2533:1975, Standard Atmosphere, section 2.3 and Table 4; U.S. Standard Atmosphere, 1976
// (NOAA-S/T 76-1562), eq. 23 and Table 4.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.laws;

float<K> standard_temperature(relation<lusine.bodies> world, float<m> altitude)
{
  let layers = world.atmosphere_layers;
  return fold(t = world.atmospheres.sea_level_temperature)
           lusine.laws.layer_temperature(altitude, layers.base, layers.base_temperature, layers.gradient)
         where layers.layer == 0 || layers.base <= altitude;
}

The program¶

A lattice of items, each the host's world (a relation to its body's row) and altitude, and main() binding the kernel over them. temperature(world, altitudes) runs it once for a world's altitudes, geopotential metres for the Earth as ISO 2533 tabulates them.

In [3]:

AIR = """namespace air;

lusine.version version = { { number = 1 } };

struct item { host relation<lusine.bodies> world; host float<m> altitude; float<K> temperature; };
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.temperature = lusine.laws.standard_temperature(items.world, items.altitude); }
"""
print(AIR)


def temperature(world, altitudes):
    altitudes = np.atleast_1d(np.asarray(altitudes, dtype=float))
    run = trait_run.Program(['laws/standard_temperature.trait'], AIR, {'air.items': [len(altitudes)]})
    run.store('air.items', 'world', run.row('lusine.bodies', world))
    run.store('air.items', 'altitude', altitudes)
    return run.steps('air.main', 1, [('air.items', 'temperature')])[0][0], run

AIR = """namespace air; lusine.version version = { { number = 1 } }; struct item { host relation<lusine.bodies> world; host float<m> altitude; float<K> temperature; }; lattice<item, 1> items(spacing = 1 m, boundary = wall); void main() { items.temperature = lusine.laws.standard_temperature(items.world, items.altitude); } """ print(AIR) def temperature(world, altitudes): altitudes = np.atleast_1d(np.asarray(altitudes, dtype=float)) run = trait_run.Program(['laws/standard_temperature.trait'], AIR, {'air.items': [len(altitudes)]}) run.store('air.items', 'world', run.row('lusine.bodies', world)) run.store('air.items', 'altitude', altitudes) return run.steps('air.main', 1, [('air.items', 'temperature')])[0][0], run

namespace air;

lusine.version version = { { number = 1 } };

struct item { host relation<lusine.bodies> world; host float<m> altitude; float<K> temperature; };
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.temperature = lusine.laws.standard_temperature(items.world, items.altitude); }

Reference: ISO 2533, every layer's base and points inside¶

ISO 2533's temperature (the U.S. Standard Atmosphere 1976's, Table 4, which agrees) at each of the seven layers' bases and the seventh's top, 84.852 km, and at points inside the layers and 2 km below sea level, against the kernel's, each within binary32's rounding of the standard's value (its last place is 0.001 K; a binary32 near 300 K is 3e-5 K apart).

In [4]:

ISO = {-2000: 301.15, 0: 288.15, 1000: 281.65, 5000: 255.65, 8000: 236.15, 11000: 216.65, 15000: 216.65,
       20000: 216.65, 25000: 221.65, 32000: 228.65, 40000: 251.05, 47000: 270.65, 51000: 270.65, 60000: 245.45,
       71000: 214.65, 80000: 196.65, 84852: 186.946}
hs = np.array(sorted(ISO), dtype=float)
T, _ = temperature('earth', hs)
for h, t in zip(hs, T):
    print(f'{h:8.0f} m  ISO {ISO[int(h)]:8.3f} K   kernel {t:10.5f} K  ({t - ISO[int(h)]:+.1e})')
assert np.max(np.abs(T - np.array([ISO[int(h)] for h in hs]))) < 1e-4

ISO = {-2000: 301.15, 0: 288.15, 1000: 281.65, 5000: 255.65, 8000: 236.15, 11000: 216.65, 15000: 216.65, 20000: 216.65, 25000: 221.65, 32000: 228.65, 40000: 251.05, 47000: 270.65, 51000: 270.65, 60000: 245.45, 71000: 214.65, 80000: 196.65, 84852: 186.946} hs = np.array(sorted(ISO), dtype=float) T, _ = temperature('earth', hs) for h, t in zip(hs, T): print(f'{h:8.0f} m ISO {ISO[int(h)]:8.3f} K kernel {t:10.5f} K ({t - ISO[int(h)]:+.1e})') assert np.max(np.abs(T - np.array([ISO[int(h)] for h in hs]))) < 1e-4

   -2000 m  ISO  301.150 K   kernel  301.14999 K  (-6.1e-06)
       0 m  ISO  288.150 K   kernel  288.14999 K  (-6.1e-06)
    1000 m  ISO  281.650 K   kernel  281.64999 K  (-6.1e-06)
    5000 m  ISO  255.650 K   kernel  255.64999 K  (-6.1e-06)
    8000 m  ISO  236.150 K   kernel  236.14999 K  (-6.1e-06)
   11000 m  ISO  216.650 K   kernel  216.64999 K  (-6.1e-06)
   15000 m  ISO  216.650 K   kernel  216.64999 K  (-6.1e-06)
   20000 m  ISO  216.650 K   kernel  216.64999 K  (-6.1e-06)
   25000 m  ISO  221.650 K   kernel  221.64999 K  (-6.1e-06)
   32000 m  ISO  228.650 K   kernel  228.64999 K  (-6.1e-06)
   40000 m  ISO  251.050 K   kernel  251.04999 K  (-1.2e-05)
   47000 m  ISO  270.650 K   kernel  270.64999 K  (-6.1e-06)
   51000 m  ISO  270.650 K   kernel  270.64999 K  (-6.1e-06)
   60000 m  ISO  245.450 K   kernel  245.45000 K  (-3.1e-06)
   71000 m  ISO  214.650 K   kernel  214.64999 K  (-6.1e-06)
   80000 m  ISO  196.650 K   kernel  196.64999 K  (-6.1e-06)
   84852 m  ISO  186.946 K   kernel  186.94600 K  (-8.5e-07)

Reference: bit for bit, the profile in binary32¶

Over 4,000 altitudes from 2 km below sea level to the top, on the Earth and on Mars, the kernel's temperature is a binary32 model's bit for bit: the layer chosen as the law says from the rows the program holds (the cells as binary32 stores them), and Tb + a (h - hb) rounded after each operation, as Trait's floats are. The numpy reference in float64 agrees within 2e-4 K: binary32's rounding of the altitude (an altitude near 85 km is 8 mm apart, times a gradient of 0.0028 K/m), the cells and the result.

In [5]:

import sys
sys.path.insert(0, os.path.join(root, 'kernels'))
from lusine.kernels.laws.standard_temperature import standard_temperature
for world, top in (('earth', 84852.0), ('mars', 40000.0)):
    hs = F32(np.linspace(-2000.0, top, 4001))
    T, run = temperature(world, hs)
    layers = [r for r in rows.rows(run.unit, 'lusine.atmosphere_layers') if r['key'].startswith(world + ', ')]
    model = np.empty(len(hs), dtype=F32)
    for i, h in enumerate(hs):
        r = layers[max(j for j, r in enumerate(layers) if F32(r['base']) <= h or j == 0)]
        model[i] = F32(r['base_temperature']) + F32(F32(r['gradient']) * F32(h - F32(r['base'])))
    assert np.array_equal(F32(T), model), world
    worst = float(np.max(np.abs(T - standard_temperature(hs.astype(float), body=world))))
    print(f'{world}: {len(hs)} altitudes bit for bit with the binary32 model; the float64 reference within '
          f'{worst:.1e} K')
    assert worst < 2e-4

import sys sys.path.insert(0, os.path.join(root, 'kernels')) from lusine.kernels.laws.standard_temperature import standard_temperature for world, top in (('earth', 84852.0), ('mars', 40000.0)): hs = F32(np.linspace(-2000.0, top, 4001)) T, run = temperature(world, hs) layers = [r for r in rows.rows(run.unit, 'lusine.atmosphere_layers') if r['key'].startswith(world + ', ')] model = np.empty(len(hs), dtype=F32) for i, h in enumerate(hs): r = layers[max(j for j, r in enumerate(layers) if F32(r['base']) <= h or j == 0)] model[i] = F32(r['base_temperature']) + F32(F32(r['gradient']) * F32(h - F32(r['base']))) assert np.array_equal(F32(T), model), world worst = float(np.max(np.abs(T - standard_temperature(hs.astype(float), body=world)))) print(f'{world}: {len(hs)} altitudes bit for bit with the binary32 model; the float64 reference within ' f'{worst:.1e} K') assert worst < 2e-4

earth: 4001 altitudes bit for bit with the binary32 model; the float64 reference within 2.1e-05 K
mars: 4001 altitudes bit for bit with the binary32 model; the float64 reference within 2.1e-05 K

Invariant: one line a layer, meeting at the bases¶

Inside each of the Earth's layers the profile's slope is that layer's gradient (6.5 K a kilometre falling in the troposphere, none in the tropopause), and at every base the layers meet: the profile is continuous. Mars's fit steps by 0.95 K at 7 km, kept as its source gives it (lusine.atmospheres).

In [6]:

earth = [r for r in rows.rows(run.unit, 'lusine.atmosphere_layers') if r['key'].startswith('earth, ')]
bases = [r['base'] for r in earth]
for r, top in zip(earth, bases[1:] + [84852.0]):
    (lo, hi), _ = temperature('earth', [r['base'] + 500.0, top - 500.0])
    slope = (hi - lo) / (top - r['base'] - 1000.0)
    assert abs(slope - r['gradient']) < 1e-6, (r['key'], slope)
    print(f"{r['key']:10s} from {r['base']:7.0f} m: {slope * 1000:+.4f} K/km, the layer's {r['gradient'] * 1000:+.4f}")
below, _ = temperature('earth', [b - 0.01 for b in bases[1:]])
above, _ = temperature('earth', bases[1:])
assert np.max(np.abs(below - above)) < 1e-3
line = np.linspace(-2000.0, 84852.0, 800)
T, _ = temperature('earth', line)
M, _ = temperature('mars', np.linspace(0.0, 40000.0, 400))
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(T, line / 1000, label='Earth (ISO 2533)')
ax.plot(M, np.linspace(0.0, 40.0, 400), label='Mars')
ax.plot([ISO[h] for h in sorted(ISO)], np.array(sorted(ISO)) / 1000, 'o', ms=3, label='ISO 2533')
ax.set_xlabel('temperature (K)')
ax.set_ylabel('altitude (km)')
ax.legend()
plt.show()

earth = [r for r in rows.rows(run.unit, 'lusine.atmosphere_layers') if r['key'].startswith('earth, ')] bases = [r['base'] for r in earth] for r, top in zip(earth, bases[1:] + [84852.0]): (lo, hi), _ = temperature('earth', [r['base'] + 500.0, top - 500.0]) slope = (hi - lo) / (top - r['base'] - 1000.0) assert abs(slope - r['gradient']) < 1e-6, (r['key'], slope) print(f"{r['key']:10s} from {r['base']:7.0f} m: {slope * 1000:+.4f} K/km, the layer's {r['gradient'] * 1000:+.4f}") below, _ = temperature('earth', [b - 0.01 for b in bases[1:]]) above, _ = temperature('earth', bases[1:]) assert np.max(np.abs(below - above)) < 1e-3 line = np.linspace(-2000.0, 84852.0, 800) T, _ = temperature('earth', line) M, _ = temperature('mars', np.linspace(0.0, 40000.0, 400)) fig, ax = plt.subplots(figsize=(7, 4)) ax.plot(T, line / 1000, label='Earth (ISO 2533)') ax.plot(M, np.linspace(0.0, 40.0, 400), label='Mars') ax.plot([ISO[h] for h in sorted(ISO)], np.array(sorted(ISO)) / 1000, 'o', ms=3, label='ISO 2533') ax.set_xlabel('temperature (K)') ax.set_ylabel('altitude (km)') ax.legend() plt.show()

earth, 0   from       0 m: -6.5000 K/km, the layer's -6.5000
earth, 1   from   11000 m: +0.0000 K/km, the layer's +0.0000
earth, 2   from   20000 m: +1.0000 K/km, the layer's +1.0000
earth, 3   from   32000 m: +2.8000 K/km, the layer's +2.8000
earth, 4   from   47000 m: +0.0000 K/km, the layer's +0.0000
earth, 5   from   51000 m: -2.8000 K/km, the layer's -2.8000
earth, 6   from   71000 m: -2.0000 K/km, the layer's -2.0000

Limit: sea level, and below it¶

At zero altitude the kernel gives the world's sea-level temperature (lusine.atmospheres) exactly, and below it the lowest layer goes on: 2 km down on the Earth is 13 K warmer.

In [7]:

for world in ('earth', 'mars'):
    (t0, t_down), run = temperature(world, [0.0, -2000.0])
    air = rows.row(run.unit, 'lusine.atmospheres', world)
    assert t0 == air['sea_level_temperature']
    print(f"{world}: {t0} K at its datum, {air['sea_level_temperature']} K its row's; {t_down:.2f} K 2 km below")
assert abs(temperature('earth', [-2000.0])[0][0] - 301.15) < 1e-4

for world in ('earth', 'mars'): (t0, t_down), run = temperature(world, [0.0, -2000.0]) air = rows.row(run.unit, 'lusine.atmospheres', world) assert t0 == air['sea_level_temperature'] print(f"{world}: {t0} K at its datum, {air['sea_level_temperature']} K its row's; {t_down:.2f} K 2 km below") assert abs(temperature('earth', [-2000.0])[0][0] - 301.15) < 1e-4

earth: 288.1499938964844 K at its datum, 288.1499938964844 K its row's; 301.15 K 2 km below
mars: 242.14999389648438 K at its datum, 242.14999389648438 K its row's; 244.15 K 2 km below