Eclipse solver

// -------------------------------------------------------------------------- //
// BUILT IN MATH FUNCTIONS AND GENERAL TRIG FUNCTIONS ----------------------- //
// -------------------------------------------------------------------------- //

const { PI, abs, sqrt, sign, sin, cos, tan, asin, acos, atan } = Math;

const TAU = PI * 2;
const D90 = PI / 2;
const D30 = PI / 6;

const toRad = (deg) => deg / 180 * PI;
const toDeg = (rad) => rad / PI * 180;

const oneWayAngleDiff = (a, b) => {
    return a < b ? a + TAU - b : a - b;
};

// -------------------------------------------------------------------------- //
// CONSTANTS ---------------------------------------------------------------- //
// -------------------------------------------------------------------------- //

const SEC = 1000;
const MIN = SEC * 60;
const HOUR = MIN * 60;
const KM = 1;
const MILE = 1.609344 * KM;
const AU = 150e6 * KM;
const KMPH = KM / HOUR;
const MPH = MILE / HOUR;

const moonRadius = 1737.4 * KM;
const earthRadius = 6371 * KM;
const sunRadius = (1392700 / 2) * KM;

const overrideMoonDist = 359829 * KM;
const overrideSunDist = 149820729 * KM;

// -------------------------------------------------------------------------- //
// AUXILIAR FUNCTIONS ------------------------------------------------------- //
// -------------------------------------------------------------------------- //

