Retrograde Mars: why it stops and loops backward through the stars (Python)

Play: watch Mars for two years

Mars drifts east among the stars by about half a degree a day, close to the ecliptic, the line the Sun follows through the stars in a year and the plane of Earth's orbit seen from inside. Every 26 months it slows, stops, runs west for about ten weeks, and stops again, next from 10 January 2027. The game shows the real bright stars along the ecliptic, east to the left as when you face south, Mars with its track of the past 200 days, and beside them the two orbits seen from above, paused on 1 December 2026. Play runs time, − and + set the speed in days per second (below zero it runs backward), and a drag across the sky moves time by hand. Watch the track bend back above Regulus while Earth catches up with Mars. The readouts give the date, Mars's motion in degrees per day, its distance in AU, and its magnitude, where smaller is brighter; the motion passes through zero and falls to −0.40°/day on 19 February.

Download the game

Show code

// game.js: watch Mars move against the stars along the ecliptic, and the two orbits from above.
// Units: days, AU, degrees. Every position is in the ecliptic frame of J2000, which is also the
// frame of the star positions, so the stars stand still. The notebook writes the orbital elements,
// the stars, and the magnitude law into #data.
"use strict";

const DATA = JSON.parse(document.getElementById("data").textContent);
const DEG = Math.PI / 180;
const NIGHT = "#05070d";             // the sky, in both themes
const STARLIGHT = "#e8ecf4", LABEL = "#c9d0e3", MARS_SKY = "#efd5ae", GUIDE = "#8e97ab";
const WEST = "#e3826b";              // the accent on the night sky: the westward run, and only that
const DAYS_PER_PX = 1;               // a drag across the sky moves time by this much per pixel

// ---- the state: everything that changes, in one object
const DEFAULTS = {t: DATA.t_start, speed: DATA.speed_start, running: false};
const state = {...DEFAULTS};

// ---- where a planet is at time t (days from J2000): elements, Kepler's equation, three rotations
const wrap180 = (x) => ((x + 180) % 360 + 360) % 360 - 180;

function kepler(M, e) {              // solve E - e sin E = M by Newton steps
  let E = M + e * Math.sin(M);
  for (let k = 0; k < 10; k++) {
    const dE = (E - e * Math.sin(E) - M) / (1 - e * Math.cos(E));
    E -= dE;
    if (Math.abs(dE) < 1e-12) break;
  }
  return E;
}

function helio(el, t) {              // heliocentric position in AU
  const T = t / 36525, at = (k) => el[k][0] + el[k][1] * T;   // value plus rate per century
  const a = at("a"), e = at("e"), I = at("I") * DEG, Om = at("Omega") * DEG, varpi = at("varpi") * DEG;
  const M = wrap180(at("L") - at("varpi")) * DEG, E = kepler(M, e), w = varpi - Om;
  const xp = a * (Math.cos(E) - e), yp = a * Math.sqrt(1 - e * e) * Math.sin(E);   // in the orbit's plane
  const x1 = xp * Math.cos(w) - yp * Math.sin(w), y1 = xp * Math.sin(w) + yp * Math.cos(w);  // turn by ω
  const y2 = y1 * Math.cos(I), z = y1 * Math.sin(I);                                        // tilt by I
  return [x1 * Math.cos(Om) - y2 * Math.sin(Om), x1 * Math.sin(Om) + y2 * Math.cos(Om), z];  // turn by Ω
}

function seen(t) {                   // Mars from Earth: longitude, latitude, distance, and both radii
  const p = helio(DATA.elements.mars, t), q = helio(DATA.elements.earth, t);
  const d = [p[0] - q[0], p[1] - q[1], p[2] - q[2]], delta = Math.hypot(...d);
  return {lam: Math.atan2(d[1], d[0]) / DEG, beta: Math.asin(d[2] / delta) / DEG, delta,
    r: Math.hypot(...p), R: Math.hypot(...q), mars: p, earth: q};
}

function magnitude(s) {              // V = H + 5 log10(r Δ) + phase term, α the Sun-Mars-Earth angle
  const alpha = Math.acos((s.r ** 2 + s.delta ** 2 - s.R ** 2) / (2 * s.r * s.delta)) / DEG;
  const [c1, c2] = DATA.phase;
  return DATA.h + 5 * Math.log10(s.r * s.delta) + c1 * alpha + c2 * alpha ** 2;
}

