validate-eclipse-data.mjs

import fs from 'node:fs';
import vm from 'node:vm';

const context = { window: {} };
vm.runInNewContext(fs.readFileSync(new URL('./eclipse-data.js', import.meta.url), 'utf8'), context);
const eclipses = context.window.ECLIPSE_DATA;
const radians = Math.PI / 180;

function clockHours(value) {
  const [hours, minutes, seconds] = value.split(':').map(Number);
  return hours + minutes / 60 + seconds / 3600;
}

function rowHours(row) {
  return Number.isFinite(row[9]) ? row[9] / 3600 : clockHours(row[0]);
}

function state(eclipse, utcHours) {
  const elements = eclipse.elements;
  let time = utcHours - elements[1] + elements[4] / 3600;
  while (time > 12) time -= 24;
  while (time < -12) time += 24;
  const x = ((elements[8] * time + elements[7]) * time + elements[6]) * time + elements[5];
  const y = ((elements[12] * time + elements[11]) * time + elements[10]) * time + elements[9];
  const declination = ((elements[15] * time + elements[14]) * time + elements[13]) * radians;
  return {
    x,
    y,
    sinDeclination: Math.sin(declination),
    cosDeclination: Math.cos(declination),
    hourAngle: ((elements[18] * time + elements[17]) * time + elements[16]) * radians,
    umbralRadius: (elements[24] * time + elements[23]) * time + elements[22]
  };
}

function margin(eclipse, shadow, point) {
  const elements = eclipse.elements;
  const latitude = point[0] * radians;
  const longitude = point[1] * radians;
  const geocentricLatitude = Math.atan(.99664719 * Math.tan(latitude));
  const rhoSin = .99664719 * Math.sin(geocentricLatitude);
  const rhoCos = Math.cos(geocentricLatitude);
  const localHourAngle = shadow.hourAngle + longitude - elements[4] / 13713.44;
  const xi = rhoCos * Math.sin(localHourAngle);
  const eta = rhoSin * shadow.cosDeclination - rhoCos * Math.cos(localHourAngle) * shadow.sinDeclination;
  const zeta = rhoSin * shadow.sinDeclination + rhoCos * Math.cos(localHourAngle) * shadow.cosDeclination;
  const localRadius = shadow.umbralRadius - zeta * elements[26];
  return Math.hypot(shadow.x - xi, shadow.y - eta) + localRadius;
}

function solarAltitude(eclipse, shadow, point) {
  const latitude = point[0] * radians;
  const localHourAngle = shadow.hourAngle + point[1] * radians - eclipse.elements[4] / 13713.44;
  return Math.asin(shadow.sinDeclination * Math.sin(latitude) + shadow.cosDeclination * Math.cos(latitude) * Math.cos(localHourAngle)) / radians;
}

function destination(center, bearing, distanceKm) {
  const angularDistance = distanceKm / 6371.0088;
  const latitude = center[0] * radians;
  const longitude = center[1] * radians;
  const direction = bearing * radians;
  const resultLatitude = Math.asin(Math.sin(latitude) * Math.cos(angularDistance) + Math.cos(latitude) * Math.sin(angularDistance) * Math.cos(direction));
  const resultLongitude = longitude + Math.atan2(Math.sin(direction) * Math.sin(angularDistance) * Math.cos(latitude), Math.cos(angularDistance) - Math.sin(latitude) * Math.sin(resultLatitude));
  let longitudeDegrees = resultLongitude / radians;
  while (longitudeDegrees - center[1] > 180) longitudeDegrees -= 360;
  while (longitudeDegrees - center[1] < -180) longitudeDegrees += 360;
  return [resultLatitude / radians, longitudeDegrees];
}

function outline(eclipse, row) {
  const shadow = state(eclipse, rowHours(row));
  const centerAltitude = solarAltitude(eclipse, shadow, row[1]);
  const points = [];
  for (let bearing = 0; bearing < 360; bearing += 3.75) {
    let inside = 0;
    let outside = Math.max(40, row[4]);
    while (margin(eclipse, shadow, destination(row[1], bearing, outside)) < 0 && outside < 5000) outside *= 1.6;
    if (outside >= 5000 && margin(eclipse, shadow, destination(row[1], bearing, outside)) < 0) continue;
    for (let iteration = 0; iteration < 20; iteration++) {
      const midpoint = (inside + outside) / 2;
      if (margin(eclipse, shadow, destination(row[1], bearing, midpoint)) < 0) inside = midpoint;
      else outside = midpoint;
    }
    let boundaryDistance = (inside + outside) / 2;
    let boundaryPoint = destination(row[1], bearing, boundaryDistance);
    if (solarAltitude(eclipse, shadow, boundaryPoint) < 0) {
      if (centerAltitude <= 0) {
        boundaryPoint = row[1];
      } else {
        let visibleDistance = 0;
        let hiddenDistance = boundaryDistance;
        for (let iteration = 0; iteration < 24; iteration++) {
          const midpoint = (visibleDistance + hiddenDistance) / 2;
          if (solarAltitude(eclipse, shadow, destination(row[1], bearing, midpoint)) >= 0) visibleDistance = midpoint;
          else hiddenDistance = midpoint;
        }
        boundaryPoint = destination(row[1], bearing, (visibleDistance + hiddenDistance) / 2);
      }
    }
    points.push(boundaryPoint);
  }
  return points;
}

const failures = [];
let maximumBoundaryResidual = 0;
let maximumOutlineResidual = 0;
let completeMidpointOutlines = 0;
for (const eclipse of eclipses) {
  if (eclipse.elements?.length !== 27) failures.push(`${eclipse.id}: invalid Besselian elements`);
  if (eclipse.timeline.length < 2) failures.push(`${eclipse.id}: insufficient timeline rows`);
  if (eclipse.limits?.length !== 2) failures.push(`${eclipse.id}: missing terminal path limits`);
  for (const limit of eclipse.limits || []) {
    for (const point of [limit.north, limit.south, limit.center]) {
      if (!point || point.length !== 2 || point.some(value => !Number.isFinite(value))) failures.push(`${eclipse.id}: invalid terminal limit coordinate`);
    }
  }
  for (let index = 1; index < eclipse.timeline.length; index++) {
    if (rowHours(eclipse.timeline[index]) <= rowHours(eclipse.timeline[index - 1])) failures.push(`${eclipse.id}: non-increasing timeline time`);
    if (!(eclipse.timeline[index][6] > 0)) failures.push(`${eclipse.id}: non-positive speed`);
  }
  for (const row of eclipse.timeline) {
    const shadow = state(eclipse, rowHours(row));
    for (const point of [row[7], row[8]]) {
      if (point) maximumBoundaryResidual = Math.max(maximumBoundaryResidual, Math.abs(margin(eclipse, shadow, point)));
    }
    if (margin(eclipse, shadow, row[1]) >= 0) failures.push(`${eclipse.id} ${row[0]}: centerline outside umbra`);
  }
  const middle = eclipse.timeline[Math.floor(eclipse.timeline.length / 2)];
  const footprint = outline(eclipse, middle);
  if (footprint.length !== 96) failures.push(`${eclipse.id}: incomplete midpoint outline (${footprint.length}/96)`);
  else completeMidpointOutlines++;
  const midpointState = state(eclipse, rowHours(middle));
  for (const point of footprint) {
    const coneResidual = Math.abs(margin(eclipse, midpointState, point));
    const horizonResidual = Math.abs(solarAltitude(eclipse, midpointState, point)) * radians;
    maximumOutlineResidual = Math.max(maximumOutlineResidual, Math.min(coneResidual, horizonResidual));
    if (margin(eclipse, midpointState, point) > 1e-6) failures.push(`${eclipse.id}: outline point outside umbra`);
    if (solarAltitude(eclipse, midpointState, point) < -1e-5) failures.push(`${eclipse.id}: outline point below horizon`);
  }
}

if (failures.length) {
  console.error(failures.join('\n'));
  process.exitCode = 1;
} else {
  console.log(JSON.stringify({
    eclipses: eclipses.length,
    first: eclipses[0].id,
    last: eclipses.at(-1).id,
    maximumBoundaryResidual,
    maximumOutlineResidual,
    completeMidpointOutlines: `${completeMidpointOutlines}/${eclipses.length}`
  }, null, 2));
}