source ⟩ hypertext ⟩ public ⟩ latbs ⟩ tools ⟩ kepler.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 planets = {
	hermes: {
		semiMajorAxis: [57.9090829e9, 0, 0, 0],
		eccentricity: [205.63175e-3, 20.407e-6, -28.3e-9, -180e-12],
		meanLon: [4.40260885, 2.608790314157e3, -93.5e-9, 30e-12],
		inclination: [122.2601e-3, -103.875e-6, 14e-9, 750e-12],
		ascendingNodeLon: [843.53321e-3, -2.189039e-3, -1.542e-6, -3.49e-9],
		periapsisLon: [1.3518643, 2.772705e-3, -234.2e-9, -100e-12]
	},
	venus: {
		semiMajorAxis: [108.208601e9, 0, 0, 0],
		eccentricity: [6.77192e-3, -47.765e-6, 98.1e-9, 460e-12],
		meanLon: [3.1761467, 1.02132855462e3, 28.8e-9, -30e-12],
		inclination: [59.24803e-3, -14.95e-6, -566.2e-9, 150e-12],
		ascendingNodeLon: [1.3383171, -4.8522492e-3, -2.4883e-6, -2.86e-9],
		periapsisLon: [2.29621979, 85.078e-6, -24.1671e-6, -99.4e-9]
	},
	earth: {
		semiMajorAxis: [149.598023e9, 0, 0, 0],
		eccentricity: [16.70863e-3, -42.037e-6, -126.7e-9, 140e-12],
		meanLon: [1.75347046, 628.307584999, -99.1e-9, -200e-12],
		inclination: [0, 227.849e-6, -162e-9, -590e-12],
		ascendingNodeLon: [3.05211269, -4.207828e-3, 743.86e-9, 20e-12],
		periapsisLon: [1.79659565, 5.62983e-3, 2.5829e-6, -680e-12]
	},
	mars: {
		semiMajorAxis: [227.9391852e9, 0, 0, 0],
		eccentricity: [93.40065e-3, 90.484e-6, -80.6e-9, -250e-12],
		meanLon: [6.20347612, 334.06124267, 45.7e-9, -50e-12],
		inclination: [32.28381e-3, -142.2e-6, -393.6e-9, -510e-12],
		ascendingNodeLon: [864.95189e-3, -5.149158e-3, -11.178e-6, -34.28e-9],
		periapsisLon: [5.86535773, 7.747544e-3, -3.0217e-6, 9.04e-9]
	},
	ceres: {
		semiMajorAxis: [413730212640.28864, 0, 0, 0],
		eccentricity: [0.07957631994408416, 0, 0, 0],
		meanLon: [2.7716297474987073, 136.610452404061, 0, 0],
		inclination: [0.18479348168482485, 0, 0, 0],
		ascendingNodeLon: [1.4006202828577676, 0, 0, 0],
		periapsisLon: [2.679942342337361, 0, 0, 0]
	},
	jupiter: {
		semiMajorAxis: [778.2983622e9, 28.62e3, 0, 0],
		eccentricity: [48.49793e-3, 163.225e-6, -471.4e-9, -2.01e-9],
		meanLon: [599.54711e-3, 52.9690962649, -1.484e-6, 280e-12],
		inclination: [22.7463e-3, -34.692e-6, 579.4e-9, 1.7e-9],
		ascendingNodeLon: [1.75343468, 3.084402e-3, 15.83e-6, 126.9e-9],
		periapsisLon: [250.12675e-3, 3.761549e-3, 12.603e-6, -78.28e-9]
	},
	saturn: {
		semiMajorAxis: [1.42939407e12, -319.99e3, 600, 0],
		eccentricity: [55.54814e-3, -346.641e-6, -643.6e-9, 3.4e-9],
		meanLon: [874.01628e-3, 21.3299104958, 3.6659e-6, -800e-12],
		inclination: [43.43913e-3, 44.53e-6, -856.3e-9, 300e-12],
		ascendingNodeLon: [1.98383727, -4.479775e-3, -3.2112e-6, 8.38e-9],
		periapsisLon: [1.6241552, 9.888015e-3, 9.2241e-6, 85.73e-9]
	},
	ouranos: {
		semiMajorAxis: [2.875038609e12, -5.57e3, 150, 0],
		eccentricity: [46.38122e-3, -27.293e-6, 78.9e-9, 240e-12],
		meanLon: [5.48129387, 7.478159856, -84.8e-9, 100e-12],
		inclination: [13.4948e-3, -29.442e-6, 60.9e-9, 280e-12],
		ascendingNodeLon: [1.2916476, 1.29404e-3, 7.0754e-6, 2.08e-9],
		periapsisLon: [3.01951195, 1.55895e-3, -1.653e-6, 7.23e-9]
	},
	neptune: {
		semiMajorAxis: [4.5044497616e12, -24.88e3, 100, 0],
		eccentricity: [9.45575e-3, 6.033e-6, 0, -50e-12],
		meanLon: [5.31188628, 3.813303564, 10e-9, -30e-12],
		inclination: [30.89151e-3, 3.937e-6, 4e-9, 0],
		ascendingNodeLon: [2.3000657, -107.6e-6, -38.2e-9, -1.4e-9],
		periapsisLon: [839.85725e-3, 509.402e-6, 1.328e-6, 0]
	},
	pluto: {
		semiMajorAxis: [5.907150229e12, 672.818e6, 0, 0],
		eccentricity: [248.85238e-3, 60.16e-6, 0, 0],
		meanLon: [4.17073215760049, 2.53387649603146, -220.386913439529e-6, 0],
		inclination: [299.167630594609e-3, 87.4409955249159e-9, 0, 0],
		ascendingNodeLon: [1.92512748403772, -141.368353285962e-6, 0, 0],
		periapsisLon: [3.91123094727827, -169.092210322191e-6, 0, 0]
	},
	persephone: {
		semiMajorAxis: [10.1955736e12, 0, 0, 0],
		eccentricity: [432.3204e-3, 0, 0, 0],
		meanLon: [6.65283224592425, 1.1167140784747, 0, 0],
		inclination: [763.491122238973e-3, 0, 0, 0],
		ascendingNodeLon: [629.666832471709e-3, 0, 0, 0],
		periapsisLon: [3.26292179962947, 0, 0, 0]
	},
	haumea: {
		semiMajorAxis: [6433531050866.602, 0, 0, 0],
		eccentricity: [0.1957748188564788, 0, 0, 0],
		meanLon: [9.63366795275941, 2.2278486427835054, 0, 0],
		inclination: [0.4923295578119718, 0, 0, 0],
		ascendingNodeLon: [2.1257637339034994, 0, 0, 0],
		periapsisLon: [6.330058331275721, 0, 0, 0]
	},
	makemake: {
		semiMajorAxis: [6808301094571.345, 0, 0, 0],
		eccentricity: [0.1604249998705246, 0, 0, 0],
		meanLon: [8.993863839985694, 2.046452426537862, 0, 0],
		inclination: [0.5067093316791631, 0, 0, 0],
		ascendingNodeLon: [1.3835036678328163, 0, 0, 0],
		periapsisLon: [6.568448083737754, 0, 0, 0]
	},
	sedna: {
		semiMajorAxis: [82210355196446.73, 0, 0, 0],
		eccentricity: [0.8612979286929314, 0, 0, 0],
		meanLon: [7.912841105614632, 0.048771936225330555, 0, 0],
		inclination: [0.20814651761955802, 0, 0, 0],
		ascendingNodeLon: [2.5216293860468753, 0, 0, 0],
		periapsisLon: [7.949774840138166, 0, 0, 0]
	}
};

