connor/resolved
1
0 Fork 0
Code Issues 31 Pull Requests Actions Packages Projects Releases Wiki Activity

Added stars at L1

Browse Source
This commit is contained in:
connor committed 2026-10-16 14:43:50 -04:00
1 parent e0381ccdd7
commit 970551b4e1
26 files changed
+2118 -55
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

scripts/fetch-stars.py
Executable
+399
View File
@@ -0,0 +1,399 @@
#!/usr/bin/env python3
"""Downloads the Hipparcos, Tycho-2 and Gaia DR3 star catalogs and writes the renderer's
star file.
Sources:
- Hipparcos main catalog, ESA 1997 (CDS I/239, hip_main.dat): 118,218 entries,
complete to about V 7.3, with Johnson V and B-V. Positions are ICRS at epoch
J1991.25. The 2007 new reduction (I/311) improves the astrometry by about a
milliarcsecond, far below any pixel here, but carries no V, so the original is used.
- Tycho-2, Hog et al. 2000 (CDS I/259, tyc2.dat.00-19): 2.54 million stars to about
VT 11.5 (99 % complete to V 11.0), with Tycho BT and VT magnitudes. Mean positions
are ICRS at epoch J2000.
- Gaia DR3, Gaia Collaboration 2023 (gaiadr3.gaia_source_lite, from the ESA Gaia
archive's TAP service, https://gea.esac.esa.int/archive/): every star brighter than
G 15 by default (--gaia), 36.9 million, with its G magnitude and BP - RP color.
Positions are ICRS at epoch J2016.0. Gaia's photometry saturates for the brightest
stars, which Hipparcos and Tycho-2 cover.
Every star is taken from one catalog: Hipparcos first, then Tycho-2 for the stars it
does not have (by Hipparcos number), then Gaia for the stars neither has (by Gaia's own
cross-matches to them, gaiadr3.hipparcos2_best_neighbour and
gaiadr3.tycho2tdsc_merge_best_neighbour). Magnitudes are converted to Johnson V and
B - V:
- Tycho: V = VT - 0.090 (BT - VT), B - V = 0.850 (BT - VT) (ESA 1997, vol. 1, 1.3,
appendix 4). A star with no BT keeps V = VT and no color.
- Gaia: V from G and BP - RP by Riello et al. 2021 (A&A 649, A3, table C.2),
G - V = -0.02704 + 0.01424 c - 0.2156 c^2 + 0.01426 c^3 for c = BP - RP in -0.5 to
2.75, 0.03 mag scatter. B - V from BP - RP by a cubic fitted here to the Hipparcos
stars Gaia matched (printed when the script runs). A star with no BP - RP gets
c = 0.82, the median of stars near G 15, and no B - V.
Positions are carried to J2000 by their proper motion. The few Tycho-2 stars with no
mean position (pflag X) use their observed position, at about epoch 1991.5, with no
proper motion. Parallax is not kept: it moves every star by under 1 arcsecond.
Hipparcos and Tycho-2 come from the CDS archive (https://cdsarc.cds.unistra.fr) and are
checked by sha256. Gaia comes from about 2,000 queries of the archive, each a range of
source_id (a patch of sky), 4 at a time, split in half when a reply reaches the
service's 50,000-row limit; a query's results are kept, so a run can be resumed, but
they are not pinned by checksum. At G 15 that is about 2.3 GB of text and takes about
an hour. Everything goes to data/stars/source/; the output is data/stars/stars.bin,
about 600 MB at G 15, both gitignored. Needs numpy and pandas.
The output's format, little-endian:
8 bytes b"RSTARS02"
u32 number of stars
u32 0
then per star, brightest first, 16 bytes:
u32 right ascension, ICRS, epoch J2000, in units of 2 pi / 2^32
i32 declination, ICRS, epoch J2000, in units of pi / 2^32
i16 proper motion in right ascension times cos(declination), mas/yr
i16 proper motion in declination, mas/yr
u16 Johnson V magnitude plus 5, in thousandths
i16 Johnson B - V in thousandths, or -32768 if unknown
Usage: scripts/fetch-stars.py [--gaia G] [--faint V]
--gaia G take Gaia stars brighter than G (default 15; 0 for none)
--faint V keep stars to V (default: all those the catalogs give)
"""
import argparse
import concurrent.futures
import gzip
import hashlib
import io
import math
import os
import sys
import time
import urllib.parse
import urllib.request
import numpy as np
import pandas as pd
CDS = "https://cdsarc.cds.unistra.fr/ftp/I"
GAIA_TAP = "https://gea.esac.esa.int/tap-server/tap/sync"
# The service's limit on rows per synchronous query.
GAIA_ROWS = 50_000
# Gaia source_id carries the HEALPix level-12 cell in its high bits.
CELLS = 12 * 4**12
CELL = 2**35
FILES = {
"hip_main.dat": ("239/hip_main.dat", "58ceabb104d647160d9437ce6e513a02a036bb4ad9f8879a5a22fd52943616e0"),
"tyc2.dat.00.gz": ("259/tyc2.dat.00.gz", "1b224570fc6eb151984ce106fdf797728649d1ea77aaf3effd2d3444cfac6df6"),
"tyc2.dat.01.gz": ("259/tyc2.dat.01.gz", "36221f4a5cdd9c5009d64299d2e6409f8d4cdf924321768780165bd3e2a2a99a"),
"tyc2.dat.02.gz": ("259/tyc2.dat.02.gz", "e2ea4eeb7d204dde70ac54d172b1c429540c0d23d62ec0d7500561adc44cc57b"),
"tyc2.dat.03.gz": ("259/tyc2.dat.03.gz", "29e0cb57cba7e1651455efcb528276b0c28e1fe83f47f8ebe963e48d8b180afb"),
"tyc2.dat.04.gz": ("259/tyc2.dat.04.gz", "d04868329b470aa320f0fe5e9c525be619a7a228b364433e5f114f2178946d9b"),
"tyc2.dat.05.gz": ("259/tyc2.dat.05.gz", "51f08772197dcf0cd26ff098c866036340b797f71ee04be22712517c10c73c4e"),
"tyc2.dat.06.gz": ("259/tyc2.dat.06.gz", "f55503a66abce0e11d1a4ccbcde1f4c97a6e3e769de85e005f57ba963086ab03"),
"tyc2.dat.07.gz": ("259/tyc2.dat.07.gz", "bc15b5b1f6308477360fbb93f6a98d9029601cab68973202be0f8764c127f524"),
"tyc2.dat.08.gz": ("259/tyc2.dat.08.gz", "315ac34c678bc2e0f568b43b6b2d7c3d7c9cf9c089b51115d18b872864e50426"),
"tyc2.dat.09.gz": ("259/tyc2.dat.09.gz", "d66c671fa29aad10bc5fd697f68c918bb774c1ffc3c7db9c20234895268390d5"),
"tyc2.dat.10.gz": ("259/tyc2.dat.10.gz", "d86529df819ebdf3f4b8892510973ee10e1db06930814e014dd467193e740d52"),
"tyc2.dat.11.gz": ("259/tyc2.dat.11.gz", "fc0c3203b3d2787da54c43d9963e363c8a3926576c383b101a15350a7ca23e9c"),
"tyc2.dat.12.gz": ("259/tyc2.dat.12.gz", "68fd9d7e7353d52dea0de043cc62b45c04d234104eb05c063dc23ea7b6759576"),
"tyc2.dat.13.gz": ("259/tyc2.dat.13.gz", "a08ecc092e6d134742c0f594c56c5c2e02b8fd5e76cffb1779877f2fcefac3fc"),
"tyc2.dat.14.gz": ("259/tyc2.dat.14.gz", "2846e8cc489795a4d8160d2d5b64ff0e9a00598adba35a52a4f8ccf4c5e38b3e"),
"tyc2.dat.15.gz": ("259/tyc2.dat.15.gz", "74784a3656382090a70d778820335f4a781509dbe3fcb9e92d18cf96b2c46c72"),
"tyc2.dat.16.gz": ("259/tyc2.dat.16.gz", "320560e3d551cc40fa1f54ef7133709bf4bd45efa6688f808be7c09be5544a4d"),
"tyc2.dat.17.gz": ("259/tyc2.dat.17.gz", "a8154e940aa0a0e8f31d91b4bd9fd56e6849da42cddb13c7dda772439b06991a"),
"tyc2.dat.18.gz": ("259/tyc2.dat.18.gz", "5e951a1a51df7956f205b83cbfa5501d357c843640a253c1e6d3917bebe7d928"),
"tyc2.dat.19.gz": ("259/tyc2.dat.19.gz", "f59605a38116f517a31a7dbdee3469c077658f2f40b8afe5da2aeb832eaee3dd"),
}
# A star's fields as arrays: right ascension and declination (radians, J2000), proper
# motion (mas/yr), V and B - V (NaN if unknown).
FIELDS = ("ra", "dec", "pm_ra", "pm_dec", "v", "bv")
def fetch(source):
os.makedirs(source, exist_ok=True)
for name, (path, digest) in FILES.items():
target = os.path.join(source, name)
if not os.path.exists(target):
print(f"fetching {name}")
urllib.request.urlretrieve(f"{CDS}/{path}", target + ".part")
os.replace(target + ".part", target)
h = hashlib.sha256()
with open(target, "rb") as f:
for block in iter(lambda: f.read(1 << 20), b""):
h.update(block)
if h.hexdigest() != digest:
sys.exit(f"{target}: sha256 {h.hexdigest()}, expected {digest}")
def number(line, first, last):
"""The value in 1-based columns first to last, or None if blank."""
text = line[first - 1:last].strip()
return float(text) if text else None
def propagate(ra, dec, pm_ra, pm_dec, years):
"""Moves positions (radians) by their proper motion (mas/yr, right ascension times
cos declination) over `years`, along the sphere."""
scale = np.radians(1 / 3.6e6) * years
ca, sa, cd, sd = np.cos(ra), np.sin(ra), np.cos(dec), np.sin(dec)
x = cd * ca + scale * (pm_ra * -sa + pm_dec * -sd * ca)
y = cd * sa + scale * (pm_ra * ca + pm_dec * -sd * sa)
z = sd + scale * (pm_dec * cd)
norm = np.sqrt(x * x + y * y + z * z)
return np.mod(np.arctan2(y, x), 2 * np.pi), np.arcsin(z / norm)
def table(rows):
columns = list(zip(*rows)) if rows else [[]] * 6
return {k: np.array(c, dtype=float) for k, c in zip(FIELDS, columns)}
def hipparcos(path):
rows, numbers = [], set()
with open(path, encoding="ascii") as f:
for line in f:
v, ra, dec = number(line, 42, 46), number(line, 52, 63), number(line, 65, 76)
if v is None or ra is None or dec is None:
continue
bv = number(line, 246, 251)
rows.append((math.radians(ra), math.radians(dec), number(line, 88, 95) or 0.0,
number(line, 97, 104) or 0.0, v, math.nan if bv is None else bv))
numbers.add(int(line[8:14]))
stars = table(rows)
stars["ra"], stars["dec"] = propagate(
stars["ra"], stars["dec"], stars["pm_ra"], stars["pm_dec"], 2000.0 - 1991.25)
return stars, numbers
def hipparcos_colors(path):
"""B - V by Hipparcos number."""
colors = {}
with open(path, encoding="ascii") as f:
for line in f:
bv = number(line, 246, 251)
if bv is not None:
colors[int(line[8:14])] = bv
return colors
def tycho2(paths, hip_numbers):
rows, skipped = [], {"hipparcos": 0, "no VT": 0}
for path in paths:
with gzip.open(path, "rt", encoding="ascii") as f:
for line in f:
hip = line[142:148].strip()
if hip and int(hip) in hip_numbers:
skipped["hipparcos"] += 1
continue
bt, vt = number(line, 111, 116), number(line, 124, 129)
if vt is None:
skipped["no VT"] += 1
continue
if line[13] == "X":
ra, dec = number(line, 153, 164), number(line, 166, 177)
pm_ra = pm_dec = 0.0
else:
ra, dec = number(line, 16, 27), number(line, 29, 40)
pm_ra, pm_dec = number(line, 42, 48) or 0.0, number(line, 50, 56) or 0.0
if bt is None:
v, bv = vt, math.nan
else:
v, bv = vt - 0.090 * (bt - vt), 0.850 * (bt - vt)
rows.append((math.radians(ra), math.radians(dec), pm_ra, pm_dec, v, bv))
return table(rows), skipped
def gaia_query(query):
"""The rows of an ADQL query as a DataFrame, retrying a few times."""
data = urllib.parse.urlencode({"REQUEST": "doQuery", "LANG": "ADQL", "FORMAT": "csv",
"QUERY": query}).encode()
for attempt in range(6):
try:
with urllib.request.urlopen(GAIA_TAP, data=data, timeout=600) as reply:
return pd.read_csv(io.BytesIO(reply.read()))
except Exception as e: # network or service errors: wait and retry
if attempt == 5:
raise
print(f" retrying after {e}")
time.sleep(10 * (attempt + 1))
def gaia_range(directory, faint, lo, hi):
"""Gaia stars brighter than G `faint` with source_id in [lo, hi), from the cache in
`directory` or the archive, split until each query is within the row limit."""
path = os.path.join(directory, f"{lo}-{hi}.csv")
split = os.path.join(directory, f"{lo}-{hi}.split")
mid = (lo + hi) // 2
if os.path.exists(path):
return [path]
if os.path.exists(split):
return gaia_range(directory, faint, lo, mid) + gaia_range(directory, faint, mid, hi)
query = (
"SELECT g.ra, g.dec, g.pmra, g.pmdec, g.phot_g_mean_mag, g.bp_rp, "
"t.source_id AS tycho, h.source_id AS hip "
"FROM gaiadr3.gaia_source_lite AS g "
"LEFT OUTER JOIN gaiadr3.tycho2tdsc_merge_best_neighbour AS t ON g.source_id = t.source_id "
"LEFT OUTER JOIN gaiadr3.hipparcos2_best_neighbour AS h ON g.source_id = h.source_id "
f"WHERE g.phot_g_mean_mag < {faint} AND g.source_id >= {lo} AND g.source_id < {hi}"
)
rows = gaia_query(query)
if len(rows) >= GAIA_ROWS:
# Too many for one reply: split, and remember it, so a rerun does not ask again.
open(split, "w").close()
return gaia_range(directory, faint, lo, mid) + gaia_range(directory, faint, mid, hi)
rows.to_csv(path + ".part", index=False)
os.replace(path + ".part", path)
return [path]
def gaia(source, faint):
directory = os.path.join(source, f"gaia-g{faint:g}")
os.makedirs(directory, exist_ok=True)
chunks = 2048
edges = [CELLS * i // chunks * CELL for i in range(chunks + 1)]
paths = []
with concurrent.futures.ThreadPoolExecutor(4) as pool:
jobs = [pool.submit(gaia_range, directory, faint, lo, hi) for lo, hi in zip(edges, edges[1:])]
for done, job in enumerate(concurrent.futures.as_completed(jobs), 1):
paths += job.result()
if done % 100 == 0:
print(f" Gaia: {done} of {chunks} ranges")
return pd.concat((pd.read_csv(p) for p in sorted(paths)), ignore_index=True)
def gaia_color_fit(source, hip_main):
"""A cubic giving B - V from BP - RP, fitted to the Hipparcos stars Gaia matched, and
the scatter of the fit."""
path = os.path.join(source, "gaia-hipparcos.csv")
if not os.path.exists(path):
parts = []
edges = [CELLS * i // 8 * CELL for i in range(9)]
for lo, hi in zip(edges, edges[1:]):
parts.append(gaia_query(
"SELECT h.original_ext_source_id AS hip, g.bp_rp "
"FROM gaiadr3.hipparcos2_best_neighbour AS h "
"JOIN gaiadr3.gaia_source_lite AS g ON h.source_id = g.source_id "
f"WHERE g.source_id >= {lo} AND g.source_id < {hi}"
))
pd.concat(parts).to_csv(path, index=False)
matched = pd.read_csv(path).dropna()
colors = hipparcos_colors(hip_main)
c = matched["bp_rp"].to_numpy()
bv = np.array([colors.get(int(h), np.nan) for h in matched["hip"]])
ok = np.isfinite(bv) & (c > -0.5) & (c < 3.5)
coefficients = np.polyfit(c[ok], bv[ok], 3)
residual = bv[ok] - np.polyval(coefficients, c[ok])
return coefficients, np.median(np.abs(residual)) * 1.4826, ok.sum()
def gaia_stars(rows, coefficients):
"""Gaia rows that neither Hipparcos nor Tycho-2 has, as V and B - V at J2000."""
rows = rows[rows["tycho"].isna() & rows["hip"].isna()]
c = rows["bp_rp"].to_numpy()
known = np.isfinite(c)
c = np.where(known, c, 0.82)
cc = np.clip(c, -0.5, 2.75)
g_minus_v = -0.02704 + 0.01424 * cc - 0.2156 * cc**2 + 0.01426 * cc**3
pm_ra = np.nan_to_num(rows["pmra"].to_numpy())
pm_dec = np.nan_to_num(rows["pmdec"].to_numpy())
ra, dec = propagate(np.radians(rows["ra"].to_numpy()), np.radians(rows["dec"].to_numpy()),
pm_ra, pm_dec, 2000.0 - 2016.0)
return {
"ra": ra, "dec": dec, "pm_ra": pm_ra, "pm_dec": pm_dec,
"v": rows["phot_g_mean_mag"].to_numpy() - g_minus_v,
"bv": np.where(known, np.polyval(coefficients, np.clip(c, -0.5, 3.5)), np.nan),
}
def unit(stars):
cd = np.cos(stars["dec"])
return np.stack([cd * np.cos(stars["ra"]), cd * np.sin(stars["ra"]), np.sin(stars["dec"])], 1)
def drop_near(stars, others, radius):
"""`stars` without those within `radius` radians of any of `others`, and how many
were dropped: cells of the radius on a side, searched with their neighbors."""
key = lambda c: ((c[:, 0] + (1 << 20)) << 42) | ((c[:, 1] + (1 << 20)) << 21) | (c[:, 2] + (1 << 20))
points = unit(stars)
keys = key(np.floor(points / radius).astype(np.int64))
order = np.argsort(keys)
# A key no cell has ends the list, so a walk past the last star stops there.
sorted_keys = np.append(keys[order], np.iinfo(np.int64).max)
near = unit(others)
cells = np.floor(near / radius).astype(np.int64)
drop = np.zeros(len(points), dtype=bool)
for offset in np.ndindex(3, 3, 3):
k = key(cells + np.array(offset) - 1)
i = np.searchsorted(sorted_keys, k)
# Every star in the cell: walk forward while the key repeats.
while True:
hit = sorted_keys[i] == k
if not hit.any():
break
j = order[i[hit]]
close = np.linalg.norm(points[j] - near[hit], axis=1) < radius
drop[j[close]] = True
i = i + hit
return {f: stars[f][~drop] for f in FIELDS}, int(drop.sum())
def write(path, stars, faint):
keep = stars["v"] <= faint
order = np.argsort(stars["v"][keep], kind="stable")
s = {k: stars[k][keep][order] for k in FIELDS}
record = np.dtype([("ra", "<u4"), ("dec", "<i4"), ("pm_ra", "<i2"), ("pm_dec", "<i2"),
("v", "<u2"), ("bv", "<i2")])
out = np.zeros(len(s["v"]), dtype=record)
out["ra"] = np.round(np.mod(s["ra"], 2 * np.pi) / (2 * np.pi) * 2**32).astype(np.uint64) % 2**32
out["dec"] = np.clip(np.round(s["dec"] / np.pi * 2**32), -(2**31), 2**31 - 1)
out["pm_ra"] = np.clip(np.round(s["pm_ra"]), -32767, 32767)
out["pm_dec"] = np.clip(np.round(s["pm_dec"]), -32767, 32767)
out["v"] = np.clip(np.round((s["v"] + 5.0) * 1000), 0, 65535)
out["bv"] = np.where(np.isfinite(s["bv"]), np.clip(np.round(s["bv"] * 1000), -32767, 32767), -32768)
with open(path + ".part", "wb") as f:
f.write(b"RSTARS02" + np.array([len(out), 0], dtype="<u4").tobytes())
f.write(out.tobytes())
os.replace(path + ".part", path)
return len(out), s["v"][-1] if len(out) else math.nan
def main():
parser = argparse.ArgumentParser(description=__doc__.split("\n")[0])
parser.add_argument("--gaia", type=float, default=15.0, help="Gaia stars brighter than G")
parser.add_argument("--faint", type=float, default=math.inf, help="faintest V kept")
args = parser.parse_args()
root = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "data", "stars")
source = os.path.join(root, "source")
fetch(source)
hip_main = os.path.join(source, "hip_main.dat")
hip, numbers = hipparcos(hip_main)
tyc, skipped = tycho2(
sorted(os.path.join(source, n) for n in FILES if n.startswith("tyc2")), numbers
)
print(f"Hipparcos: {len(hip['v'])} stars; Tycho-2: {len(tyc['v'])} more "
f"({skipped['hipparcos']} in Hipparcos, {skipped['no VT']} without VT)")
parts = [hip, tyc]
if args.gaia > 0:
coefficients, scatter, used = gaia_color_fit(source, hip_main)
print(f"B - V from BP - RP: cubic {np.round(coefficients, 5).tolist()} "
f"(highest power first), {scatter:.3f} mag scatter over {used} Hipparcos stars")
rows = gaia(source, args.gaia)
extra = gaia_stars(rows, coefficients)
print(f"Gaia DR3 to G {args.gaia:g}: {len(rows)} stars, {len(extra['v'])} not in "
f"Hipparcos or Tycho-2 by Gaia's cross-match")
# Gaia's cross-match misses about 1 % of Tycho-2: drop Gaia stars within 1
# arcsecond of one, which Tycho's and Hipparcos's photometry, unable to resolve
# pairs that close, already includes.
others = {k: np.concatenate([hip[k], tyc[k]]) for k in FIELDS}
extra, dropped = drop_near(extra, others, np.radians(1 / 3600))
print(f" {dropped} more within 1 arcsecond of a Hipparcos or Tycho-2 star")
parts.append(extra)
stars = {k: np.concatenate([p[k] for p in parts]) for k in FIELDS}
out = os.path.join(root, "stars.bin")
count, faintest = write(out, stars, args.faint)
print(f"{count} stars to V {faintest:.2f} in {out}")
if __name__ == "__main__":
main()
scripts/make-example-scenes.py
+140 -2
View File
@@ -1,8 +1,10 @@
#!/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):
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.
from geostationary height, and scenes for the stars of camera tier C1: a star field,
examples/scenes/stars-pleiades.toml, and a telescope on one geostationary satellite
imaging another along the arc, examples/scenes/geo-telescope-*.toml.
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
@@ -447,6 +449,135 @@ attitude = {fmt(quaternion(attitude), 9)}
"""
# Earth's orbital speed, m/s, and the geostationary orbital speed.
EARTH_SPEED = 29.78e3
GEO_SPEED = 3074.7
# The telescope's longitude on the geostationary arc, degrees: looking along the arc
# from here points at right ascension 103 degrees, where the Milky Way crosses the
# equator in Monoceros, the densest stretch of the equator (190 stars per square degree
# to V 12, against 50 where it crosses in Aquila, behind the dust of the Aquila Rift).
GEO_LONGITUDE = 13.0
# The ecliptic's north pole in the inertial (equatorial J2000) frame.
ECLIPTIC_POLE = [0.0, -math.sin(math.radians(23.4393)), math.cos(math.radians(23.4393))]
def earth_velocity(sun):
"""Earth's velocity around the Sun when the Sun is in direction `sun`: along the
ecliptic, a quarter turn ahead of the Sun's direction from Earth."""
return [EARTH_SPEED * x for x in norm(cross(ECLIPTIC_POLE, [-x for x in sun]))]
def stars_pleiades():
# The Pleiades, right ascension 56.75 and declination 24.12 degrees, from
# geostationary orbit, with the Sun behind the camera.
ra, dec = math.radians(56.75), math.radians(24.12)
look = [math.cos(dec) * math.cos(ra), math.cos(dec) * math.sin(ra), math.sin(dec)]
camera = [GEO, 0.0, 0.0]
sun = [-x for x in look]
fmt = lambda v, d=2: "[" + ", ".join(f"{x:.{d}f}" for x in v) + "]"
return f"""# A star field: the Pleiades through the example camera (12 by 7 degrees), from
# geostationary orbit, long-exposed at 5 s. Stars are drawn at lighting L1 with camera
# C1, from the catalog scripts/fetch-stars.py writes. Alcyone, the brightest at V 2.9,
# gives about 200,000 electrons, ten times the full well; a V 8 star about 2,000 and a
# V 11 star about 130. Generated by scripts/make-example-scenes.py.
#
# resolved render --scene examples/scenes/stars-pleiades.toml --lighting L1 --camera C1 --pretty --preview pleiades.png
sun = {fmt([AU * x for x in sun], 0)}
moon = [0, -384400000, 0]
earth_velocity = {fmt(earth_velocity(sun))}
[camera]
exposure = 5
[observer]
position = {fmt(camera)}
look_at = {fmt([camera[i] + 1e9 * look[i] for i in range(3)])}
up = [0, 0, 1]
velocity = [0.00, {GEO_SPEED:.2f}, 0.00]
# At lighting L1, Earth's surface shows July (from scripts/fetch-earth.sh).
[appearance]
earth_month = 7
"""
def geo_telescope(range_m, phase_deg, f_number, exposure, description):
# The camera at GEO_LONGITUDE on the geostationary arc, the target ahead along it.
file, scale, center = "ssl1300.glb", 0.1, MODELS[3][3]
lon = math.radians(GEO_LONGITUDE)
camera = [GEO * math.cos(lon), GEO * math.sin(lon), 0.0]
angle = lon + range_m / GEO
spot = [GEO * math.cos(angle), GEO * math.sin(angle), 0.0]
velocity = [-GEO_SPEED * math.sin(lon), GEO_SPEED * math.cos(lon), 0.0]
look = norm([spot[i] - camera[i] for i in range(3)])
image_up = [0.0, 0.0, 1.0]
attitude = triad([1, 1, 1], [0, 0, 1], [-x for x in look], image_up)
offset = mul(attitude, center)
position = [spot[i] - offset[i] for i in range(3)]
# The Sun `phase_deg` from the camera's direction as the target sees it, turned
# about the image's up, 10 degrees out of the equator so the target is not shadowed
# by Earth at any phase.
back = [-x for x in look]
side = cross(image_up, back)
p, t = math.radians(phase_deg), math.radians(10.0)
sun = norm([
math.cos(t) * (math.cos(p) * back[i] + math.sin(p) * side[i]) + math.sin(t) * image_up[i]
for i in range(3)
])
fmt = lambda v, d=2: "[" + ", ".join(f"{x:.{d}f}" for x in v) + "]"
heading = textwrap.fill(description, 86, initial_indent="# ", subsequent_indent="# ")
return f"""{heading}
# A 1 m, f/{f_number:.0f} telescope ({100 / f_number:g} cm aperture, 5.5 urad per pixel, 0.6 by 0.34
# degrees) on one geostationary satellite, imaging an SSL-1300 {range_m / 1000:.0f} km ahead along
# the arc, lit at a phase of {phase_deg} degrees, against the stars toward the Milky Way in
# Monoceros. Stars are drawn at lighting L1 with camera C1 (scripts/fetch-stars.py).
# Generated by scripts/make-example-scenes.py; the model comes from
# scripts/fetch-models.sh.
#
# resolved render --scene <this file> --lighting L1 --camera C1 --pretty --preview geo.png
sun = {fmt([AU * x for x in sun], 0)}
moon = [0, -384400000, 0]
earth_velocity = {fmt(earth_velocity(sun))}
[camera]
focal_length = 1.0
f_number = {f_number:.1f}
exposure = {exposure}
[observer]
position = {fmt(camera)}
look_at = {fmt(spot)}
up = [0, 0, 1]
velocity = {fmt(velocity)}
# At lighting L1, Earth's surface shows July (from scripts/fetch-earth.sh), and lights
# the target's Earth-facing side.
[appearance]
earth_month = 7
[[target]]
model = "../models/nasa/{file}"
scale = {scale}
position = {fmt(position, 4)}
attitude = {fmt(quaternion(attitude), 9)}
"""
GEO_TELESCOPE = [
("geo-telescope-100km", 100e3, 30, 8.0, 6e-4,
"GEO to GEO at 100 km: the target about 55 px across, lit nearly face on, exposed "
"short of saturating it, with only the brighter stars behind it."),
("geo-telescope-100km-side", 100e3, 90, 8.0, 6e-4,
"GEO to GEO at 100 km, lit from the side: half the target in shadow."),
("geo-telescope-1000km", 1000e3, 30, 4.0, 0.2,
"GEO to GEO at 1000 km, exposed 0.2 s at f/4 as a surveillance camera would: the "
"target a few pixels across and far past saturation, in a field of stars to the "
"catalog's limit, about V 15, which gives a few DN."),
]
out = pathlib.Path(__file__).resolve().parent.parent / "examples" / "scenes"
for model in MODELS:
path = out / f"{model[0]}.toml"
@@ -458,3 +589,10 @@ print(f"wrote {path}")
path = out / "moonrise.toml"
path.write_text(moonrise())
print(f"wrote {path}")
path = out / "stars-pleiades.toml"
path.write_text(stars_pleiades())
print(f"wrote {path}")
for name, range_m, phase, f_number, exposure, description in GEO_TELESCOPE:
path = out / f"{name}.toml"
path.write_text(geo_telescope(range_m, phase, f_number, exposure, description))
print(f"wrote {path}")