// Grey line — the day/night terminator, and the twilight band around it where
// HF propagation briefly does things it does at no other time.
//
// The maths is the standard low-precision solar position (Meeus, chapter 25, as
// used by NOAA's solar calculator): good to well under a degree, which is far
// finer than the band being drawn. No ephemeris, no network, no dependency —
// this has to work in a field shack with no internet.
//
// Everything here is pure so it can be tested against known positions; the
// drawing lives in the map component.

const DEG = Math.PI / 180;

// julianDay for a date, including the fraction of the day.
export function julianDay(d: Date): number {
	return d.getTime() / 86400000 + 2440587.5;
}

// Sun's declination and Greenwich hour angle, both in degrees.
//
// Returned together because they are computed from the same intermediates and
// every caller needs both — splitting them into two exported functions would
// double the work at each call for the sake of a tidier signature.
export function sunPosition(date: Date): { dec: number; gha: number } {
	const jd = julianDay(date);
	const n = jd - 2451545.0; // days from J2000.0

	// Mean longitude and mean anomaly of the Sun.
	const L = (280.460 + 0.9856474 * n) % 360;
	const g = ((357.528 + 0.9856003 * n) % 360) * DEG;

	// Ecliptic longitude: the equation of centre applied to the mean longitude.
	const lambda = (L + 1.915 * Math.sin(g) + 0.020 * Math.sin(2 * g)) * DEG;

	// Obliquity of the ecliptic, slowly decreasing.
	const eps = (23.439 - 0.0000004 * n) * DEG;

	const dec = Math.asin(Math.sin(eps) * Math.sin(lambda)) / DEG;

	// Right ascension, kept in the same quadrant as the ecliptic longitude —
	// atan2 alone lands 180° out for half the year, which mirrors the whole
	// terminator about the prime meridian. Six months of a plausible-looking
	// wrong answer is exactly the kind of bug that survives review.
	let ra = Math.atan2(Math.cos(eps) * Math.sin(lambda), Math.cos(lambda)) / DEG;
	if (ra < 0) ra += 360;

	// Greenwich mean sidereal time, in degrees.
	const gmst = (280.46061837 + 360.98564736629 * n) % 360;

	let gha = gmst - ra;
	gha = ((gha % 360) + 360) % 360;
	return { dec, gha };
}

// terminatorLatitude returns the latitude, at a given longitude, where the Sun
// sits at `altitude` degrees above (or below) the horizon.
//
// altitude 0 is the geometric terminator; -6 is civil twilight, the outer edge
// of the grey line as operators use the term.
//
// Returns null where no such latitude exists at that longitude — inside the
// polar day or polar night, where the Sun never crosses that altitude at all.
export function terminatorLatitude(lonDeg: number, dec: number, gha: number, altitude = 0): number | null {
	const ha = (gha + lonDeg) * DEG; // local hour angle
	const d = dec * DEG;
	const h = altitude * DEG;

	// From sin(alt) = sin(dec)sin(lat) + cos(dec)cos(lat)cos(ha), solved for lat:
	//   A·sin(lat) + B·cos(lat) = sin(alt),  A = sin(dec), B = cos(dec)cos(ha)
	//
	// That has TWO solutions per meridian, and picking the wrong one is not a
	// small error — it mirrors the answer about the equator. New York at 02:00
	// local in December came out as daylight because the branch chosen put the
	// terminator at +63° where it belongs at −63°: a shaded map that looks
	// entirely plausible and is exactly wrong.
	//
	// So anchor on the geometric terminator, where the answer is unambiguous and
	// classical, and take whichever branch is nearest it. At altitude 0 that
	// reproduces it exactly; for the twilight band it follows it by continuity,
	// which is the physical requirement — the twilight curve is a small offset
	// from the terminator, never on the other side of the planet.
	const A = Math.sin(d);
	const B = Math.cos(d) * Math.cos(ha);
	const R = Math.hypot(A, B);
	if (R === 0) return null;
	const s = Math.sin(h) / R;
	if (s < -1 || s > 1) return null; // the Sun never reaches this altitude here

	// The classical terminator. Division by zero is deliberate and correct here:
	// at an equinox tan(dec) → 0 and the terminator becomes a meridian, which is
	// exactly what atan(±Infinity) = ±90° expresses.
	const t = Math.tan(d);
	if (t === 0 && Math.cos(ha) === 0) return null; // degenerate: no crossing to name
	const lat0 = Math.atan(-Math.cos(ha) / t) / DEG;

	const phi = Math.atan2(B, A); // phase of the A·sin(lat) + B·cos(lat) sum
	const fold = (x: number) => {
		let v = ((x + 180) % 360 + 360) % 360 - 180;
		if (v > 90) v = 180 - v;
		if (v < -90) v = -180 - v;
		return v;
	};
	const c1 = fold((Math.asin(s) - phi) / DEG);
	const c2 = fold((Math.PI - Math.asin(s) - phi) / DEG);
	return Math.abs(c1 - lat0) <= Math.abs(c2 - lat0) ? c1 : c2;
}

// nightPolygon returns a ring covering the part of the world where the Sun is
// below `altitude`, ready for L.polygon (as [lat, lon] pairs).
//
// The ring runs along the terminator from west to east and closes along
// whichever pole is in darkness — which pole that is flips with the season, and
// getting it wrong shades the lit half of the planet.
export function nightPolygon(date: Date, altitude = 0, stepDeg = 2): [number, number][] {
	const { dec, gha } = sunPosition(date);
	const ring: [number, number][] = [];

	for (let lon = -180; lon <= 180; lon += stepDeg) {
		const lat = terminatorLatitude(lon, dec, gha, altitude);
		// Inside a polar day/night there is no crossing; clamp to the pole so the
		// ring stays closed rather than tearing open across the map.
		ring.push([lat === null ? (dec > 0 ? -90 : 90) : lat, lon]);
	}

	// Northern summer (positive declination) → the SOUTH pole is dark.
	const pole = dec > 0 ? -90 : 90;
	ring.push([pole, 180]);
	ring.push([pole, -180]);
	return ring;
}