const preferredNumbers = [1, 1.2, 1.5, 1.75, 2, 2.5, 3, 4, 5, 6, 7, 8, 9, 10];
const roundToPreferredNumber = (num, mode = 0, prefList = preferredNumbers) => {
	if (num == 0) return num;

	const exponent = 10 ** Math.floor(Math.log10(num));
	const mantissa = num / exponent;

	if (mode == 1) {
		return exponent * prefList.find(el => el >= mantissa);
	}
	if (mode == -1) {
		return mantissa == 1
			? num
			: exponent * prefList[prefList.findIndex(el => el > mantissa) - 1];
	}

	const differences = prefList.map(el => Math.abs(el - mantissa));
	return exponent * prefList[differences.indexOf(Math.min(...differences))];
};

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;
};

const getLightTravelTime = (r1, r2) =>
	Math.hypot(r2.x - r1.x, r2.y - r1.y, r2.z - r1.z) / LIGHTSECOND;

const getNormalisedVector = vec => {
	let output = { xyz: vec };
	output.abs = Math.hypot(...vec);
	output.unit = vec.map(n => n / output.abs);

	return output;
};

const cross = (a, b) => [
	a[1] * b[2] - a[2] * b[1],
	a[2] * b[0] - a[0] * b[2],
	a[0] * b[1] - a[1] * b[0]
];

