The Mathematics of a Striking Moonrise

A first-principles derivation of solar and lunar position, rise/set times, and the geometry that produces a perfect moonrise — with interactive visualizers.
mathematics
astronomy
celestial-mechanics
Published

May 5, 2026

雪と雪今宵師走の名月か

Yuki to yuki koyoi shiwasu no meigetsu ka

It snows and snows.

Tonight’s the full moon of December,

One wonders

-Matsuo Basho, Winter 1684

Every full moon is a gift which commands a certain immediate attention and yet even amongst full moons, individual differences make themselves apparent. The full moon always marks the start of night, but does so unequally throughout the year.

There is a particular kind of moonrise that commands attention: an enormous, brilliant disk cresting the horizon just as twilight fades into true night. It looks perfectly full, the sky is fully dark, contrast is at its peak. To get this, the moon must rise roughly one hour after sunset.

But on any given calendar date, when does the moon actually rise relative to the sun, and what phase is it then? The answers are seasonal, not constant — and the reason is geometric. This post derives the answer from first principles, building Sun and Moon position, sidereal time, and the rise–set equation from scratch, ending in four interactive diagrams you can drive yourself.

This post was inspired by the spring full moon on May 1st 2026, the Flower moon. My hope is the reader leaves with an enhanced appreciation of the difference between the waxing and waning gibbous.

The default location for this post is Washington, D.C. (\(\varphi = 38.9°\,\text{N}\), \(L = -77.04°\)) but every visualizer accepts your own coordinates.

1. Celestial mechanics from first principles

Reference frames and the rise equation

Astronomy is bookkeeping in three frames, each more local than the last:

  1. Ecliptic frame — defined by Earth’s orbital plane around the Sun. Coordinates: ecliptic longitude \(\lambda\) (measured from the vernal equinox, eastward) and ecliptic latitude \(\beta\). The Sun, by definition, has \(\beta_\odot = 0\).

  2. Equatorial frame — defined by Earth’s rotation axis, tilted from the ecliptic pole by the obliquity \(\varepsilon \approx 23.44°\). Coordinates: right ascension \(\alpha\) and declination \(\delta\). Stars are nearly fixed in this frame.

  3. Horizontal (topocentric) frame — defined by the observer’s local vertical. Coordinates: azimuth \(A\) (from north, clockwise) and altitude \(a\) (above horizon). A body’s altitude depends on the observer’s latitude \(\varphi\), the body’s \(\delta\), and the local hour angle \(H = \theta - \alpha\), where \(\theta\) is local sidereal time.

A body rises when its altitude \(a\) crosses a small negative threshold \(h_0\) from below; the corresponding hour angle satisfies

\[ \cos H_0 = \frac{\sin h_0 - \sin\varphi\sin\delta}{\cos\varphi\cos\delta} \tag{1.3} \]

This is the central equation of the post. Everything else is figuring out \(\delta\) and \(\alpha\) for the Sun and the Moon at the time we care about.

The transformations between the three frames are rotations of unit vectors:

\[ \begin{pmatrix}\cos\delta\cos\alpha\\ \cos\delta\sin\alpha\\ \sin\delta\end{pmatrix} = R_x(-\varepsilon) \begin{pmatrix}\cos\beta\cos\lambda\\ \cos\beta\sin\lambda\\ \sin\beta\end{pmatrix} \]

\[ \begin{pmatrix}\cos a \cos A\\ -\cos a \sin A\\ \sin a\end{pmatrix} = R_y\!\left(\tfrac{\pi}{2}-\varphi\right) \begin{pmatrix}\cos\delta\cos H\\ \cos\delta\sin H\\ \sin\delta\end{pmatrix} \]

Expanding the \(z\)-component of the second product gives the altitude:

\[ \sin a = \sin\varphi\sin\delta + \cos\varphi\cos\delta\cos H \tag{1.1} \]

and the \(x,y\) components together give azimuth:

\[ \tan A = \frac{-\sin H}{\tan\delta\cos\varphi - \sin\varphi\cos H} \tag{1.2} \]

Setting \(a = h_0\) in (1.1) and solving for \(\cos H\) produces (1.3) above. Two solutions \(\pm H_0\) correspond to rise (eastern, \(-H_0\)) and set (western, \(+H_0\)).

The Sun’s apparent motion — Kepler’s equation

Earth’s orbit is a slight ellipse with eccentricity \(e \approx 0.0167\). Earth moves faster at perihelion (early January) and slower at aphelion (early July), so the Sun’s apparent ecliptic longitude does not advance uniformly in time. The equation of center \(C\) is the correction from uniform “mean longitude” \(L_0\) to the actual ecliptic longitude:

\[ \lambda_\odot = L_0 + C \approx L_0 + 2e\sin M + \tfrac{5}{4}e^2\sin 2M + \cdots \]

with \(M\) the mean anomaly. Adding aberration and the dominant nutation term gives a model accurate to \(\sim 0.01°\) from 1900–2100 — well inside our \(\pm 2\) min rise-time target.

Kepler’s equation links the uniformly-advancing mean anomaly \(M\) to the geometric eccentric anomaly \(E\):

\[ M = E - e \sin E \tag{2.1} \]

The true anomaly \(\nu\) — the actual angle swept from perihelion — relates to \(E\) via

\[ \tan\tfrac{\nu}{2} = \sqrt{\tfrac{1+e}{1-e}}\,\tan\tfrac{E}{2}. \tag{2.2} \]

Series-expanding for small \(e\) yields

\[ C \equiv \nu - M = 2e\sin M + \tfrac{5}{4}e^2\sin 2M + O(e^3). \tag{2.3} \]

Following Meeus Astronomical Algorithms ch. 25, with \(T\) in Julian centuries from J2000.0 (JD 2451545.0):

\[ \begin{aligned} L_0 &= 280.46646° + 36000.76983°\,T + 0.0003032°\,T^2 \\ M &= 357.52911° + 35999.05029°\,T - 0.0001537°\,T^2 \\ e &= 0.016708634 - 0.000042037\,T \\ C &= (1.914602° - 0.004817°\,T)\sin M + (0.019993° - 0.000101°\,T)\sin 2M + 0.000289°\sin 3M \\ \lambda_\odot &= L_0 + C - 0.00569° - 0.00478°\sin\Omega \end{aligned} \]

The last two corrections are aberration and the dominant nutation term; \(\Omega = 125.04° - 1934.136°\,T\) is the longitude of the Moon’s ascending node. With \(\beta_\odot \equiv 0\), the equatorial conversion is

\[ \alpha_\odot = \arctan\!\big(\cos\varepsilon\sin\lambda_\odot,\;\cos\lambda_\odot\big), \qquad \sin\delta_\odot = \sin\varepsilon\sin\lambda_\odot. \tag{2.4} \]

Why winter inverts the moonrise lag

A perfectly full moon sits at \(\lambda_M = \lambda_\odot + 180°\). In December, the Sun is at \(\lambda_\odot \approx 270°\) and dips to \(\delta \approx -23.4°\) — far south, short days. The full moon is then at \(\lambda_M \approx 90°\) with \(\delta_M \approx +23.4°\) — far north, long time above the horizon. From mid-northern latitudes a high-declination body has \(H_0\) close to \(\pi\) (it spends almost the whole sidereal day up), so it rises early relative to the Sun’s setting. The “winter paradox” in the table below is just (1.3) being applied honestly.

