Surface gravity: GM over the radius squared¶
The acceleration of free fall at a body's surface, g = GM / R^2, from its gravitational parameter and its radius (lusine.bodies' gm and radius): Newton's gravitation of a spherically symmetric body, its rotation left out (a body's weight at its equator is less by omega^2 R, 0.034 m/s2 on the Earth).
The Earth's is 9.798 m/s2; the standard gravity, 9.80665 m/s2, is the convention lusine.conventions.g_n, not this. Source: Murray & Dermott, Solar System Dynamics, 1.2; Newton, Principia, book I, prop. 71. Checked against a reference in numpy, which has its own tests.
- Source
- Murray & Dermott, Solar System Dynamics, 1.2; Newton, Principia, book I, prop. 71.
- Proof
- reference, invariant.
surface_gravity.trait
namespace lusine.sky;
float<m/s2> surface_gravity(float<m3/s2> gm, float<m> radius)
{
return gm / (radius * radius);
}
The acceleration of free fall at a body's surface, g = GM/R², from its row of lusine.bodies, its gravitational parameter and its radius; the body spherically symmetric and its rotation left out. lusine.bodies stores no surface gravity its GM and radius give: this kernel derives it (Mars's is kept by hand while no pinned source gives its radius). The Earth's is 9.798 m/s², not the standard gravity 9.80665 m/s², which is a convention, lusine.conventions.g_n.
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.library_sources import parse_rows
EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32
from lusine.library_sources import bodies
The kernel¶
The law as the library states it.
print(trait_run.source('sky/surface_gravity.trait'))
// lusine.sky.surface_gravity: the acceleration of free fall at a body's surface, g = GM / R^2, from its gravitational
// parameter and its radius (lusine.bodies' gm and radius): Newton's gravitation of a spherically symmetric body, its
// rotation left out (a body's weight at its equator is less by omega^2 R, 0.034 m/s2 on the Earth). The Earth's is
// 9.798 m/s2; the standard gravity, 9.80665 m/s2, is the convention lusine.conventions.g_n, not this.
// Source: Murray & Dermott, Solar System Dynamics, 1.2; Newton, Principia, book I, prop. 71.
// Checked against a reference in numpy, which has its own tests.
namespace lusine.sky;
float<m/s2> surface_gravity(float<m3/s2> gm, float<m> radius)
{
return gm / (radius * radius);
}
The program¶
Each section runs the kernel over a lattice of items, one a case.
print(trait_run.law_program('sky/surface_gravity.trait', 'surface_gravity'))
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host float<m3/s2> gm;
host float<m> radius;
float<m/s2> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.sky.surface_gravity(items.gm, items.radius); }
Reference: the Sun, the Earth and the Moon, from their rows¶
Each body whose row holds a GM and a radius (lusine.bodies), held to GM/R² of the two, to binary32's rounding: the Earth's 9.798 m/s² from the IAU's nominal GM and equatorial radius, the Moon's 1.624 m/s² from its mean radius, the Sun's 274.2 m/s². The Earth's is 0.008 m/s² short of the standard gravity: g_n is a convention agreed in 1901, near the gravity at 45° latitude with the rotation in it, and not the Earth's GM/R².
rows = bodies.rows()
names = [b for b in rows if 'gm' in rows[b] and 'radius' in rows[b]]
gm = np.array([rows[b]['gm'] for b in names])
r = np.array([rows[b]['radius'] for b in names])
g = trait_run.law('sky/surface_gravity.trait', 'surface_gravity', gm=gm, radius=r)
want = np.float32(gm).astype(float) / np.float32(r).astype(float) ** 2
assert np.all(np.abs(g / want - 1) <= 3 * EPS), np.abs(g / want - 1).max()
at = dict(zip(names, g))
assert abs(at['earth'] - 9.798) < 1e-3 and abs(at['moon'] - 1.624) < 1e-3 and abs(at['sun'] - 274.2) < 0.1, at
import re
g_n = float(re.search(r'float<m/s2> g_n = ([0-9.]+) m/s2;', trait_run.source('conventions.trait')).group(1))
assert 0.008 < g_n - at['earth'] < 0.009
for b in names:
print(f'{b:6s} g = {at[b]:9.4f} m/s2')
print(f"the Earth's GM / R^2 is {g_n - at['earth']:.4f} m/s2 short of lusine.conventions.g_n, {g_n}")
sun g = 274.2001 m/s2 earth g = 9.7984 m/s2 moon g = 1.6242 m/s2 the Earth's GM / R^2 is 0.0083 m/s2 short of lusine.conventions.g_n, 9.80665
Invariant: the inverse square, and the rotation left out¶
Twice the radius gives a quarter of the gravity, bit for bit (a power of two scales without rounding). The kernel leaves the rotation out: at the Earth's equator a body's weight is less by ω²R, 0.034 m/s², taken from the Earth's row's rotation rate and radius, which a game that wants it applies as a centrifugal term of its own.
gm = np.float32(10 ** np.random.default_rng(326).uniform(10, 20, 32)).astype(float)
r = np.float32(10 ** np.random.default_rng(327).uniform(5, 9, 32)).astype(float)
g1 = trait_run.law('sky/surface_gravity.trait', 'surface_gravity', gm=gm, radius=r)
g2 = trait_run.law('sky/surface_gravity.trait', 'surface_gravity', gm=gm, radius=2 * r)
assert np.array_equal(g2, g1 / 4)
earth = rows['earth']
spin = earth['rotation_rate'] ** 2 * earth['radius']
assert abs(spin - 0.0339) < 1e-4
print(f'omega^2 R at the equator: {spin:.4f} m/s2')
omega^2 R at the equator: 0.0339 m/s2