// Using Izzo’s algorithm, which i don’t fully understand myself: https://arxiv.org/pdf/1403.2705.pdf
const lambert = {
	flightTime: (x, M, λ, λ2) => {
		const oneMinusX2 = 1 - x ** 2;
		const y = Math.sqrt(1 - λ2 * oneMinusX2);
		const ψ = Math.acos(x * y + λ * oneMinusX2);

		if (Math.abs(x - 1) < 0.01) {
			const η = y - λ * x;
			const S1 = (1 - λ - x * η) / 2;

			// calculate hypergeometric function
			const ε = 1e-7;
			let Sj = 1,
				Cj = 1,
				idx = 0;
			while (Math.abs(Cj) > ε && idx < 15) {
				Cj =
					(((Cj * (idx + 3) * (idx + 1)) / (idx + 2.5)) * S1) /
					(idx + 1);
				Sj += Cj;
				idx++;
			}

			return 2 * η * (λ + (η ** 2 * Sj) / 3);
		}

		return (
			((ψ + M * Math.PI) / Math.sqrt(Math.abs(oneMinusX2)) - x + λ * y) /
			oneMinusX2
		);
	},
	dTdx: (T, x, λ, λ2) => {
		const y2 = 1 - λ2 * (1 - x ** 2);
		const y = Math.sqrt(y2);

		let dT = [null, 0, 0, 0];
		dT[1] = (3 * T * x - 2 + (2 * λ2 * λ * x) / y) / (1 - x ** 2);
		dT[2] =
			(3 * T + 5 * x * dT[1] + 2 * (1 - λ2) * ((λ2 * λ) / (y2 * y))) /
			(1 - x ** 2);
		dT[3] =
			(7 * x * dT[2] +
				8 * dT[1] -
				(6 * λ2 * λ2 * λ * (1 - λ2) * x) / (y2 * y2 * y)) /
			(1 - x ** 2);

		return dT;
	},
	householder: (T, x0, M, λ, λ2) => {
		const ε = 1e-7;
		let err = 1;

		let idx = 0;
		while (err > ε && idx < 15) {
			const tof = lambert.flightTime(x0, M, λ, λ2);
			const dtof = lambert.dTdx(tof, x0, λ, λ2);
			const Δ = tof - T;

			const xNew =
				x0 -
				(Δ * (dtof[1] ** 2 - (Δ * dtof[2]) / 2)) /
					(dtof[1] * (dtof[1] ** 2 - Δ * dtof[2]) +
						(dtof[3] * Δ ** 2) / 6);
			err = Math.abs(x0 - xNew);
			x0 = xNew;

			idx++;
		}

		return x0;
	},
	findX: (λ, λ2, T) => {
		const ε = 1e-7;

		let maxM = Math.floor(T / Math.PI);
		const T00 = Math.acos(λ) + λ * Math.sqrt(1 - λ2);
		let T0 = T00 * maxM + Math.PI;
		let T1 = ((1 - λ2 * λ) * 2) / 3;

		if (T < T0 && maxM > 0 && maxM <= MAX_MULTIREVOLUTIONS) {
			let minT = T0,
				xOld = 0,
				xNew = 0,
				err = 1,
				idx = 0;
			while (err > ε && idx < 15) {
				const dT = lambert.dTdx(minT, xOld, λ, λ2);
				if (dT[1] != 0) {
					xNew =
						xOld -
						(dT[1] * dT[2]) / (dT[2] ** 2 - (dT[1] * dT[3]) / 2);
				}
				err = Math.abs(xNew - xOld);
				minT = lambert.flightTime(xNew, maxM, λ, λ2);
				xOld = xNew;
				idx++;
			}

			if (minT > T) maxM -= 1;
		}

		maxM = Math.min(maxM, MAX_MULTIREVOLUTIONS);

		let solutions = [];

		const x0 =
			T >= T00
				? (T00 / T) ** (2 / 3) - 1
				: T < T1
				? 2.5 * ((T1 * (T1 - T)) / (T * (1 - λ2 * λ2 * λ))) + 1
				: (T00 / T) ** Math.log2(T1 / T00);
		solutions.push(lambert.householder(T, x0, 0, λ, λ2));

		while (maxM > 0) {
			const x0lTemp = (((maxM + 1) * Math.PI) / (8 * T)) ** (2 / 3);
			const x0rTemp = ((8 * T) / (maxM * Math.PI)) ** (2 / 3);

			const x0l = (x0lTemp - 1) / (x0lTemp + 1);
			const x0r = (x0rTemp - 1) / (x0rTemp + 1);

			solutions.push(
				lambert.householder(T, x0l, maxM, λ, λ2),
				lambert.householder(T, x0r, maxM, λ, λ2)
			);

			maxM--;
		}

		return solutions;
	},
	solve: (r1Raw, r2Raw, time, μ = GRAV_PARAM) => {
		const r1 = getNormalisedVector(r1Raw.xyz),
			r2 = getNormalisedVector(r2Raw.xyz);

		if (time <= 0) return null;

		const c = getNormalisedVector(
			r2.xyz.map((_, i) => r2.xyz[i] - r1.xyz[i])
		);
		const s = (r1.abs + r2.abs + c.abs) / 2;

		const λ2 = 1 - c.abs / s;
		let λ = Math.sqrt(λ2);
		let ît1, ît2;
		const îh = getNormalisedVector(cross(r1.unit, r2.unit)).unit;

		if (r1.xyz[0] * r2.xyz[1] - r1.xyz[1] * r2.xyz[0] < 0) {
			λ *= -1;
			ît1 = cross(r1.unit, îh);
			ît2 = cross(r2.unit, îh);
		} else {
			ît1 = cross(îh, r1.unit);
			ît2 = cross(îh, r2.unit);
		}

		const T = time * Math.sqrt((2 * μ) / s ** 3);
		const γ = Math.sqrt((μ * s) / 2);
		const ρ = (r1.abs - r2.abs) / c.abs;
		const σ = Math.sqrt(1 - ρ ** 2);

		const solutions = lambert.findX(λ, λ2, T).map(x => {
			const y = Math.sqrt(1 - λ2 * (1 - x ** 2));
			const λy_x = λ * y - x,
				ρλyx = ρ * (λ * y + x),
				γσyλx = γ * σ * (y + λ * x);

			const absVr1 = (γ * (λy_x - ρλyx)) / r1.abs;
			const absVr2 = (-γ * (λy_x + ρλyx)) / r2.abs;
			const absVt1 = γσyλx / r1.abs;
			const absVt2 = γσyλx / r2.abs;

			const vr1 = r1.unit.map(n => n * absVr1);
			const vr2 = r2.unit.map(n => n * absVr2);
			const vt1 = ît1.map(n => n * absVt1);
			const vt2 = ît2.map(n => n * absVt2);

			return [
				vr1.map((_, i) => vr1[i] + vt1[i]),
				vr2.map((_, i) => vr2[i] + vt2[i])
			];
		});

		return solutions;
	}
};