This can be visualized by imaging a sphere, say a naval orange, and dividing it in half with a mark, one half will be the hemisphere that is visible to a person on the ground. Then place a ribbon around the orange along this mark and spin it so the center of the ribbon is no longer following the center of the mark. From here, the edges of the ribbon represent the angle of the sun during the solstices, and the center of the ribbon will represent the equinox. One edge of the ribbon will intersect the hemisphere more, that is the summer solstice, the one that is less is the winter solstice. The movement of the moon will be the opposite edge of the ribbon relative to whatever edge is representing the sun. This is why the lunar day is long when the solar day is short and vice versa.

For an excellent visual aid, consult this link and imagine the yellow band as the ribbon on the orange in the above analogy.

https://andrewmarsh.com/apps/staging/sunpath3d.html

The Moon — five fundamental arguments

The Moon’s motion is harder because the Sun perturbs its orbit substantially. Every modern lunar theory (Brown, ELP-2000, etc.) parametrizes by five Delaunay-style fundamental arguments, all linear in \(T\):

Symbol Name Linear-in-\(T\) expression (degrees)
\(L'\) Moon’s mean longitude \(218.31644 + 481267.88123\,T - 0.00158\,T^2\)
\(D\) Mean elongation from Sun \(297.85019 + 445267.11140\,T - 0.00188\,T^2\)
\(M\) Sun’s mean anomaly \(357.52911 + 35999.05029\,T\)
\(M'\) Moon’s mean anomaly \(134.96340 + 477198.86751\,T + 0.00874\,T^2\)
\(F\) Argument of latitude (from ascending node) \(93.27210 + 483202.01752\,T - 0.00365\,T^2\)

Geometrically: \(L'\) tracks the Moon along the ecliptic; \(D\) measures the synodic phase (\(D = 0\) is new moon, \(D = 180°\) is full); \(M'\) encodes orbital eccentricity; \(F\) measures distance from the ascending node and so drives ecliptic latitude.

The Moon’s true ecliptic longitude and latitude are trigonometric series in these arguments. Truncated to the dozen-or-so leading terms, the model is accurate to \(\sim 0.1°\), translating to under \(\pm 2\) min of rise-time error at \(\varphi \approx 39°\).

The full ELP-2000 series has thousands of terms; we keep Meeus Table 47.A truncated. For longitude (degrees):

\[ \begin{aligned} \lambda_M - L' \approx\;& 6.289\sin M' + 1.274\sin(2D - M') + 0.658\sin 2D + 0.214\sin 2M' \\ & -\,0.186\sin M - 0.114\sin 2F + 0.059\sin(2D - 2M') + 0.057\sin(2D - M - M') \\ & +\,0.053\sin(2D + M') + 0.046\sin(2D - M) - 0.041\sin(M - M') \\ & -\,0.035\sin D - 0.030\sin(M + M') + \cdots \end{aligned} \]

For latitude, the dominant inclination term is \(5.128°\) — the Moon’s orbital inclination \(i\). Higher-order terms come from the perturbed motion of the node:

\[ \beta_M \approx 5.128\sin F + 0.281\sin(M' + F) + 0.278\sin(M' - F) + 0.173\sin(2D - F) + \cdots \]

The Earth–Moon distance is dominated by the elliptical term \(-20905\cos M'\) km, oscillating between perigee (\(\sim 363{,}000\) km) and apogee (\(\sim 405{,}000\) km). Conversion to equatorial uses the full transform (with \(\beta \neq 0\)):

\[ \sin\delta_M = \sin\beta_M\cos\varepsilon + \cos\beta_M\sin\varepsilon\sin\lambda_M, \tag{3.1} \] \[ \alpha_M = \arctan\!\big(\sin\lambda_M\cos\varepsilon - \tan\beta_M\sin\varepsilon,\; \cos\lambda_M\big). \tag{3.2} \]

Rise time — refraction, parallax, and sidereal iteration

A body’s geometric center is not at altitude \(0\) when its limb touches the horizon. Atmospheric refraction lifts the image by \(\sim 34'\); the Sun’s and Moon’s apparent semidiameters add another \(\sim 16'\); lunar horizontal parallax of \(\sim 57'\) shifts the Moon further. Net effect:

\[ h_0^{\odot} = -50' \qquad h_0^{M} = 0.7275\,\pi_M - 34' \approx +7'. \tag{4.1} \]

These are not optional — together they shift moonrise by \(5{-}10\) min, comparable to the seasonal effect we are studying. The hour angle is then \(H = \theta - \alpha\), with \(\theta\) the local sidereal time. Because the Moon’s RA drifts \(\sim 0.5°/\text{hr}\), “rise time” is a fixed point and must be iterated.

Refraction and parallax. Three corrections matter:

  • Atmospheric refraction at the horizon raises the apparent position by \(\rho \approx 34'\).
  • The Sun’s semidiameter \(s_\odot \approx 16'\). We see the Sun “rise” when its upper limb appears, so the geometric center must be \(-s_\odot\) below.
  • Lunar horizontal parallax \(\pi_M = \arcsin(R_\oplus / \Delta) \approx 57'\) raises apparent altitude. Its semidiameter is \(\sim 15'-17'\). The factor \(0.7275\) in (4.1) comes from the difference between geometric horizontal parallax and the contribution to the center’s altitude (Meeus 15.1).

Greenwich Mean Sidereal Time (Meeus 12.4), in degrees:

\[ \theta_G = 280.46062 + 360.98564737\,(JD - 2451545) + 0.000388\,T^2 - T^3/38710000 \tag{4.2} \]

Local sidereal time \(\theta = \theta_G + L\) for east-positive longitude \(L\).

Iteration. Starting from noon UT of the desired calendar date,

  1. Compute \(\alpha(t), \delta(t)\) from the lunar series above.
  2. Compute \(H_0\) from (1.3) using (4.1) for \(h_0\).
  3. Compute current \(H = \theta(t) - \alpha(t)\), normalized to \((-\pi, \pi]\).
  4. Update \(t \leftarrow t + (\Delta H)/(2\pi/0.9972696\,\text{days})\) where \(\Delta H = -H_0 - H\), normalized.

Two to four iterations converge to under a second. Setting and transit use \(+H_0\) and \(0\).

The harvest moon effect

The daily lag in moonrise — how much later the moon rises tonight than yesterday — is set by how fast the moon’s right ascension grows (about \(13°/\text{day}\)) and by how that RA increase projects onto the eastern horizon at the observer’s latitude.

If the Moon’s path on the sky met the eastern horizon perpendicularly, a \(13°\) RA shift would translate directly into \(13°/15° = 52\,\text{min}\) of delay. But near the autumnal equinox at northern mid-latitudes, the ecliptic crosses the eastern horizon at a shallow angle in the early evening: the Moon’s daily march along the ecliptic projects mostly along the horizon, not down into it. The lag shrinks to as little as \(20-25\) min for several consecutive nights — the harvest moon effect.

If \(\psi\) is the angle the ecliptic makes with the local horizon at the moment of moonrise, the daily lag is approximately

\[ \Delta t_\text{daily} \approx \frac{13°/15°}{1 + (\dot\alpha/\dot\lambda - 1)} \cdot \sin\psi \approx \frac{52\,\text{min}}{1} \cdot \sin\psi. \tag{5.1} \]

A shallow \(\psi \approx 25°\) gives \(\sim 22\) min, while a steep \(\psi \approx 70°\) (spring evenings) gives \(\sim 49\) min. The orbital inclination of \(5.14°\) adds a secondary modulation depending on where the ascending node sits — when the node is near the vernal equinox (“major standstill”), the Moon’s declination range is the full \(\pm 28.6°\) rather than \(\pm 18.3°\), exaggerating both effects.

2. The compute engine

The next four visualizers all share one OJS module implementing the Sun position, Moon position, and rise/set iteration from §1. Code is folded; expand to inspect or modify.

Code 
DEG = Math.PI / 180
RAD = 180 / Math.PI
J2000 = 2451545.0

julianCenturies = (jd) => (jd - J2000) / 36525

// Calendar date → JD (Meeus 7.1). Month 1–12, day can be fractional.
jdFromDate = function(year, month, day) {
  let y = year, m = month;
  if (m <= 2) { y -= 1; m += 12; }
  const a = Math.floor(y / 100);
  const b = 2 - a + Math.floor(a / 4);
  return Math.floor(365.25*(y + 4716)) + Math.floor(30.6001*(m + 1)) + day + b - 1524.5;
}

obliquity = function(T) {
  // Mean obliquity, Meeus 22.2 (low-precision)
  return (23.4392911 - 0.0130042*T - 1.64e-7*T*T + 5.04e-7*T*T*T) * DEG;
}

normRad = function(x) {
  let y = x % (2*Math.PI);
  if (y < 0) y += 2*Math.PI;
  return y;
}

sunPosition = function(jd) {
  const T = julianCenturies(jd);
  const L0 = (280.46646 + 36000.76983*T + 0.0003032*T*T) * DEG;
  const M  = (357.52911 + 35999.05029*T - 0.0001537*T*T) * DEG;
  const C  = ((1.914602 - 0.004817*T - 0.000014*T*T)*Math.sin(M)
            + (0.019993 - 0.000101*T)*Math.sin(2*M)
            + 0.000289*Math.sin(3*M)) * DEG;
  const trueLong = L0 + C;
  const omega = (125.04 - 1934.136*T) * DEG;
  const lambda = trueLong - (0.00569 + 0.00478*Math.sin(omega)) * DEG;
  const eps = obliquity(T);
  const alpha = Math.atan2(Math.cos(eps)*Math.sin(lambda), Math.cos(lambda));
  const delta = Math.asin(Math.sin(eps)*Math.sin(lambda));
  return { lambda: normRad(lambda), beta: 0, alpha: normRad(alpha), delta, eps };
}

moonPosition = function(jd) {
  const T = julianCenturies(jd);
  const Lp = ((218.3164477 + 481267.88123421*T - 0.0015786*T*T) * DEG);
  const D  = ((297.8501921 + 445267.1114034*T - 0.0018819*T*T) * DEG);
  const M  = ((357.5291092 + 35999.0502909*T) * DEG);
  const Mp = ((134.9633964 + 477198.8675055*T + 0.0087414*T*T) * DEG);
  const F  = (( 93.2720950 + 483202.0175233*T - 0.0036539*T*T) * DEG);
  const dLambda = (
      6.288774 * Math.sin(Mp)
    + 1.274027 * Math.sin(2*D - Mp)
    + 0.658314 * Math.sin(2*D)
    + 0.213618 * Math.sin(2*Mp)
    - 0.185116 * Math.sin(M)
    - 0.114332 * Math.sin(2*F)
    + 0.058793 * Math.sin(2*D - 2*Mp)
    + 0.057066 * Math.sin(2*D - M - Mp)
    + 0.053322 * Math.sin(2*D + Mp)
    + 0.045758 * Math.sin(2*D - M)
    - 0.040923 * Math.sin(M - Mp)
    - 0.034720 * Math.sin(D)
    - 0.030383 * Math.sin(M + Mp)
    + 0.015327 * Math.sin(2*D - 2*F)
    - 0.012528 * Math.sin(Mp + 2*F)
    + 0.010980 * Math.sin(Mp - 2*F)
  ) * DEG;
  const dBeta = (
      5.128122 * Math.sin(F)
    + 0.280602 * Math.sin(Mp + F)
    + 0.277693 * Math.sin(Mp - F)
    + 0.173237 * Math.sin(2*D - F)
    + 0.055413 * Math.sin(2*D + F - Mp)
    + 0.046271 * Math.sin(2*D - F - Mp)
    + 0.032573 * Math.sin(2*D + F)
    + 0.017198 * Math.sin(2*Mp + F)
    + 0.009266 * Math.sin(2*D + Mp - F)
    + 0.008822 * Math.sin(2*Mp - F)
  ) * DEG;
  const dist = 385000.56
    - 20905.355 * Math.cos(Mp)
    -  3699.111 * Math.cos(2*D - Mp)
    -  2955.968 * Math.cos(2*D)
    -   569.925 * Math.cos(2*Mp)
    +    48.888 * Math.cos(M)
    -   246.158 * Math.cos(2*Mp - 2*D);
  const lambda = Lp + dLambda;
  const beta = dBeta;
  const eps = obliquity(T);
  const sinD = Math.sin(beta)*Math.cos(eps) + Math.cos(beta)*Math.sin(eps)*Math.sin(lambda);
  const delta = Math.asin(sinD);
  const alpha = Math.atan2(
    Math.sin(lambda)*Math.cos(eps) - Math.tan(beta)*Math.sin(eps),
    Math.cos(lambda)
  );
  return { lambda: normRad(lambda), beta, alpha: normRad(alpha), delta, distance: dist };
}

gmst = function(jd) {
  const T = julianCenturies(jd);
  const d = jd - J2000;
  let theta = 280.46061837 + 360.98564736629*d + 0.000387933*T*T - T*T*T/38710000;
  theta = ((theta % 360) + 360) % 360;
  return theta * DEG;
}

lstFor = (jd, lonEastDeg) => normRad(gmst(jd) + lonEastDeg*DEG)

altAz = function(latRad, dec, alpha, jd, lonEastDeg) {
  const H = lstFor(jd, lonEastDeg) - alpha;
  const sinAlt = Math.sin(latRad)*Math.sin(dec) + Math.cos(latRad)*Math.cos(dec)*Math.cos(H);
  const alt = Math.asin(Math.max(-1, Math.min(1, sinAlt)));
  const y = -Math.sin(H);
  const x = Math.tan(dec)*Math.cos(latRad) - Math.sin(latRad)*Math.cos(H);
  let az = Math.atan2(y, x);
  if (az < 0) az += 2*Math.PI;
  return { alt, az, H };
}

sunH0 = () => -0.8333 * DEG
moonH0 = (m) => {
  const parallax = Math.asin(6378.14 / m.distance);
  return 0.7275*parallax - 0.5667*DEG;
}

// Find rise / transit / set on the calendar date whose 0h UT is jd0.
// Sign: -1 for rise, 0 for transit, +1 for set.
findEvent = function(bodyFn, h0Fn, jd0, latRad, lonEastDeg, sign) {
  let t = jd0 + 0.5; // start at noon UT
  for (let i = 0; i < 8; i++) {
    const body = bodyFn(t);
    const h0 = h0Fn(body);
    const cosH0 = (Math.sin(h0) - Math.sin(latRad)*Math.sin(body.delta))
                / (Math.cos(latRad)*Math.cos(body.delta));
    if (cosH0 >  1 + 1e-9) return { jd: null, reason: 'never-rises' };
    if (cosH0 < -1 - 1e-9) return { jd: null, reason: 'never-sets' };
    const H0 = Math.acos(Math.max(-1, Math.min(1, cosH0)));
    const desiredH = sign === 0 ? 0 : (sign < 0 ? -H0 : +H0);
    let currentH = lstFor(t, lonEastDeg) - body.alpha;
    while (currentH >  Math.PI) currentH -= 2*Math.PI;
    while (currentH < -Math.PI) currentH += 2*Math.PI;
    let dH = desiredH - currentH;
    while (dH >  Math.PI) dH -= 2*Math.PI;
    while (dH < -Math.PI) dH += 2*Math.PI;
    const dt = dH / (2*Math.PI / 0.9972695663); // sidereal day
    t += dt;
    if (Math.abs(dt) < 1e-5) break;
  }
  // Confirm event lies within [jd0, jd0+1]
  if (t < jd0 - 0.05 || t > jd0 + 1.05) {
    // try alternative branch
  }
  return { jd: t, reason: 'ok' };
}

riseSet = function(bodyFn, h0Fn, jd0, latRad, lonEastDeg) {
  return {
    rise:    findEvent(bodyFn, h0Fn, jd0, latRad, lonEastDeg, -1),
    transit: findEvent(bodyFn, h0Fn, jd0, latRad, lonEastDeg,  0),
    set:     findEvent(bodyFn, h0Fn, jd0, latRad, lonEastDeg, +1)
  };
}

// Phase angle and illumination, Meeus 48
illumination = function(jd) {
  const s = sunPosition(jd);
  const m = moonPosition(jd);
  const psi = Math.acos(
    Math.sin(s.delta)*Math.sin(m.delta) +
    Math.cos(s.delta)*Math.cos(m.delta)*Math.cos(s.alpha - m.alpha)
  );
  const R = 1.496e8;
  const Δ = m.distance;
  const i = Math.atan2(R*Math.sin(psi), Δ - R*Math.cos(psi));
  // Synodic phase angle (0 = new, π = full), signed by D
  const D = normRad(m.lambda - s.lambda);
  return { phaseAngle: i, illuminated: (1 + Math.cos(i))/2, elongation: D };
}

// Sub-solar / sub-lunar geographic point (lat, lon east) at given JD
subPoint = function(jd, body) {
  const b = body(jd);
  const gast = gmst(jd); // ignore equation of equinoxes (tiny)
  let lon = b.alpha - gast;
  while (lon >  Math.PI) lon -= 2*Math.PI;
  while (lon < -Math.PI) lon += 2*Math.PI;
  return { lat: b.delta, lon };
}

mod360 = x => ((x % 360) + 360) % 360

// Format a JD (UT) as local clock time HH:MM given a TZ offset in hours.
fmtLocal = function(jd, tzHours) {
  if (jd == null) return "—";
  const local = jd + tzHours/24;
  const frac = local - Math.floor(local) - 0.5;
  let hours24 = (frac * 24 + 24) % 24;
  const h = Math.floor(hours24);
  const m = Math.round((hours24 - h)*60);
  return `${String(h).padStart(2,'0')}:${String(m).padStart(2,'0')}`;
}

3. Orbital geometry & horizon timeline

Set your location, date, and time of day. Type values directly into the number boxes, drag sliders, or hit Animate day-of-year to sweep through the calendar. The first diagram shows the Earth–Moon–Sun geometry from above the ecliptic; the second is a 24-hour timeline of when the Sun and Moon are above the horizon at your latitude.

Code 
todayObj = new Date()
yearToday = todayObj.getUTCFullYear()
todayDateStr = todayObj.toISOString().slice(0,10)
maxDateStr = (() => {
  const d = new Date(todayObj.getTime() + 2*365.25*86400000);
  return d.toISOString().slice(0,10);
})()
Code 
viewof loc = Inputs.form({
  lat: Inputs.number([-89.9, 89.9], {value: 38.9,  step: 0.001, label: "Latitude (°N)"}),
  lon: Inputs.number([-180, 180],  {value: -77.04, step: 0.001, label: "Longitude (°E, west negative)"}),
  tz:  Inputs.number([-12, 14],    {value: -5,    step: 1,     label: "Timezone offset (hours from UTC)"})
})
viewof picker = Inputs.form({
  date: Inputs.date({value: todayDateStr, min: todayDateStr, max: maxDateStr, label: "Date"}),
  hour: Inputs.range([0, 24], {value: 19, step: 0.05, label: "Local clock hour"})
})
viewof animState = Inputs.form({
  animate: Inputs.toggle({value: false, label: "▶ Animate day-of-year"}),
  speed:   Inputs.range([1, 60], {value: 10, step: 1, label: "Animation speed (days/sec)"})
})
Code 
pickerDateObj = (typeof picker.date === 'string')
  ? new Date(picker.date + "T00:00:00Z")
  : picker.date
pickerYear = pickerDateObj.getUTCFullYear()
pickerDoy = {
  const start = Date.UTC(pickerYear, 0, 1);
  return Math.floor((pickerDateObj.getTime() - start)/86400000) + 1;
}

// Reactive day-of-year: while animating, advance with a generator; otherwise pass through.
doy = animState.animate
  ? Generators.observe(notify => {
      let d = pickerDoy;
      notify(d);
      const id = setInterval(() => {
        const isLeap = (pickerYear % 4 === 0 && pickerYear % 100 !== 0) || pickerYear % 400 === 0;
        const maxDoy = isLeap ? 366 : 365;
        d = (d % maxDoy) + 1;
        notify(d);
      }, 1000 / Math.max(1, animState.speed));
      return () => clearInterval(id);
    })
  : pickerDoy

// Derived JD at 0h UT of (year, doy), and at chosen local clock time
jd0 = jdFromDate(pickerYear, 1, doy)
jdNow = jd0 + (picker.hour - loc.tz)/24
sun = sunPosition(jdNow)
moon = moonPosition(jdNow)
ill = illumination(jdNow)
events = ({
  sun:  riseSet(sunPosition,  sunH0,  jd0, loc.lat*DEG, loc.lon),
  moon: riseSet(moonPosition, moonH0, jd0, loc.lat*DEG, loc.lon)
})

// Convenience: compute the human-readable date for the current animated doy
currentDateLabel = {
  const ms = Date.UTC(pickerYear, 0, 1) + (doy - 1)*86400000;
  return new Date(ms).toISOString().slice(0,10);
}
Code 
// ===== Visualizer A: overhead geometric diagram =====
{
  const W = 960, H = 540, cx = W/2, cy = H/2;
  const Rsun = 230, Rmoon = 120;

  const svg = d3.create("svg")
    .attr("viewBox", [0, 0, W, H])
    .style("background", "#0a0a1a")
    .style("border-radius", "8px")
    .style("max-width", "100%")
    .style("height", "auto");

  // labels for ecliptic cardinal directions (vernal equinox = +x)
  const dir = [["♈ 0°", Rsun+22, 0],["♋ 90°", 0, Rsun+22],["♎ 180°", -(Rsun+30), 0],["♑ 270°", 0, -(Rsun+22)]];
  for (const [t,dx,dy] of dir) {
    svg.append("text").attr("x", cx+dx).attr("y", cy-dy).text(t)
      .attr("fill","#666").attr("font-size","10px").attr("text-anchor","middle");
  }

  // Sun's ecliptic position
  const sx = cx + Rsun*Math.cos(sun.lambda);
  const sy = cy - Rsun*Math.sin(sun.lambda); // y inverted on screen
  // Moon ecliptic + small radial bump for β
  const r_moon_eff = Rmoon + 50*Math.sin(moon.beta); // visualize β
  const mx = cx + r_moon_eff*Math.cos(moon.lambda);
  const my = cy - r_moon_eff*Math.sin(moon.lambda);

  // d3.arc angle to a screen vector (dx, dy): 0 = up, increases clockwise.
  // For screen vector (dx, dy): arcAngle = atan2(dx, -dy).
  function screenToArcAngle(dx, dy) { return Math.atan2(dx, -dy); }

  // Earth (with a real day/night terminator: lit hemisphere faces the sun)
  const earth = svg.append("g").attr("transform", `translate(${cx},${cy})`);
  earth.append("circle").attr("r", 20).attr("fill", "#234").attr("stroke","#4b9cd3");
  const earthSunArc = screenToArcAngle(sx - cx, sy - cy);
  const litEarth = d3.arc()({innerRadius: 0, outerRadius: 20,
    startAngle: earthSunArc - Math.PI/2,
    endAngle:   earthSunArc + Math.PI/2});
  earth.append("path").attr("d", litEarth).attr("fill", "#4b9cd3").attr("opacity", 0.9);
  // Observer pin: place on Earth's surface at the angle pointing to the observer's instantaneous
  // local meridian — i.e. towards screen direction matching local sidereal time.
  const obsArcAngle = screenToArcAngle(
    Math.cos(lstFor(jdNow, loc.lon)),
    -Math.sin(lstFor(jdNow, loc.lon))
  );
  earth.append("circle")
    .attr("cx", 22*Math.sin(obsArcAngle))
    .attr("cy", -22*Math.cos(obsArcAngle))
    .attr("r", 2.5).attr("fill","#ff6").attr("stroke","#000").attr("stroke-width",0.5);

  // Sun
  svg.append("circle").attr("cx", sx).attr("cy", sy).attr("r", 13).attr("fill","#ffcc33")
    .attr("filter","drop-shadow(0 0 10px #ffcc33)");

  // Moon orbit ring
  svg.append("circle").attr("cx", cx).attr("cy", cy).attr("r", Rmoon)
    .attr("fill","none").attr("stroke","rgba(255,255,255,0.12)").attr("stroke-dasharray","3,4");

  // Moon
  const moonG = svg.append("g").attr("transform", `translate(${mx},${my})`);
  moonG.append("circle").attr("r", 8).attr("fill", "#ddd");
  // Lit hemisphere faces the sun; dark hemisphere is opposite.
  const moonSunArc = screenToArcAngle(sx - mx, sy - my);
  const darkMoon = d3.arc()({innerRadius:0, outerRadius:8,
    startAngle: moonSunArc + Math.PI/2,
    endAngle:   moonSunArc + 3*Math.PI/2});
  moonG.append("path").attr("d", darkMoon).attr("fill","rgba(0,0,0,0.82)");

  // Annotations
  svg.append("text").attr("x", 12).attr("y", 22).attr("fill","#ccc").attr("font-size","12px")
    .text(`Sun λ = ${(sun.lambda*RAD).toFixed(1)}°, δ = ${(sun.delta*RAD).toFixed(2)}°`);
  svg.append("text").attr("x", 12).attr("y", 40).attr("fill","#ccc").attr("font-size","12px")
    .text(`Moon λ = ${(moon.lambda*RAD).toFixed(1)}°, β = ${(moon.beta*RAD).toFixed(2)}°, Δ = ${moon.distance.toFixed(0)} km`);
  svg.append("text").attr("x", 12).attr("y", 58).attr("fill","#ccc").attr("font-size","12px")
    .text(`Phase angle = ${(ill.phaseAngle*RAD).toFixed(1)}°  →  illuminated ${(ill.illuminated*100).toFixed(1)}%`);

  return svg.node();
}
Code 
// ===== Visualizer B: 24-h horizon timeline =====
{
  const W = 960, H = 200, m = {top: 36, right: 24, bottom: 36, left: 36};
  const x = d3.scaleLinear().domain([0, 24]).range([m.left, W - m.right]);
  const svg = d3.create("svg").attr("viewBox", [0,0,W,H])
    .style("max-width","100%").style("height","auto");

  const barTop = 80, barH = 36;

  // Background bar (night)
  svg.append("rect").attr("x", m.left).attr("y", barTop).attr("width", W-m.left-m.right).attr("height", barH)
    .attr("fill", "#10103a");

  function localHours(ev) {
    if (!ev.jd) return null;
    const local = ev.jd + loc.tz/24;
    return ((local - Math.floor(local) - 0.5) * 24 + 24) % 24;
  }
  const sr = localHours(events.sun.rise);
  const ss = localHours(events.sun.set);
  const mr = localHours(events.moon.rise);
  const ms = localHours(events.moon.set);

  // Daylight band
  if (sr != null && ss != null) {
    if (sr < ss) {
      svg.append("rect").attr("x", x(sr)).attr("y", barTop).attr("width", x(ss)-x(sr))
        .attr("height", barH).attr("fill","#ffd866");
    } else {
      svg.append("rect").attr("x", x(0)).attr("y", barTop).attr("width", x(ss)-x(0))
        .attr("height", barH).attr("fill","#ffd866");
      svg.append("rect").attr("x", x(sr)).attr("y", barTop).attr("width", x(24)-x(sr))
        .attr("height", barH).attr("fill","#ffd866");
    }
  }
  // Moon band
  function drawMoonBand(a, b) {
    svg.append("rect").attr("x", x(a)).attr("y", barTop+8).attr("width", x(b)-x(a))
      .attr("height", barH-16).attr("fill", "rgba(220,220,220,0.75)").attr("rx",3);
  }
  if (mr != null && ms != null) {
    if (mr < ms) drawMoonBand(mr, ms);
    else { drawMoonBand(mr, 24); drawMoonBand(0, ms); }
  }

  // Tick at sunset+1h
  if (ss != null) {
    const target = (ss + 1) % 24;
    svg.append("line").attr("x1", x(target)).attr("x2", x(target))
      .attr("y1", barTop-6).attr("y2", barTop+barH+6).attr("stroke","#e44").attr("stroke-dasharray","3,3");
    svg.append("text").attr("x", x(target)).attr("y", barTop-10).attr("fill","#e44")
      .attr("font-size","11px").attr("text-anchor","middle").text("sunset + 1 h");
  }

  // Axis
  svg.append("g").attr("transform",`translate(0,${barTop+barH})`)
    .call(d3.axisBottom(x).ticks(12).tickFormat(d => d+"h"))
    .attr("color","#666").attr("font-size","11px");

  // Title
  svg.append("text").attr("x", W/2).attr("y", 24).attr("text-anchor","middle")
    .attr("font-size","16px").attr("fill","#333").attr("font-weight","bold")
    .text(`Sun & Moon above horizon — ${currentDateLabel}`);

  // Inline readout
  const lag = (mr != null && ss != null)
    ? (((mr - ss) + 24) % 24) * 60
    : null;
  svg.append("text").attr("x", W/2).attr("y", H-6).attr("text-anchor","middle")
    .attr("font-size","11px").attr("fill","#444")
    .text(`Sunset ${ss==null?'—':fmtLocal(events.sun.set.jd, loc.tz)}, `
        + `Moonrise ${mr==null?'—':fmtLocal(events.moon.rise.jd, loc.tz)}, `
        + `Δ = ${lag==null?'—':lag.toFixed(0)+' min'}, illum ${(ill.illuminated*100).toFixed(0)}%`);

  return svg.node();
}

4. Two views of the same sky

The next two diagrams show the same instant from two perspectives. Drag the time-of-day slider above; both diagrams update.

The stereographic plot (left) is what you would see looking up: zenith at center, horizon as the bounding circle, cardinal directions on the rim. The Sun and Moon are drawn at their true altitude and azimuth; their shaded sides face away from each other. Colored faint trails show each body’s path over the past 24 hours.

The wireframe globe (right) is the same instant viewed from outside Earth. The latitude/longitude wireframe rotates with the planet. The day/night terminator is drawn as the great circle 90° from the sub-solar point. The sub-solar point ☉, sub-lunar point ☾, and your location pin are marked. The inset reports whether each body is currently above your horizon.

Code 
{
  const W = 540, H = 540, cx = W/2, cy = H/2, R = 240;
  const svg = d3.create("svg").attr("viewBox", [0,0,W,H])
    .style("background","radial-gradient(circle at center, #001 0%, #000 70%)")
    .style("border-radius","8px").style("max-width","100%").style("height","auto");

  // Stereographic from nadir: r = R * tan((90° - alt)/2). Returns null when the
  // point is below the horizon so callers can break their path there.
  function project(alt, az) {
    if (alt < -0.001) return null;
    const r = R * Math.tan((Math.PI/2 - alt)/2);
    // az measured from N clockwise; on screen, x = r sin(az), y = -r cos(az)
    return [cx + r*Math.sin(az), cy - r*Math.cos(az)];
  }

  // Path generator that BREAKS at null entries instead of connecting across them.
  // Also breaks when consecutive points are far apart on screen (azimuth wrap-around).
  const lineGen = d3.line()
    .defined(d => d != null)
    .x(d => d[0]).y(d => d[1]);

  // Horizon ring + altitude circles
  for (const a of [30, 60]) {
    const r = R * Math.tan((Math.PI/2 - a*DEG)/2);
    svg.append("circle").attr("cx",cx).attr("cy",cy).attr("r",r)
      .attr("fill","none").attr("stroke","#234").attr("stroke-dasharray","2,3");
  }
  svg.append("circle").attr("cx",cx).attr("cy",cy).attr("r",R)
    .attr("fill","none").attr("stroke","#445");

  // Cardinal labels
  const card = [["N",0],["E",Math.PI/2],["S",Math.PI],["W",3*Math.PI/2]];
  for (const [t,a] of card) {
    svg.append("text").attr("x", cx + (R+12)*Math.sin(a))
      .attr("y", cy - (R+12)*Math.cos(a) + 4)
      .attr("text-anchor","middle").attr("fill","#88a").attr("font-size","12px").text(t);
  }

  // Celestial equator + ecliptic + body trails, sampled.
  // We post-process: any pair of consecutive non-null points that are unrealistically
  // far apart on screen (a wrap-around chord) gets a null inserted between them so
  // .defined() breaks the path there.
  function drawGreatCircle(samples, color) {
    const broken = [];
    for (let i = 0; i < samples.length; i++) {
      const p = samples[i];
      const prev = broken[broken.length - 1];
      if (p && prev) {
        const dx = p[0] - prev[0], dy = p[1] - prev[1];
        if (dx*dx + dy*dy > (R*0.5)*(R*0.5)) broken.push(null);
      }
      broken.push(p);
    }
    svg.append("path").attr("d", lineGen(broken)).attr("fill","none")
      .attr("stroke",color).attr("stroke-width",0.6).attr("opacity",0.55);
  }
  // Equator: δ = 0, vary α
  const eqPts = [];
  for (let a=0; a<2*Math.PI; a+=Math.PI/120) {
    const H = lstFor(jdNow, loc.lon) - a;
    const sinAlt = Math.cos(loc.lat*DEG)*Math.cos(H);
    const alt = Math.asin(sinAlt);
    const az = Math.atan2(-Math.sin(H), -Math.sin(loc.lat*DEG)*Math.cos(H));
    eqPts.push(project(alt, (az+2*Math.PI)%(2*Math.PI)));
  }
  drawGreatCircle(eqPts, "#5a8");

  // Ecliptic
  const eps0 = obliquity(julianCenturies(jdNow));
  const ecPts = [];
  for (let l=0; l<2*Math.PI; l+=Math.PI/120) {
    const dec = Math.asin(Math.sin(eps0)*Math.sin(l));
    const ra = Math.atan2(Math.cos(eps0)*Math.sin(l), Math.cos(l));
    const H = lstFor(jdNow, loc.lon) - ra;
    const sinAlt = Math.sin(loc.lat*DEG)*Math.sin(dec)+Math.cos(loc.lat*DEG)*Math.cos(dec)*Math.cos(H);
    const alt = Math.asin(Math.max(-1, Math.min(1, sinAlt)));
    const y = -Math.sin(H);
    const x = Math.tan(dec)*Math.cos(loc.lat*DEG) - Math.sin(loc.lat*DEG)*Math.cos(H);
    let az = Math.atan2(y,x); if (az<0) az += 2*Math.PI;
    ecPts.push(project(alt, az));
  }
  drawGreatCircle(ecPts, "#a85");

  // 24-hour trails for Sun & Moon — keep nulls in-place so the path breaks
  // when the body dips below the horizon.
  function trail(bodyFn, color) {
    const pts = [];
    for (let dh=-24; dh<=0; dh+=0.25) {
      const t = jdNow + dh/24;
      const b = bodyFn(t);
      const aa = altAz(loc.lat*DEG, b.delta, b.alpha, t, loc.lon);
      pts.push(project(aa.alt, aa.az));
    }
    drawGreatCircle(pts, color);
  }
  trail(sunPosition, "#ffcc33");
  trail(moonPosition, "#cce");

  // Current Sun & Moon
  const sunAA = altAz(loc.lat*DEG, sun.delta, sun.alpha, jdNow, loc.lon);
  const moonAA = altAz(loc.lat*DEG, moon.delta, moon.alpha, jdNow, loc.lon);
  const sp = project(sunAA.alt, sunAA.az);
  const mp = project(moonAA.alt, moonAA.az);
  if (sp) {
    svg.append("circle").attr("cx",sp[0]).attr("cy",sp[1]).attr("r",7)
      .attr("fill","#ffcc33").attr("filter","drop-shadow(0 0 6px #ffcc33)");
  }
  if (mp) {
    const g = svg.append("g").attr("transform",`translate(${mp[0]},${mp[1]})`);
    g.append("circle").attr("r",6).attr("fill","#eee");
    // shade based on phase
    const k = ill.illuminated;
    g.append("path").attr("d", d3.arc()({innerRadius:0, outerRadius:6,
      startAngle: 0, endAngle: Math.PI*(1-k)*2}))
      .attr("fill","rgba(0,0,0,0.7)");
  }

  svg.append("text").attr("x", 8).attr("y", 16).attr("fill","#bbb").attr("font-size","11px")
    .text(`zenith view · ${currentDateLabel} · ${picker.hour.toFixed(2)}h local`);
  svg.append("text").attr("x", 8).attr("y", H-8).attr("fill","#5a8").attr("font-size","10px")
    .text("— equator");
  svg.append("text").attr("x", 90).attr("y", H-8).attr("fill","#a85").attr("font-size","10px")
    .text("— ecliptic");
  return svg.node();
}
{
  const W = 540, H = 540;
  const svg = d3.create("svg").attr("viewBox", [0,0,W,H])
    .style("background","#001").style("border-radius","8px")
    .style("max-width","100%").style("height","auto");

  // Center the globe on the observer's longitude rotated by Earth's spin
  const rot = [-loc.lon, -loc.lat*0.6, 0];
  const proj = d3.geoOrthographic().scale(W*0.45).translate([W/2, H/2]).rotate(rot).clipAngle(90);
  const path = d3.geoPath(proj);

  // Globe sphere
  svg.append("path").datum({type:"Sphere"}).attr("d", path)
    .attr("fill","#001a33").attr("stroke","#345");

  // Wireframe graticule (15° spacing)
  const graticule = d3.geoGraticule().step([15,15]);
  svg.append("path").datum(graticule).attr("d", path)
    .attr("fill","none").attr("stroke","#356").attr("stroke-width",0.4);
  // Equator highlighted
  svg.append("path").datum(d3.geoGraticule().step([180,180]).extentMinor([[-180,0],[180,0]]))
    .attr("d", path).attr("fill","none").attr("stroke","#588").attr("stroke-width",0.7);

  // Sub-solar / sub-lunar
  const subSun  = subPoint(jdNow, sunPosition);
  const subMoon = subPoint(jdNow, moonPosition);

  // Night side: great circle 90° from sub-solar, hemisphere AWAY from sun is dark
  const antiSun = [(subSun.lon*RAD + 180 + 540) % 360 - 180, -subSun.lat*RAD];
  const night = d3.geoCircle().radius(90).center(antiSun)();
  svg.append("path").datum(night).attr("d", path)
    .attr("fill","rgba(0,0,30,0.55)").attr("stroke","#446").attr("stroke-width",0.5);

  // Sub-solar point
  svg.append("path").datum({type:"Point", coordinates:[subSun.lon*RAD, subSun.lat*RAD]})
    .attr("d", path.pointRadius(6))
    .attr("fill","#ffcc33").attr("stroke","#aa6");
  // Sub-lunar point
  svg.append("path").datum({type:"Point", coordinates:[subMoon.lon*RAD, subMoon.lat*RAD]})
    .attr("d", path.pointRadius(5))
    .attr("fill","#eee").attr("stroke","#888");
  // Observer
  svg.append("path").datum({type:"Point", coordinates:[loc.lon, loc.lat]})
    .attr("d", path.pointRadius(3.5))
    .attr("fill","#f44").attr("stroke","#fff").attr("stroke-width",0.5);

  // Inset readout
  const sunAlt = altAz(loc.lat*DEG, sun.delta, sun.alpha, jdNow, loc.lon).alt*RAD;
  const moonAlt = altAz(loc.lat*DEG, moon.delta, moon.alpha, jdNow, loc.lon).alt*RAD;
  svg.append("text").attr("x", 8).attr("y", 16).attr("fill","#bbb").attr("font-size","11px")
    .text(`Sun alt ${sunAlt.toFixed(1)}°, Moon alt ${moonAlt.toFixed(1)}°`);
  svg.append("text").attr("x", 8).attr("y", H-8).attr("fill","#888").attr("font-size","10px")
    .text(`☉ sub-solar (${(subSun.lat*RAD).toFixed(1)}°, ${(subSun.lon*RAD).toFixed(1)}°)`);

  return svg.node();
}

5. Full moon calendar & seasonal pattern

A year of full moon timings

Every row below is a real moonrise–sunset offset for the next 13 lunations at your selected location. Negative values mean the moon rises before sunset.

Code 
fullMoonsTable = {
  // Find approximate full moons over the next ~13 months by scanning for elongation = π
  const start = jdFromDate(yearToday, todayObj.getUTCMonth()+1, todayObj.getUTCDate());
  const rows = [];
  let prev = null;
  let lastFound = -1e9;
  for (let d = 0; d < 400; d += 0.05) {
    const jd = start + d;
    const elong = normRad(moonPosition(jd).lambda - sunPosition(jd).lambda);
    const dev = elong - Math.PI;
    if (prev != null && prev * dev < 0 && Math.abs(prev - dev) < 0.5 && (jd - lastFound) > 20) {
      // refine by linear interpolation
      const jdFull = jd - 0.05 * dev/(dev - prev);
      lastFound = jdFull;
      // sunset on the full-moon's local calendar date
      const localJd = jdFull + loc.tz/24;
      const dateInt = Math.floor(localJd - 0.5) + 0.5; // local 0h UT-equivalent
      const jd0 = dateInt - loc.tz/24; // back to UT 0h of that local day
      const ev = riseSet(sunPosition,  sunH0,  jd0, loc.lat*DEG, loc.lon);
      const em = riseSet(moonPosition, moonH0, jd0, loc.lat*DEG, loc.lon);
      const dt = (ev.set.jd != null && em.rise.jd != null)
        ? (em.rise.jd - ev.set.jd) * 1440
        : null;
      // age past full at "1h after sunset"
      let age = null, illumAtTarget = null;
      if (ev.set.jd != null) {
        const tTarget = ev.set.jd + 1/24;
        age = (tTarget - jdFull); // days
        illumAtTarget = illumination(tTarget).illuminated;
      }
      const dObj = new Date((jdFull + loc.tz/24 - 2440587.5)*86400000);
      rows.push({
        date: dObj.toISOString().slice(0,10),
        sunset:    fmtLocal(ev.set.jd, loc.tz),
        moonrise:  fmtLocal(em.rise.jd, loc.tz),
        deltaMin:  dt==null ? "—" : dt.toFixed(1),
        ageDays:   age==null ? "—" : age.toFixed(2),
        illum1h:   illumAtTarget==null ? "—" : (illumAtTarget*100).toFixed(1)+"%"
      });
    }
    prev = dev;
  }
  return rows;
}

Inputs.table(fullMoonsTable, {
  columns: ["date","sunset","moonrise","deltaMin","ageDays","illum1h"],
  header:  {date:"Full moon (local)", sunset:"Sunset", moonrise:"Moonrise",
            deltaMin:"Δ (min)", ageDays:"Age at sunset+1h (d)", illum1h:"Illum at sunset+1h"}
})

The seasonal pattern is now reproducible: in late autumn at northern mid-latitudes, the moonrise–sunset offset goes negative (full moon already up at sunset), so a viewer wanting a “rising-just-as-night-falls” moon must actually catch one \(1{-}2\) days past full, when illumination has fallen to \(\sim 95\%\).

Moonrise solver

Given any calendar date, what lunar phase gives a moonrise exactly \(1\) hour after sunset? Solve for the elongation \(D = \lambda_M - \lambda_\odot\) that makes the moonrise event align with sunset \(+ 1\) h. We treat the true lunar series of §1 as known and adjust a phase offset \(\Delta D\) added to the natural elongation that day; in practice we just bisect on \(\Delta D \in (-\pi, \pi]\).

Code 
solverResult = {
  const target = events.sun.set.jd;
  if (target == null) return md`*Sun never sets at this latitude/date — no solution.*`;
  const targetT = target + 1/24;
  // Find the lunar elongation D (lambda_M - lambda_sun) at which a body with
  // alpha_M, delta_M produced by adjusting only longitude by ΔD rises at targetT.
  // Bisection on ΔD: residual = (LST(rise) - alpha_M_shifted) + H0
  function residual(dD) {
    // Build a synthetic moon: take the moon's β, but replace longitude with sun.lambda + π + dD.
    const lam = sunPosition(targetT).lambda + Math.PI + dD;
    const eps = obliquity(julianCenturies(targetT));
    const synth = {
      delta: Math.asin(Math.sin(moon.beta)*Math.cos(eps) + Math.cos(moon.beta)*Math.sin(eps)*Math.sin(lam)),
      alpha: normRad(Math.atan2(Math.sin(lam)*Math.cos(eps) - Math.tan(moon.beta)*Math.sin(eps), Math.cos(lam))),
      distance: moon.distance, beta: moon.beta, lambda: lam
    };
    const h0 = moonH0(synth);
    const cosH0 = (Math.sin(h0) - Math.sin(loc.lat*DEG)*Math.sin(synth.delta))
                 / (Math.cos(loc.lat*DEG)*Math.cos(synth.delta));
    if (cosH0 > 1 || cosH0 < -1) return NaN;
    const H0 = Math.acos(cosH0);
    let H = lstFor(targetT, loc.lon) - synth.alpha; // current hour angle at target time
    while (H >  Math.PI) H -= 2*Math.PI;
    while (H < -Math.PI) H += 2*Math.PI;
    return H + H0; // zero when target is exactly at moonrise
  }
  let lo = -Math.PI*0.45, hi = Math.PI*0.45;
  let fLo = residual(lo), fHi = residual(hi);
  // expand search
  if (!(fLo*fHi < 0)) {
    // brute search
    let best = {dD:0, f:Infinity}, last = residual(-Math.PI*0.99);
    for (let i=-99; i<=99; i++) {
      const x = i/100*Math.PI;
      const f = residual(x);
      if (Number.isFinite(f) && Math.abs(f) < Math.abs(best.f)) best = {dD:x, f};
      if (Number.isFinite(f) && Number.isFinite(last) && f*last < 0) { lo = x - Math.PI/100; hi = x; break; }
      last = f;
    }
    fLo = residual(lo); fHi = residual(hi);
  }
  for (let i=0; i<60 && fLo*fHi < 0; i++) {
    const mid = (lo+hi)/2, fm = residual(mid);
    if (fm*fLo < 0) { hi = mid; fHi = fm; } else { lo = mid; fLo = fm; }
  }
  const dD = (lo+hi)/2;
  // True synodic phase at the rise event: sun.lambda + π + dD - sun.lambda = π + dD
  const phase = Math.PI + dD;             // 0 = new, π = full
  const illumPct = (1 - Math.cos(phase))/2 * 100;
  const ageDaysFromFull = dD * 29.5306 / (2*Math.PI);
  return md`
**For ${currentDateLabel} at (${loc.lat.toFixed(2)}°, ${loc.lon.toFixed(2)}°):**
- Sunset: \`${fmtLocal(events.sun.set.jd, loc.tz)}\`
- Required elongation $D$: **${(phase*RAD).toFixed(2)}°** (full = 180°)
- Equivalent age: **${ageDaysFromFull.toFixed(2)} days** ${ageDaysFromFull > 0 ? "past full" : "before full"}
- Illumination at moonrise: **${illumPct.toFixed(1)}%**
`;
}

Seasonal phase pattern

The “striking moonrise” — the moon rising one hour after sunset — sweeps through the lunar phases over the course of a year. Below is a single plot of the required lunar phase, computed for every day of the next two years at your selected location. The horizontal axis is calendar date; the vertical axis is the synodic phase angle \(D = \lambda_M - \lambda_\odot\) (in degrees) needed for moonrise to land exactly at sunset \(+1\) h.

Code 
seasonalCurve = {
  // For each day of the next ~2 years, find the elongation D such that
  // moonrise occurs at sunset+1h. Reuse the bisection from §10.
  const phi = loc.lat*DEG;
  const startJd = jdFromDate(yearToday, todayObj.getUTCMonth()+1, todayObj.getUTCDate());
  const rows = [];
  for (let k = 0; k < 730; k += 3) {
    const jd0_k = startJd + k;
    const ev = riseSet(sunPosition, sunH0, jd0_k, phi, loc.lon);
    if (ev.set.jd == null) continue;
    const target = ev.set.jd + 1/24;
    const sunAtTarget = sunPosition(target);
    function residual(dD) {
      const lam = sunAtTarget.lambda + Math.PI + dD;
      const eps = obliquity(julianCenturies(target));
      const beta = 0; // approximate: ignore moon's β for the seasonal sweep
      const dec = Math.asin(Math.sin(eps)*Math.sin(lam));
      const ra = normRad(Math.atan2(Math.cos(eps)*Math.sin(lam), Math.cos(lam)));
      const cosH0 = (Math.sin(0.125*DEG) - Math.sin(phi)*Math.sin(dec))
                   / (Math.cos(phi)*Math.cos(dec));
      if (cosH0 > 1 || cosH0 < -1) return NaN;
      const H0 = Math.acos(cosH0);
      let H = lstFor(target, loc.lon) - ra;
      while (H >  Math.PI) H -= 2*Math.PI;
      while (H < -Math.PI) H += 2*Math.PI;
      return H + H0;
    }
    // bisection on dD ∈ (-π/2, π/2)
    let lo = -Math.PI*0.4, hi = Math.PI*0.4;
    let fLo = residual(lo), fHi = residual(hi);
    if (!(Number.isFinite(fLo) && Number.isFinite(fHi) && fLo*fHi < 0)) continue;
    for (let i=0; i<50; i++) {
      const mid = (lo+hi)/2, fm = residual(mid);
      if (!Number.isFinite(fm)) break;
      if (fm*fLo < 0) { hi = mid; fHi = fm; } else { lo = mid; fLo = fm; }
    }
    const dD = (lo+hi)/2;
    const D_deg = (Math.PI + dD) * RAD;       // 0 = new, 180 = full
    const illum = (1 - Math.cos(Math.PI + dD))/2 * 100;
    rows.push({
      date: new Date((jd0_k + 0.5 - 2440587.5)*86400000),
      phaseDeg: D_deg,
      illum
    });
  }
  return rows;
}

Plot.plot({
  width: 900, height: 320,
  marginLeft: 60, marginRight: 60,
  x: { type: "utc", label: "Date" },
  y: { label: "Required phase angle D (°) — 180 = full", domain: [150, 210] },
  marks: [
    Plot.ruleY([180], { stroke: "#888", strokeDasharray: "3,3" }),
    Plot.text([{date: seasonalCurve[0]?.date, label: "exact full"}],
      {x: "date", y: () => 180.6, text: "label", fill: "#888", fontSize: 11, dx: 6}),
    Plot.line(seasonalCurve, { x: "date", y: "phaseDeg", stroke: "#B35A44", strokeWidth: 2 }),
    Plot.dot(seasonalCurve.filter((_,i)=>i%5===0), {
      x: "date", y: "phaseDeg",
      fill: d => d.phaseDeg > 180 ? "#B35A44" : "#1B2E26",
      r: 2.5
    })
  ]
})

How to read the plot. The horizontal line at \(D = 180°\) is “exact full.” The curve sits above the line when the striking moon must be past full (waning gibbous) — i.e., a full moon would have already risen earlier than sunset+1h, so we need to wait a day or two for the moon to drop into position. The curve sits below the line when the striking moon must be just before full (waxing gibbous) — a “true” full moon would not yet have cleared the horizon.

For mid-northern observers (D.C., \(\sim 39°\text{N}\)), the pattern is:

Season (full-moon date roughly) Curve relative to \(180°\) Phase you actually see at sunset+1h
Late spring / early summer \(D < 180°\) (a few °) Near-full waxing gibbous; full moon rises after sunset, so at sunset+1h it has just cleared the horizon
Midsummer very close to \(180°\) Essentially full
Late summer → autumnal equinox \(D\) rises above \(180°\) Slightly waning gibbous; the full moon rose before sunset (harvest-moon geometry)
Late autumn / early winter \(D \approx 185°-195°\) — peak deviation Waning gibbous \(\sim 92\)–\(97\%\) — the full moon was already up well before sunset; you must wait \(1\)–\(2\) days past full to catch a moonrise an hour after sunset
Late winter / early spring \(D\) returns toward \(180°\) Back to near-full waxing gibbous

This corrects a common intuition: it is late autumn / early winter — not summer — when the striking moon visible 1 hour after sunset is a clearly waning gibbous, because that is when the geometry pushes the full moon’s natural rise time earliest relative to sunset. The “harvest moon” of September is the gateway into this regime; by November the deviation peaks. Conversely, summer full moons rise late (south-of-equator declination, short time above horizon) so a sunset+1h viewing catches them essentially at full.

In the southern hemisphere the pattern flips by six months. At the equator the curve is nearly flat at \(D \approx 180°\) year-round.