weather-station/server/celestial.mjs

const RAD = Math.PI / 180;
const DEG = 180 / Math.PI;

function jdn(date) {
	return (date instanceof Date ? date : new Date(date)) / 86400000 + 2440587.5;
}

function jdnToDate(jd) {
	return new Date((jd - 2440587.5) * 86400000);
}

function sunPosition(lat, lon, date = new Date()) {
	const d = jdn(date) - 2451545.0;
	const L = (280.46 + 0.9856474 * d) % 360;
	const g = (357.528 + 0.9856003 * d) % 360;
	const lam = L + 1.915 * Math.sin(g * RAD) + 0.02 * Math.sin(2 * g * RAD);
	const eps = 23.439 - 0.0000004 * d;
	const sinL = Math.sin(lam * RAD);
	const ra = Math.atan2(Math.cos(eps * RAD) * sinL, Math.cos(lam * RAD)) * DEG;
	const dec = Math.asin(Math.sin(eps * RAD) * sinL) * DEG;
	const UT = date.getUTCHours() + date.getUTCMinutes() / 60 + date.getUTCSeconds() / 3600;
	const GMST = (6.697375 + 0.0657098242 * d + UT) % 24;
	const LMST = (GMST + lon / 15) % 24;
	const ha = LMST * 15 - ra;
	const elev = Math.asin(
		Math.sin(lat * RAD) * Math.sin(dec * RAD) +
		Math.cos(lat * RAD) * Math.cos(dec * RAD) * Math.cos(ha * RAD)
	) * DEG;
	const az = Math.atan2(
		-Math.sin(ha * RAD),
		Math.tan(dec * RAD) * Math.cos(lat * RAD) - Math.sin(lat * RAD) * Math.cos(ha * RAD)
	) * DEG;
	return {
		sun_elevation: Math.round(elev * 10) / 10,
		sun_azimuth: Math.round(((az + 360) % 360) * 10) / 10,
	};
}

function sunriseSunset(lat, lon, date = new Date()) {
	const d = Math.floor(jdn(date)) - 2451545;
	const noon = 2451545 + 0.0009 + ((-lon) / 360) + Math.round(d - (-lon) / 360);
	const M = (357.5291 + 0.98560028 * (noon - 2451545)) % 360;
	const C = 1.9148 * Math.sin(M * RAD) + 0.02 * Math.sin(2 * M * RAD) + 0.0003 * Math.sin(3 * M * RAD);
	const lam = (M + C + 180 + 102.9372) % 360;
	const jnoon = noon + 0.0053 * Math.sin(M * RAD) - 0.0069 * Math.sin(2 * lam * RAD);
	const dec = Math.asin(Math.sin(23.4397 * RAD) * Math.sin(lam * RAD)) * DEG;
	const cosH = (Math.sin(-0.8333 * RAD) - Math.sin(lat * RAD) * Math.sin(dec * RAD)) / (Math.cos(lat * RAD) * Math.cos(dec * RAD));

	if(Math.abs(cosH) > 1) {
		return {sunrise: null, sunset: null, daytime: cosH < -1};
	}

	const H = Math.acos(cosH) * DEG;
	const rise = jdnToDate(jnoon - H / 360);
	const set = jdnToDate(jnoon + H / 360);
	const noon_d = jdnToDate(jnoon);

	// Round to nearest second, drop milliseconds
	const fmt = d => new Date(Math.round(d.getTime() / 1000) * 1000).toISOString();

	return {
		sunrise: fmt(rise),
		sunset: fmt(set),
		daylight: Math.round((set - rise) / 36000) / 100,
		daytime: cosH < -1
	};
}

// Smooth sunspot number approximation from solar flux F10.7
// SSN ≈ 1.61 * F10.7 - 63.7  (linear regression, valid for F10.7 > 70)
export function ssnFromFlux(f107) {
	if(!f107) return null;
	return Math.max(0, Math.round(1.61 * f107 - 63.7));
}

