resolved/scripts/make-example-scenes.py

#!/usr/bin/env python3
"""Writes examples/scenes/<model>.toml for the NASA models fetch-models.sh downloads.

Every scene has the same composition: the spacecraft in an isometric view (the
camera looks along its body (+1, +1, +1) diagonal, body +Z toward image up) with
Earth's limb crossing the image 2 degrees above center, Earth below and space
above. The camera distance makes the model's bounding sphere span about 5.6
degrees, half the example camera's width. The Sun direction is the one, on a
3-degree grid, that best lights the three visible body faces and Earth near the
limb, without putting the spacecraft in Earth's shadow.

Model centers and radii were measured from the decoded files at the listed
scales (the bounding box of every vertex); remeasure if a model changes.
"""

import math
import pathlib
import textwrap

RE = 6378137.0
AU = 1.495978707e11
LEO = 7.0e6
GEO = 42_164_000.0

# name, file, scale, bounding-box center (model frame, scaled, m), radius (m),
# orbit radius (m), description
MODELS = [
    ("landsat8", "landsat8.glb", 1.0, (0.0001, -3.6797, 3.6108), 8.3027, RE + 705e3,
     "Landsat 8 (85k triangles) in its 705 km orbit"),
    ("icesat2", "icesat2.glb", 1.0, (1.8979, -0.0013, 4.1380), 5.5159, RE + 496e3,
     "ICESat-2 (752k triangles) in its 496 km orbit; the heavy model for benchmarks"),
    ("ostm", "ostm.glb", 3.5, (0.0, -0.1750, 0.0), 4.3988, RE + 1336e3,
     "OSTM / Jason-2 (65k triangles) in its 1336 km orbit"),
    ("ssl1300", "ssl1300.glb", 0.1, (0.0, 0.5135, 0.0), 14.3366, GEO,
     "an SSL-1300 communications satellite (182k triangles) in geostationary orbit"),
    ("ace", "ace.glb", 0.00485, (0.0, 0.0, 0.0), 4.0114, LEO,
     "the Advanced Composition Explorer (26k triangles). The real ACE sits at the "
     "Sun-Earth L1 point; here it is put in low orbit"),
]


def norm(v):
    n = math.sqrt(sum(x * x for x in v))
    return [x / n for x in v]


def cross(a, b):
    return [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]]


def dot(a, b):
    return sum(x * y for x, y in zip(a, b))


def mul(m, v):
    return [dot(row, v) for row in m]


def transpose(m):
    return [[m[j][i] for j in range(3)] for i in range(3)]


def matmul(a, b):
    return [[sum(a[i][k] * b[k][j] for k in range(3)) for j in range(3)] for i in range(3)]


def triad(b1, b2, i1, i2):
    """The rotation taking body b1 onto inertial i1, with b2 as close to i2 as it can."""
    def frame(v1, v2):
        t1 = norm(v1)
        t2 = norm(cross(v1, v2))
        t3 = cross(t1, t2)
        return [[t1[i], t2[i], t3[i]] for i in range(3)]
    return matmul(frame(i1, i2), transpose(frame(b1, b2)))


def quaternion(r):
    """[w, x, y, z] of rotation matrix r, with w >= 0."""
    tr = r[0][0] + r[1][1] + r[2][2]
    if tr > 0:
        s = math.sqrt(tr + 1) * 2
        q = [0.25 * s, (r[2][1] - r[1][2]) / s, (r[0][2] - r[2][0]) / s, (r[1][0] - r[0][1]) / s]
    else:
        i = max(range(3), key=lambda k: r[k][k])
        j, k = (i + 1) % 3, (i + 2) % 3
        s = math.sqrt(1 + r[i][i] - r[j][j] - r[k][k]) * 2
        v = [0.0, 0.0, 0.0]
        v[i] = 0.25 * s
        v[j] = (r[j][i] + r[i][j]) / s
        v[k] = (r[k][i] + r[i][k]) / s
        q = [(r[k][j] - r[j][k]) / s] + v
    return [-x for x in q] if q[0] < 0 else q


def scene(name, file, scale, center, radius, orbit, description):
    limb = math.asin(RE / orbit)
    theta = limb - math.radians(2.0)
    zenith = [1.0, 0.0, 0.0]
    look = [orbit, 0.0, 0.0]
    boresight = [-math.cos(theta), math.sin(theta), 0.0]
    attitude = triad([1, 1, 1], [0, 0, 1], [-x for x in boresight], zenith)
    distance = radius / math.tan(math.radians(2.8))
    camera = [look[i] - distance * boresight[i] for i in range(3)]
    # Place the model so its bounding-box center is at the look point.
    offset = mul(attitude, center)
    position = [look[i] - offset[i] for i in range(3)]

    to_limb = [-math.cos(limb), math.sin(limb), 0.0]
    reach = math.sqrt(orbit * orbit - RE * RE)
    limb_normal = norm([look[i] + reach * to_limb[i] for i in range(3)])
    best = None
    for az in range(0, 360, 3):
        for el in range(-87, 90, 3):
            a, e = math.radians(az), math.radians(el)
            sun = [math.cos(e) * math.cos(a), math.cos(e) * math.sin(a), math.sin(e)]
            if dot(sun, [-1.0, 0.0, 0.0]) > math.cos(limb):
                continue  # behind Earth: the spacecraft would be in shadow
            faces = [dot(mul(attitude, n), sun) for n in ([1, 0, 0], [0, 1, 0], [0, 0, 1])]
            score = min(min(faces), dot(limb_normal, sun))
            if best is None or score > best[0]:
                best = (score, sun)
    sun = best[1]

    fmt = lambda v, d=2: "[" + ", ".join(f"{x:.{d}f}" for x in v) + "]"
    q = quaternion(attitude)
    header = textwrap.fill(
        f"{description[0].upper() + description[1:]}, in an isometric view with Earth's limb "
        "behind it. Generated by scripts/make-example-scenes.py; the model comes from "
        "scripts/fetch-models.sh.",
        width=80, initial_indent="# ", subsequent_indent="# ", break_on_hyphens=False)
    return f"""{header}
#
#   resolved render --scene examples/scenes/{name}.toml --preview {name}.png

sun = {fmt([AU * x for x in sun], 0)}
moon = [0, 0, -3.844e8]

[observer]
position = {fmt(camera)}
look_at = {fmt(look)}
up = [1, 0, 0]

[[target]]
model = "../models/nasa/{file}"
scale = {scale}
position = {fmt(position, 4)}
attitude = {fmt(q, 9)}
"""


out = pathlib.Path(__file__).resolve().parent.parent / "examples" / "scenes"
for model in MODELS:
    path = out / f"{model[0]}.toml"
    path.write_text(scene(*model))
    print(f"wrote {path}")