// Computes the light travel time between two planets for a date and the week around it
const getCommTimeBetween = (planet1, planet2, date) => {
	const dateOffset = getJ2000Offset(date);
	const week = [-3, -2, -1, 0, 1, 2, 3].map(x => x + dateOffset);

	let times = {
		daily: null,
		weekly: null
	};
	times.weekly = week.map(d => {
		const r1 = getElements(planet1, d).r;
		const r2 = getElements(planet2, d).r;

		return getLightTravelTime(r1, r2);
	});
	times.daily = times.weekly[3];

	return times;
};

const findLeastDeltaV = (planet1, planet2, departureJ2000, arrivalJ2000) => {
	const r1 = getElements(planet1, departureJ2000);
	const r2 = getElements(planet2, arrivalJ2000);
	const tripTime = (arrivalJ2000 - departureJ2000) * 86400;

	const getDeltaV = solution =>
		Math.hypot(...solution[0].map((el, idx) => r1.velocity.xyz[idx] - el)) +
		Math.hypot(...solution[1].map((el, idx) => r2.velocity.xyz[idx] - el));

	const solutions = lambert.solve(r1.r, r2.r, tripTime);

	return solutions == null ? Infinity : Math.min(...solutions.map(getDeltaV));
};

// WIP - Solar System orrery model

