tower / files / src / world / celestial.ts
celestial.ts
| 1 | /** |
| 2 | * The real sky behind Tower: which way the stars lie from its axis at a given instant. |
| 3 | * See design.md §1.6 and decisions.md D-075. |
| 4 | * |
| 5 | * Tower's axis points at the sun (§1.6), so the sky it sees is fixed by one number, the |
| 6 | * sun's ecliptic longitude — the far cap looks at the antisolar point, and that walks |
| 7 | * round the ecliptic once a year. Roll about the axis is not fixed by anything outside: |
| 8 | * inside it is the spin, and outside this code puts **ecliptic north up** (scene +Y), |
| 9 | * which is the one choice that is square to the axis on every day of the year. |
| 10 | * |
| 11 | * Scene frame, as everywhere in the 3D views: +Z runs from the sun cap to the far cap, |
| 12 | * so the sun is at −Z. The matrix here takes a direction in that frame to J2000 |
| 13 | * equatorial coordinates, the frame the star catalogue is written in. |
| 14 | * |
| 15 | * Nothing in here knows about Three.js, so it is tested directly (test/celestial.test.ts). |
| 16 | */ |
| 17 | |
| 18 | const RAD = Math.PI / 180 |
| 19 | /** Obliquity of the ecliptic at J2000. It drifts 47″ a century, which nobody can see. */ |
| 20 | const OBLIQUITY = 23.4392911 * RAD |
| 21 | /** General precession in ecliptic longitude, per Julian year. */ |
| 22 | const PRECESSION_PER_YEAR = (50.29 / 3600) * RAD |
| 23 | const J2000_MS = Date.UTC(2000, 0, 1, 12) |
| 24 | const DAY_MS = 86_400_000 |
| 25 | |
| 26 | /** |
| 27 | * The sun's apparent ecliptic longitude at an instant, radians, of date. The |
| 28 | * Astronomical Almanac's low-precision formula: good to 0.01° between 1950 and 2050, |
| 29 | * a fiftieth of the sun's own width. |
| 30 | */ |
| 31 | export function sunLongitude(at: Date): number { |
| 32 | const n = (at.getTime() - J2000_MS) / DAY_MS |
| 33 | const L = 280.46 + 0.9856474 * n |
| 34 | const g = (357.528 + 0.9856003 * n) * RAD |
| 35 | const lambda = (L + 1.915 * Math.sin(g) + 0.02 * Math.sin(2 * g)) * RAD |
| 36 | return ((lambda % (2 * Math.PI)) + 2 * Math.PI) % (2 * Math.PI) |
| 37 | } |
| 38 | |
| 39 | /** |
| 40 | * Earth's Moon's phase at an instant: how much of its disc is lit, 0 at new and 1 at full, |
| 41 | * and whether it is waxing. Tower's moon keeps it (D-088), as its sun keeps Pisa's days. |
| 42 | * |
| 43 | * Meeus, *Astronomical Algorithms* ch. 48, low precision: the phase angle from the mean |
| 44 | * elongation and the two mean anomalies, with the six largest periodic terms. Good to about |
| 45 | * a hundredth of the disc, which is finer than an eye can judge a moon. |
| 46 | */ |
| 47 | export function moonPhase(at: Date): { fraction: number; waxing: boolean } { |
| 48 | const T = (at.getTime() - J2000_MS) / DAY_MS / 36525 |
| 49 | const D = (297.8501921 + 445267.1114034 * T) * RAD |
| 50 | const M = (357.5291092 + 35999.0502909 * T) * RAD |
| 51 | const Mp = (134.9633964 + 477198.8675055 * T) * RAD |
| 52 | const i = Math.PI - D |
| 53 | - 6.289 * RAD * Math.sin(Mp) + 2.1 * RAD * Math.sin(M) - 1.274 * RAD * Math.sin(2 * D - Mp) |
| 54 | - 0.658 * RAD * Math.sin(2 * D) - 0.214 * RAD * Math.sin(2 * Mp) - 0.11 * RAD * Math.sin(D) |
| 55 | // Waxing while the moon's elongation from the sun runs from 0° to 180°. |
| 56 | return { fraction: (1 + Math.cos(i)) / 2, waxing: Math.sin(D) > 0 } |
| 57 | } |
| 58 | |
| 59 | /** |
| 60 | * Scene direction → J2000 equatorial direction at an instant, as a row-major 3×3. |
| 61 | * |
| 62 | * Built in the ecliptic: the axis points away from the sun, and ecliptic north is up. |
| 63 | * The sun's longitude is of date and the catalogue's are of J2000, so the sun is moved |
| 64 | * back by the precession since then rather than every star forward by it — a third of |
| 65 | * a degree by 2026, which is the width of the moon. |
| 66 | */ |
| 67 | export function skyFrame(at: Date): number[] { |
| 68 | const years = (at.getTime() - J2000_MS) / (365.25 * DAY_MS) |
| 69 | const lam = sunLongitude(at) - PRECESSION_PER_YEAR * years |
| 70 | // The scene's axes written in J2000 ecliptic coordinates. X = Y × Z, so the frame is |
| 71 | // right-handed like the scene's and the sky is not drawn mirrored. |
| 72 | const X = [Math.sin(lam), -Math.cos(lam), 0] |
| 73 | const Y = [0, 0, 1] |
| 74 | const Z = [-Math.cos(lam), -Math.sin(lam), 0] |
| 75 | const c = Math.cos(OBLIQUITY) |
| 76 | const s = Math.sin(OBLIQUITY) |
| 77 | // Ecliptic → equatorial is a turn about the equinox direction by the obliquity. |
| 78 | const toEq = (v: number[]) => [v[0]!, v[1]! * c - v[2]! * s, v[1]! * s + v[2]! * c] |
| 79 | const cols = [toEq(X), toEq(Y), toEq(Z)] |
| 80 | return [0, 1, 2].flatMap((row) => cols.map((col) => col[row]!)) |
| 81 | } |
| 82 | |
| 83 | /** A J2000 unit vector from right ascension and declination, degrees. */ |
| 84 | export function equatorial(raDeg: number, decDeg: number): [number, number, number] { |
| 85 | const a = raDeg * RAD |
| 86 | const d = decDeg * RAD |
| 87 | return [Math.cos(d) * Math.cos(a), Math.cos(d) * Math.sin(a), Math.sin(d)] |
| 88 | } |
| 89 | |
| 90 | // ------------------------------------------------------------------ the bake |
| 91 | |
| 92 | /** |
| 93 | * Cells along a cube face's edge. A cell is a third of a degree across at the middle |
| 94 | * of a face and a seventh at its corners, where 9,000 stars 2° apart on average almost |
| 95 | * never share one. |
| 96 | */ |
| 97 | export const STAR_FACE_CELLS = 256 |
| 98 | /** A face's grid with one cell of border all round (see `bakeStars`). */ |
| 99 | export const STAR_FACE_SIZE = STAR_FACE_CELLS + 2 |
| 100 | /** The faintest star in the catalogue is about this; its flux is 1. */ |
| 101 | const MAG_LIMIT = 6.5 |
| 102 | |
| 103 | /** |
| 104 | * Which cube face a direction falls on, and where: face = 2·axis + (negative ? 1 : 0), |
| 105 | * and (u, v) in −1..1 across it. Each axis takes the next two round as its u and v — |
| 106 | * x: (y, z), y: (z, x), z: (x, y) — and STAR_GLSL uses the same rule, so the two agree |
| 107 | * on every face without either needing the other's vectors. |
| 108 | */ |
| 109 | function project(d: readonly number[], axis: number): [number, number] { |
| 110 | const n = d[axis]! |
| 111 | return [d[(axis + 1) % 3]! / Math.abs(n), d[(axis + 2) % 3]! / Math.abs(n)] |
| 112 | } |
| 113 | |
| 114 | /** |
| 115 | * The catalogue as a texture the sky shader can look stars up in: six cube faces, each a |
| 116 | * grid of cells, each cell holding at most one star as (u, v, flux, B−V) — its position on |
| 117 | * that face, its brightness relative to the faintest the eye can see, and its colour. |
| 118 | * |
| 119 | * A shader cannot loop over 9,000 stars per pixel, but it can look in the nine cells |
| 120 | * round the one it is in. At a face's edge some of those nine are on the next face, so |
| 121 | * every grid has a border of one cell holding the neighbouring faces' stars projected |
| 122 | * onto this one — a star near an edge is baked twice, once in each face's own terms, and |
| 123 | * a pixel only ever reads the face it is on, so it is never counted twice. |
| 124 | * |
| 125 | * Two stars in one cell are a double the eye would not split: their light is added, |
| 126 | * and the brighter one's place and colour kept. |
| 127 | * |
| 128 | * `stars` is the catalogue's flat [raDeg, decDeg, vmag, bv, …]. Returns RGBA float |
| 129 | * texels, STAR_FACE_SIZE wide and 6 × STAR_FACE_SIZE tall, face by face. |
| 130 | */ |
| 131 | export function bakeStars(stars: readonly number[]): Float32Array { |
| 132 | const S = STAR_FACE_SIZE |
| 133 | const out = new Float32Array(S * S * 6 * 4) |
| 134 | const brightest = new Float32Array(S * S * 6) |
| 135 | for (let i = 0; i + 3 < stars.length; i += 4) { |
| 136 | const d = equatorial(stars[i]!, stars[i + 1]!) |
| 137 | const flux = 10 ** (-0.4 * (stars[i + 2]! - MAG_LIMIT)) |
| 138 | for (let axis = 0; axis < 3; axis++) { |
| 139 | if (d[axis] === 0) continue |
| 140 | const face = 2 * axis + (d[axis]! < 0 ? 1 : 0) |
| 141 | const [u, v] = project(d, axis) |
| 142 | const cu = Math.floor(((u + 1) / 2) * STAR_FACE_CELLS) + 1 |
| 143 | const cv = Math.floor(((v + 1) / 2) * STAR_FACE_CELLS) + 1 |
| 144 | if (cu < 0 || cu >= S || cv < 0 || cv >= S) continue |
| 145 | const cell = (face * S + cv) * S + cu |
| 146 | const o = cell * 4 |
| 147 | if (flux > brightest[cell]!) { |
| 148 | brightest[cell] = flux |
| 149 | out[o] = u |
| 150 | out[o + 1] = v |
| 151 | out[o + 3] = stars[i + 3]! |
| 152 | } |
| 153 | out[o + 2] = out[o + 2]! + flux |
| 154 | } |
| 155 | } |
| 156 | return out |
| 157 | } |