Gravity at an altitude: the inverse square of the distance from the centre¶

The acceleration of free fall at an altitude above a body's surface, g = GM / (R + h)^2, from its gravitational parameter, its radius (lusine.bodies' gm and radius) and the altitude: the inverse square of the distance from its centre, lusine.sky.surface_gravity at h = 0.

The body is spherically symmetric and its rotation left out, as there. 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, limit, invariant.

gravity_at.trait

namespace lusine.sky;

float<m/s2> gravity_at(float<m3/s2> gm, float<m> radius, float<m> altitude)
{
  let distance = radius + altitude;
  return gm / (distance * distance);
}

The acceleration of free fall at an altitude h above a body's surface, g = GM/(R + h)², from its row of lusine.bodies; lusine.sky.surface_gravity at h = 0.

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.library_sources import parse_rows
EPS = float(np.finfo(np.float32).eps)   # Trait's float is binary32
from lusine.library_sources import bodies

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.

In [2]:

print(trait_run.source('sky/gravity_at.trait'))

print(trait_run.source('sky/gravity_at.trait'))

// lusine.sky.gravity_at: the acceleration of free fall at an altitude above a body's surface, g = GM / (R + h)^2, from
// its gravitational parameter, its radius (lusine.bodies' gm and radius) and the altitude: the inverse square of the
// distance from its centre, lusine.sky.surface_gravity at h = 0. The body is spherically symmetric and its rotation
// left out, as there.
// 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> gravity_at(float<m3/s2> gm, float<m> radius, float<m> altitude)
{
  let distance = radius + altitude;
  return gm / (distance * distance);
}

The program¶

Each section runs the kernel over a lattice of items, one a case.

In [3]:

print(trait_run.law_program('sky/gravity_at.trait', 'gravity_at'))

print(trait_run.law_program('sky/gravity_at.trait', 'gravity_at'))

namespace probe;

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

struct item {
  host float<m3/s2> gm;
  host float<m> radius;
  host float<m> altitude;
  float<m/s2> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.sky.gravity_at(items.gm, items.radius, items.altitude); }

Limit: at the surface it is the surface gravity¶

At h = 0, R + h is R exactly, so the kernel gives surface_gravity's bits for the Sun, the Earth and the Moon.

In [4]:

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])
at0 = trait_run.law('sky/gravity_at.trait', 'gravity_at', gm=gm, radius=r, altitude=0.0)
surface = trait_run.law('sky/surface_gravity.trait', 'surface_gravity', gm=gm, radius=r)
assert np.array_equal(at0, surface)
print(dict(zip(names, at0)))

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]) at0 = trait_run.law('sky/gravity_at.trait', 'gravity_at', gm=gm, radius=r, altitude=0.0) surface = trait_run.law('sky/surface_gravity.trait', 'surface_gravity', gm=gm, radius=r) assert np.array_equal(at0, surface) print(dict(zip(names, at0)))

{'sun': np.float64, 'earth': np.float64, 'moon': np.float64}

Reference: the Earth's gravity in orbit¶

From the Earth's row: 8.68 m/s² at the ISS's 400 km, 89% of the surface's (an astronaut floats because the station falls, not because gravity is gone), and 0.224 m/s² at geostationary altitude, 35 786 km; each held to GM/(R + h)², to binary32's rounding, and g(R + h)² is GM at every altitude to 1000 Earth radii.

In [5]:

earth = rows['earth']
gm, r = earth['gm'], earth['radius']
hs = np.array([0.0, 400e3, 35786e3])
g = trait_run.law('sky/gravity_at.trait', 'gravity_at', gm=gm, radius=r, altitude=hs)
assert abs(g[1] - 8.68) < 0.01 and abs(g[2] - 0.2242) < 1e-3, g
print(f'ISS {g[1]:.3f} m/s2 ({g[1] / g[0]:.1%} of the surface), geostationary {g[2]:.4f} m/s2')
hs = np.float32(np.geomspace(1.0, 1000 * r, 64)).astype(float)
g = trait_run.law('sky/gravity_at.trait', 'gravity_at', gm=gm, radius=r, altitude=hs)
d = np.float32(np.float32(r) + np.float32(hs)).astype(float)
assert np.all(np.abs(g * d * d / np.float32(gm) - 1) <= 4 * EPS)
assert np.all(np.diff(g) < 0)
fig, ax = plt.subplots(figsize=(7, 3.2))
for b in ('earth', 'moon'):
    h = np.linspace(0.0, 4 * rows[b]['radius'], 201)
    ax.plot(h / 1e3, trait_run.law('sky/gravity_at.trait', 'gravity_at', gm=rows[b]['gm'], radius=rows[b]['radius'], altitude=h), label=b)
ax.set_xlabel('altitude (km)')
ax.set_ylabel('g (m/s2)')
ax.legend()
plt.show()

earth = rows['earth'] gm, r = earth['gm'], earth['radius'] hs = np.array([0.0, 400e3, 35786e3]) g = trait_run.law('sky/gravity_at.trait', 'gravity_at', gm=gm, radius=r, altitude=hs) assert abs(g[1] - 8.68) < 0.01 and abs(g[2] - 0.2242) < 1e-3, g print(f'ISS {g[1]:.3f} m/s2 ({g[1] / g[0]:.1%} of the surface), geostationary {g[2]:.4f} m/s2') hs = np.float32(np.geomspace(1.0, 1000 * r, 64)).astype(float) g = trait_run.law('sky/gravity_at.trait', 'gravity_at', gm=gm, radius=r, altitude=hs) d = np.float32(np.float32(r) + np.float32(hs)).astype(float) assert np.all(np.abs(g * d * d / np.float32(gm) - 1) <= 4 * EPS) assert np.all(np.diff(g) < 0) fig, ax = plt.subplots(figsize=(7, 3.2)) for b in ('earth', 'moon'): h = np.linspace(0.0, 4 * rows[b]['radius'], 201) ax.plot(h / 1e3, trait_run.law('sky/gravity_at.trait', 'gravity_at', gm=rows[b]['gm'], radius=rows[b]['radius'], altitude=h), label=b) ax.set_xlabel('altitude (km)') ax.set_ylabel('g (m/s2)') ax.legend() plt.show()

ISS 8.676 m/s2 (88.5% of the surface), geostationary 0.2242 m/s2

Invariant: one radius up, a quarter¶

At h = R the distance from the centre doubles and the gravity is a quarter of the surface's, bit for bit.

In [6]:

g = trait_run.law('sky/gravity_at.trait', 'gravity_at', gm=gm, radius=r, altitude=r)
assert g[0] == at0[names.index('earth')] / 4
print(g)

g = trait_run.law('sky/gravity_at.trait', 'gravity_at', gm=gm, radius=r, altitude=r) assert g[0] == at0[names.index('earth')] / 4 print(g)

[2.4495995]