// 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;
};

// Presentational stuff

const tr = {
	lang: {
		current: "en",
		pl: new Intl.PluralRules("en-GB")
	},
	en: {
		secondUnits: {
			short: ["d ", "h ", "m ", "s"],
			long: [
				{
					one: " day, ",
					other: " days, "
				},
				{
					one: " hour, ",
					other: " hours, "
				},
				{
					one: " minute, ",
					other: " minutes, "
				},
				{
					one: " second",
					other: " seconds"
				}
			]
		},
		porkchop: {
			transitTime: "Days spent in transit",
			departure: "Departure date",
			deltaV: "Required Δv"
		},
			planetNames: {
				hermes: "Hermes (Mercury)",
				venus: "Venus",
				earth: "Earth",
				mars: "Mars",
				ceres: "Ceres",
				jupiter: "Jupiter",
				saturn: "Saturn",
				ouranos: "Ouranos (Uranus)",
				neptune: "Neptune",
				pluto: "Pluto",
				persephone: "Persephone (Eris)",
				haumea: "Haumea",
				makemake: "Makemake",
				sedna: "Sedna"
			}
	},
	nl: {
		secondUnits: {
			short: ["d ", "u ", "m ", "s"],
			long: [
				{
					one: " dag, ",
					other: " dagen, "
				},
				{
					one: " uur, ",
					other: " uur, "
				},
				{
					one: " minuten, ",
					other: " minuten, "
				},
				{
					one: " seconde",
					other: " seconden"
				}
			]
		},
		porkchop: {
			transitTime: "Dagen onderweg",
			departure: "Vertrekdatum",
			deltaV: "Δv nodig"
		},
			planetNames: {
				hermes: "Hermes (Mercurius)",
				venus: "Venus",
				earth: "Aarde",
				mars: "Mars",
				ceres: "Ceres",
				jupiter: "Jupiter",
				saturn: "Saturnus",
				ouranos: "Ouranos (Uranus)",
				neptune: "Neptunus",
				pluto: "Pluto",
				persephone: "Persephone (Eris)",
				haumea: "Haumea",
				makemake: "Makemake",
				sedna: "Sedna"
			}
	}
};

const formatSeconds = (secs, format = "short") => {
	const roundSecs = Math.round(secs);

	const daysHoursMinsSecs = [
		Math.floor(roundSecs / 86400),
		Math.floor(roundSecs / 3600) % 24,
		Math.floor(roundSecs / 60) % 60,
		roundSecs % 60
	];

	format = format in tr[tr.lang.current].secondUnits ? format : "short";
	const outputStrings = daysHoursMinsSecs.map(
		(x, idx) =>
			`${x}${
				format == "short"
					? tr[tr.lang.current].secondUnits.short[idx]
					: tr[tr.lang.current].secondUnits[format][idx][
							tr.lang.pl.select(x)
					  ]
			}`
	);

	return outputStrings
		.slice(daysHoursMinsSecs.findIndex(x => x > 0))
		.join("");
};

