source ⟩ hypertext ⟩ public ⟩ sandbox ⟩ orrery_calculations.js

const AU = 1.495978707e11; // metres
const LIGHTSECOND = 299792458; // metres

const CENTURY = 36525; // Julian days
const J2000_EPOCH = new Date("2000-01-01T12:00:00Z");
const DEGREE = Math.PI / 180;

const GRAV_PARAM = 1.32712440042e20; //m³·s⁻²

const MAX_MULTIREVOLUTIONS = 5;

const getJ2000Offset = date => (date - J2000_EPOCH) / 86400000;

const adjustElements = (orbit, centuriesSinceEpoch) =>
	Object.fromEntries(
		Object.entries(orbit).map(param => [
			param[0],
			param[1]
				.map((x, idx) => (idx ? x * centuriesSinceEpoch ** idx : x))
				.reduce((a, b) => a + b, 0)
		])
	);

const findEccentricAnomaly = (ecc, meanAnom) => {
	const ε = 1e-6;
	const f = eccAnom => eccAnom - ecc * Math.sin(eccAnom) - meanAnom;
	const f_ = eccAnom => 1 - ecc * Math.cos(eccAnom);

	let estEccAnom = ecc > 0.8 ? Math.PI : meanAnom;
	while (Math.abs(f(estEccAnom)) > ε) {
		estEccAnom -= f(estEccAnom) / f_(estEccAnom);
	}
	return estEccAnom;
};

const getElements = (planet, daysSinceJ2000) => {
	let orbit = adjustElements(planet, daysSinceJ2000 / CENTURY);

	orbit.periapsisArg = orbit.periapsisLon - orbit.ascendingNodeLon;

	orbit.meanAnom = (orbit.meanLon - orbit.periapsisLon) % (Math.PI * 2);

	orbit.eccAnom = findEccentricAnomaly(orbit.eccentricity, orbit.meanAnom);

	orbit.trueAnom =
		2 *
		Math.atan(
			Math.tan(orbit.eccAnom / 2) *
				Math.sqrt((1 + orbit.eccentricity) / (1 - orbit.eccentricity))
		);

	const distance =
		orbit.semiMajorAxis *
		(1 - orbit.eccentricity * Math.cos(orbit.eccAnom));

	const posVec = {
		x: distance * Math.cos(orbit.trueAnom),
		y: distance * Math.sin(orbit.trueAnom)
	};
	const velVec = {
		x:
			(-1 *
				Math.sin(orbit.eccAnom) *
				Math.sqrt(GRAV_PARAM * orbit.semiMajorAxis)) /
			distance,
		y:
			(Math.cos(orbit.eccAnom) *
				Math.sqrt(1 - orbit.eccentricity ** 2) *
				Math.sqrt(GRAV_PARAM * orbit.semiMajorAxis)) /
			distance
	};

	// calculate Cartesian coördinates
	// i = inclination, ω = arg. of periapsis, Ω = lon. of asc. node
	const si = Math.sin(orbit.inclination),
		ci = Math.cos(orbit.inclination),
		sω = Math.sin(orbit.periapsisArg),
		cω = Math.cos(orbit.periapsisArg),
		sΩ = Math.sin(orbit.ascendingNodeLon),
		cΩ = Math.cos(orbit.ascendingNodeLon);

	orbit.r = {
		x:
			posVec.x * (cω * cΩ - sω * ci * sΩ) -
			posVec.y * (sω * cΩ + cω * ci * sΩ),
		y:
			posVec.x * (cω * sΩ + sω * ci * cΩ) +
			posVec.y * (cω * ci * cΩ - sω * sΩ),
		z: posVec.x * sω * si + posVec.y * cω * si
	};
	orbit.r.xyz = [orbit.r.x, orbit.r.y, orbit.r.z];

	orbit.velocity = {
		x:
			velVec.x * (cω * cΩ - sω * ci * sΩ) -
			velVec.y * (sω * cΩ + cω * ci * sΩ),
		y:
			velVec.x * (cω * sΩ + sω * ci * cΩ) +
			velVec.y * (cω * ci * cΩ - sω * sΩ),
		z: velVec.x * sω * si + velVec.y * cω * si
	};
	orbit.velocity.xyz = [orbit.velocity.x, orbit.velocity.y, orbit.velocity.z];

	return orbit;
};

// Returns a 121-length array of orbital information for the past year of a planet,
// with element [0] being for the requested date

const getOrbitTrail = (planet, date, detail = 120) => {
	const dateOffset = getJ2000Offset(date);
	const radiansPerCentury = planet.meanLon[1];
	const daysPerIncrement =
		Math.min((CENTURY * 2 * Math.PI) / (detail * radiansPerCentury), CENTURY / 10);

	let offset = dateOffset;
	let trailData = [];
	for (let idx = detail; idx >= 0; idx--) {
		trailData.push(getElements(planet, offset).r);
		offset -= daysPerIncrement;
	}
	return trailData;
};

export { getJ2000Offset, getElements, getOrbitTrail }