const motion = (t) => wrap180(seen(t + 0.5).lam - seen(t - 0.5).lam);   // degrees per day

function center(t) {                 // Mars's longitude averaged over ±60 days, so the view drifts slowly
  let c = 0, s = 0;
  for (let k = -6; k <= 6; k++) {
    const lam = seen(t + 10 * k).lam * DEG;
    c += Math.cos(lam);
    s += Math.sin(lam);
  }
  return Math.atan2(s, c) / DEG;
}

// ---- one tick: real seconds times days per second, capped so a hidden tab does not jump
function tick(dtReal) {
  if (!state.running) return;
  const t = state.t + Math.min(dtReal, 0.05) * DATA.speeds[state.speed];
  state.t = Math.min(DATA.t_max, Math.max(DATA.t_min, t));
  if (state.t !== t) state.running = false;   // stop at the ends of the elements' range
}

// ---- drawing
const $ = (id) => document.getElementById(id);
const color = (name) => getComputedStyle(document.documentElement).getPropertyValue(`--${name}`).trim();
const FONT = "14px system-ui, sans-serif";
const views = {sky: {id: "sky"}, orbits: {id: "orbits"}};
const radius = (V) => Math.max(0.6, 1 + 0.55 * (4.5 - V));      // px; the same law for stars and Mars
const opacity = (V) => Math.min(1, 0.35 + 0.15 * (6 - V));
const MONTHS = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"];

function dateText(t) {
  const d = new Date(Date.UTC(2000, 0, 1, 12) + t * 86400000);
  return `${d.getUTCDate()} ${MONTHS[d.getUTCMonth()]} ${d.getUTCFullYear()}`;
}

function fit(v) {                    // size a canvas to its box, sharp on any screen
  const c = $(v.id), box = c.parentElement.getBoundingClientRect(), dpr = window.devicePixelRatio || 1;
  c.width = Math.round(box.width * dpr);
  c.height = Math.round(box.height * dpr);
  v.ctx = c.getContext("2d");
  v.ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
  Object.assign(v, {w: box.width, h: box.height});
}

function layout() {
  const W = document.documentElement.clientWidth, H = document.documentElement.clientHeight;
  document.body.className = H > W ? "upright" : "wide";
  fit(views.sky);
  fit(views.orbits);
  render();
}

function disk(ctx, x, y, r, fill) {
  ctx.beginPath();
  ctx.arc(x, y, r, 0, 2 * Math.PI);
  ctx.fillStyle = fill;
  ctx.fill();
}

function label(ctx, text, x, y, style, align = "left") {
  ctx.font = FONT;
  ctx.textAlign = align;
  ctx.fillStyle = style;
  ctx.fillText(text, x, y);
}

function path(ctx, pts, style, width, dash = []) {
  ctx.beginPath();
  pts.forEach((p, k) => (k ? ctx.lineTo(p[0], p[1]) : ctx.moveTo(p[0], p[1])));
  ctx.strokeStyle = style;
  ctx.lineWidth = width;
  ctx.setLineDash(dash);
  ctx.stroke();
  ctx.setLineDash([]);
}

