resolved/scripts/make-example-scenes.py

#!/usr/bin/env python3
"""Writes examples/scenes/<model>.toml for the NASA models fetch-models.sh downloads,
and two scenes with Earth, the Moon and a spacecraft together (below):
examples/scenes/earthrise.toml from lunar orbit and examples/scenes/moonrise.toml
from geostationary height.

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.

Earth is turned so the ground behind the spacecraft, at the image center, is
25 N 33 E (Egypt, the Nile and the Red Sea), with the camera looking north
toward the Mediterranean on the limb. The Sun is the direction, on a 2-degree
grid, that puts it at least 58 degrees above that ground (cos i >= 0.85) and
best lights the spacecraft's three visible faces, without putting the
spacecraft in Earth's shadow. With the Sun that high over the ground, one face
is lit at a slant.

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
# The ground at the image center, geocentric latitude and longitude in degrees.
GROUND = (25.0, 33.0)
# The least cos(incidence) the Sun may have on that ground.
GROUND_SUN = 0.85
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)]

    # The ground behind the spacecraft at the image center: where the boresight,
    # continued past the spacecraft, meets Earth.
    b = dot(look, boresight)
    reach = -b - math.sqrt(b * b - (dot(look, look) - RE * RE))
    ground = norm([look[i] + reach * boresight[i] for i in range(3)])
    # Turn Earth so GROUND is there, with north along the camera's view across it.
    lat, lon = (math.radians(x) for x in GROUND)
    fixed = [math.cos(lat) * math.cos(lon), math.cos(lat) * math.sin(lon), math.sin(lat)]
    north = norm([boresight[i] - dot(boresight, ground) * ground[i] for i in range(3)])
    earth = triad(fixed, [0, 0, 1], ground, north)

    best = None
    for az in range(0, 360, 2):
        for el in range(-88, 90, 2):
            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
            if dot(ground, sun) < GROUND_SUN:
                continue
            faces = [dot(mul(attitude, n), sun) for n in ([1, 0, 0], [0, 1, 0], [0, 0, 1])]
            if best is None or min(faces) > best[0]:
                best = (min(faces), 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
#   resolved render --scene examples/scenes/{name}.toml --lighting L1 --preview {name}.png

sun = {fmt([AU * x for x in sun], 0)}
moon = [0, 0, -3.844e8]
# Earth-fixed to inertial: 25 N 33 E behind the spacecraft, looking north.
earth_attitude = {fmt(quaternion(earth), 9)}

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

# At lighting L1, Earth's surface shows July (from scripts/fetch-earth.sh).
[appearance]
earth_month = 7

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


# Earthrise: from 1000 km above the Moon, Earth just above the curved lunar horizon
# through the example camera's 50 mm lens, Landsat 8 standing in for a lunar orbiter
# in the foreground.
RM = 1737400.0
MOON_DISTANCE = 3.844e8
EARTHRISE_ALTITUDE = 1000e3
EARTHRISE_FOCAL = 0.05
EARTHRISE_GROUND = (15.0, 25.0)  # Africa, geocentric degrees, on Earth's lit half
# (degrees above the lunar limb, degrees right of Earth) of the image center and of
# the satellite's center, Earth's center above the limb, and the satellite's distance.
EARTHRISE_CENTER = (2.0, -3.0)
EARTHRISE_EARTH = 2.6
EARTHRISE_SATELLITE = (2.4, -6.0)
EARTHRISE_RANGE = 500.0
EARTHRISE_EXPOSURE = 2e-4  # s


def earthrise():
    file, scale, center = "landsat8.glb", 1.0, MODELS[0][3]
    moon = [0.0, MOON_DISTANCE, 0.0]
    to_earth = norm([-x for x in moon])
    orbit = RM + EARTHRISE_ALTITUDE
    dip = math.acos(RM / orbit)

    # The observer over the Moon, tilted from the Earth-facing point along the lunar
    # equator (toward +X; toward +Z it would look over the pole, where the map has no
    # data) until Earth's center, as seen from there, stands EARTHRISE_EARTH above the
    # limb.
    want = math.radians(EARTHRISE_EARTH) - dip
    theta = want
    for _ in range(20):
        u = norm([math.sin(theta) * to_earth[i] + math.cos(theta) * (i == 0) for i in range(3)])
        camera = [moon[i] + orbit * u[i] for i in range(3)]
        earth_dir = norm([-x for x in camera])
        theta += want - math.asin(dot(earth_dir, u))
    up = u
    ahead = norm([earth_dir[i] - dot(earth_dir, up) * up[i] for i in range(3)])
    right = cross(ahead, up)

    def direction(above_limb, rightward):
        el, az = math.radians(above_limb) - dip, math.radians(rightward)
        return norm([math.cos(el) * math.cos(az) * ahead[i] + math.sin(el) * up[i]
                     + math.cos(el) * math.sin(az) * right[i] for i in range(3)])

    boresight = direction(*EARTHRISE_CENTER)
    look = [camera[i] + MOON_DISTANCE * boresight[i] for i in range(3)]
    to_satellite = direction(*EARTHRISE_SATELLITE)
    attitude = triad([1, 1, 1], [0, 0, 1], [-x for x in to_satellite], up)
    spot = [camera[i] + EARTHRISE_RANGE * to_satellite[i] for i in range(3)]
    offset = mul(attitude, center)
    position = [spot[i] - offset[i] for i in range(3)]

    # The lunar ground at the limb toward Earth.
    reach = math.sqrt(orbit ** 2 - RM ** 2)
    limb = direction(0.0, 0.0)
    ground = norm([camera[i] + reach * limb[i] - moon[i] for i in range(3)])

    # The Sun: Earth and the Moon are the subject, so the direction that best lights
    # both, the product of Earth's lit fraction as seen and the Sun's height over the
    # lunar ground at the limb, with the satellite out of the Moon's shadow.
    to_camera = norm(camera)  # from Earth's center
    best = None
    for az in range(0, 360, 2):
        for el in range(-88, 90, 2):
            a, e = math.radians(az), math.radians(el)
            sun = [math.cos(e) * math.cos(a), math.cos(e) * math.sin(a), math.sin(e)]
            phase = math.degrees(math.acos(dot(sun, to_camera)))
            if dot(sun, [-x for x in up]) > math.cos(math.asin(RM / orbit)):
                continue
            score = (1 + math.cos(math.radians(phase))) / 2 * max(dot(ground, sun), 0.0)
            if best is None or score > best[0]:
                best = (score, sun, phase)
    _, sun, phase = best

    # Earth turned with its axis along image up, and EARTHRISE_GROUND as near as that
    # allows to midway across its lit half as seen.
    across = norm([sun[i] - dot(sun, to_camera) * to_camera[i] for i in range(3)])
    shown = norm([to_camera[i] + across[i] for i in range(3)])
    north = norm([up[i] - dot(up, to_camera) * to_camera[i] for i in range(3)])
    lat, lon = (math.radians(x) for x in EARTHRISE_GROUND)
    fixed = [math.cos(lat) * math.cos(lon), math.cos(lat) * math.sin(lon), math.sin(lat)]
    earth = triad([0, 0, 1], fixed, north, shown)

    fmt = lambda v, d=2: "[" + ", ".join(f"{x:.{d}f}" for x in v) + "]"
    return f"""# Earthrise from lunar orbit: 1000 km above the Moon, Earth partly lit (phase