// Page stuff

const showCommTime = () => {
	const times = getCommTimeBetween(
		planets[$("#commtimecalc-planet1").value],
		planets[$("#commtimecalc-planet2").value],
		new Date($("#commtimecalc-date").value + "Z")
	);

	$("#commtimecalc-daily-result").innerHTML = formatSeconds(
		times.daily,
		"long"
	);

	const weekdays = $$("#commtimecalc-weekly-table td b");

	times.weekly.forEach((time, idx) => {
		weekdays[idx].innerHTML = formatSeconds(time, "short");
	});
};

let porkchopX = null;
let porkchopY = null;
let porkchopLegend = null;
let porkchopChart = null;

const showPorkchop = () => {
	if ($("#porkchop-planet1").value == $("#porkchop-planet2").value) {
		return;
	}

	const svg = $("#porkchop-plot");
	const svgWidth = svg.clientWidth;
	const svgHeight = svg.clientHeight;
	const chart = d3.select("#porkchop-plot");

	const p1 = planets[$("#porkchop-planet1").value];
	const p2 = planets[$("#porkchop-planet2").value];
	const departureYear = $("#porkchop-departure").value;

	const hohmannTransferTime =
		(Math.PI *
			Math.sqrt(
				(p1.semiMajorAxis[0] + p2.semiMajorAxis[0]) ** 3 /
					(8 * GRAV_PARAM)
			)) /
		86400;
	const travelTimeRange = [
		roundToPreferredNumber(hohmannTransferTime / 2.5),
		roundToPreferredNumber(hohmannTransferTime * 2.5)
	];

	const xGrid = d3.range(
		getJ2000Offset(new Date(`${departureYear}-01-01T12:00Z`)),
		getJ2000Offset(new Date(`${departureYear}-12-31T12:00Z`)),
		1
	);
	const yGrid = d3.range(
		...travelTimeRange,
		roundToPreferredNumber(
			(travelTimeRange[1] - travelTimeRange[0]) / 150,
			0,
			[1, 2, 5, 10]
		)
	);
	const grid = yGrid.flatMap(y =>
		xGrid.map(x => findLeastDeltaV(p1, p2, x, x + y))
	);

	const xAxis = d3
		.axisBottom(
			d3
				.scaleUtc()
				.domain([
					new Date("2020-01-01T12:00Z"),
					new Date("2020-12-31T12:00Z")
				])
				.range([70, svgWidth - 80])
		)
		.tickFormat(d3.utcFormat("%b"));
	const yAxis = d3.axisLeft(
		d3
			.scaleLinear()
			.domain([yGrid[0], yGrid.at(-1)])
			.range([svgHeight - 60, 20])
	);

	const minDeltaV = Math.min(
		...grid.filter(x => x !== null && x !== Infinity && x > 1)
	);
	const scaleMin = roundToPreferredNumber(minDeltaV, -1);
	const scaleMax = roundToPreferredNumber(minDeltaV * 2);
	const scale = d3.scaleSequential([scaleMin, scaleMax], d3.interpolateTurbo);

	const legendAxis = d3
		.axisRight(
			d3.scaleLinear().domain([scaleMax, scaleMin]).range([25, 175])
		)
		.ticks(5)
		.tickFormat(n =>
			d3
				.format("~s")(n)
				.replace(/([0-9])([a-zA-Zµ])?$/, "$1 $2m/s")
		);

	const rectWidth = (svgWidth - 150) / xGrid.length;
	const rectHeight = (svgHeight - 80) / yGrid.length;

	if (porkchopChart === null) {
		porkchopChart = chart.append("g");
	}

	porkchopChart
		.selectAll("rect")
		.data(grid)
		.join("rect")
		.attr("fill", d =>
			d > scaleMax || !(d > 1) ? "transparent" : scale(d)
		)
		.attr("width", rectWidth)
		.attr("height", rectHeight)
		.attr("x", (_, idx) => 60 + (idx % xGrid.length) * rectWidth)
		.attr(
			"y",
			(_, idx) =>
				svgHeight -
				60 -
				(1 + Math.floor(idx / xGrid.length)) * rectHeight
		);

	if (porkchopX === null) {
		porkchopX = chart
			.append("g")
			.attr("transform", `translate(0,${svgHeight - 40})`)
			.call(xAxis)
			.attr("font-size", "12");

		chart
			.append("text")
			.attr("class", "x-axis-label")
			.attr("font-size", "12")
			.attr("x", svgWidth / 2 - 5)
			.attr("y", svgHeight - 10)
			.text(tr[tr.lang.current].porkchop.departure);
	} else {
		porkchopX.transition().duration(1000).call(xAxis);
	}

	if (porkchopY === null) {
		porkchopY = chart
			.append("g")
			.attr("transform", `translate(50, 0)`)
			.call(yAxis)
			.attr("font-size", "12");

		chart
			.append("text")
			.attr("class", "y-axis-label")
			.attr("font-size", "12")
			.attr("x", 10)
			.attr("y", svgHeight / 2 - 20)
			.text(tr[tr.lang.current].porkchop.transitTime);
	} else {
		porkchopY.transition().duration(1000).call(yAxis);
	}

	if (porkchopLegend === null) {
		const legendCanvas = document.createElement("canvas");
		legendCanvas.width = 1;
		legendCanvas.height = 256;
		const legendCtx = legendCanvas.getContext("2d");
		for (let i = 0; i < 256; i++) {
			legendCtx.fillStyle = d3.interpolateTurbo(1 - i / 255);
			legendCtx.fillRect(0, i, 1, 1);
		}

		chart
			.append("image")
			.attr("x", svgWidth - 75)
			.attr("y", 25)
			.attr("width", 15)
			.attr("height", 150)
			.attr("preserveAspectRatio", "none")
			.attr("xlink:href", legendCanvas.toDataURL());

		porkchopLegend = chart
			.append("g")
			.attr("transform", `translate(${svgWidth - 60}, 0)`)
			.call(legendAxis)
			.attr("font-size", "12");

		chart
			.append("text")
			.attr("class", "legend-label")
			.attr("font-size", "12")
			.attr("x", svgWidth - 75)
			.attr("y", 15)
			.text(tr[tr.lang.current].porkchop.deltaV);
	} else {
		porkchopLegend.transition().duration(1000).call(legendAxis);
	}
};