// The sky: at least 36° of ecliptic longitude and 20° of latitude, east to the left as when you face south
function drawSky(v, now) {
  const {ctx, w, h} = v, k = Math.min(h / 20, w / 36), lc = center(state.t);
  const at = (lam, beta) => [w / 2 + wrap180(lc - lam) * k, h / 2 - beta * k];
  ctx.fillStyle = NIGHT;
  ctx.fillRect(0, 0, w, h);
  for (const [lam, beta, V] of DATA.stars) {
    const [x, y] = at(lam, beta);
    if (x < -6 || x > w + 6 || y < -6 || y > h + 6) continue;
    ctx.globalAlpha = opacity(V);
    disk(ctx, x, y, radius(V), STARLIGHT);
  }
  ctx.globalAlpha = 1;
  path(ctx, [[0, h / 2], [w, h / 2]], GUIDE, 1, [6, 6]);       // the ecliptic
  for (let lam = 0; lam < 360; lam += 10) {                    // longitude marks along the bottom
    const [x] = at(lam, 0);
    if (x < 20 || x > w - 20) continue;
    path(ctx, [[x, h - 22], [x, h - 16]], GUIDE, 1);
    label(ctx, `${lam}°`, x, h - 4, GUIDE, "center");
  }
  for (const s of DATA.names) {
    const [x, y] = at(s.lam, s.beta);
    if (x > 0 && x < w) label(ctx, s.name, x + 7, y - 5, LABEL);   // above right, clear of the ecliptic
  }
  label(ctx, "← east", 8, 18, LABEL);
  label(ctx, "west →", w - 8, 18, LABEL, "right");
  // the track of the past 200 days, then Mars itself, its size from its magnitude
  const trail = [];
  for (let j = 0; j < DATA.trail_points; j++) {
    const tj = state.t - DATA.trail_days * (1 - j / (DATA.trail_points - 1)), s = seen(tj);
    trail.push([...at(s.lam, s.beta), motion(tj) < 0]);
  }
  ctx.globalAlpha = 0.6;
  path(ctx, trail, MARS_SKY, 1.5);
  for (const [x, y] of trail) disk(ctx, x, y, 1.5, MARS_SKY);
  ctx.globalAlpha = 1;
  for (let j = 1; j < trail.length; j++)                       // the westward stretches on top, in the accent
    if (trail[j][2] && trail[j - 1][2]) path(ctx, [trail[j - 1], trail[j]], WEST, 2.5);
  for (const [x, y, west] of trail) if (west) disk(ctx, x, y, 2.2, WEST);
  const [x, y] = at(now.lam, now.beta);
  disk(ctx, x, y, radius(now.V), MARS_SKY);
  label(ctx, "Mars", x + radius(now.V) + 5, y - 6, MARS_SKY);
}

// The orbits from above, from the north side of the ecliptic: Sun, Earth, Mars, and the line of sight
function drawOrbits(v, now, c) {
  const {ctx, w, h} = v, k = (Math.min(w, h) / 2 - 22) / 1.68, cx = w / 2, cy = h / 2;   // room for names
  const at = (p) => [cx + p[0] * k, cy - p[1] * k];
  ctx.clearRect(0, 0, w, h);
  path(ctx, ORBITS.earth.map(at), c.muted, 1);
  path(ctx, ORBITS.mars.map(at), c.muted, 1);
  const e = at(now.earth), m = at(now.mars), len = Math.hypot(m[0] - e[0], m[1] - e[1]);
  const far = [e[0] + (m[0] - e[0]) / len * 2 * Math.max(w, h), e[1] + (m[1] - e[1]) / len * 2 * Math.max(w, h)];
  path(ctx, [e, far], c.muted, 1.2, [4, 4]);                  // the line of sight, out to the stars
  disk(ctx, cx, cy, 6, c.ink);
  disk(ctx, e[0], e[1], 5, c.second);
  disk(ctx, m[0], m[1], 5, c.ink);
  label(ctx, "Sun", cx, cy + 22, c.ink, "center");
  // Earth's name toward the Sun, Mars's outward and to one side of the line of sight: neither crosses a line
  const unit = (x, y) => { const n = Math.hypot(x, y) || 1; return [x / n, y / n]; };
  const los = unit(m[0] - e[0], m[1] - e[1]), out = unit(m[0] - cx, m[1] - cy);
  const side = -los[1] * out[0] + los[0] * out[1] < 0 ? -1 : 1;
  tag(v, "Earth", e, unit(cx - e[0], cy - e[1]), c.second);
  tag(v, "Mars", m, unit(out[0] - los[1] * side, out[1] + los[0] * side), c.ink);
}

function tag(v, text, p, u, style) {  // a name offset from a dot in direction u, kept inside the view
  const {ctx, w, h} = v;
  ctx.font = FONT;
  const tw = ctx.measureText(text).width;
  const d = 10 + Math.abs(u[0]) * tw / 2 + Math.abs(u[1]) * 8;
  const x = Math.min(w - tw / 2 - 2, Math.max(tw / 2 + 2, p[0] + u[0] * d));
  const y = Math.min(h - 9, Math.max(9, p[1] + u[1] * d));
  ctx.textBaseline = "middle";
  label(ctx, text, x, y, style, "center");
  ctx.textBaseline = "alphabetic";
}

let ORBITS = null;                   // one revolution of each planet, drawn as its orbit
function orbits() {
  const ring = (el, period) => Array.from({length: 121}, (_, j) => helio(el, DATA.t_start + period * j / 120));
  ORBITS = {earth: ring(DATA.elements.earth, 365.256), mars: ring(DATA.elements.mars, 686.98)};
}

