Advertisement
onlinebacon

Eclipse solver

Apr 13th, 2024
908
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. // -------------------------------------------------------------------------- //
  2. // BUILT IN MATH FUNCTIONS AND GENERAL TRIG FUNCTIONS ----------------------- //
  3. // -------------------------------------------------------------------------- //
  4.  
  5. const { PI, abs, sqrt, sign, sin, cos, tan, asin, acos, atan } = Math;
  6.  
  7. const TAU = PI * 2;
  8. const D90 = PI / 2;
  9. const D30 = PI / 6;
  10.  
  11. const toRad = (deg) => deg / 180 * PI;
  12. const toDeg = (rad) => rad / PI * 180;
  13.  
  14. const oneWayAngleDiff = (a, b) => {
  15.     return a < b ? a + TAU - b : a - b;
  16. };
  17.  
  18. // -------------------------------------------------------------------------- //
  19. // CONSTANTS ---------------------------------------------------------------- //
  20. // -------------------------------------------------------------------------- //
  21.  
  22. const SEC = 1000;
  23. const MIN = SEC * 60;
  24. const HOUR = MIN * 60;
  25. const KM = 1;
  26. const MILE = 1.609344 * KM;
  27. const AU = 150e6 * KM;
  28. const KMPH = KM / HOUR;
  29. const MPH = MILE / HOUR;
  30.  
  31. const moonRadius = 1737.4 * KM;
  32. const earthRadius = 6371 * KM;
  33. const sunRadius = (1392700 / 2) * KM;
  34.  
  35. const overrideMoonDist = 359829 * KM;
  36. const overrideSunDist = 149820729 * KM;
  37.  
  38. // -------------------------------------------------------------------------- //
  39. // AUXILIAR FUNCTIONS ------------------------------------------------------- //
  40. // -------------------------------------------------------------------------- //
  41.  
  42. const parseAngle = (string) => {
  43.     const sepRegex = /^\s*[°'"]\s*|\s+/;
  44.     const numRegex = /^\d+(\.\d+)?/;
  45.     const unitMap = {
  46.         '°': 1,
  47.         "'": 1 / 60,
  48.         '"': 1 / 60 / 60,
  49.     };
  50.     string = string.trim();
  51.     let sign = 1;
  52.     if (string.startsWith('-')) {
  53.         sign = -1;
  54.         string = string.substring(1).trim();
  55.     }
  56.     let sum = 0;
  57.     let unit = 1;
  58.     while (string.length > 0) {
  59.         const num = string.match(numRegex)?.[0];
  60.         if (!num) {
  61.             return NaN;
  62.         }
  63.         string = string.substring(num.length);
  64.         const sep = string.match(sepRegex)?.[0];
  65.         if (sep) {
  66.             string = string.substring(sep.length);
  67.             const short = sep.trim();
  68.             if (short !== '') {
  69.                 unit = unitMap[short];
  70.             }
  71.         } else if (string !== '') {
  72.             return NaN;
  73.         }
  74.         sum += num * unit;
  75.         unit /= 60;
  76.     }
  77.     return toRad(sum * sign);
  78. };
  79.  
  80. const round = (value, n = 0) => {
  81.     return Number(value.toFixed(n));
  82. };
  83.  
  84. const strAngle = (angle, pos = '', neg = '-') => {
  85.     const tSec = round(toDeg(abs(angle)) * 3600, 1);
  86.     const sec = round(tSec % 60, 1);
  87.     const tMin = round((tSec - sec)/60);
  88.     const min = round(tMin % 60);
  89.     const deg = round((tMin - min)/60);
  90.     const sign = angle < 0 ? neg : pos;
  91.     return `${sign}${deg}°${min}'${sec}" / ${round(toDeg(angle), 6).toString()}°`;
  92. };
  93.  
  94. const numberWithNDigits = (num, n, suffix = '') => {
  95.     if (abs(num) >= (10 ** n)) {
  96.         return Math.round(num).toString() + suffix;
  97.     }
  98.     return Number(num.toPrecision(n)) + suffix;
  99. };
  100.  
  101. const strDist = (dist, n = 5) => {
  102.     const km = numberWithNDigits(dist/KM, n, ' km');
  103.     const mi = numberWithNDigits(dist/MILE, n, ' mi');
  104.     return  `${km} / ${mi}`;
  105. };
  106.  
  107. const strSpeed = (speed, n = 5) => {
  108.     const kmph = numberWithNDigits(speed / KMPH, n, ' km/h');
  109.     const mph = numberWithNDigits(speed / MPH, n, ' mph');
  110.     return `${kmph} / ${mph}`;
  111. };
  112.  
  113. const strTime = (timestamp) => {
  114.     const date = new Date(timestamp).toISOString();
  115.     return date.replace(/^.*T/, '').replace('Z', ' (UTC)');
  116. };
  117.  
  118. const strDuration = (duration) => {
  119.     const tSec = round(duration / SEC);
  120.     const sec = tSec % 60;
  121.     const tMin = (tSec - sec)/60;
  122.     const min = tMin % 60;
  123.     const hrs = (tMin - min)/60;
  124.     const strSec = sec.toString().padStart(2, '0');
  125.     const strMin = min.toString().padStart(2, '0');
  126.     const strHrs = hrs.toString().padStart(2, '0');
  127.     return `${strHrs}h${strMin}m${strSec}s`;
  128. };
  129.  
  130. // -------------------------------------------------------------------------- //
  131. // DATA PARSING ------------------------------------------------------------- //
  132. // -------------------------------------------------------------------------- //
  133.  
  134. const table = [];
  135.  
  136. const addAlmanacData = (text) => {
  137.     let date = null;
  138.     const lines = text.trim().split(/\n/);
  139.     for (let line of lines) {
  140.         line = line.trim();
  141.         if (line.startsWith('#')) {
  142.             continue;
  143.         }
  144.         if (/^date:/i.test(line)) {
  145.             const strDate = line.replace(/^date:/i, '').trim();
  146.             const isoDate = strDate + 'T00:00:00Z';
  147.             date = new Date(isoDate).getTime();
  148.             continue;
  149.         }
  150.         const [ hour, ...angles ] = line.trim().split(/\s*\|\s*/);
  151.         const time = date + hour * HOUR;
  152.         const [ sunGHA, sunDec, moonGHA, moonDec, moonHP ] = angles.map(parseAngle);
  153.         const obj = { time, sunGHA, sunDec, moonGHA, moonDec, moonHP };
  154.         table.push(obj);
  155.     }
  156. };
  157.  
  158. const interpolate = (time, field) => {
  159.     for (let i=1; i<table.length; i++) {
  160.         const bef = table[i - 1];
  161.         if (time < bef.time) {
  162.             continue;
  163.         }
  164.         const aft = table[i];
  165.         if (time > aft.time) {
  166.             continue;
  167.         }
  168.         const t = (time - bef.time) / (aft.time - bef.time);
  169.         const a = bef[field];
  170.         const b = aft[field];
  171.         return a + (b - a) * t;
  172.     }
  173.     throw new Error(`time ${new Date(time).toISOString()} is out of bounds`);
  174. };
  175.  
  176. // -------------------------------------------------------------------------- //
  177. // VECTOR OPERATIONS (LINEAR ALGEBRA) --------------------------------------- //
  178. // -------------------------------------------------------------------------- //
  179.  
  180. const subVec = ([ ax, ay, az ], [ bx, by, bz ]) => {
  181.     return [ ax - bx, ay - by, az - bz ];
  182. };
  183.  
  184. const sumVec = ([ ax, ay, az ], [ bx, by, bz ]) => {
  185.     return [ ax + bx, ay + by, az + bz ];
  186. };
  187.  
  188. const normalizeVec = ([ x, y, z ]) => {
  189.     const len = sqrt(x*x + y*y + z*z);
  190.     return [ x/len, y/len, z/len ];
  191. };
  192.  
  193. const dotProd = (a, b) => {
  194.     const [ ax, ay, az ] = a;
  195.     const [ bx, by, bz ] = b;
  196.     return (ax * bx) + (ay * by) + (az * bz);
  197. };
  198.  
  199. const scaleVec = (vec, val) => {
  200.     const [ x, y, z ] = vec;
  201.     return [ x*val, y*val, z*val ];
  202. };
  203.  
  204. const vecLen = ([ x, y, z ]) => {
  205.     return sqrt(x*x + y*y + z*z);
  206. };
  207.  
  208. const rotVecY = ([ x, y, z ], angle) => {
  209.     const s = sin(angle);
  210.     const c = cos(angle);
  211.     return [ x*c - z*s, y, z*c + x*s ];
  212. };
  213.  
  214. const rotVecX = ([ x, y, z ], angle) => {
  215.     const s = sin(angle);
  216.     const c = cos(angle);
  217.     return [ x, y*c + z*s, z*c - y*s ];
  218. };
  219.  
  220. const rotVecZ = ([ x, y, z ], angle) => {
  221.     const s = sin(angle);
  222.     const c = cos(angle);
  223.     return [ x*c + y*s, y*c - x*s, z ];
  224. };
  225.  
  226. const vecIsValid = ([ x, y, z ]) => {
  227.     return !(isNaN(x) || isNaN(y) || isNaN(z));
  228. };
  229.  
  230. // -------------------------------------------------------------------------- //
  231. // LAT/LONG MATH ------------------------------------------------------------ //
  232. // -------------------------------------------------------------------------- //
  233.  
  234. // Geocentric
  235. const GC = {
  236.     vecToGP: (vec) => {
  237.         const [ x, y, z ] = normalizeVec(vec);
  238.         const lat = asin(z);
  239.         const len = sqrt(x*x + y*y);
  240.         const lon = len === 0 ? 0 : acos(x/len) * (y < 0 ? -1 : 1);
  241.         return [ lat, lon ];
  242.     },
  243.     gpToUnitVec: (gp) => {
  244.         const [ lat, lon ] = gp;
  245.         const x = + cos(lat) * cos(lon);
  246.         const y = + cos(lat) * sin(lon);
  247.         const z = + sin(lat);
  248.         return [ x, y, z ];
  249.     },
  250.     gpDistToVec: (gp, dist) => {
  251.         const vec = GC.gpToUnitVec(gp);
  252.         return scaleVec(vec, dist);
  253.     },
  254.     calcAngle: (gp1, gp2) => {
  255.         const aVec = GC.gpToUnitVec(gp1);
  256.         const bVec = GC.gpToUnitVec(gp2);
  257.         return acos(dotProd(aVec, bVec));
  258.     },
  259.     calcDistance: (gp1, gp2) => {
  260.         return GC.calcAngle(gp1, gp2) * earthRadius;
  261.     },
  262.     calcAzimuth: (gp1, gp2) => {
  263.         const [ lat1, lon1 ] = gp1;
  264.         const [ lat2, lon2 ] = gp2;
  265.         const dlon = lon2 - lon1;
  266.         const x = cos(lat2) * cos(dlon);
  267.         const y = cos(lat2) * sin(dlon);
  268.         const z = sin(lat2);
  269.         const newZ = z*cos(lat1) - x*sin(lat1);
  270.         const east = acos(newZ/sqrt(newZ ** 2 + y ** 2));
  271.         return y >= 0 ? east : TAU - east;
  272.     },
  273. };
  274.  
  275. // -------------------------------------------------------------------------- //
  276. // ECLIPSE MATH ------------------------------------------------------------- //
  277. // -------------------------------------------------------------------------- //
  278.  
  279. const ghaToLon = (gha) => {
  280.     return (TAU*1.5 - gha) % TAU - PI;
  281. };
  282.  
  283. const getMoonGhaDec = (time) => {
  284.     const gha = interpolate(time, 'moonGHA');
  285.     const dec = interpolate(time, 'moonDec');
  286.     return [ gha, dec ];
  287. };
  288.  
  289. const getSunGhaDec = (time) => {
  290.     const gha = interpolate(time, 'sunGHA');
  291.     const dec = interpolate(time, 'sunDec');
  292.     return [ gha, dec ];
  293. };
  294.  
  295. const hpToDist = (hp) => {
  296.     const dist = earthRadius / tan(hp);
  297.     return dist;
  298. };
  299.  
  300. const findGHAMeridianMatch = () => {
  301.     for (let i=1; i<table.length; ++i) {
  302.         const bef = table[i - 1];
  303.         const aft = table[i];
  304.         const sunApproaching = oneWayAngleDiff(bef.moonGHA, bef.sunGHA) < D30;
  305.         const sunPasses = oneWayAngleDiff(aft.sunGHA, aft.moonGHA) < D30;
  306.         if (!sunApproaching || !sunPasses) {
  307.             continue;
  308.         }
  309.         const dt = aft.time - bef.time;
  310.         const dSun = aft.sunGHA - bef.sunGHA;
  311.         const dMoon = aft.moonGHA - bef.moonGHA;
  312.         const t = (bef.moonGHA - bef.sunGHA) / (dSun - dMoon);
  313.         const time = bef.time + dt * t;
  314.         return time;
  315.     }
  316.     return null;
  317. };
  318.  
  319. const ghaDecDistToVec = (gha, dec, dist) => {
  320.     const gp = [ dec, ghaToLon(gha) ];
  321.     return GC.gpDistToVec(gp, dist);
  322. };
  323.  
  324. const getMoonDist = (time) => {
  325.     return overrideMoonDist ?? hpToDist(interpolate(time, 'moonHP'));
  326. };
  327.  
  328. const getMoonVec = (time) => {
  329.     const [ gha, dec ] = getMoonGhaDec(time);
  330.     const dist = getMoonDist(time);
  331.     return ghaDecDistToVec(gha, dec, dist);
  332. };
  333.  
  334. const getSunDist = (time) => {
  335.     return overrideSunDist ?? AU;
  336. };
  337.  
  338. const getSunVec = (time) => {
  339.     const [ gha, dec ] = getSunGhaDec(time);
  340.     const dist = getSunDist(time);
  341.     return ghaDecDistToVec(gha, dec, dist);
  342. };
  343.  
  344. const getCenterShadowRay = (moonVec, sunVec) => {
  345.     const dir = normalizeVec(subVec(moonVec, sunVec));
  346.     const origin = moonVec;
  347.     return { origin, dir };
  348. };
  349.  
  350. const rayEarthIntersection = ({ origin, dir }) => {
  351.     const midT = dotProd(subVec([ 0, 0, 0 ], origin), dir);
  352.     if (midT < 0) {
  353.         return [ NaN, NaN, NaN ];
  354.     }
  355.     const midVec = sumVec(origin, scaleVec(dir, midT));
  356.     const centerDist = vecLen(midVec);
  357.     if (centerDist > earthRadius) {
  358.         return [ NaN, NaN, NaN ];
  359.     }
  360.     const tDiff = sqrt(earthRadius**2 - centerDist**2);
  361.     const t = midT - tDiff;
  362.     return sumVec(origin, scaleVec(dir, t));
  363. };
  364.  
  365. const calcUmbraAngle = (moonSunDist) => {
  366.     return asin((sunRadius - moonRadius)/moonSunDist);
  367. };
  368.  
  369. const calcPenumbraAngle = (moonSunDist) => {
  370.     return - asin((sunRadius + moonRadius)/moonSunDist);
  371. };
  372.  
  373. const buildShadowEdgeRay = (moonVec, sunVec, azimuth, calcAngle) => {
  374.     const moonSunDist = vecLen(subVec(moonVec, sunVec));
  375.     const angle = calcAngle(moonSunDist);
  376.    
  377.     let origin = [ 0, 0, moonRadius ];
  378.     let dir = [ -1, 0, 0 ];
  379.  
  380.     origin = rotVecY(origin, angle);
  381.     dir = rotVecY(dir, angle)
  382.  
  383.     origin = rotVecX(origin, azimuth);
  384.     dir = rotVecX(dir, azimuth)
  385.  
  386.     const [ lat, lon ] = GC.vecToGP(subVec(sunVec, moonVec));
  387.  
  388.     origin = rotVecY(origin, lat);
  389.     dir = rotVecY(dir, lat);
  390.  
  391.     origin = rotVecZ(origin, -lon);
  392.     dir = rotVecZ(dir, -lon);
  393.  
  394.     origin = sumVec(origin, moonVec);
  395.  
  396.     return { origin, dir };
  397. };
  398.  
  399. const getShadowEdgeCoord = (moonVec, sunVec, azimuth, calcAngle) => {
  400.     return GC.vecToGP(
  401.         rayEarthIntersection(
  402.             buildShadowEdgeRay(moonVec, sunVec, azimuth, calcAngle),
  403.         ),
  404.     );
  405. };
  406.  
  407. const calcShadowSize = (moonVec, sunVec, azimuth, calcAngle) => {
  408.     const az1 = azimuth;
  409.     const az2 = (azimuth + PI) % TAU;
  410.     const a = getShadowEdgeCoord(moonVec, sunVec, az1, calcAngle);
  411.     const b = getShadowEdgeCoord(moonVec, sunVec, az2, calcAngle);
  412.     return GC.calcDistance(a, b);
  413. };
  414.  
  415. const calcVecs = (time) => {
  416.     const moonVec = getMoonVec(time);
  417.     const [ sunGHA, sunDec ] = getSunGhaDec(time);
  418.     const sunGP = [ sunDec, ghaToLon(sunGHA) ];
  419.     const sunDist = getSunDist(time);
  420.     const sunVec = ghaDecDistToVec(sunGHA, sunDec, sunDist);
  421.     const locVec = rayEarthIntersection(getCenterShadowRay(moonVec, sunVec));
  422.     const loc = GC.vecToGP(locVec);
  423.     return { sunGP, moonVec, sunVec, loc, locVec };
  424. };
  425.  
  426. const hasTotality = (locVec, time) => {
  427.     const { sunVec, moonVec } = calcVecs(time);
  428.     const moonAngRad = calcAngularRadius(locVec, moonVec, moonRadius);
  429.     const sunAngRad = calcAngularRadius(locVec, sunVec, sunRadius);
  430.     const sunDir = normalizeVec(subVec(sunVec, locVec));
  431.     const moonDir = normalizeVec(subVec(moonVec, locVec));
  432.     const prod = dotProd(moonDir, sunDir);
  433.     if (prod >= 1) {
  434.         return true;
  435.     }
  436.     const angDist = acos(dotProd(moonDir, sunDir));
  437.     const gap = moonAngRad - sunAngRad;
  438.     return angDist <= gap;
  439. };
  440.  
  441. const calcAngularRadius = (locVec, bodyVec, bodyRad) => {
  442.     const dist = vecLen(subVec(locVec, bodyVec));
  443.     return asin(bodyRad/dist);
  444. };
  445.  
  446. const findTotalityEdgeAt = (locVec, startTime, endTime) => {
  447.     let at = startTime;
  448.     let av = hasTotality(locVec, at);
  449.     let bt = endTime;
  450.     let bv = hasTotality(locVec, bt);
  451.     if (av === bv) {
  452.         return null;
  453.     }
  454.     while (bt - at > 1) {
  455.         const mt = (at + bt)/2;
  456.         const mv = hasTotality(locVec, mt);
  457.         if (mv === av) {
  458.             at = mt;
  459.             av = mv;
  460.         } else {
  461.             bt = mt;
  462.             bv = mv;
  463.         }
  464.     }
  465.     return (at + bt)/2;
  466. };
  467.  
  468. // -------------------------------------------------------------------------- //
  469. // MAIN --------------------------------------------------------------------- //
  470. // -------------------------------------------------------------------------- //
  471.  
  472. addAlmanacData(`
  473.     Date: 2024-04-08
  474.     # Hour | Sun GHA    | Sun Dec  | Moon GHA   | Moon Dec | Moon HP
  475.     #------|------------|----------|------------|----------|---------
  476.       18   |  89° 35.5' | 7° 35.2' |  89° 54.4' | 7° 48.9' | 60.9'
  477.       19   | 104° 35.6' | 7° 36.2' | 104° 23.2' | 8° 06.3' | 60.9'
  478. `);
  479.  
  480. const peakTime = findGHAMeridianMatch();
  481. const { sunGP, moonVec, sunVec, loc, locVec } = calcVecs(peakTime);
  482. const [ lat, lon ] = loc;
  483. const umbraSizeNS = calcShadowSize(moonVec, sunVec,   0, calcUmbraAngle);
  484. const umbraSizeEW = calcShadowSize(moonVec, sunVec, D90, calcUmbraAngle);
  485. const penumbraSizeNS = calcShadowSize(moonVec, sunVec,   0, calcPenumbraAngle);
  486. const penumbraSizeEW = calcShadowSize(moonVec, sunVec, D90, calcPenumbraAngle);
  487. const sunAlt = D90 - GC.calcAngle(sunGP, loc);
  488.  
  489. const dtTime = 1 * SEC;
  490. const nextTime = peakTime + dtTime;
  491. const nextLoc = calcVecs(nextTime).loc;
  492. const dtDist = GC.calcDistance(loc, nextLoc);
  493. const speed = dtDist/dtTime;
  494. const dir = GC.calcAzimuth(loc, nextLoc);
  495. const pathWidth = calcShadowSize(moonVec, sunVec, dir + D90, calcUmbraAngle);
  496. const startTime = findTotalityEdgeAt(locVec, table[0].time, peakTime);
  497. const endTime = findTotalityEdgeAt(locVec, peakTime, table.at(-1).time);
  498. const duration = endTime - startTime;
  499. const sunAngRad = calcAngularRadius(locVec, sunVec, sunRadius);
  500. const moonAngRad = calcAngularRadius(locVec, moonVec, moonRadius);
  501.  
  502. console.log(' 1. Latitude:', strAngle(lat, 'N ', 'S '));
  503. console.log(' 2. Longitude:', strAngle(lon, 'E ', 'W '));
  504. console.log(' 3. Time:', strTime(peakTime));
  505. console.log(' 4. Sun Alt.:', strAngle(sunAlt));
  506. console.log(' 5. Moon:Sun ratio:', (moonAngRad/sunAngRad).toFixed(4));
  507. console.log(' 6. Umbra E-W size:', strDist(umbraSizeEW));
  508. console.log(' 7. Umbra N-S size:', strDist(umbraSizeNS));
  509. console.log(' 8. Path angle:', strAngle(dir));
  510. console.log(' 9. Path width:', strDist(pathWidth));
  511. console.log('10. Ground speed:', strSpeed(speed));
  512. console.log('11. Duration:', strDuration(duration));
  513. console.log('12. Penumbra E-W size:', strDist(penumbraSizeEW));
  514. console.log('12. Penumbra N-S size:', strDist(penumbraSizeNS));
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement