Local forecasting
This commit is contained in:
237
server/src/celestial.mjs
Normal file
237
server/src/celestial.mjs
Normal file
@@ -0,0 +1,237 @@
|
||||
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))))
|
||||
}
|
||||
Reference in New Issue
Block a user