function render() {
  const c = Object.fromEntries(["bg", "fg", "ink", "accent", "second", "muted"].map((n) => [n, color(n)]));
  const now = seen(state.t);
  now.V = magnitude(now);
  drawSky(views.sky, now);
  drawOrbits(views.orbits, now, c);
  const rate = motion(state.t);
  $("date").textContent = dateText(state.t);
  $("motion").classList.toggle("acc", rate < 0);           // the accent only while Mars runs west
  $("rate").textContent = `${rate < 0 ? "−" : "+"}${Math.abs(rate).toFixed(2)}°/day ${rate < 0 ? "westward" : "eastward"}`;
  $("delta").textContent = `${now.delta.toFixed(3)} AU`;
  $("mag").textContent = `${now.V < 0 ? "−" : "+"}${Math.abs(now.V).toFixed(1)}`;
  $("speed-val").textContent = `${String(DATA.speeds[state.speed]).replace("-", "−")} days/s`;
  $("go").textContent = state.running ? "Pause" : "Play";
}

// ---- controls, keys, and the loop
function set(partial) {
  Object.assign(state, partial);
  render();
}
const nudge = (sign) => set({speed: Math.min(DATA.speeds.length - 1, Math.max(0, state.speed + sign))});
const go = () => set({running: !state.running});
const reset = () => set({...DEFAULTS});

let drag = null;                     // a horizontal drag on the sky moves time and pauses the clock
$("sky").addEventListener("pointerdown", (e) => {
  $("sky").setPointerCapture(e.pointerId);
  drag = {x: e.clientX, t: state.t};
  set({running: false});
});
$("sky").addEventListener("pointermove", (e) => {
  if (!drag) return;
  const t = drag.t + (e.clientX - drag.x) * DAYS_PER_PX;
  set({t: Math.min(DATA.t_max, Math.max(DATA.t_min, t))});
});
$("sky").addEventListener("pointerup", () => { drag = null; });
$("sky").addEventListener("pointercancel", () => { drag = null; });
$("go").addEventListener("click", go);
$("speed-minus").addEventListener("click", () => nudge(-1));
$("speed-plus").addEventListener("click", () => nudge(1));
$("reset").addEventListener("click", reset);
document.addEventListener("keydown", (e) => {
  const actions = {" ": go, "-": () => nudge(-1), "=": () => nudge(1), "+": () => nudge(1), r: reset};
  const act = actions[e.key.toLowerCase()];
  if (act && !(e.target instanceof HTMLButtonElement && e.key === " ")) {
    e.preventDefault();
    act();
  }
});

let last = null;
function frame(now) {                // called by the browser once per screen refresh
  const dt = last === null ? 0 : (now - last) / 1000;
  last = now;
  tick(dt);
  render();
  requestAnimationFrame(frame);
}

orbits();
window.addEventListener("resize", layout);
new ResizeObserver(layout).observe($("orbitbox"));   // the boxes also change when the readouts rewrap
window.matchMedia("(prefers-color-scheme: dark)").addEventListener("change", render);
layout();
requestAnimationFrame(frame);
window.game = {state, set, render};

Mars stops and moves backward

Press Play and from mid-December the dot slows. It stands still on 10 January 2027, runs west past the longitude of Regulus, stands still again on 1 April, and moves on east. In the view from above, the middle of the backward run is the opposition of 19 February, when the Sun, Earth, and Mars line up. The hidden cell holds a Python copy of the game's model.

Show code

import gzip
import json
from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt
import pooch
from scipy.optimize import brentq, newton

plt.rcParams.update({
    "figure.figsize": (7, 3.6), "figure.dpi": 110,
    "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25,
    "font.size": 11, "lines.linewidth": 1.8,
})
INK, ACCENT, SECOND, MUTED = "#1f2a44", "#c8553d", "#2a7f9e", "#8a8f98"
NIGHT, STARLIGHT, LABEL, MARS_SKY = "#05070d", "#e8ecf4", "#c9d0e3", "#efd5ae"   # the game's sky, stars, names, and Mars

# Keplerian elements for approximate positions of the major planets, valid 1800 to 2050:
# E. M. Standish, JPL Solar System Dynamics, Table 1. J2000 ecliptic and equinox.
# Each entry: (value at J2000, rate per Julian century); a in AU, angles in degrees.
ELEMENTS = {
    "earth": dict(a=(1.00000261, 0.00000562), e=(0.01671123, -0.00004392), I=(-0.00001531, -0.01294668),
                  L=(100.46457166, 35999.37244981), varpi=(102.93768193, 0.32327364), Omega=(0.0, 0.0)),
    "mars": dict(a=(1.52371034, 0.00001847), e=(0.09339410, 0.00007882), I=(1.84969142, -0.00813131),
                 L=(-4.55343205, 19140.30268499), varpi=(-23.94362959, 0.44441088), Omega=(49.55953891, -0.29257343)),
}
H_MARS, PHASE = -1.601, (0.02267, -0.0001302)   # Mallama and Hilton (2018): V = H + 5 log10(r Δ) + c1 α + c2 α²
J2000 = np.datetime64("2000-01-01T12:00")
DEG = np.pi / 180

def day(when):
    """Days from J2000 (noon on 1 January 2000) for a date such as '2027-02-19'."""
    return (np.datetime64(when) - J2000) / np.timedelta64(1, "D")

def date(t):
    return str(J2000 + np.timedelta64(int(round(t * 1440)), "m"))[:10]

def wrap(x):
    """An angle in degrees, brought into [-180, 180)."""
    return (np.asarray(x) + 180) % 360 - 180

def elements(body, t):
    T = np.asarray(t, dtype=float) / 36525
    return {k: v0 + v1 * T for k, (v0, v1) in ELEMENTS[body].items()}

def kepler(M, e):
    """Eccentric anomaly E from E - e sin E = M, in radians, by Newton's method."""
    return newton(lambda E: E - e * np.sin(E) - M, M + e * np.sin(M),
                  fprime=lambda E: 1 - e * np.cos(E), tol=1e-12, maxiter=50)

def helio(body, t):
    """Heliocentric ecliptic position in AU, shape (3,) + shape of t."""
    el = elements(body, t)
    a, e = el["a"], el["e"]
    I, Om, w = el["I"] * DEG, el["Omega"] * DEG, (el["varpi"] - el["Omega"]) * DEG
    M = wrap(el["L"] - el["varpi"]) * DEG
    E = kepler(M, e)
    xp, yp = a * (np.cos(E) - e), a * np.sqrt(1 - e**2) * np.sin(E)
    x1, y1 = xp * np.cos(w) - yp * np.sin(w), xp * np.sin(w) + yp * np.cos(w)
    y2, z = y1 * np.cos(I), y1 * np.sin(I)
    return np.array([x1 * np.cos(Om) - y2 * np.sin(Om), x1 * np.sin(Om) + y2 * np.cos(Om), z])

def geocentric(t, target="mars", observer="earth"):
    """Longitude and latitude in degrees, distance Δ in AU, of target seen from observer."""
    d = helio(target, t) - helio(observer, t)
    delta = np.linalg.norm(d, axis=0)
    return np.degrees(np.arctan2(d[1], d[0])) % 360, np.degrees(np.arcsin(d[2] / delta)), delta

def daily_motion(t, **bodies):
    """dλ/dt in degrees per day, a central difference over half a day either side."""
    return wrap(geocentric(t + 0.5, **bodies)[0] - geocentric(t - 0.5, **bodies)[0])

def magnitude(t):
    """Mars's V magnitude; α is the Sun-Mars-Earth angle in degrees."""
    r, R = np.linalg.norm(helio("mars", t), axis=0), np.linalg.norm(helio("earth", t), axis=0)
    delta = geocentric(t)[2]
    alpha = np.degrees(np.arccos((r**2 + delta**2 - R**2) / (2 * r * delta)))
    return H_MARS + 5 * np.log10(r * delta) + PHASE[0] * alpha + PHASE[1] * alpha**2

def roots(f, t0, t1, ok=lambda a, b: True):
    """Every root of f between t0 and t1: a sign change from one day to the next, refined by brentq."""
    t = np.arange(t0, t1, 1.0)
    y = f(t)
    return [brentq(f, t[i], t[i + 1]) for i in np.nonzero(np.sign(y[:-1]) != np.sign(y[1:]))[0]
            if ok(y[i], y[i + 1])]

def stations(t0, t1, **bodies):
    return roots(lambda t: daily_motion(t, **bodies), t0, t1)

def oppositions(t0, t1):
    """Mars 180° from the Sun in longitude; the Sun's longitude is that of minus Earth's position."""
    def f(t):
        e = helio("earth", t)
        return wrap(geocentric(t)[0] - np.degrees(np.arctan2(-e[1], -e[0])) - 180)
    return roots(f, t0, t1, ok=lambda a, b: abs(a) < 90 and abs(b) < 90)   # skip the jump at conjunction

def closest(t_opp):
    """Closest approach near an opposition: the root of dΔ/dt."""
    rate = lambda t: geocentric(t + 0.5)[2] - geocentric(t - 0.5)[2]
    return brentq(rate, t_opp - 30, t_opp + 30)

def crosses_itself(lam, beta):
    """True if the track (unwrapped λ, β) has two segments that cross."""
    p, q = np.column_stack([lam[:-1], beta[:-1]]), np.column_stack([lam[1:], beta[1:]])
    d = q - p
    cross = lambda u, v: u[..., 0] * v[..., 1] - u[..., 1] * v[..., 0]
    denom = cross(d[:, None], d[None, :])
    with np.errstate(divide="ignore", invalid="ignore"):
        s = cross(p[None, :] - p[:, None], d[None, :]) / denom
        u = cross(p[None, :] - p[:, None], d[:, None]) / denom
    hit = (s > 0) & (s < 1) & (u > 0) & (u < 1)
    return bool(np.triu(hit, k=2).any())

def track(t_mid, half=160):
    t = np.arange(t_mid - half, t_mid + half + 1, 1.0)
    lam, beta, _ = geocentric(t)
    return t, np.degrees(np.unwrap(np.radians(lam))), beta

# The Yale Bright Star Catalogue, 5th revised edition, CDS catalog V/50: positions for J2000.
path = pooch.retrieve(
    "https://cdsarc.cds.unistra.fr/ftp/V/50/catalog.gz",
    known_hash="sha256:3dc44b1e90be8fbe5bcc7656032560f51275f985c7e3f783c9028e1838ec7bed",
)
with gzip.open(path, "rt", encoding="ascii") as f:
    records = [line.ljust(197) for line in f.read().splitlines()]
rows = [r for r in records if r[75:83].strip()]

def column(first, last):
    """Bytes first to last (1-based, as in the CDS ReadMe) of every row, as floats."""
    return np.array([r[first - 1:last].strip() or np.nan for r in rows], dtype=float)

hr = column(1, 4).astype(int)
ra = np.radians(15 * (column(76, 77) + column(78, 79) / 60 + column(80, 83) / 3600))
dec_sign = np.where([r[83] == "-" for r in rows], -1, 1)
dec = np.radians(dec_sign * (column(85, 86) + column(87, 88) / 60 + column(89, 90) / 3600))
vmag = column(103, 107)
EPS = 23.4392911 * DEG                           # the tilt of the ecliptic against the equator, J2000
star_lam = np.degrees(np.arctan2(np.sin(ra) * np.cos(EPS) + np.tan(dec) * np.sin(EPS), np.cos(ra))) % 360
star_beta = np.degrees(np.arcsin(np.sin(dec) * np.cos(EPS) - np.cos(dec) * np.sin(EPS) * np.sin(ra)))
near = np.abs(star_beta) <= 12
NAMES = {1457: "Aldebaran", 2990: "Pollux", 3982: "Regulus", 5056: "Spica", 6134: "Antares"}
regulus = hr == 3982

def radius(V):
    """Dot radius in CSS px, the game's law for stars and Mars alike."""
    return np.maximum(0.6, 1 + 0.55 * (4.5 - np.asarray(V)))

def opacity(V):
    return np.minimum(1, 0.35 + 0.15 * (6 - np.asarray(V)))

def draw_sky(ax, lam0, lam1, b0, b1, names=True, zoom=1):
    """The game's sky between two longitudes and two latitudes, east (larger λ) to the left; zoom scales the dots."""
    ax.set_facecolor(NIGHT)
    inside = near & (star_lam > lam0 - 1) & (star_lam < lam1 + 1)
    ax.scatter(star_lam[inside], star_beta[inside], s=(2 * 0.75 * zoom * radius(vmag[inside])) ** 2, lw=0,
               color=STARLIGHT, alpha=opacity(vmag[inside]), zorder=1)
    ax.axhline(0, color=MUTED, lw=1, ls="--", zorder=0)
    if names:
        for n, name in NAMES.items():
            k = hr == n
            if lam0 < star_lam[k][0] < lam1:
                ax.text(star_lam[k][0] - 0.6, star_beta[k][0] - 0.6, name, color=LABEL, va="top")
    ax.set(xlim=(lam1, lam0), ylim=(b0, b1), aspect="equal")
    ax.grid(False)