let orrery = {
	chart: null,
	graticule: null,
	orbits: null,
	paths: null,
	circles: null
};
const planetColours = {
	hermes: "#876",
	venus: "#c70",
	earth: "#295",
	mars: "#c32",
	jupiter: "#c78",
	saturn: "#da4",
	ouranos: "#2bb",
	neptune: "#26b",
	pluto: "#72a",
	persephone: "#b18",
	ceres: "#cdb",
	haumea: "#bdc",
	makemake: "#bcd",
	sedna: "#cbd",
};

const showOrrery = settings => {
	const svg = $("#orrery-plot");
	const svgHeight = svg.clientHeight;

	const detailLevel = 120;
	// Zoom will be from 3e11 to 3e13 metres
	const zoom = settings?.zoom ?? 1.5e13;
	const perspective = settings?.perspective ?? 0;

	const today = settings?.today ?? new Date("2558-05-26T12:00Z");

	const orbits = [
		"ceres",
		"haumea",
		"makemake",
		"sedna",
		"hermes",
		"venus",
		"earth",
		"mars",
		"jupiter",
		"saturn",
		"ouranos",
		"neptune",
		"pluto",
		"persephone"
	].map(el => ({
		planet: el,
		orbit: getOrbitTrail(planets[el], today, detailLevel)
	}));

	const scale = d3.scaleLinear([-zoom / 2, zoom / 2], [0, svgHeight]);

	// takes in XYZ coördinates and spits out nice XY ones
	const project = xyz => {
		return {
			x: scale(xyz[0]),
			y: scale(
				-xyz[1] * Math.cos(perspective) - xyz[2] * Math.sin(perspective)
			)
		};
	};

	if (orrery.chart === null) {
		orrery.chart = d3.select("#orrery-chart");

		orrery.chart
			.append("circle")
			.attr("cx", project([0, 0, 0]).x)
			.attr("cy", project([0, 0, 0]).y)
			.attr("r", 7)
			.attr("fill", "#ff8")
			.attr("stroke", "#880")
			.attr("stroke-width", "1px")
			.attr("z-index", "0");
	}

	const line = d3.line(
		d => project(d.xyz).x,
		d => project(d.xyz).y
	);

	orrery.orbits = d3
		.select("#orrery-orbits")
		.selectAll("g")
		.data(orbits)
		.join("g")
		.attr("color", d => planetColours[d.planet]);

	orrery.paths = orrery.orbits
		.selectAll("path")
		.data(d =>
			d.orbit
				.map((el, idx, array) => [el, array[idx + 1] || null])
				.slice(0, detailLevel - 1)
		)
		.join("path")
		.attr("d", d => line(d))
		.attr("opacity", (d, idx) => 1 - idx / detailLevel)
		.attr("stroke", "currentColor")
		.attr("stroke-width", "2px")
		.attr("fill", "none");

	orrery.circles = d3
		.select("#orrery-circles")
		.selectAll("circle")
		.data(orbits)
		.join("circle")
		.attr("cx", d => Math.round(project(d.orbit[0].xyz).x))
		.attr("cy", d => Math.round(project(d.orbit[0].xyz).y))
		.attr("r", 5.5)
		.attr("fill", d => planetColours[d.planet])
		.html(d => `<title>${tr[tr.lang.current].planetNames[d.planet]}</title>`);
};

