In The Epiphany of Gliese 581, a group of explorers search the mortal remains of a dead superintelligence. The expedition begins in Beta Pictoris—today an unremarkable blue star; in the story, a posthuman Dyson swarm civilization of quadrillions—then passes through Gliese 581, and ends where it started.
I wanted to build a timeline of the story, and, because this is hard science fiction, this means doing real math on real astronomical data. I ended up writing a small framework for doing astronomical calculations. In Common Lisp, for old time’s sake.
The code is in this repo. The rest of this post is a walkthrough of the astronomy framework, followed by the storyspecific code.
The Problem
Building a timeline requires knowing travel times, which requires knowing the distances between the stars. Which requires knowing the positions of the stars^{1}. So I had to download the HYG database, which has all the information I need.
But. The characters don’t travel in a straightline ballistic trajectory. They’re digital people, so they can travel at the speed of light by sending their mental states over an interstellar communications network. And, because lasers decohere with distance, each “jump” is limited to a relatively short distance. So the fastest route from A to B is not a straight line on a rocket, but a zigzagging trajectory on an optical link. To find the fastest network route I wrote an implementation of Dijkstra’s algorithm and ran it over a graph of stars, the edges linking all stars whose distance is less than the laser cutoff distance.
An added constraint: Gliese 581 is not on the network, because it was inhabited by a taciturn superintelligence. So I had to write even more code to find the star closest to Gliese 581, from which the characters complete the last leg of the journey on a fusion rocket.
Distances
The most widely used unit of length in astronomy is the parsec, some weird trigonometry thing about parallax. It’s just as geocentric as lightyears but harder to intuit. We will want to convert to light years to present the results.
We could just use floats, but I have a peculiar malady where I enjoy writing CLOS class definitions too much. Also, if I don’t have distinct types for the two units, I will confuse them and ruin all subsequent calculations.
(defclass lightyears ()
((value :reader value
:initarg :value
:type real
:documentation "The underlying value."))
(:documentation "Represents distance in light years."))
(defclass parsecs ()
((value :reader value
:initarg :value
:type real
:documentation "The underlying value."))
(:documentation "Represents distance in parsecs."))
(defun makeparsecs (value)
"Create an instance of PARSECS from a numeric value."
(makeinstance 'parsecs :value value))
We can easily convert between the two units by multiplying by a constant:
\[\begin{align*} pc &= ly \times 0.306601 \\ ly &= pc \times 3.26156 \end{align*}\]In Common Lisp these can be implemented as:
(defun lightyearstoparsecs (ly)
"Convert the given distance in light years to parsecs."
(makeparsecs (* (value ly) 0.306601)))
(defun parsecstolightyears (pc)
"Convert the given distance in parsecs to light years."
(makeinstance 'lightyears :value (* (value pc) 3.26156)))
For example:
CLUSER> (makeparsecs 5.0)
#<PARSECS 5.0pc>
CLUSER> (parsecstolightyears *)
#<LIGHTYEARS 16.3ly>
CLUSER> (lightyearstoparsecs *)
#<PARSECS 5.0pc>
Star Positions
Star positions are typically given in equatorial coordinates, a spherical coordinate system where positions in space are represented by two angles (called right ascension (RA) and declination (DEC)), and a radius (the star’s distance from the Sun). For historical reasons, right ascension is reported in hoursminutesseconds format, and declination is reported in degreesminutesseconds.
Luckily, the HYG database contains the Cartesian (X, Y, Z) coordinates^{2} as well. This makes it much easier to calculate distances and, later, to make star maps (see the appendix for dealing with equatorial coordinates).
The cartesianposition
class represents a Cartesian triple:
(defclass cartesianposition ()
((x :reader x
:initarg :x
:type parsecs
:documentation "The X coordinate in parsecs.")
(y :reader y
:initarg :y
:type parsecs
:documentation "The Y coordinate in parsecs.")
(z :reader z
:initarg :z
:type parsecs
:documentation "The Z coordinate in parsecs."))
(:documentation "A position in Cartesian (X, Y, Z) coordinates."))
And the euclideandistance
function calculates the distance between two points:
(defun euclideandistance (p1 p2)
"Calculate the Euclidean distance between two Cartesian coordinates.
Returns a value in parsecs."
(withslots ((x1 x) (y1 y) (z1 z)) p1
(withslots ((x2 x) (y2 y) (z2 z)) p2
(let ((x1 (value x1))
(y1 (value y1))
(z1 (value z1))
(x2 (value x2))
(y2 (value y2))
(z2 (value z2)))
(makeparsecs (sqrt (+ (expt ( x1 x2) 2)
(expt ( y1 y2) 2)
(expt ( z1 z2) 2))))))))
We can use this like this:
CLUSER> (defparameter a
(makeinstance 'cartesianposition
:x (makeparsecs 3.4)
:y (makeparsecs 6.7)
:z (makeparsecs 1.2)))
#<CARTESIANPOSITION X=3.4pc Y=6.7pc Z=1.2pc>
CLUSER> (defparameter b
(makeinstance 'cartesianposition
:x (makeparsecs 9.1)
:y (makeparsecs 4.3)
:z (makeparsecs 7.2)))
#<CARTESIANPOSITION X=9.1pc Y=4.3pc Z=7.2pc>
CLUSER> (euclideandistance a b)
#<PARSECS 13.8pc>
Parsing the HYG Database
The HYG database is a straightforward CSV, so we can parse it easily. We use the
parsenumber
library to parse floats (the alternative is readfromstring
,
which is heavy and insecure).
First, a class to store the database:
(defclass hygdatabase ()
((stars :reader databasestars
:initarg :stars
:type (vector star)
:documentation "The star vector."))
(:documentation "The inmemory HYG database."))
(defun starcount (db)
(length (databasestars db)))
And another to represent the star data we care about:
(defclass star ()
((id :reader starid
:initarg :id
:type integer
:documentation "The star's ID in the HYG database.")
(proper :reader starproper
:initarg :proper
:type (or null string)
:documentation "The star's proper name, if known.")
(hip :reader starhip
:initarg :hip
:type (or null string)
:documentation "The star's name in the Hipparcos catalog, if known.")
(hd :reader starhd
:initarg :hd
:type (or null string)
:documentation "The star's name in the Henry Draper catalog, if known.")
(gliese :reader stargliese
:initarg :gliese
:type (or null string)
:documentation "The star's name in the the third edition of the Gliese Catalog of Nearby Stars, if known.")
(bayer :reader starbayer
:initarg :bayer
:type (or null string)
:documentation "The star's BayerFlamsteed designation, if known.")
(distance :reader stardistance
:initarg :distance
:type parsecs
:documentation "The star's distance from the Sun in parsecs.")
(cartesianposition :reader starcartesianposition
:initarg :cartesianposition
:documentation "The star's Cartesian (X, Y, Z) position."))
(:documentation "Represents a star from the HYG database."))
(defun starname (star)
"A star's name. The following are tried in order: proper name, Bayer
designation, Gliese name, HIP name, HD name. If the star doesn't have any names,
returns '?'."
(withslots (proper bayer gliese hip) star
(or proper
bayer
gliese
(if hip
(concatenate 'string "HIP " hip)
"?"))))
And the parsing code is very straightforward:
(defun loadhygdatabase (pathname)
"Load the HYG database from a CSV file."
(withopenfile (stream pathname :direction :input)
;; Discard the header.
(readline stream nil)
(let ((stars (makearray 0 :adjustable t :elementtype 'star :fillpointer 0)))
(loop for line = (readline stream nil)
while line
do
(let ((columns (uiop:splitstring line :separator ",")))
(let ((star (parsestar columns)))
(vectorpushextend star stars))))
(makeinstance 'hygdatabase :stars stars))))
(defun parsestar (cells)
"Parse a star from a row (a list of CSV cells)."
(destructuringbind (id hip hd hr gliese bayer proper ra dec dist prma prdec rv mag absmag spect ci x y z &rest etc) cells
(declare (ignore hr ra dec prma prdec rv mag absmag spect ci etc))
(makeinstance 'star
:id (parseinteger id)
:proper (stringornil proper)
:hip (stringornil hip)
:hd (stringornil hd)
:gliese (stringornil gliese)
:bayer (stringornil bayer)
:distance (makeparsecs (parsenumber:parserealnumber dist))
:cartesianposition (makeinstance 'cartesianposition
:x (parseparsecs x)
:y (parseparsecs y)
:z (parseparsecs z)))))
(defun stringornil (str)
(if (string= str "") nil str))
(defun parseparsecs (str)
(makeparsecs (parsenumber:parserealnumber str)))
And, finally, some code to find stars by name:
(defun findstarbyname (db name)
(loop for star across (databasestars db) do
(when (string= name (starname star))
(returnfrom findstarbyname star)))
nil)
Now we can query the database:
CLUSER> (defparameter db (loadhygdatabase #p"hygdata_v3.csv"))
DB
CLUSER> (defparameter star (findstarbyname db "Gl 581"))
#<STAR Gl 581>
CLUSER> (starname star)
"Gl 581"
CLUSER> (stardistance star)
#<PARSECS 6.2pc>
CLUSER> (starcartesianposition star)
#<CARTESIANPOSITION X=4.0pc Y=4.7pc Z=.8pc>
Nearest Stars
Since Gliese 581 is not on the network, we need to find the star closest to it. First, this function returns all stars within a given radius of a position:
(defun findstarswithinradius (db pos radius)
"Given an HYG database, a position in Cartesian coordinates, and a radius in
parsecs, return a vector of all the stars that are within the radius from that
position."
(let ((stars (makearray 0 :adjustable t :elementtype 'star :fillpointer 0)))
(loop for star across (databasestars db) do
(let ((starpos (starcartesianposition star)))
(when (< (value (euclideandistance pos starpos))
(value radius))
(vectorpushextend star stars))))
stars))
Now, the driver code. First, we load the database:
(defparameter +db+
(loadhygdatabase #p"hygdata_v3.csv"))
(assert (= (starcount +db+) 119614))
Then we find all the stars within 3 parsecs of Gliese 581, sort them, and print them out:
(defparameter +g581+ (findstarbyname +db+ "Gl 581"))
(let ((stars (findstarswithinradius +db+
(starcartesianposition +g581+)
(makeparsecs 3.0))))
;; Sort stars by distance.
(flet ((dist (star)
;; Distance to Gliese 581.
(stareuclideandistance +g581+ star)))
(let ((sorted (sort stars
#'(lambda (stara starb)
(< (value (dist stara))
(value (dist starb)))))))
;; Print the list of stars.
(format t "Ten stars closest to Gliese 581:~%~%")
(format t "~8@A ~12@A~%" "Dist" "Star")
(format t "~%")
; subseq because the first star is Gl581 itself, and because we only want the top 10
(loop for star across (subseq sorted 1 11) do
(format t "~6,2fly ~12@A~%"
(value (parsecstolightyears (dist star)))
(starname star)))
(let ((star (elt sorted 1)))
(format t "~%The star closest to Gliese 581 is ~A at ~0,2fly"
(starname star)
(value (parsecstolightyears (dist star))))))))
(defun stareuclideandistance (a b)
"The Euclidean distance between two stars in parsecs."
(euclideandistance (starcartesianposition a) (starcartesianposition b)))
The output here is:
Ten stars closest to Gliese 581:
Dist Star

4.25ly Gl 555
5.15ly Gl 570B
5.17ly Gl 570A
7.43ly NN 3877
7.87ly Gl 563.2A
8.17ly Gl 644B
8.18ly Gl 644C
8.21ly Gl 628
8.34ly Gl 644A
8.83ly Gl 643
The star closest to Gliese 581 is Gl 555 at 4.25ly
Now we have to find a route from Beta Pictoris to Gliese 555.
The Network Route
The shortest path between two points is the Euclidean distance. But rockets are slow, and light is fast. And the characters, being uploads, have the option of travelling by simply sending copies of their mind over a communications network. Optical transceivers have a limited range, but even then, an indirect, zigzagging trajectory at the speed of light is faster than a straightline trajectory on a spacecraft moving at maybe 10% of the speed of light.
So: what’s the shortest path between two stars, when each hop is limited by technology and economics to some constant distance? To find the answer, I just implemented Dijkstra’s algorithm and ran it over the graph of stars and network connections.
The edge
and graph
classes represent a graph where edges have a cost:
(defclass edge ()
((start :reader edgestart
:initarg :start
:type integer
:documentation "The ID of the start node.")
(end :reader edgeend
:initarg :end
:type integer
:documentation "The ID of the end node.")
(cost :reader edgecost
:initarg :cost
:type number
:documentation "The cost of traversing this edge."))
(:documentation "Represents an edge in a graph."))
(defmethod initializeinstance :after ((edge edge) &key)
(withslots (start end cost) edge
;; Check edges are not degenerate.
(when (= start end)
(error "Degenerate edge: start and end IDs are the same: ~A" start))
;; Check cost is nonnegative.
(when (< cost 0.0)
(error "Edge cost is negative: ~A" cost))))
(defclass graph ()
((vertices :accessor graphvertices
:initarg :vertices
:type (vector integer)
:documentation "A vector of vertex IDs.")
(edges :accessor graphedges
:initarg :edges
:type (vector edge)
:documentation "A vector of edge objects."))
(:documentation "Represents a graph."))
This convenience function builds a graph object from the vector of edges:
(defun makegraphfromedges (edges)
"Construct a graph from a vector of edges."
(let ((vertextable (makehashtable :test 'equal)) ; table of seen IDs
(vertices (makearray 0 :adjustable t ; vertex accumulator
:elementtype 'integer
:fillpointer 0)))
(loop for edge across edges do
(withslots (start end) edge
(unless (gethash start vertextable)
(vectorpushextend start vertices)
(setf (gethash start vertextable) t))
(unless (gethash end vertextable)
(vectorpushextend end vertices)
(setf (gethash end vertextable) t))))
(makeinstance 'graph :vertices vertices
:edges edges)))
Given a graph whose edges have a cost, Dijkstra’s algorithm finds the cheapest path between given start and end nodes. When cost represents distance, this is the shortest path.
At the highest level: Dijkstra is breadthfirst beam search on a tree rooted at the start vertex.
More detailed: Dijkstra begin with the start vertex, and performs breadthfirst search by building up a tree through the neighbouring nodes. Each leaf node knows its cost, that is, the sum of the edge costs in the path from the root to the leaf. In normal breadthfirst search, the search order is just the iteration order of the node’s children. Dijkstra picks the next node that minimizes cost. That’s the beam part of beam search.
Even more detailed, here’s the full pseudocode:
Dijkstra’s Algorithm
 Inputs:
 $G : \text{Graph}$
 $V_i: \text{Vertex}$: the start vertex.
 $V_f: \text{Vertex}$: the end vertex.
 Let:
 $D: \text{Map}[\text{Vertex}, \mathbb{R}]$ is the table of distances from $V_i$ to every other vertex. Initially, we set $D[V_i] := 0$ and $D[v] := +\infty, \forall v \in G$ for all other vertices.
 $L : \text{Map}[\text{Vertex}, \text{Option}[\text{Vertex}]]$ is the previous links table, which keeps track of the path we build while the algorithm runs. It maps a vertex to the previous vertex in the path. Initially, $L[v] := \text{NIL}, \forall v \in G$.
 $Q : \text{Queue}[\text{Vertex}]$ is a priority queue of vertices ordered by $D[v]$. This is initialized to contain every vertex in $G$. This supports one operation, $\text{pop}$, which takes the vertex $v$ with the minimum value of $D[v]$, removes it from $Q$, and returns it.
 While $Q$ is nonempty:
 Let $u := \text{pop}(Q)$.
 If $u = V_f$:
 Break out of the loop (found the target).
 If $D[u] = +\infty$:
 Failed: there is no path from $V_i$ to $V_f$.
 Else:
 For each pair $(v, c)$ in the neighbours of $u$:
 Let $d := c + D[u]$.
 If $d < D[v]$:
 $D[v] := d$.
 $L[v] := u$.
 For each pair $(v, c)$ in the neighbours of $u$:
 Let $P: \text{List}[\text{Vertex}] = ()$.
 Let $l := V_f$.
 While $L[l] \neq \text{NIL}$:
 Append $l$ to $P$.
 $l := P[l]$.
 Reverse $P$ and return it.
In Common Lisp we can realize this as follows:
(defun dijkstra (graph source destination)
"Find a path from the source node to the destination in the given graph.
GRAPH is an instance of GRAPH. SOURCE and DESTINATION are integer vertex IDs.
Returns a vector of integer vertex IDs."
;; Table of distances.
(let ((dist (makehashtable :test 'equal)))
(loop for vertex across (graphvertices graph) do
(setf (gethash vertex dist) doublefloatpositiveinfinity))
(setf (gethash source dist) 0.0)
;; Table of previous nodes.
(let ((previous (makehashtable :test 'equal)))
(loop for vertex across (graphvertices graph) do
(setf (gethash vertex previous) nil))
;; Table of neighbor costs.
(let ((neighbors (makeneighbors graph)))
;; Queue.
(let ((q (makehashtable :test 'equal)))
(loop for vertex across (graphvertices graph) do
(setf (gethash vertex q) t))
(loop while (> (hashtablecount q) 0) do
(let ((u (popmin q dist)))
(when (or (= u destination)
(= (gethash u dist) doublefloatpositiveinfinity))
(return))
;; Adjust neighbor distances.
(let ((neighbors (gethash u neighbors)))
(loop for v being the hashkeys of neighbors do
(let ((cost (gethash v neighbors)))
(let ((alt (+ cost (gethash u dist))))
(when (< alt (gethash v dist))
(setf (gethash v dist) alt)
(setf (gethash v previous) u))))))))
;; Build path
(buildpath destination previous))))))
Using the following auxiliary functions:
(defun buildpath (destination previous)
"Given a vertex ID, and a hash table from vertex IDs to vertex IDs,
follow the path starting from DESTINATION through the hash table, and
return a vector of vertex IDs from the last visited node to the DESTINATION."
(let ((path (makearray 0 :adjustable t :elementtype 'integer :fillpointer 0)))
(let ((u destination))
(loop while (gethash u previous) do
(vectorpushextend u path)
(setf u (gethash u previous)))
(vectorpushextend u path)
(reverse path))))
(defun makeneighbors (graph)
"Construct the neighbors map. This is a hash table from vertex IDs to a hash
table of vertex IDs to costs. That is:
Map[ID, Map[ID, Float]]"
(withslots (vertices edges) graph
(let ((neighbors (makehashtable :test 'equal)))
;; Initialize.
(loop for vertex across vertices do
(setf (gethash vertex neighbors) (makehashtable :test 'equal)))
;; Fill.
(loop for edge across edges do
(withslots (start end cost) edge
(setf (gethash end (gethash start neighbors)) cost)
(setf (gethash start (gethash end neighbors)) cost)))
;; Return
neighbors)))
(defun popmin (queue dist)
"The worst priority queue implementation ever written."
(let ((minid nil)
(mindist nil))
(loop for vertex being the hashkeys of queue do
(when (or (null minid)
(< (gethash vertex dist) mindist))
(setf minid vertex)
(setf mindist (gethash vertex dist))))
(assert (not (null minid)))
(assert (not (null mindist)))
(remhash minid queue)
minid))
Now we can build the star graph. This function takes a vector of stars, and creates a graph where the nodes are star IDs, and two nodes are connected if the distance between two stars is less than the threshold:
(defun makegraph (stars dist)
"Create a graph from a vector of stars. Only create edges between
stars that are less than DIST parsecs apart."
(let ((n (length stars)) ; Number of stars.
(edges (makearray 0 :adjustable t :elementtype 'edge :fillpointer 0))) ; edge accumulator
;; Do the Cartesian product.
(loop for i from 0 to ( n 1) do
(loop for j from (+ i 1) to ( n 1) do
(let ((a (elt stars i))
(b (elt stars j)))
(assert (not (= (starid a) (starid b))))
;; Find the distance between the two stars.
(let ((d (euclideandistance (starcartesianposition a)
(starcartesianposition b))))
;; If the distance is within the radius, add the edge to the graph.
(when (<= (value d) (value dist))
(vectorpushextend (makeinstance 'edge
:start (starid a)
:end (starid b)
:cost (value d))
edges))))))
;; Construct the graph from the edges.
(makegraphfromedges edges)))
The HYG database has over 100k stars, but all the stars in the story are very close to the Sun. So we can pare down the graph considerably by dropping every star over 70 light years (22 parsecs) from the Sun:
(defparameter +stars+
(removeif #'(lambda (star)
(> (value (stardistance star)) 22.0))
(copyseq (databasestars +db+))))
Interstellar laser links reach ~16 light years (5 parsecs).
(defparameter +laserlimit+ (makeparsecs 5.0))
(This number is completely arbitrary, and arguably too high.)
Now we build the star graph:
(defparameter +graph+
(makegraph +stars+ +laserlimit+))
(format t "Star graph has ~A vertices and ~A edges.~%~%"
(length (graphvertices +graph+))
(length (graphedges +graph+)))
And run Dijkstra’s algorithm to get the shortest network route and print it out:
(defparameter +path+
(loop for id across (dijkstra +graph+ (starid +bpic+) (starid +g555+))
collecting (find id (databasestars +db+) :key #'starid)))
(format t "Network route has ~A jumps.~%~%" (1 (length +path+)))
;;; Print the network route.
(format t "~12@A ~12@A ~10@A~%" "Start" "End" "Dist")
(format t "~%")
(loop for (a b) on +path+ by #'cdr while b do
(format t "~12@A ~12@A ~8,2fly~%"
(starname a)
(starname b)
(value (parsecstolightyears (stareuclideandistance a b)))))
The output is:
Star graph has 2333 vertices and 32266 edges.
Network route has 7 jumps.
Start End Dist

Bet Pic Gl 238 14.08ly
Gl 238 HIP 27887 12.11ly
HIP 27887 HIP 31293 16.27ly
HIP 31293 HIP 31292 1.10ly
HIP 31292 HIP 58910 15.46ly
HIP 58910 Gl 563.2A 12.99ly
Gl 563.2A Gl 555 6.24ly
Now let’s add up the distances, and compare it to the Euclidean distance:
(format t "Distance from Beta Pictoris to Gliese 555: ~,2fly~%~%"
(value (parsecstolightyears (stareuclideandistance +bpic+ +g555+))))
(let ((length 0.0))
(loop for (a b) on +path+ by #'cdr while b do
(incf length (value (stareuclideandistance a b))))
(format t "Total network route length: ~,2fly~%"
(value (parsecstolightyears (makeparsecs length)))))
This outputs:
Distance from Beta Pictoris to Gliese 555: 70.80ly
Total network route length: 78.25ly
A journey on the network is eight years longer than a straightline journey, but fusion rockets are limited to ~10% of the speed of light, so the straightline ballistic trajectory would take 700 years.
Star Maps
I wanted to make a map of all the stars mentioned in the story, as well as a map showing the network route. Part to verify that the calculations were actually correct, part to have a better understanding of the story’s geography: a picture is easier to understand that a table of numbers.
The simplest way to do this—without writing my own graphics code—is to use a 3D scatter plot. There aren’t any good plotting libraries for Common Lisp, so I used Python’s matplotlib, and used a CSV to transfer the data.
Plotting is always tedious, but happily I was able to use ChatGPT to write most of the plotting code. I asked it to generate a simple example of a scatterplot with labeled points—that is, stars. Then I modified the presentation a bit, changing colours and font sizes, but the static plots have the problem that the perspective makes it hard to know where places really are.
So I showed ChatGPT my plotting code, and asked it to rewrite it to create an animated version where the entire plot is rotated about the vertical axis. It rewrote the script, preserving bitforbit identical output, and added animation support. I got an inscrutable error, showed it to ChatGPT, and it suggested a fix. I’m really happy with this approach.
I’ll start with the output first. Here’s the animated map of all place names mentioned in the story^{3}:
The map uses the inuniverse names: Ctesiphon is Beta Pictoris, Wepwawet is Gliese 555, Ararat is Gliese 570A, Tigranes is HD 35650.
Zoomed in around Gliese 581 for clarity:
And this is the network route from Ctesiphon to Wepwawet:
To transfer data from Common Lisp to Python, I dumped the stars to different
CSVs, one for each plot, with X,Y,Z,Label
as the columns. I won’t post the
code since it’s very ordinary.
The main body of the Python plotting code (sans CSV parsing, etc.) is:
def plot_stars(input_path, output_path, route=False):
"""
Plot the stars from the CSV in `input_path`, writes an MP4
animated scatterplot to `output_path`.
If `route=True`, draws lines between adjacent points.
"""
# Data.
x: list[float] = []
y: list[float] = []
z: list[float] = []
labels: list[str] = []
# Parse the CSV.
# snipped
# Create a figure and a 3D Axes.
dpi = 600
fig = plt.figure(figsize=(2,2), dpi=dpi)
ax = fig.add_subplot(111, projection='3d')
# Create the scatterplot.
scatter = ax.scatter(x, y, z, c='b', marker='.', s=1, alpha=0.5)
# Add labels to each star.
text = [
ax.text(x[i], y[i], z[i] + 0.1, label, fontsize=2) for i, label in enumerate(labels)
]
# Draw the impulses.
impulses = [
ax.plot([x[i], x[i]], [y[i], y[i]], [0, z[i]], ':', c='k', linewidth=0.2) for i in range(len(x))
]
# Are we plotting a route? If so, draw the lines between the stars.
if route:
lines = [
ax.plot(
[x[i], x[i+1]],
[y[i], y[i+1]],
[z[i], z[i+1]],
color='r',
linewidth=0.1
)
for i in range(len(labels)1)
]
else:
lines = []
# Plot the origin plane.
GAP = 1
X, Y = np.meshgrid(
np.linspace(min(x)  GAP, max(x) + GAP, 10),
np.linspace(min(y)  GAP, max(y) + GAP, 10)
)
Z = np.zeros_like(X)
origin_plane = ax.plot_wireframe(X, Y, Z, rstride=1, cstride=1, linewidths=0.1)
# Hide the axes and grid planes.
plt.axis("off")
# Define the animation function.
def animate(i):
ax.view_init(elev=30, azim=i)
# Create the animation object.
anim = animation.FuncAnimation(fig, animate, frames=360, interval=20, blit=False)
# Save the animation as an MP4 file.
anim.save(output_path, fps=30, extra_args=['vcodec', 'libx264'], dpi=dpi)
Conclusion
All of this had very little consequence to the story: a few throwaway sentences. But I see it as a matter of respect for the audience. If it’s hard SF, it’s hard SF.
Appendix: Equatorial Coordinates
This section describes to how accurately represent positions in equatorial coordinates, and how to convert them to Cartesian coordinates for ease of use.
None of this is strictly necessary, because the HYG database has star positions both in equatorial and Cartesian coordinates, but I had to solve this problem, for two reasons. The first is that software is crystallized and verified understanding: you know you understand something when you can write a computer program that concretizes that understanding. The second is that when software solves a problem, it solves it in perpetuity.
Representation
A position in equatorial coordinates is a right ascension (RA, an angle, in hoursminutesseconds), a declination (DEC, an angle, in degreesminutesseconds), and a distance from the Sun in parsecs.
The equatorialposition
class represents this:
(defclass equatorialposition ()
((rightascension :reader rightascension
:initarg :rightascension
:type hmsdegrees
:documentation "The right ascension in HMS.")
(declination :reader declination
:initarg :declination
:type dmsdegrees
:documentation "The declination in DMS.")
(distance :reader distance
:initarg :distance
:type parsecs
:documentation "The distance in parsecs."))
(:documentation "A position in equatorial (RA, DEC, DIST) coordinates."))
We could convert the right ascension and declination to angles eagerly, but I’d rather represent things explicitly, and add conversion functions as needed.
We have three representations of angles: decimal degrees, hoursminutesseconds, and degreesminutesseconds. Each corresponds to a CLOS class:
(defclass decimaldegrees ()
((value :reader value
:initarg :value
:type real))
(:documentation "Represents an angle in decimal degrees."))
(defclass hmsdegrees ()
((hours :reader hours
:initarg :hours
:type real)
(minutes :reader minutes
:initarg :minutes
:type real)
(seconds :reader seconds
:initarg :seconds
:type real))
(:documentation "Represents an angle in HMS (hoursminutesseconds) format."))
(defclass dmsdegrees ()
((degrees :reader degrees
:initarg :degrees
:type real)
(minutes :reader minutes
:initarg :minutes
:type real)
(seconds :reader seconds
:initarg :seconds
:type real))
(:documentation "Represents an angle in DMS (degreesminutesseconds) format."))
We use initializeinstance
methods to verify that all the values are in
range^{4}:
(defmethod initializeinstance :after ((d decimaldegrees) &key)
(let ((d (value d)))
(unless (and (> d 180.0) (<= d 180.0))
(error "Value out of range for decimal degrees: ~A" d))))
(defmethod initializeinstance :after ((angle hmsdegrees) &key)
(withslots (hours minutes seconds) angle
(unless (and (>= hours 0.0) (< hours 24.0))
(error "Invalid value for hours in hmsdegrees: ~a" hours))
(unless (and (>= minutes 0.0) (< minutes 60.0))
(error "Invalid value for minutes in hmsdegrees: ~a" minutes))
(unless (and (>= seconds 0.0) (< seconds 60.0))
(error "Invalid value for seconds in hmsdegrees: ~a" seconds))))
(defmethod initializeinstance :after ((angle dmsdegrees) &key)
(withslots (degrees minutes seconds) angle
(unless (and (> degrees 90.0) (<= degrees 90.0))
(error "Invalid value for degrees in hmsdegrees: ~a" degrees))
(unless (and (>= minutes 0.0) (< minutes 60.0))
(error "Invalid value for minutes in hmsdegrees: ~a" minutes))
(unless (and (>= seconds 0.0) (< seconds 60.0))
(error "Invalid value for seconds in hmsdegrees: ~a" seconds))))
Angle Conversion
To convert from HMS and DMS to decimal, we have these functions:
(defun hmstodecimal (hms)
"Convert an HMS (hoursminutesseconds) angle to decimal degrees."
(withslots (hours minutes seconds) hms
(let ((d (+ (* 15.0 hours)
(/ (* 15.0 minutes) 60.0)
(/ (* 15.0 seconds) 3600.0))))
(makeinstance 'decimaldegrees :value d))))
(defun dmstodecimal (dms)
"Convert a DMS (degreesminutesseconds) angle to decimal degrees."
(withslots (degrees minutes seconds) dms
(let ((d (* (sign degrees)
(+ (abs degrees)
(/ minutes 60.0)
(/ seconds 3600.0)))))
(makeinstance 'decimaldegrees :value d))))
(defun sign (x)
(cond ((< x 0.0)
1)
((= x 0.0)
0)
(t
1.0)))
And we can use these like so:
CLUSER> (makeinstance 'hmsdegrees :hours 7.2 :minutes 2.24 :seconds 1.42)
#<HMSDEGREES 7.2h2.2m1.4s>
CLUSER> (hmstodecimal *)
#<DECIMALDEGREES 108.6°>
CLUSER> (makeinstance 'dmsdegrees :degrees 51.2 :minutes 2.24 :seconds 1.42)
#<DMSDEGREES 51.2°2.2m1.4s>
CLUSER> (dmstodecimal *)
#<DECIMALDEGREES 51.2°>
Implementing decimaltohms
and decimaltodms
is left as an exercise to the
reader.
Equatorial to Cartesian
I got the formula from the Atomic Rockets website:
(defun equatorialtocartesian (pos)
"Convert a position from equatorial to cartesian coordinates."
(withslots (rightascension declination distance) pos
(let ((φ (value (hmstodecimal rightascension)))
(θ (value (dmstodecimal declination)))
(ρ (value distance)))
(let ((rvect (* ρ (cosr θ))))
(let ((x (* rvect (cosr φ)))
(y (* rvect (sinr φ)))
(z (* ρ (sinr θ))))
(makeinstance 'cartesianposition
:x (makeparsecs x)
:y (makeparsecs y)
:z (makeparsecs z)))))))
(defun rad (x)
(* x 0.0174532925))
(defun sinr (x) (sin (rad x)))
(defun cosr (x) (cos (rad x)))
We can use this like so:
CLUSER> (defparameter tauceti
(makeinstance 'equatorialposition
:rightascension (makeinstance 'hmsdegrees :hours 1 :minutes 41 :seconds 45)
:declination (makeinstance 'dmsdegrees :degrees 16.0 :minutes 12.0 :seconds 0.0)
:distance (makeparsecs 3.61)))
#<EQUATORIALPOSITION RA=1.0h41.0m45.0s DEC=16.0°12.0m.0s D=3.6pc>
CLUSER> (equatorialtocartesian tauceti)
#<CARTESIANPOSITION X=3.1pc Y=1.5pc Z=1.0pc>
Footnotes

What about drift? The HYG database has the velocity vectors of the stars, but it doesn’t matter. In the short run, stars move too slowly for their motion to affect travel times appreciably. In the long run, it’s too hard to solve: nbody with very inaccurate data spanning tens of kiloyears. The story doesn’t take place far enough into the future for this to matter, so all coordinates are J2000. ↩

The Cartesian coordinate system used by the HYG database has:

$+X$ towards the vernal point in J2000 (RA 0h, DEC 0°, in the constellation Pisces).

$+Y$ towards RA 6h, DEC 0° (in Monoceros).

$+Z$ towards the north celestial pole.
The $Z=0$ plane is the plane of the ecliptic. ↩


This is a simple check that the math is correct. Note that Gliese 581 is close to the plane of the ecliptic, while Ctesiphon (Beta Pictoris) is far to galactic south. Gliese 581 has a declination of 7°, while Beta Pictoris is 51°. If you look at a map of the constellations, Libra (where Gliese 581 is) is just off the ecliptic, while Pictor is far to the south. They are also in roughly opposite directions in the sky. ↩

Much of this repetitive errorchecking code was written by ChatGPT, but I had to anneal the inequalities to get it right. ↩