def draw_orbits(ax, dates):
    """The orbits from above (north of the ecliptic), Earth and Mars at the dates, lines of sight extended."""
    for body, period in [("earth", 365.256), ("mars", 686.980)]:
        x, y, _ = helio(body, dates[0] + np.linspace(0, period, 241))
        ax.plot(x, y, color=MUTED, lw=1)
    ax.plot(0, 0, "o", ms=8, color=INK)
    for j, t in enumerate(dates):
        alpha = 0.3 + 0.7 * j / (len(dates) - 1)
        e, m = helio("earth", t), helio("mars", t)
        far = e + 3.2 * (m - e) / np.linalg.norm(m - e)
        ax.plot([e[0], far[0]], [e[1], far[1]], color=INK, lw=1, alpha=0.6 * alpha)
        ax.plot(*e[:2], "o", ms=6, color=SECOND, alpha=alpha)
        ax.plot(*m[:2], "o", ms=6, color=INK, alpha=alpha)
    ax.set(xlim=(-2.7, 1.9), ylim=(-1.9, 1.9), aspect="equal", xlabel="x / AU", ylabel="y / AU")

t_opp_2027 = oppositions(day("2027-01-01"), day("2027-06-01"))[0]
print(f"catalog.gz, {Path(path).stat().st_size / 1024:.0f} kB: {near.sum()} of {len(rows)} stars lie within 12° of the ecliptic; "
      f"Mars is at opposition on {date(t_opp_2027)}")

catalog.gz, 560 kB: 1706 of 9096 stars lie within 12° of the ecliptic; Mars is at opposition on 2027-02-19

The game measures the sky in ecliptic coordinates. Longitude λ is the angle along the ecliptic, eastward from the vernal equinox, where the Sun crosses the celestial equator in March. Latitude β is the angle north or south of the ecliptic. A star keeps its λ and β; a planet does not. The line of sight to Mars is the difference of two arrows from the Sun, \(\mathbf d = \mathbf r_\text{Mars} - \mathbf r_\text{Earth}\). With \(\mathbf d = (x, y, z)\),

\[\lambda = \operatorname{atan2}(y, x), \qquad \beta = \arcsin\frac{z}{\Delta}, \qquad \Delta = |\mathbf d| .\]

The still shows November to May, a dot every five days, and the view from above at seven dates.

Show code

st_2027 = stations(day("2026-11-01"), day("2027-06-01"))
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3.4), gridspec_kw=dict(width_ratios=[2.5, 1]))
draw_sky(ax1, 134, 166, -2.5, 7.5)
t = np.arange(day("2026-11-01"), day("2027-06-01"), 5.0)
lam, beta, _ = geocentric(t)
back = daily_motion(t) < 0
ax1.plot(lam[~back], beta[~back], "o", ms=4, mfc="none", mec=MARS_SKY, mew=1, zorder=2)    # eastward: rings
ax1.plot(lam[back], beta[back], "o", ms=4.5, color=ACCENT, zorder=3)                          # westward: filled
for ts in st_2027:
    l, b, _ = geocentric(ts)
    ax1.annotate(date(ts), (l, b + 0.4), xytext=(l, b + 2.3), color=LABEL, ha="center",
                 arrowprops=dict(arrowstyle="-", color=LABEL, lw=0.8))
ax1.set(xlabel="ecliptic longitude λ / degrees", ylabel="β / degrees")
draw_orbits(ax2, np.linspace(day("2026-12-01"), day("2027-04-20"), 7))
e_last, m_last = helio("earth", day("2027-04-20")), helio("mars", day("2027-04-20"))
ax2.text(0.2, 0, "Sun", color=INK, va="center")
ax2.text(e_last[0] + 0.18, e_last[1] + 0.02, "Earth", color=SECOND, va="center")   # inside its orbit
ax2.text(m_last[0] - 0.12, m_last[1] - 0.1, "Mars", color=INK, ha="right", va="top")
fig.tight_layout()
plt.show()