# {phase:.0f} degrees) just above the curved lunar horizon, and Landsat 8 standing in
# for a lunar orbiter {EARTHRISE_RANGE:.0f} m off in the foreground. Generated by
# scripts/make-example-scenes.py; the model comes from scripts/fetch-models.sh.
#
#   resolved render --scene examples/scenes/earthrise.toml --lighting L1 --preview earthrise.png
#
# At lighting L1 Earth shows Africa on its lit half, the Moon's ground its reflectance
# map, and the satellite is lit by the Sun and, faintly, by Earth and the Moon. The
# Moon has no atmosphere, so its horizon is a sharp edge, as in the Apollo photographs.

sun = {fmt([AU * x for x in sun], 0)}
moon = {fmt(moon)}
# Earth-fixed to inertial: north up, 15 N 25 E on the lit side.
earth_attitude = {fmt(quaternion(earth), 9)}

# Exposed so Earth, the Moon and the satellite stay short of saturation.
[camera]
focal_length = {EARTHRISE_FOCAL}
exposure = {EARTHRISE_EXPOSURE}

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

# At lighting L1, Earth's surface shows July (from scripts/fetch-earth.sh).
[appearance]
earth_month = 7

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


# Moonrise from geostationary height: looking toward Earth's limb through a 100 mm lens,
# Australia filling the lower half of the frame, the Moon rising over the limb, and the
# SSL-1300 in the foreground.
MOONRISE_ORBIT = GEO
MOONRISE_FOCAL = 0.1
MOONRISE_GROUND = (-25.0, 134.0)  # central Australia, geocentric degrees
# The image center's height above the limb (degrees). Then, as (degrees right of the
# image center, degrees above the limb at that column): the ground at the bottom of the
# frame, the Moon's center and the satellite's; and the satellite's distance (m).
MOONRISE_CENTER = 0.15
MOONRISE_SWATH = (0.0, -1.4)
MOONRISE_MOON = (1.7, 0.2)
MOONRISE_SATELLITE = (-1.6, 0.8)
MOONRISE_RANGE = 1100.0
MOONRISE_EXPOSURE = 2e-4  # s