const parseAngle = (string) => {
    const sepRegex = /^\s*[°'"]\s*|\s+/;
    const numRegex = /^\d+(\.\d+)?/;
    const unitMap = {
        '°': 1,
        "'": 1 / 60,
        '"': 1 / 60 / 60,
    };
    string = string.trim();
    let sign = 1;
    if (string.startsWith('-')) {
        sign = -1;
        string = string.substring(1).trim();
    }
    let sum = 0;
    let unit = 1;
    while (string.length > 0) {
        const num = string.match(numRegex)?.[0];
        if (!num) {
            return NaN;
        }
        string = string.substring(num.length);
        const sep = string.match(sepRegex)?.[0];
        if (sep) {
            string = string.substring(sep.length);
            const short = sep.trim();
            if (short !== '') {
                unit = unitMap[short];
            }
        } else if (string !== '') {
            return NaN;
        }
        sum += num * unit;
        unit /= 60;
    }
    return toRad(sum * sign);
};

const round = (value, n = 0) => {
    return Number(value.toFixed(n));
};

const strAngle = (angle, pos = '', neg = '-') => {
    const tSec = round(toDeg(abs(angle)) * 3600, 1);
    const sec = round(tSec % 60, 1);
    const tMin = round((tSec - sec)/60);
    const min = round(tMin % 60);
    const deg = round((tMin - min)/60);
    const sign = angle < 0 ? neg : pos;
    return `${sign}${deg}°${min}'${sec}" / ${round(toDeg(angle), 6).toString()}°`;
};

const numberWithNDigits = (num, n, suffix = '') => {
    if (abs(num) >= (10 ** n)) {
        return Math.round(num).toString() + suffix;
    }
    return Number(num.toPrecision(n)) + suffix;
};

const strDist = (dist, n = 5) => {
    const km = numberWithNDigits(dist/KM, n, ' km');
    const mi = numberWithNDigits(dist/MILE, n, ' mi');
    return  `${km} / ${mi}`;
};

const strSpeed = (speed, n = 5) => {
    const kmph = numberWithNDigits(speed / KMPH, n, ' km/h');
    const mph = numberWithNDigits(speed / MPH, n, ' mph');
    return `${kmph} / ${mph}`;
};

const strTime = (timestamp) => {
    const date = new Date(timestamp).toISOString();
    return date.replace(/^.*T/, '').replace('Z', ' (UTC)');
};

const strDuration = (duration) => {
    const tSec = round(duration / SEC);
    const sec = tSec % 60;
    const tMin = (tSec - sec)/60;
    const min = tMin % 60;
    const hrs = (tMin - min)/60;
    const strSec = sec.toString().padStart(2, '0');
    const strMin = min.toString().padStart(2, '0');
    const strHrs = hrs.toString().padStart(2, '0');
    return `${strHrs}h${strMin}m${strSec}s`;
};

// -------------------------------------------------------------------------- //
// DATA PARSING ------------------------------------------------------------- //
// -------------------------------------------------------------------------- //

const table = [];

const addAlmanacData = (text) => {
    let date = null;
    const lines = text.trim().split(/\n/);
    for (let line of lines) {
        line = line.trim();
        if (line.startsWith('#')) {
            continue;
        }
        if (/^date:/i.test(line)) {
            const strDate = line.replace(/^date:/i, '').trim();
            const isoDate = strDate + 'T00:00:00Z';
            date = new Date(isoDate).getTime();
            continue;
        }
        const [ hour, ...angles ] = line.trim().split(/\s*\|\s*/);
        const time = date + hour * HOUR;
        const [ sunGHA, sunDec, moonGHA, moonDec, moonHP ] = angles.map(parseAngle);
        const obj = { time, sunGHA, sunDec, moonGHA, moonDec, moonHP };
        table.push(obj);
    }
};

const interpolate = (time, field) => {
    for (let i=1; i<table.length; i++) {
        const bef = table[i - 1];
        if (time < bef.time) {
            continue;
        }
        const aft = table[i];
        if (time > aft.time) {
            continue;
        }
        const t = (time - bef.time) / (aft.time - bef.time);
        const a = bef[field];
        const b = aft[field];
        return a + (b - a) * t;
    }
    throw new Error(`time ${new Date(time).toISOString()} is out of bounds`);
};

// -------------------------------------------------------------------------- //
// VECTOR OPERATIONS (LINEAR ALGEBRA) --------------------------------------- //
// -------------------------------------------------------------------------- //

const subVec = ([ ax, ay, az ], [ bx, by, bz ]) => {
    return [ ax - bx, ay - by, az - bz ];
};

const sumVec = ([ ax, ay, az ], [ bx, by, bz ]) => {
    return [ ax + bx, ay + by, az + bz ];
};

const normalizeVec = ([ x, y, z ]) => {
    const len = sqrt(x*x + y*y + z*z);
    return [ x/len, y/len, z/len ];
};

const dotProd = (a, b) => {
    const [ ax, ay, az ] = a;
    const [ bx, by, bz ] = b;
    return (ax * bx) + (ay * by) + (az * bz);
};

const scaleVec = (vec, val) => {
    const [ x, y, z ] = vec;
    return [ x*val, y*val, z*val ];
};

const vecLen = ([ x, y, z ]) => {
    return sqrt(x*x + y*y + z*z);
};

const rotVecY = ([ x, y, z ], angle) => {
    const s = sin(angle);
    const c = cos(angle);
    return [ x*c - z*s, y, z*c + x*s ];
};

const rotVecX = ([ x, y, z ], angle) => {
    const s = sin(angle);
    const c = cos(angle);
    return [ x, y*c + z*s, z*c - y*s ];
};

const rotVecZ = ([ x, y, z ], angle) => {
    const s = sin(angle);
    const c = cos(angle);
    return [ x*c + y*s, y*c - x*s, z ];
};

const vecIsValid = ([ x, y, z ]) => {
    return !(isNaN(x) || isNaN(y) || isNaN(z));
};

// -------------------------------------------------------------------------- //
// LAT/LONG MATH ------------------------------------------------------------ //
// -------------------------------------------------------------------------- //

// Geocentric
const GC = {
    vecToGP: (vec) => {
        const [ x, y, z ] = normalizeVec(vec);
        const lat = asin(z);
        const len = sqrt(x*x + y*y);
        const lon = len === 0 ? 0 : acos(x/len) * (y < 0 ? -1 : 1);
        return [ lat, lon ];
    },
    gpToUnitVec: (gp) => {
        const [ lat, lon ] = gp;
        const x = + cos(lat) * cos(lon);
        const y = + cos(lat) * sin(lon);
        const z = + sin(lat);
        return [ x, y, z ];
    },
    gpDistToVec: (gp, dist) => {
        const vec = GC.gpToUnitVec(gp);
        return scaleVec(vec, dist);
    },
    calcAngle: (gp1, gp2) => {
        const aVec = GC.gpToUnitVec(gp1);
        const bVec = GC.gpToUnitVec(gp2);
        return acos(dotProd(aVec, bVec));
    },
    calcDistance: (gp1, gp2) => {
        return GC.calcAngle(gp1, gp2) * earthRadius;
    },
    calcAzimuth: (gp1, gp2) => {
        const [ lat1, lon1 ] = gp1;
        const [ lat2, lon2 ] = gp2;
        const dlon = lon2 - lon1;
        const x = cos(lat2) * cos(dlon);
        const y = cos(lat2) * sin(dlon);
        const z = sin(lat2);
        const newZ = z*cos(lat1) - x*sin(lat1);
        const east = acos(newZ/sqrt(newZ ** 2 + y ** 2));
        return y >= 0 ? east : TAU - east;
    },
};

// -------------------------------------------------------------------------- //
// ECLIPSE MATH ------------------------------------------------------------- //
// -------------------------------------------------------------------------- //

const ghaToLon = (gha) => {
    return (TAU*1.5 - gha) % TAU - PI;
};

const getMoonGhaDec = (time) => {
    const gha = interpolate(time, 'moonGHA');
    const dec = interpolate(time, 'moonDec');
    return [ gha, dec ];
};

const getSunGhaDec = (time) => {
    const gha = interpolate(time, 'sunGHA');
    const dec = interpolate(time, 'sunDec');
    return [ gha, dec ];
};

const hpToDist = (hp) => {
    const dist = earthRadius / tan(hp);
    return dist;
};

const findGHAMeridianMatch = () => {
    for (let i=1; i<table.length; ++i) {
        const bef = table[i - 1];
        const aft = table[i];
        const sunApproaching = oneWayAngleDiff(bef.moonGHA, bef.sunGHA) < D30;
        const sunPasses = oneWayAngleDiff(aft.sunGHA, aft.moonGHA) < D30;
        if (!sunApproaching || !sunPasses) {
            continue;
        }
        const dt = aft.time - bef.time;
        const dSun = aft.sunGHA - bef.sunGHA;
        const dMoon = aft.moonGHA - bef.moonGHA;
        const t = (bef.moonGHA - bef.sunGHA) / (dSun - dMoon);
        const time = bef.time + dt * t;
        return time;
    }
    return null;
};

const ghaDecDistToVec = (gha, dec, dist) => {
    const gp = [ dec, ghaToLon(gha) ];
    return GC.gpDistToVec(gp, dist);
};

const getMoonDist = (time) => {
    return overrideMoonDist ?? hpToDist(interpolate(time, 'moonHP'));
};

const getMoonVec = (time) => {
    const [ gha, dec ] = getMoonGhaDec(time);
    const dist = getMoonDist(time);
    return ghaDecDistToVec(gha, dec, dist);
};

const getSunDist = (time) => {
    return overrideSunDist ?? AU;
};

const getSunVec = (time) => {
    const [ gha, dec ] = getSunGhaDec(time);
    const dist = getSunDist(time);
    return ghaDecDistToVec(gha, dec, dist);
};

const getCenterShadowRay = (moonVec, sunVec) => {
    const dir = normalizeVec(subVec(moonVec, sunVec));
    const origin = moonVec;
    return { origin, dir };
};

const rayEarthIntersection = ({ origin, dir }) => {
    const midT = dotProd(subVec([ 0, 0, 0 ], origin), dir);
    if (midT < 0) {
        return [ NaN, NaN, NaN ];
    }
    const midVec = sumVec(origin, scaleVec(dir, midT));
    const centerDist = vecLen(midVec);
    if (centerDist > earthRadius) {
        return [ NaN, NaN, NaN ];
    }
    const tDiff = sqrt(earthRadius**2 - centerDist**2);
    const t = midT - tDiff;
    return sumVec(origin, scaleVec(dir, t));
};

const calcUmbraAngle = (moonSunDist) => {
    return asin((sunRadius - moonRadius)/moonSunDist);
};

const calcPenumbraAngle = (moonSunDist) => {
    return - asin((sunRadius + moonRadius)/moonSunDist);
};

const buildShadowEdgeRay = (moonVec, sunVec, azimuth, calcAngle) => {
    const moonSunDist = vecLen(subVec(moonVec, sunVec));
    const angle = calcAngle(moonSunDist);

    let origin = [ 0, 0, moonRadius ];
    let dir = [ -1, 0, 0 ];

    origin = rotVecY(origin, angle);
    dir = rotVecY(dir, angle)

    origin = rotVecX(origin, azimuth);
    dir = rotVecX(dir, azimuth)

    const [ lat, lon ] = GC.vecToGP(subVec(sunVec, moonVec));

    origin = rotVecY(origin, lat);
    dir = rotVecY(dir, lat);

    origin = rotVecZ(origin, -lon);
    dir = rotVecZ(dir, -lon);

    origin = sumVec(origin, moonVec);

    return { origin, dir };
};

const getShadowEdgeCoord = (moonVec, sunVec, azimuth, calcAngle) => {
    return GC.vecToGP(
        rayEarthIntersection(
            buildShadowEdgeRay(moonVec, sunVec, azimuth, calcAngle),
        ),
    );
};

const calcShadowSize = (moonVec, sunVec, azimuth, calcAngle) => {
    const az1 = azimuth;
    const az2 = (azimuth + PI) % TAU;
    const a = getShadowEdgeCoord(moonVec, sunVec, az1, calcAngle);
    const b = getShadowEdgeCoord(moonVec, sunVec, az2, calcAngle);
    return GC.calcDistance(a, b);
};

const calcVecs = (time) => {
    const moonVec = getMoonVec(time);
    const [ sunGHA, sunDec ] = getSunGhaDec(time);
    const sunGP = [ sunDec, ghaToLon(sunGHA) ];
    const sunDist = getSunDist(time);
    const sunVec = ghaDecDistToVec(sunGHA, sunDec, sunDist);
    const locVec = rayEarthIntersection(getCenterShadowRay(moonVec, sunVec));
    const loc = GC.vecToGP(locVec);
    return { sunGP, moonVec, sunVec, loc, locVec };
};

const hasTotality = (locVec, time) => {
    const { sunVec, moonVec } = calcVecs(time);
    const moonAngRad = calcAngularRadius(locVec, moonVec, moonRadius);
    const sunAngRad = calcAngularRadius(locVec, sunVec, sunRadius);
    const sunDir = normalizeVec(subVec(sunVec, locVec));
    const moonDir = normalizeVec(subVec(moonVec, locVec));
    const prod = dotProd(moonDir, sunDir);
    if (prod >= 1) {
        return true;
    }
    const angDist = acos(dotProd(moonDir, sunDir));
    const gap = moonAngRad - sunAngRad;
    return angDist <= gap;
};

const calcAngularRadius = (locVec, bodyVec, bodyRad) => {
    const dist = vecLen(subVec(locVec, bodyVec));
    return asin(bodyRad/dist);
};

const findTotalityEdgeAt = (locVec, startTime, endTime) => {
    let at = startTime;
    let av = hasTotality(locVec, at);
    let bt = endTime;
    let bv = hasTotality(locVec, bt);
    if (av === bv) {
        return null;
    }
    while (bt - at > 1) {
        const mt = (at + bt)/2;
        const mv = hasTotality(locVec, mt);
        if (mv === av) {
            at = mt;
            av = mv;
        } else {
            bt = mt;
            bv = mv;
        }
    }
    return (at + bt)/2;
};

// -------------------------------------------------------------------------- //
// MAIN --------------------------------------------------------------------- //
// -------------------------------------------------------------------------- //

addAlmanacData(`
    Date: 2024-04-08
    # Hour | Sun GHA    | Sun Dec  | Moon GHA   | Moon Dec | Moon HP
    #------|------------|----------|------------|----------|---------
      18   |  89° 35.5' | 7° 35.2' |  89° 54.4' | 7° 48.9' | 60.9'
      19   | 104° 35.6' | 7° 36.2' | 104° 23.2' | 8° 06.3' | 60.9'
`);

const peakTime = findGHAMeridianMatch();
const { sunGP, moonVec, sunVec, loc, locVec } = calcVecs(peakTime);
const [ lat, lon ] = loc;
const umbraSizeNS = calcShadowSize(moonVec, sunVec,   0, calcUmbraAngle);
const umbraSizeEW = calcShadowSize(moonVec, sunVec, D90, calcUmbraAngle);
const penumbraSizeNS = calcShadowSize(moonVec, sunVec,   0, calcPenumbraAngle);
const penumbraSizeEW = calcShadowSize(moonVec, sunVec, D90, calcPenumbraAngle);
const sunAlt = D90 - GC.calcAngle(sunGP, loc);

const dtTime = 1 * SEC;
const nextTime = peakTime + dtTime;
const nextLoc = calcVecs(nextTime).loc;
const dtDist = GC.calcDistance(loc, nextLoc);
const speed = dtDist/dtTime;
const dir = GC.calcAzimuth(loc, nextLoc);
const pathWidth = calcShadowSize(moonVec, sunVec, dir + D90, calcUmbraAngle);
const startTime = findTotalityEdgeAt(locVec, table[0].time, peakTime);
const endTime = findTotalityEdgeAt(locVec, peakTime, table.at(-1).time);
const duration = endTime - startTime;
const sunAngRad = calcAngularRadius(locVec, sunVec, sunRadius);
const moonAngRad = calcAngularRadius(locVec, moonVec, moonRadius);

console.log(' 1. Latitude:', strAngle(lat, 'N ', 'S '));
console.log(' 2. Longitude:', strAngle(lon, 'E ', 'W '));
console.log(' 3. Time:', strTime(peakTime));
console.log(' 4. Sun Alt.:', strAngle(sunAlt));
console.log(' 5. Moon:Sun ratio:', (moonAngRad/sunAngRad).toFixed(4));
console.log(' 6. Umbra E-W size:', strDist(umbraSizeEW));
console.log(' 7. Umbra N-S size:', strDist(umbraSizeNS));
console.log(' 8. Path angle:', strAngle(dir));
console.log(' 9. Path width:', strDist(pathWidth));
console.log('10. Ground speed:', strSpeed(speed));
console.log('11. Duration:', strDuration(duration));
console.log('12. Penumbra E-W size:', strDist(penumbraSizeEW));
console.log('12. Penumbra N-S size:', strDist(penumbraSizeNS));