let shiftOrrDate, todayOrrDate, resetOrrDate;
documentReady(() => {
	showOrrery();

	const reloadOrrery = () => {
		showOrrery({
			today: new Date($("#orrery-date").value + "Z"),
			perspective:
				(Math.PI / 180) * (+$("#orrery-perspective-slider").value - 90),
			zoom: 3 * 10 ** +$("#orrery-zoom-slider").value
		});
	};

	const orreryDate = $("#orrery-date");
	const orreryDinkus = $("#orrery-perspective-dinkus");
	const zoomDinkus = $("#orrery-zoom-dinkus");

	shiftOrrDate = days => {
		console.log(orreryDate.value);
		let shifted = new Date(`${orreryDate.value}Z`);
		shifted.setDate(shifted.getDate() + days);
		orreryDate.value = shifted.toISOString().slice(0, 16);
		console.log(orreryDate.value);
		reloadOrrery();
	};
	todayOrrDate = date => {
		orreryDate.value = new Date().toISOString().slice(0, 16);
		reloadOrrery();
	};
	resetOrrDate = date => {
		orreryDate.value = "2558-05-26T12:00";
		reloadOrrery();
	};

	$("#orrery-perspective-slider").on("input", evt => {
		orreryDinkus.innerHTML = `${evt.target.value}°`;
		orreryDinkus.style.setProperty(
			"--dinkus-position",
			`${evt.target.value}`
		);
	});

	$("#orrery-zoom-slider").on("input", evt => {
		const rawZoom = 3 * 10 ** +evt.target.value;
		const roundedZoom = roundToPreferredNumber(rawZoom);

		zoomDinkus.innerHTML = d3
			.format("~s")(roundedZoom)
			.replace(/([0-9])([a-zA-Zµ])?$/, "$1 $2m");
		zoomDinkus.style.setProperty(
			"--dinkus-position",
			`${evt.target.value}`
		);
	});

	$$("#orrery-tool input").forEach(el => {
		el.on("input", reloadOrrery);
	});
});