export function sunActivityLabel(f107) {
	if(!f107) return null;
	if(f107 < 80) return 'Very Low';
	if(f107 < 100) return 'Low';
	if(f107 < 150) return 'Moderate';
	if(f107 < 200) return 'High';
	return 'Very High';
}

function moonPosition(lat, lon, date = new Date()) {
	const d = jdn(date) - 2451545;
	const L = (218.316 + 13.176396 * d) % 360;
	const M = (134.963 + 13.064993 * d) % 360;
	const F = (93.272 + 13.229350 * d) % 360;
	const lon_ = L + 6.289 * Math.sin(M * RAD);
	const b = 5.128 * Math.sin(F * RAD);
	const dec = Math.asin(
		Math.sin(b * RAD) * Math.cos(23.4397 * RAD) +
		Math.cos(b * RAD) * Math.sin(23.4397 * RAD) * Math.sin(lon_ * RAD)
	) * DEG;
	const ra = Math.atan2(
		Math.sin(lon_ * RAD) * Math.cos(23.4397 * RAD) - Math.tan(b * RAD) * Math.sin(23.4397 * RAD),
		Math.cos(lon_ * RAD)
	) * DEG;
	const UT = date.getUTCHours() + date.getUTCMinutes() / 60 + date.getUTCSeconds() / 3600;
	const GMST = (6.697375 + 0.0657098242 * d + UT) % 24;
	const LMST = (GMST + lon / 15) % 24;
	const ha = LMST * 15 - ra;
	const elev = Math.asin(
		Math.sin(lat * RAD) * Math.sin(dec * RAD) +
		Math.cos(lat * RAD) * Math.cos(dec * RAD) * Math.cos(ha * RAD)
	) * DEG;
	const az = Math.atan2(
		-Math.sin(ha * RAD),
		Math.tan(dec * RAD) * Math.cos(lat * RAD) - Math.sin(lat * RAD) * Math.cos(ha * RAD)
	) * DEG;
	return {
		moon_elevation: Math.round(elev * 10) / 10,
		moon_azimuth: Math.round(((az + 360) % 360) * 10) / 10,
	};
}

function sunMoonDistance(date = new Date()) {
	const d = jdn(date) - 2451545;
	// Sun ecliptic longitude
	const Ls = (280.46 + 0.9856474 * d) % 360;
	const gs = (357.528 + 0.9856003 * d) % 360;
	const sunLon = Ls + 1.915 * Math.sin(gs * RAD) + 0.02 * Math.sin(2 * gs * RAD);
	// Moon ecliptic longitude
	const Lm = (218.316 + 13.176396 * d) % 360;
	const Mm = (134.963 + 13.064993 * d) % 360;
	const Fm = (93.272 + 13.229350 * d) % 360;
	const moonLon = Lm + 6.289 * Math.sin(Mm * RAD);
	const moonLat = 5.128 * Math.sin(Fm * RAD);
	// Angular separation
	const dLon = (moonLon - sunLon) * RAD;
	const sep = Math.acos(
		Math.cos(moonLat * RAD) * Math.cos(dLon)
	) * DEG;
	return { sun_moon_separation: Math.round(sep * 10) / 10 };
}

function moonPhase(date = new Date()) {
	const jd = jdn(date);
	const cycle = 29.53058867;
	const known = 2451550.1;
	const phase = ((jd - known) % cycle + cycle) % cycle;
	// More accurate illumination using proper phase angle
	const illum = Math.round((1 - Math.cos(phase / cycle * 2 * Math.PI)) / 2 * 1000) / 10;
	let name;
	if(phase < 1.85) name = 'New Moon';
	else if(phase < 7.38) name = 'Waxing Crescent';
	else if(phase < 9.22) name = 'First Quarter';
	else if(phase < 14.76) name = 'Waxing Gibbous';
	else if(phase < 16.61) name = 'Full Moon';
	else if(phase < 22.15) name = 'Waning Gibbous';
	else if(phase < 23.99) name = 'Last Quarter';
	else name = 'Waning Crescent';
	return { moon_phase: name, moon_illumination: illum };
}