def moonrise():
    file, scale, center = "ssl1300.glb", 0.1, MODELS[3][3]
    up, east = [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]
    b = RE * (1 - 1 / 298.257223563)
    lat = math.radians(MOONRISE_GROUND[0])
    radius_here = RE * b / math.hypot(b * math.cos(lat), RE * math.sin(lat))
    dip = math.acos(radius_here / MOONRISE_ORBIT)
    camera = [MOONRISE_ORBIT, 0.0, 0.0]

    def hits(d, earth=None):
        """How far along d the line of sight meets Earth, or None: the WGS84 ellipsoid
        turned by `earth` (Earth-fixed to inertial), or without it a sphere of the
        ellipsoid's radius at MOONRISE_GROUND."""
        if earth is None:
            o, v, k = camera, d, 1.0
        else:
            back = transpose(earth)
            stretch = [1.0, 1.0, RE / b]
            o = [x * f for x, f in zip(mul(back, camera), stretch)]
            v = [x * f for x, f in zip(mul(back, d), stretch)]
            k = 1.0
        r = radius_here if earth is None else RE
        aa, bb, cc = dot(v, v), dot(o, v), dot(o, o) - r * r
        disc = bb * bb - aa * cc
        return (-bb - math.sqrt(disc)) / aa * k if disc >= 0 else None

    # The camera's axes: looking east and MOONRISE_CENTER above the limb, up toward the
    # local vertical.
    el = math.radians(MOONRISE_CENTER) - dip
    forward = norm([math.cos(el) * east[i] + math.sin(el) * up[i] for i in range(3)])
    image_up = norm([up[i] - dot(up, forward) * forward[i] for i in range(3)])
    image_right = cross(forward, image_up)

    def ray(x, y):
        """The line of sight x degrees right of the image center and y degrees up."""
        tx, ty = math.tan(math.radians(x)), math.tan(math.radians(y))
        return norm([forward[i] + tx * image_right[i] + ty * image_up[i] for i in range(3)])

    def limb(x, earth=None):
        """Degrees up from the image center to the limb in column x."""
        lo, hi = -10.0, 10.0  # on Earth, in space
        for _ in range(60):
            mid = (lo + hi) / 2
            lo, hi = (mid, hi) if hits(ray(x, mid), earth) is not None else (lo, mid)
        return lo

    def place(x, above, earth=None):
        return ray(x, limb(x, earth) + above)

    look = [camera[i] + MOON_DISTANCE * forward[i] for i in range(3)]
    to_satellite = place(*MOONRISE_SATELLITE)
    attitude = triad([1, 1, 1], [0, 0, 1], [-x for x in to_satellite], image_up)
    spot = [camera[i] + MOONRISE_RANGE * to_satellite[i] for i in range(3)]
    offset = mul(attitude, center)
    position = [spot[i] - offset[i] for i in range(3)]

    # The ground at the bottom of the frame, turned to MOONRISE_GROUND with north
    # along the line of sight across it (toward the limb, image up).
    swath = place(*MOONRISE_SWATH)
    ground = norm([camera[i] + hits(swath) * swath[i] for i in range(3)])
    north = norm([swath[i] - dot(swath, ground) * ground[i] for i in range(3)])
    lon = math.radians(MOONRISE_GROUND[1])
    fixed = [math.cos(lat) * math.cos(lon), math.cos(lat) * math.sin(lon), math.sin(lat)]
    earth = triad(fixed, [0, 0, 1], ground, north)

    # The Moon on its line of sight, MOON_DISTANCE from Earth's center, placed against
    # the ellipsoid's limb.
    to_moon = place(*MOONRISE_MOON, earth)
    od = dot(camera, to_moon)
    reach = -od + math.sqrt(od * od - dot(camera, camera) + MOON_DISTANCE ** 2)
    moon = [camera[i] + reach * to_moon[i] for i in range(3)]

    # The Sun: Earth and the Moon are the subject, so the direction that best lights
    # both, the product of the Sun's height over that ground and the Moon's lit
    # fraction, with the satellite and the Moon out of Earth's shadow.
    from_moon = norm([camera[i] - moon[i] for i in range(3)])
    moon_dir = norm(moon)
    best = None
    for az in range(0, 360, 2):
        for el in range(-88, 90, 2):
            a, e = math.radians(az), math.radians(el)
            sun = [math.cos(e) * math.cos(a), math.cos(e) * math.sin(a), math.sin(e)]
            phase = math.degrees(math.acos(dot(sun, from_moon)))
            if dot(sun, [-x for x in norm(camera)]) > math.cos(math.asin(RE / GEO)):
                continue
            if dot(sun, [-x for x in moon_dir]) > math.cos(math.radians(3)):
                continue
            score = (1 + math.cos(math.radians(phase))) / 2 * max(dot(ground, sun), 0.0)
            if best is None or score > best[0]:
                best = (score, sun, phase)
    _, sun, phase = best

    fmt = lambda v, d=2: "[" + ", ".join(f"{x:.{d}f}" for x in v) + "]"
    return f"""# Moonrise from geostationary height: looking toward Earth's limb through a 100 mm
# lens, Australia filling the lower half of the frame, the Moon (phase {phase:.0f} degrees)
# rising over the limb, and an SSL-1300 communications satellite {MOONRISE_RANGE / 1000:.1f} km off in
# the foreground. Generated by scripts/make-example-scenes.py; the model comes from
# scripts/fetch-models.sh.
#
#   resolved render --scene examples/scenes/moonrise.toml --lighting L1 --preview moonrise.png
#
# At lighting L1 Earth's surface and the Moon's come from their maps. There is no
# atmosphere yet, so the limb is a sharp edge rather than a thin glowing band.

sun = {fmt([AU * x for x in sun], 0)}
moon = {fmt(moon)}
# Earth-fixed to inertial: 25 S 134 E at the bottom of the frame, looking north.
earth_attitude = {fmt(quaternion(earth), 9)}

# Exposed so Earth, the Moon and the satellite stay short of saturation.
[camera]
focal_length = {MOONRISE_FOCAL}
exposure = {MOONRISE_EXPOSURE}

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

# At lighting L1, Earth's surface shows July (from scripts/fetch-earth.sh).
[appearance]
earth_month = 7

[[target]]
model = "../models/nasa/{file}"
scale = {scale}
position = {fmt(position, 4)}
attitude = {fmt(quaternion(attitude), 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}")
path = out / "earthrise.toml"
path.write_text(earthrise())
print(f"wrote {path}")
path = out / "moonrise.toml"
path.write_text(moonrise())
print(f"wrote {path}")