Line of sight#
Eclipse shadows and satellite visibility from a triangle mesh and a bounding-volume tree
A satellite is in eclipse when the straight line from it to the Sun hits the Earth. Two satellites see each other when the line between them does not. Both are the same question: does this ray hit this mesh before it ends?
This page answers it with hawk.math.ray_tri and hawk.geometry. The Earth is a mesh of 5120 triangles, build_bvh puts a tree over it, and closest_hit walks the tree for every lane. Everything is in single precision and in units of the Earth’s radius, so every coordinate is near 1 and the arithmetic keeps its digits.
Colab or Kaggle: run this first cell (it installs the packages and downloads the helper files this notebook uses). Elsewhere it only checks the imports.
import sys, pathlib
sys.path.insert(0, str(pathlib.Path.cwd().parent / "_shared"))
from nb_helpers import gpu_available, to_numpy, xp
The mesh and the tree#
An icosahedron split four times and pushed out to the unit sphere gives 5120 triangles. A triangle’s corners are on the sphere, so the mesh sits a little inside it: the faces sit up to 0.11 % below the sphere, at their lowest. build_bvh splits the triangles at the median of the longest axis until a leaf holds four, and keeps the result as flat tables.
import numpy as np
import eagle, eagle.plan
import hawk
import hawk.math as m
from hawk import Index, Mutable, Scalar, Table, Vector
from hawk.geometry import build_bvh, closest_hit
def icosphere(levels):
t = (1 + 5 ** 0.5) / 2
v = [(-1, t, 0), (1, t, 0), (-1, -t, 0), (1, -t, 0), (0, -1, t), (0, 1, t), (0, -1, -t), (0, 1, -t),
(t, 0, -1), (t, 0, 1), (-t, 0, -1), (-t, 0, 1)]
f = [(0, 11, 5), (0, 5, 1), (0, 1, 7), (0, 7, 10), (0, 10, 11), (1, 5, 9), (5, 11, 4), (11, 10, 2), (10, 7, 6),
(7, 1, 8), (3, 9, 4), (3, 4, 2), (3, 2, 6), (3, 6, 8), (3, 8, 9), (4, 9, 5), (2, 4, 11), (6, 2, 10),
(8, 6, 7), (9, 8, 1)]
v = [np.array(p, float) / np.linalg.norm(p) for p in v]
for _ in range(levels):
cache, nf = {}, []
def mid(a, b):
key = (min(a, b), max(a, b))
if key not in cache:
p = v[a] + v[b]
v.append(p / np.linalg.norm(p))
cache[key] = len(v) - 1
return cache[key]
for a, b, c in f:
ab, bc, ca = mid(a, b), mid(b, c), mid(c, a)
nf += [(a, ab, ca), (b, bc, ab), (c, ca, bc), (ab, bc, ca)]
f = nf
v = np.array(v)
return v[np.array(f)]
earth = icosphere(4)
bvh = build_bvh(earth, leaf=4)
print(f"{len(earth)} triangles, {bvh.n_nodes} nodes, {bvh.depth} levels deep, leaves of at most {bvh.leaf}")
5120 triangles, 4095 nodes, 11 levels deep, leaves of at most 4
One lane: a ray against the tree#
A kernel declares the tree’s seven tables as planes and hands them to closest_hit through bvh.tables(...). The result is [hit, t, u, v, leaf position]: the nearest hit along o + t d for 0 <= t <= tmax, or hit = 0 when there is none. The stack that remembers the nodes still to visit holds depth levels (32 by default); a tree deeper than that is refused when the kernel is traced, rather than traversed with a stack that drops nodes.
The shadow kernel looks from the satellite towards the Sun, d = sun, with no limit on t. The visibility kernel looks from one satellite to the other with d = p2 - p1 and tmax = 1, so the line ends at the second satellite.
@hawk.kernel
def shadow(pos: Vector[3], sun: Vector[3], y: Mutable[Vector[5]],
lo: Table[Scalar], hi: Table[Scalar], left: Table[Index], right: Table[Index],
first: Table[Index], size: Table[Index], tri: Table[Scalar]):
y = closest_hit(bvh.tables(lo, hi, left, right, first, size, tri), pos, sun)
@hawk.kernel
def visible(p1: Vector[3], p2: Vector[3], y: Mutable[Vector[5]],
lo: Table[Scalar], hi: Table[Scalar], left: Table[Index], right: Table[Index],
first: Table[Index], size: Table[Index], tri: Table[Scalar]):
y = closest_hit(bvh.tables(lo, hi, left, right, first, size, tri), p1, p2 - p1, 1.0)
DEVICE = "cuda" if gpu_available() else "host"
print("running on", DEVICE)
def launch(kernel, **planes):
"""Run ``kernel`` over its lanes in single precision, on the GPU when there is one."""
arrays = {k: xp.asarray(np.ascontiguousarray(v, dtype=np.int32 if k in ("left", "right", "first", "size")
else np.float32)) for k, v in {**planes, **bvh.planes(np.float32)}.items()}
n = arrays["pos" if "pos" in arrays else "p1"].shape[-1]
arrays["y"] = xp.zeros((5, n), dtype=xp.float32)
plan = eagle.plan.auto([kernel], targets=(DEVICE,), scalar_type="float32")[0]
if DEVICE == "host":
plan.run(**arrays)
else:
graph = eagle.Graph()
graph.add(plan.bind(**arrays), name="los")
graph.build().launch()
xp.cuda.runtime.deviceSynchronize()
return to_numpy(arrays["y"])
running on cuda
An orbit in and out of the shadow#
The satellite flies a circle of radius 1.1 Earth radii (about 640 km up) tilted 30 degrees from the Sun’s direction. Each of 720 points on the orbit is one lane. For a circular orbit the fraction of the orbit in the cylindrical shadow is acos(sqrt(r^2 - 1) / (r cos(beta))) / pi, with beta the angle between the Sun and the orbit plane, and the mesh should agree with it to the mesh’s own 0.11 % of height.
R_ORBIT, BETA = 1.1, np.radians(30.0)
n_orbit = 720
phase = 2 * np.pi * (np.arange(n_orbit) + 0.5) / n_orbit
sun_dir = np.array([1.0, 0.0, 0.0])
# the orbit plane contains the y axis and is tilted BETA out of the Sun direction
e1 = np.array([np.cos(BETA), 0.0, np.sin(BETA)])
e2 = np.array([0.0, 1.0, 0.0])
orbit = R_ORBIT * (np.outer(np.cos(phase), e1) + np.outer(np.sin(phase), e2)) # (720, 3)
sun_planes = dict(pos=orbit.T, sun=np.tile(sun_dir, (n_orbit, 1)).T)
result = launch(shadow, **sun_planes)
in_shadow = result[0] > 0.5
fraction = in_shadow.mean()
beta_sun = np.arcsin(abs(np.cross(e1, e2) @ sun_dir)) # the Sun's angle to the orbit plane
expected = np.arccos(np.sqrt(R_ORBIT ** 2 - 1) / (R_ORBIT * np.cos(beta_sun))) / np.pi
assert abs(fraction - expected) < 0.01, (fraction, expected)
print(f"{in_shadow.sum()} of {n_orbit} points are in shadow: {fraction:.3f} of the orbit; the cylindrical-shadow formula gives {expected:.3f}")
244 of 720 points are in shadow: 0.339 of the orbit; the cylindrical-shadow formula gives 0.340
The same answer by brute force#
The tree should change how fast the answer comes, never the answer. ray_tri against all 5120 triangles for every lane, in NumPy and in single precision, must find the same shadow and the same distance to the Earth’s surface.
def brute_force(o, d, tmax=np.inf):
# nearest hit by testing every triangle, in float32 like the kernel; returns hit, t
o, d = o.astype(np.float32), d.astype(np.float32)
best = np.full(len(o), tmax, np.float32)
found = np.zeros(len(o), bool)
for a, b, c in earth.astype(np.float32):
e1, e2 = b - a, c - a
p = np.cross(d, e2)
det = (e1 * p).sum(axis=1)
ok = det != 0
inv = np.where(ok, np.float32(1) / np.where(ok, det, np.float32(1)), np.float32(0))
s = o - a
u = inv * (s * p).sum(axis=1)
q = np.cross(s, e1)
v = inv * (d * q).sum(axis=1)
t = inv * (e2 * q).sum(axis=1)
hit = ok & (u >= 0) & (v >= 0) & (u + v <= 1) & (t >= 0) & (t <= best)
better = hit & (t < best)
best = np.where(better, t, best)
found |= better
return found, np.where(found, best, 0.0)
found, t_ref = brute_force(orbit, np.tile(sun_dir, (n_orbit, 1)))
assert np.array_equal(found, in_shadow)
assert np.allclose(result[1], t_ref, rtol=2e-5, atol=1e-6)
print(f"the tree and the {len(earth)}-triangle brute force agree on all {n_orbit} lanes: the same {found.sum()} shadows, distances equal to {np.abs(result[1] - t_ref).max():.1e}")
the tree and the 5120-triangle brute force agree on all 720 lanes: the same 244 shadows, distances equal to 6.0e-08
Who sees whom#
Eight satellites sit on random circular orbits at 1.1 to 1.4 Earth radii, and each pair is one lane: 28 lines of sight. A pair is blocked when the line between them dips below the surface, which for a sphere is when the segment passes closer to the centre than 1. The mesh is a little smaller than the sphere, so the two agree except for a line that grazes the surface within 0.11 %.
rng = np.random.default_rng(12)
radii = rng.uniform(1.1, 1.4, 8)
axes = rng.normal(size=(8, 3))
axes /= np.linalg.norm(axes, axis=1, keepdims=True)
sats = np.array([r * a for r, a in zip(radii, axes)])
pairs = np.array([(i, j) for i in range(8) for j in range(i + 1, 8)])
p1, p2 = sats[pairs[:, 0]], sats[pairs[:, 1]]
seen = launch(visible, p1=p1.T, p2=p2.T)
blocked = seen[0] > 0.5
# the closest the segment p1 -> p2 gets to the centre of the Earth
d = p2 - p1
s = np.clip(-(p1 * d).sum(axis=1) / (d * d).sum(axis=1), 0, 1)
closest = np.linalg.norm(p1 + s[:, None] * d, axis=1)
sphere_blocked = closest < 1.0
margin = np.abs(closest - 1.0) < 2e-3
assert np.array_equal(blocked[~margin], sphere_blocked[~margin])
brute, _ = brute_force(p1, d, 1.0)
assert np.array_equal(blocked, brute)
print(f"{blocked.sum()} of {len(pairs)} lines are blocked by the Earth; the sphere says {sphere_blocked.sum()} "
f"({margin.sum()} graze it within 0.002); the mesh brute force agrees with the tree on every line")
20 of 28 lines are blocked by the Earth; the sphere says 20 (0 graze it within 0.002); the mesh brute force agrees with the tree on every line
Depth#
The stack holds 32 levels by default and this tree needs 11. Asking for fewer than the tree has is refused when the kernel is traced:
try:
@hawk.kernel
def too_shallow(pos: Vector[3], sun: Vector[3], y: Mutable[Vector[5]],
lo: Table[Scalar], hi: Table[Scalar], left: Table[Index], right: Table[Index],
first: Table[Index], size: Table[Index], tri: Table[Scalar]):
y = closest_hit(bvh.tables(lo, hi, left, right, first, size, tri), pos, sun, depth=4)
except hawk.HawkError as err:
print(str(err)[:160])
closest_hit: the tree is 11 levels deep, past depth=4; pass depth=11 or more (a shorter stack would drop nodes without saying)
What just happened#
ray_triis one hawk op overaether::geometry::rayTriangle;build_bvhandclosest_hitbuild and walk a tree of such tests, and the whole traversal is ordinary kernel code, so it ran on the host and, when there was a GPU, on the device, from the same source.The mesh’s shadow fraction matched the cylindrical-shadow formula, and every lane matched a 5120-triangle brute force written separately in NumPy.
The pairwise lines of sight matched the sphere’s closest-approach test except where a line grazes within the mesh’s own 0.11 %.
Single precision was enough because the scene is in Earth radii. The same scene in metres puts
7e6in every coordinate, where single precision keeps only half a metre, and a facet a metre across stops being resolved.A stack shorter than the tree is an error at trace time, never a silently wrong answer.