function nextMoonEvents(date = new Date()) {
	const cycle = 29.53058867;
	const jd = jdn(date);
	const known = 2451550.1;
	const phase = ((jd - known) % cycle + cycle) % cycle;
	const toNew = (cycle - phase) % cycle;
	const toFull = phase < 14.76 ? 14.76 - phase : cycle - phase + 14.76;
	return {
		moon_new: jdnToDate(jd + toNew).toISOString(),
		moon_full: jdnToDate(jd + toFull).toISOString(),
	};
}

function moonriseMoonset(lat, lon, date = new Date()) {
	const base = new Date(date);
	base.setUTCHours(0, 0, 0, 0);
	let prev = null;
	let moonrise = null;
	let moonset = null;

	// Scan 48h in 10min steps — covers edge cases where moon rises/sets next UTC day
	for(let m = 0; m <= 2880; m += 10) {
		const t = new Date(base.getTime() + m * 60000);
		const d = jdn(t) - 2451545;
		const L = (218.316 + 13.176396 * d) % 360;
		const M = (134.963 + 13.064993 * d) % 360;
		const F = (93.272 + 13.229350 * d) % 360;
		const lon_ = L + 6.289 * Math.sin(M * RAD);
		const b = 5.128 * Math.sin(F * RAD);
		const dec = Math.asin(
			Math.sin(b * RAD) * Math.cos(23.4397 * RAD) +
			Math.cos(b * RAD) * Math.sin(23.4397 * RAD) * Math.sin(lon_ * RAD)
		) * DEG;
		const GMST = (6.697375 + 0.0657098242 * d + (t.getUTCHours() + t.getUTCMinutes() / 60)) % 24;
		const LMST = (GMST + lon / 15) % 24;
		const ra = Math.atan2(
			Math.sin(lon_ * RAD) * Math.cos(23.4397 * RAD) - Math.tan(b * RAD) * Math.sin(23.4397 * RAD),
			Math.cos(lon_ * RAD)
		) * DEG;
		const ha = LMST * 15 - ra;
		const elev = Math.asin(
			Math.sin(lat * RAD) * Math.sin(dec * RAD) +
			Math.cos(lat * RAD) * Math.cos(dec * RAD) * Math.cos(ha * RAD)
		) * DEG;

		if(prev !== null) {
			if(prev < 0 && elev >= 0 && !moonrise) moonrise = new Date(t.getTime() - 5 * 60000).toISOString();
			if(prev >= 0 && elev < 0 && !moonset) moonset = new Date(t.getTime() - 5 * 60000).toISOString();
		}
		prev = elev;
		if(moonrise && moonset) break;
	}

	return {moonrise: moonrise, moonset: moonset};
}

function nextSolsticeEquinox(date = new Date()) {
	const y = date.getFullYear();
	return {
		summer_solstice: new Date(Date.UTC(y, 5, 21)),
		winter_solstice: new Date(Date.UTC(y, 11, 21)),
		vernal_equinox: new Date(Date.UTC(y, 2, 20)),
		autumnal_equinox: new Date(Date.UTC(y + 1, 2, 20)),
	};
}

export function getCelestialCurrent(lat, lon, date = new Date()) {
	return {
		...sunPosition(lat, lon, date),
		...sunriseSunset(lat, lon, date),
		...moonPhase(date),
		...moonPosition(lat, lon, date),
		...sunMoonDistance(date),
		...nextMoonEvents(date),
		...moonriseMoonset(lat, lon, date),
	};
}

export function getCelestialForecast(lat, lon, hours) {
	return hours.map((data) => Object.assign(data, getCelestialCurrent(lat, lon, new Date(data.time))))
}