celestial.js
184 lines · in natura-src-1.0.0.zip as natura/celestial.js · a module · its entry in the module reference · its test: celestial-test.js
// celestial.js -- Natura CELESTIAL (SYSTEMS.md, ALGORITHMS.md): the abstract clock made PHYSICAL.
// A hierarchy of orbiting bodies (a star at the origin, a planet that spins on a tilted axis, a moon
// that circles the planet) whose real geometry GENERATES what the middle-scale systems consume: the
// sun's altitude and the daylight at the world's latitude, the seasons from axial tilt, the moon's
// phase, and the tides (semidiurnal, with spring/neap from sun-moon alignment). It reaches UP to the
// cosmic scale -- the body list is an orrery you can add planets and moons to -- and DOWN into the
// climate: weather and vegetation can read celestial.daylight / .decl instead of the clock's heuristic,
// so the sun physically drives the world. Everything is ANALYTIC in the clock time (no integration
// drift, fully resumable) and deterministic. Couples to: (out) weather, vegetation, ocean/tides.
//
// A Natura system spec (init/step/fields/events). UMD (node test / window.Celestial).
(function (root, factory) {
if (typeof module !== 'undefined' && module.exports) module.exports = factory();
else root.Celestial = factory();
})(typeof self !== 'undefined' ? self : this, function () {
'use strict';
const TAU = 2 * Math.PI, DEG = Math.PI / 180;
const clamp = (v, a, b) => v < a ? a : v > b ? b : v;
const wrap = (a) => a - TAU * Math.floor(a / TAU); // -> [0, TAU)
// solve Kepler's equation M = E - e sinE for the eccentric anomaly (Newton, converges in a few iters)
function kepler(M, e) {
M = wrap(M); let E = e < 0.8 ? M : Math.PI;
for (let i = 0; i < 6; i++) { const f = E - e * Math.sin(E) - M; E -= f / (1 - e * Math.cos(E)); }
return E;
}
// make a Celestial system. config: { bodies, planet, latitude, longitude, axialTilt, rotationPeriod,
// moonInclination, sunTideRatio }. bodies default to a Sun / Earth-like planet / Moon.
function makeCelestial(config) {
config = config || {};
const YEAR = 365 * 86400, DAY = 86400;
// the orbiting bodies. each: { name, parent(name|null), a(orbit radius, AU), period(s), phase0(rad),
// e(eccentricity), inc(orbit inclination, rad), node(ascending-node longitude, rad),
// nodePeriod(node regression period, s), radius(AU, for eclipses/viz) }
const bodies = config.bodies || [
{ name: 'star', parent: null, a: 0, period: 1, phase0: 0, e: 0, radius: 0.00465 },
{ name: 'planet', parent: 'star', a: 1.0, period: YEAR, phase0: 0, e: 0.0167, radius: 4.26e-5 },
{ name: 'moon', parent: 'planet', a: 0.00257, period: 27.32 * DAY, phase0: 1.1, e: 0.0549,
inc: 5.14 * DEG, node: 0, nodePeriod: -18.6 * YEAR, radius: 1.16e-5 }
];
const byName = {}; bodies.forEach(b => { byName[b.name] = b; });
const P = {
planet: config.planet || 'planet',
moon: config.moon || 'moon',
latitude: config.latitude != null ? config.latitude : 45 * DEG, // the world's latitude on the planet
longitude: config.longitude != null ? config.longitude : 0,
axialTilt: config.axialTilt != null ? config.axialTilt : 23.44 * DEG,
rotationPeriod: config.rotationPeriod != null ? config.rotationPeriod : DAY,
sunTideRatio: config.sunTideRatio != null ? config.sunTideRatio : 0.46, // solar/lunar tide amplitude
springOffset: config.springOffset != null ? config.springOffset : 0 // year-phase of the spring equinox
};
// heliocentric-ish position (x,y in the orbital plane, z from inclination) of a body at time t
function posOf(b, t) {
if (!b.parent) return { x: 0, y: 0, z: 0 };
const parent = posOf(byName[b.parent], t);
const M = b.phase0 + TAU * t / b.period;
const e = b.e || 0;
const E = kepler(M, e);
const r = b.a * (1 - e * Math.cos(E));
const nu = Math.atan2(Math.sqrt(1 - e * e) * Math.sin(E), Math.cos(E) - e); // true anomaly
let x = r * Math.cos(nu), y = r * Math.sin(nu), z = 0;
if (b.inc) { // tilt the orbit about the node line
const node = (b.node || 0) + (b.nodePeriod ? TAU * t / b.nodePeriod : 0);
const xr = x * Math.cos(-node) - y * Math.sin(-node), yr = x * Math.sin(-node) + y * Math.cos(-node);
z = yr * Math.sin(b.inc); const yy = yr * Math.cos(b.inc);
x = xr * Math.cos(node) - yy * Math.sin(node); y = xr * Math.sin(node) + yy * Math.cos(node);
}
return { x: parent.x + x, y: parent.y + y, z: parent.z + z };
}
function compute(t) {
const planet = byName[P.planet], moon = byName[P.moon];
const pPos = posOf(planet, t);
// ecliptic longitude of the SUN as seen from the planet (opposite the planet's heliocentric dir)
const sunLon = Math.atan2(-pPos.y, -pPos.x);
// solar declination from axial tilt: the subsolar latitude oscillates over the year
const yearAngle = TAU * (t / planet.period) - TAU * P.springOffset;
const decl = P.axialTilt * Math.sin(yearAngle);
// planet rotation -> the sun's hour angle at the world's longitude (H = 0 at local noon)
const dayFrac = wrap(TAU * t / P.rotationPeriod + P.longitude) / TAU;
const H = (dayFrac - 0.5) * TAU;
const lat = P.latitude;
const sinAlt = Math.sin(lat) * Math.sin(decl) + Math.cos(lat) * Math.cos(decl) * Math.cos(H);
const sunAlt = Math.asin(clamp(sinAlt, -1, 1));
const daylight = Math.max(0, sinAlt); // 0 night .. 1 sun overhead
// day length at this latitude/season (sunrise equation), with polar day/night handled
const cosH0 = -Math.tan(lat) * Math.tan(decl);
const dayLen = cosH0 >= 1 ? 0 : cosH0 <= -1 ? 24 : 24 * Math.acos(cosH0) / Math.PI;
const insol = daylight; // top-of-atmosphere insolation proxy
// MOON: phase from the elongation between moon and sun as seen from the planet
const mPos = posOf(moon, t);
const mRel = { x: mPos.x - pPos.x, y: mPos.y - pPos.y, z: mPos.z - pPos.z };
const sRel = { x: -pPos.x, y: -pPos.y, z: -pPos.z }; // sun direction from planet
const dotms = mRel.x * sRel.x + mRel.y * sRel.y + mRel.z * sRel.z;
const lm = Math.hypot(mRel.x, mRel.y, mRel.z) || 1, ls = Math.hypot(sRel.x, sRel.y, sRel.z) || 1;
const elong = Math.acos(clamp(dotms / (lm * ls), -1, 1)); // 0 = new (toward sun), PI = full
const moonIllum = (1 - Math.cos(elong)) / 2; // illuminated fraction
const moonPhase = elong / Math.PI; // 0 new .. 1 full (magnitude)
// moon hour angle: the moon's own longitude vs the world meridian -> the semidiurnal tide
const moonLon = Math.atan2(mRel.y, mRel.x);
const Hmoon = H + (sunLon - moonLon); // moon lags/leads the sun by its elongation
const tideMoon = Math.cos(2 * Hmoon), tideSun = P.sunTideRatio * Math.cos(2 * H);
const tide = (tideMoon + tideSun) / (1 + P.sunTideRatio); // ~ -1 (low) .. 1 (high)
const springNeap = Math.abs(Math.cos(elong)); // 1 at new/full (spring), 0 at quarter (neap)
// moon ecliptic latitude (for eclipse gating): small unless the moon is near a node
const beta = Math.asin(clamp(mRel.z / lm, -1, 1));
const solarEclipse = elong < 0.04 && Math.abs(beta) < 0.009; // new moon near a node
const lunarEclipse = (Math.PI - elong) < 0.04 && Math.abs(beta) < 0.009; // full moon near a node
return {
t, decl, sunAlt, daylight, insol, dayLengthHours: dayLen,
moonPhase, moonIllum, elong, tide, springNeap, solarEclipse, lunarEclipse,
sunLon, yearAngle,
bodies: bodies.map(b => { const q = posOf(b, t); return { name: b.name, x: q.x, y: q.y, z: q.z, r: b.radius }; })
};
}
return {
name: config.name || 'celestial',
domain: 'cosmic',
scale: 'planetary',
params: P,
init(rng, params, clock) {
const s = compute(clock ? clock.t : 0);
return { cur: s, prevDecl: s.decl, prevIllum: s.moonIllum, prevElong: s.elong,
events: [], eclipseArmed: true };
},
step(st, ctx, dt) {
const t = ctx.clock.t;
const s = compute(t);
const ev = [];
// EQUINOX / SOLSTICE from the declination turning points and zero-crossings
if (st.prevDecl <= 0 && s.decl > 0) ev.push({ kind: 'equinox', which: 'spring', decl: s.decl });
if (st.prevDecl >= 0 && s.decl < 0) ev.push({ kind: 'equinox', which: 'autumn', decl: s.decl });
const dPrev = s.decl - st.prevDecl;
if (st.pDeclSlope != null && st.pDeclSlope > 0 && dPrev <= 0) ev.push({ kind: 'solstice', which: 'summer', decl: s.decl });
if (st.pDeclSlope != null && st.pDeclSlope < 0 && dPrev >= 0) ev.push({ kind: 'solstice', which: 'winter', decl: s.decl });
// NEW / FULL moon from the illuminated-fraction turning points
const dIllum = s.moonIllum - st.prevIllum;
if (st.pIllumSlope != null && st.pIllumSlope < 0 && dIllum >= 0) ev.push({ kind: 'moon', which: 'new', illum: s.moonIllum });
if (st.pIllumSlope != null && st.pIllumSlope > 0 && dIllum <= 0) ev.push({ kind: 'moon', which: 'full', illum: s.moonIllum });
// ECLIPSES (edge-triggered so one event per alignment)
if (s.solarEclipse && st.eclipseArmed) { ev.push({ kind: 'eclipse', which: 'solar' }); st.eclipseArmed = false; }
else if (s.lunarEclipse && st.eclipseArmed) { ev.push({ kind: 'eclipse', which: 'lunar' }); st.eclipseArmed = false; }
else if (!s.solarEclipse && !s.lunarEclipse) st.eclipseArmed = true;
st.pDeclSlope = dPrev; st.prevDecl = s.decl;
st.pIllumSlope = dIllum; st.prevIllum = s.moonIllum; st.prevElong = s.elong;
st.cur = s; st.events = ev;
return st;
},
// what the rest of the world reads: the sky + the tides (weather/veg use daylight/decl; ocean uses tide)
fields(st) {
const s = st.cur;
return {
daylight: s.daylight, decl: s.decl, sunAltitude: s.sunAlt, insolation: s.insol,
dayLengthHours: s.dayLengthHours,
moonPhase: s.moonPhase, moonIllum: s.moonIllum,
tide: s.tide, springNeap: s.springNeap,
season: seasonOf(s.decl, st.pDeclSlope),
bodies: s.bodies
};
},
events(st) { return st.events; }
};
// season from declination + whether it is rising or falling (N-hemisphere convention)
function seasonOf(decl, slope) {
const tilt = P.axialTilt, up = slope == null ? true : slope > 0;
if (decl > tilt * 0.5) return 'summer';
if (decl < -tilt * 0.5) return 'winter';
return up ? 'spring' : 'autumn';
}
}
return { makeCelestial };
});