Left: the sky from λ = 166° to 134°, Mars every five days from November 2026 to May 2027 above Regulus; open rings for the eastward stretches, red dots only for the westward run from 10 January to 1 April 2027. Right: Sun, Earth, and Mars from above at seven dates; the lines of sight swing back while Earth passes.

Earth is on the inside track and faster, 29.8 km/s against 24.1 km/s, so it laps Mars once per synodic period,

\[T_\text{syn} = \frac{1}{1/T_E - 1/T_M},\]

with \(T_E\) and \(T_M\) the two years. While Earth passes, the line of sight swings backward. For circles in one plane the rate at opposition is the velocity of Mars relative to Earth across the line of sight, over their distance: \(\dot\lambda = (a_M\omega_M - a_E\omega_E)/(a_M - a_E)\), with \(a\) the radii and \(\omega\) the angular speeds. The stations, the days Mars stands still, are the roots of the daily motion, the central difference of λ over half a day either side: the cell finds the days its sign flips and hands each bracket to brentq, which narrows it down to the root (root finding with scipy.optimize).

T_E, T_M = (36525 * 360 / ELEMENTS[b]["L"][1] for b in ["earth", "mars"])
T_syn = 1 / (1 / T_E - 1 / T_M)
a_E, a_M = ELEMENTS["earth"]["a"][0], ELEMENTS["mars"]["a"][0]
AU_KM = 1.495978707e8
print(f"years: Earth {T_E:.3f} d, Mars {T_M:.3f} d; synodic period {T_syn:.2f} d")
print(f"circular speeds: Earth {2 * np.pi * a_E * AU_KM / (T_E * 86400):.2f} km/s, "
      f"Mars {2 * np.pi * a_M * AU_KM / (T_M * 86400):.2f} km/s")
opp = oppositions(day("2025-01-01"), day("2036-01-01"))
print("days between oppositions 2025 to 2035:", ", ".join(f"{b - a:.0f}" for a, b in zip(opp, opp[1:])))

# circles in one plane, opposition at t = 0
w_E, w_M = 360 / T_E, 360 / T_M                  # degrees per day
lam_circ = lambda t: np.degrees(np.angle(a_M * np.exp(1j * np.radians(w_M * t)) - a_E * np.exp(1j * np.radians(w_E * t))))
rate_circ = lambda t: wrap(lam_circ(t + 0.5) - lam_circ(t - 0.5))
s1, s2 = brentq(rate_circ, -60, 0), brentq(rate_circ, 0, 60)
print(f"circles: rate at opposition {(a_M * w_M - a_E * w_E) / (a_M - a_E):+.3f}°/day "
      f"(difference quotient {rate_circ(0):+.3f}); stations at {s1:+.1f} and {s2:+.1f} d, "
      f"{s2 - s1:.1f} d and {wrap(lam_circ(s1) - lam_circ(s2)):.1f}° backward")

l1, l2 = geocentric(np.array(st_2027))[0]
print(f"2027: stations {date(st_2027[0])} and {date(st_2027[1])}, {st_2027[1] - st_2027[0]:.1f} d, "
      f"λ {l1:.2f}° to {l2:.2f}°, {l1 - l2:.2f}° backward; opposition {date(t_opp_2027)}, "
      f"rate there {daily_motion(t_opp_2027):+.3f}°/day")

years: Earth 365.256 d, Mars 686.980 d; synodic period 779.94 d
circular speeds: Earth 29.78 km/s, Mars 24.13 km/s
days between oppositions 2025 to 2035: 765, 765, 770, 786, 810
circles: rate at opposition -0.357°/day (difference quotient -0.357); stations at -36.4 and +36.4 d, 72.7 d and 15.9° backward
2027: stations 2027-01-10 and 2027-04-01, 81.0 d, λ 160.05° to 140.54°, 19.50° backward; opposition 2027-02-19, rate there -0.399°/day

The synodic period is 780 days. The circles give −0.357°/day and a backward run of 73 days over 15.9°; the real orbits give −0.399°/day and the 2027 run 19.5° west in 81 days. Circles get the mechanism right and the numbers within a fifth; the rest is eccentricity and tilt.

Play: watch Mars for two years

Mars stops and moves backward

The second half of this tutorial needs a free account