From df4dbfe53e76ff919aba95d5655553962e3ba550 Mon Sep 17 00:00:00 2001 From: Aric Camarata Date: Sun, 31 May 2026 08:48:31 -0400 Subject: [PATCH] ci: fix eslint parser devDeps, add files patterns, format src - Add @typescript-eslint/parser and @typescript-eslint/eslint-plugin as explicit devDependencies so pnpm hoists them for eslint.config.mjs - Add files: ['**/*.ts'] to eslint config entries so ESLint 10 processes TS sources instead of ignoring them - Add parserOptions.project for typed-lint rules - Run prettier --write to fix pre-existing format issues in 12 src files --- eslint.config.mjs | 11 +- package.json | 2 + pnpm-lock.yaml | 6 + src/api/index.ts | 489 ++++++++++++++++++++-------------------- src/bodies/index.ts | 180 +++++++-------- src/cli/index.ts | 178 +++++++-------- src/events/index.ts | 168 +++++++------- src/frames/index.ts | 120 +++++----- src/index.ts | 6 +- src/math/index.ts | 207 ++++++++--------- src/observer/index.ts | 118 +++++----- src/spk/index.ts | 328 +++++++++++++-------------- src/time/index.ts | 128 +++++------ src/types.ts | 286 +++++++++++------------ src/visibility/index.ts | 134 +++++------ 15 files changed, 1194 insertions(+), 1167 deletions(-) diff --git a/eslint.config.mjs b/eslint.config.mjs index 442dce0..515c70f 100644 --- a/eslint.config.mjs +++ b/eslint.config.mjs @@ -5,10 +5,17 @@ import { typescript } from '@acamarata/eslint-config'; export default [ { + files: ['**/*.ts'], plugins: { '@typescript-eslint': tsPlugin }, - languageOptions: { parser: tsParser }, + languageOptions: { + parser: tsParser, + parserOptions: { + project: true, + tsconfigRootDir: import.meta.dirname, + }, + }, }, - ...typescript, + ...typescript.map((c) => ({ ...c, files: ['**/*.ts'] })), eslintConfigPrettier, { ignores: ['dist/', 'node_modules/', 'test.mjs', 'test-cjs.cjs'], diff --git a/package.json b/package.json index 23419af..cab9aa2 100644 --- a/package.json +++ b/package.json @@ -49,6 +49,8 @@ "@acamarata/tsconfig": "^0.1.0", "@eslint/js": "^10.0.1", "@types/node": "^25.3.0", + "@typescript-eslint/eslint-plugin": "^8.0.0", + "@typescript-eslint/parser": "^8.0.0", "eslint": "^10.0.3", "eslint-config-prettier": "^10.1.8", "prettier": "^3.8.1", diff --git a/pnpm-lock.yaml b/pnpm-lock.yaml index bdee4aa..07d6eee 100644 --- a/pnpm-lock.yaml +++ b/pnpm-lock.yaml @@ -23,6 +23,12 @@ importers: '@types/node': specifier: ^25.3.0 version: 25.3.0 + '@typescript-eslint/eslint-plugin': + specifier: ^8.0.0 + version: 8.56.1(@typescript-eslint/parser@8.56.1(eslint@10.0.3)(typescript@5.9.3))(eslint@10.0.3)(typescript@5.9.3) + '@typescript-eslint/parser': + specifier: ^8.0.0 + version: 8.56.1(eslint@10.0.3)(typescript@5.9.3) c8: specifier: ^10.1.3 version: 10.1.3 diff --git a/src/api/index.ts b/src/api/index.ts index f586a71..7f9a844 100644 --- a/src/api/index.ts +++ b/src/api/index.ts @@ -28,10 +28,10 @@ import type { KernelConfig, OdehZone, Vec3, -} from '../types.js' -import { ODEH_THRESHOLDS, ODEH_DESCRIPTIONS } from '../types.js' -import { SpkKernel } from '../spk/index.js' -import { computeTimeScales, jdTTtoET, jdToDate, J2000 } from '../time/index.js' +} from "../types.js"; +import { ODEH_THRESHOLDS, ODEH_DESCRIPTIONS } from "../types.js"; +import { SpkKernel } from "../spk/index.js"; +import { computeTimeScales, jdTTtoET, jdToDate, J2000 } from "../time/index.js"; import { getMoonGeocentricState, getSunGeocentricState, @@ -39,52 +39,52 @@ import { computeCrescentWidth, getMoonSunApproximate, nearestNewMoon, -} from '../bodies/index.js' -import { geodeticToECEF, computeAzAlt } from '../observer/index.js' -import { itrsToGcrs, computeERA } from '../frames/index.js' +} from "../bodies/index.js"; +import { geodeticToECEF, computeAzAlt } from "../observer/index.js"; +import { itrsToGcrs, computeERA } from "../frames/index.js"; import { getSunMoonEvents as eventsGetSunMoonEvents, bestTimeHeuristic, bestTimeOptimized, computeObservationWindow, -} from '../events/index.js' +} from "../events/index.js"; import { computeCrescentGeometry, computeYallop, computeOdeh, buildGuidanceText, arcvMinimum, -} from '../visibility/index.js' -import { DEG2RAD } from '../math/index.js' +} from "../visibility/index.js"; +import { DEG2RAD } from "../math/index.js"; // ─── Input validation ───────────────────────────────────────────────────────── function validateDate(date: Date, label: string): void { if (!(date instanceof Date) || isNaN(date.getTime())) { - throw new RangeError(`${label}: expected a valid Date instance`) + throw new RangeError(`${label}: expected a valid Date instance`); } } function validateLatitude(lat: number, label: string): void { if (!isFinite(lat) || lat < -90 || lat > 90) { - throw new RangeError(`${label}: latitude must be a finite number in [-90, 90], got ${lat}`) + throw new RangeError(`${label}: latitude must be a finite number in [-90, 90], got ${lat}`); } } function validateLongitude(lon: number, label: string): void { if (!isFinite(lon) || lon < -180 || lon > 180) { - throw new RangeError(`${label}: longitude must be a finite number in [-180, 180], got ${lon}`) + throw new RangeError(`${label}: longitude must be a finite number in [-180, 180], got ${lon}`); } } function validateObserver(observer: Observer, label: string): void { - validateLatitude(observer.lat, label) - validateLongitude(observer.lon, label) + validateLatitude(observer.lat, label); + validateLongitude(observer.lon, label); } // ─── Module-level kernel singleton ───────────────────────────────────────────── -let activeKernel: SpkKernel | null = null +let activeKernel: SpkKernel | null = null; // ─── Cache directory resolution ──────────────────────────────────────────────── @@ -92,18 +92,18 @@ let activeKernel: SpkKernel | null = null * Resolve the platform-appropriate kernel cache directory. */ function resolveCacheDir(override?: string): string { - if (override) return override - const { platform, env } = process - if (platform === 'win32') { - return `${env['LOCALAPPDATA'] ?? env['APPDATA'] ?? 'C:\\Users\\Public\\AppData\\Local'}\\moon-sighting` + if (override) return override; + const { platform, env } = process; + if (platform === "win32") { + return `${env["LOCALAPPDATA"] ?? env["APPDATA"] ?? "C:\\Users\\Public\\AppData\\Local"}\\moon-sighting`; } - return `${env['HOME'] ?? '/tmp'}/.cache/moon-sighting` + return `${env["HOME"] ?? "/tmp"}/.cache/moon-sighting`; } // ─── Download sources ───────────────────────────────────────────────────────── -const NAIF_DE442S_URL = 'https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de442s.bsp' -const NAIF_LSK_URL = 'https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls' +const NAIF_DE442S_URL = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de442s.bsp"; +const NAIF_LSK_URL = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls"; // ─── Kernel lifecycle ───────────────────────────────────────────────────────── @@ -119,33 +119,33 @@ const NAIF_LSK_URL = 'https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/nai * @param config - Kernel source configuration. Defaults to auto-download. */ export async function initKernels(config?: KernelConfig): Promise { - const source = config?.planetary ?? { type: 'auto' as const } + const source = config?.planetary ?? { type: "auto" as const }; - let buffer: ArrayBuffer + let buffer: ArrayBuffer; - if (source.type === 'file') { - buffer = await readFileAsBuffer(source.path) - } else if (source.type === 'buffer') { - buffer = source.data - } else if (source.type === 'url') { - const res = await fetch(source.url) + if (source.type === "file") { + buffer = await readFileAsBuffer(source.path); + } else if (source.type === "buffer") { + buffer = source.data; + } else if (source.type === "url") { + const res = await fetch(source.url); if (!res.ok) - throw new Error(`Failed to fetch kernel from ${source.url}: ${res.status} ${res.statusText}`) - buffer = await res.arrayBuffer() + throw new Error(`Failed to fetch kernel from ${source.url}: ${res.status} ${res.statusText}`); + buffer = await res.arrayBuffer(); } else { // auto: download to local cache, then load - const paths = await downloadKernels(config) - buffer = await readFileAsBuffer(paths.planetaryPath) + const paths = await downloadKernels(config); + buffer = await readFileAsBuffer(paths.planetaryPath); } - activeKernel = SpkKernel.fromBuffer(buffer) + activeKernel = SpkKernel.fromBuffer(buffer); } /** Read a file into an ArrayBuffer (Node.js only). */ async function readFileAsBuffer(filePath: string): Promise { - const { readFile } = await import('node:fs/promises') - const buf = await readFile(filePath) - return buf.buffer.slice(buf.byteOffset, buf.byteOffset + buf.byteLength) as ArrayBuffer + const { readFile } = await import("node:fs/promises"); + const buf = await readFile(filePath); + return buf.buffer.slice(buf.byteOffset, buf.byteOffset + buf.byteLength) as ArrayBuffer; } /** @@ -157,59 +157,60 @@ async function readFileAsBuffer(filePath: string): Promise { * @returns Paths where kernels were saved */ export async function downloadKernels(config?: KernelConfig): Promise<{ - planetaryPath: string - leapSecondsPath: string + planetaryPath: string; + leapSecondsPath: string; }> { - const cacheDir = resolveCacheDir(config?.cacheDir) - const { mkdir, writeFile } = await import('node:fs/promises') - const { existsSync } = await import('node:fs') - const { join } = await import('node:path') + const cacheDir = resolveCacheDir(config?.cacheDir); + const { mkdir, writeFile } = await import("node:fs/promises"); + const { existsSync } = await import("node:fs"); + const { join } = await import("node:path"); - await mkdir(cacheDir, { recursive: true }) + await mkdir(cacheDir, { recursive: true }); - const planetaryPath = join(cacheDir, 'de442s.bsp') - const leapSecondsPath = join(cacheDir, 'naif0012.tls') + const planetaryPath = join(cacheDir, "de442s.bsp"); + const leapSecondsPath = join(cacheDir, "naif0012.tls"); if (!existsSync(planetaryPath)) { - process.stdout.write('Downloading de442s.bsp from NAIF... ') - const res = await fetch(NAIF_DE442S_URL) - if (!res.ok) throw new Error(`Failed to download de442s.bsp: ${res.status} ${res.statusText}`) - const buf = await res.arrayBuffer() - await writeFile(planetaryPath, Buffer.from(buf)) - console.log(`done (${(buf.byteLength / 1048576).toFixed(1)} MB)`) + process.stdout.write("Downloading de442s.bsp from NAIF... "); + const res = await fetch(NAIF_DE442S_URL); + if (!res.ok) throw new Error(`Failed to download de442s.bsp: ${res.status} ${res.statusText}`); + const buf = await res.arrayBuffer(); + await writeFile(planetaryPath, Buffer.from(buf)); + console.log(`done (${(buf.byteLength / 1048576).toFixed(1)} MB)`); if (config?.checksumOverride) { - const actual = await sha256File(planetaryPath) + const actual = await sha256File(planetaryPath); if (actual !== config.checksumOverride.toLowerCase()) { throw new Error( `de442s.bsp checksum mismatch.\n Expected: ${config.checksumOverride}\n Got: ${actual}`, - ) + ); } } } else { - console.log('de442s.bsp already cached.') + console.log("de442s.bsp already cached."); } if (!existsSync(leapSecondsPath)) { - process.stdout.write('Downloading naif0012.tls from NAIF... ') - const res = await fetch(NAIF_LSK_URL) - if (!res.ok) throw new Error(`Failed to download naif0012.tls: ${res.status} ${res.statusText}`) - const text = await res.text() - await writeFile(leapSecondsPath, text, 'utf8') - console.log('done.') + process.stdout.write("Downloading naif0012.tls from NAIF... "); + const res = await fetch(NAIF_LSK_URL); + if (!res.ok) + throw new Error(`Failed to download naif0012.tls: ${res.status} ${res.statusText}`); + const text = await res.text(); + await writeFile(leapSecondsPath, text, "utf8"); + console.log("done."); } else { - console.log('naif0012.tls already cached.') + console.log("naif0012.tls already cached."); } - return { planetaryPath, leapSecondsPath } + return { planetaryPath, leapSecondsPath }; } /** Compute the SHA-256 hex digest of a local file. */ async function sha256File(filePath: string): Promise { - const { createHash } = await import('node:crypto') - const { readFile } = await import('node:fs/promises') - const buf = await readFile(filePath) - return createHash('sha256').update(buf).digest('hex') + const { createHash } = await import("node:crypto"); + const { readFile } = await import("node:fs/promises"); + const buf = await readFile(filePath); + return createHash("sha256").update(buf).digest("hex"); } /** @@ -219,33 +220,33 @@ async function sha256File(filePath: string): Promise { * @returns { ok, errors[] } — ok is true when all checks pass */ export async function verifyKernels(config?: KernelConfig): Promise<{ - ok: boolean - errors: string[] + ok: boolean; + errors: string[]; }> { - const cacheDir = resolveCacheDir(config?.cacheDir) - const { existsSync } = await import('node:fs') - const { join } = await import('node:path') - const errors: string[] = [] + const cacheDir = resolveCacheDir(config?.cacheDir); + const { existsSync } = await import("node:fs"); + const { join } = await import("node:path"); + const errors: string[] = []; - const planetaryPath = join(cacheDir, 'de442s.bsp') - const leapSecondsPath = join(cacheDir, 'naif0012.tls') + const planetaryPath = join(cacheDir, "de442s.bsp"); + const leapSecondsPath = join(cacheDir, "naif0012.tls"); if (!existsSync(planetaryPath)) { - errors.push(`de442s.bsp not found at ${planetaryPath}. Run downloadKernels() first.`) + errors.push(`de442s.bsp not found at ${planetaryPath}. Run downloadKernels() first.`); } else if (config?.checksumOverride) { - const actual = await sha256File(planetaryPath) + const actual = await sha256File(planetaryPath); if (actual !== config.checksumOverride.toLowerCase()) { errors.push( `de442s.bsp checksum mismatch.\n Expected: ${config.checksumOverride}\n Got: ${actual}`, - ) + ); } } if (!existsSync(leapSecondsPath)) { - errors.push(`naif0012.tls not found at ${leapSecondsPath}. Run downloadKernels() first.`) + errors.push(`naif0012.tls not found at ${leapSecondsPath}. Run downloadKernels() first.`); } - return { ok: errors.length === 0, errors } + return { ok: errors.length === 0, errors }; } // ─── Kernel resolution ───────────────────────────────────────────────────────── @@ -256,25 +257,25 @@ export async function verifyKernels(config?: KernelConfig): Promise<{ */ async function resolveKernel(config?: KernelConfig): Promise { if (config?.planetary) { - const source = config.planetary - if (source.type === 'file') { - return SpkKernel.fromBuffer(await readFileAsBuffer(source.path)) - } else if (source.type === 'buffer') { - return SpkKernel.fromBuffer(source.data) - } else if (source.type === 'url') { - const res = await fetch(source.url) - if (!res.ok) throw new Error(`Failed to fetch kernel: ${res.status}`) - return SpkKernel.fromBuffer(await res.arrayBuffer()) + const source = config.planetary; + if (source.type === "file") { + return SpkKernel.fromBuffer(await readFileAsBuffer(source.path)); + } else if (source.type === "buffer") { + return SpkKernel.fromBuffer(source.data); + } else if (source.type === "url") { + const res = await fetch(source.url); + if (!res.ok) throw new Error(`Failed to fetch kernel: ${res.status}`); + return SpkKernel.fromBuffer(await res.arrayBuffer()); } } - if (activeKernel) return activeKernel + if (activeKernel) return activeKernel; // auto-init as last resort - await initKernels(config) + await initKernels(config); if (!activeKernel) - throw new Error('Kernel failed to initialize. Call initKernels() before computing.') - return activeKernel + throw new Error("Kernel failed to initialize. Call initKernels() before computing."); + return activeKernel; } // ─── Primary API ────────────────────────────────────────────────────────────── @@ -308,72 +309,72 @@ export async function getMoonSightingReport( observer: Observer, options?: SightingOptions, ): Promise { - validateDate(date, 'getMoonSightingReport') - validateObserver(observer, 'getMoonSightingReport') - const kernel = await resolveKernel(options?.kernels) + validateDate(date, "getMoonSightingReport"); + validateObserver(observer, "getMoonSightingReport"); + const kernel = await resolveKernel(options?.kernels); // Event times (sunset, moonset, twilight, rise) - const events = eventsGetSunMoonEvents(date, observer, kernel) - const { sunsetUTC, moonsetUTC } = events + const events = eventsGetSunMoonEvents(date, observer, kernel); + const { sunsetUTC, moonsetUTC } = events; if (!sunsetUTC || !moonsetUTC) { - return buildNullReport(date, observer, events, 'DE442S', false) + return buildNullReport(date, observer, events, "DE442S", false); } // Best observation time - const method = options?.bestTimeMethod ?? 'heuristic' - let bestTimeResult: { bestTimeUTC: Date; lagMinutes: number } | null = null + const method = options?.bestTimeMethod ?? "heuristic"; + let bestTimeResult: { bestTimeUTC: Date; lagMinutes: number } | null = null; - if (method === 'optimized') { - const opt = bestTimeOptimized(sunsetUTC, moonsetUTC, kernel, observer) - if (opt) bestTimeResult = { bestTimeUTC: opt.bestTimeUTC, lagMinutes: opt.lagMinutes } + if (method === "optimized") { + const opt = bestTimeOptimized(sunsetUTC, moonsetUTC, kernel, observer); + if (opt) bestTimeResult = { bestTimeUTC: opt.bestTimeUTC, lagMinutes: opt.lagMinutes }; } if (!bestTimeResult) { - bestTimeResult = bestTimeHeuristic(sunsetUTC, moonsetUTC) + bestTimeResult = bestTimeHeuristic(sunsetUTC, moonsetUTC); } if (!bestTimeResult) { - return buildNullReport(date, observer, events, 'DE442S', false) + return buildNullReport(date, observer, events, "DE442S", false); } - const { bestTimeUTC, lagMinutes } = bestTimeResult - const bestTimeWindowUTC = computeObservationWindow(bestTimeUTC) + const { bestTimeUTC, lagMinutes } = bestTimeResult; + const bestTimeWindowUTC = computeObservationWindow(bestTimeUTC); // Time scales and ephemeris time at best time - const ts = computeTimeScales(bestTimeUTC, observer.ut1utc, observer.deltaT) - const et = jdTTtoET(ts.jdTT) + const ts = computeTimeScales(bestTimeUTC, observer.ut1utc, observer.deltaT); + const et = jdTTtoET(ts.jdTT); // Body positions in GCRS (geocentric) - const moonGCRS = getMoonGeocentricState(kernel, et).position - const sunGCRS = getSunGeocentricState(kernel, et).position + const moonGCRS = getMoonGeocentricState(kernel, et).position; + const sunGCRS = getSunGeocentricState(kernel, et).position; // Observer ITRS position (km) from geodetic coordinates - const obsECEF = geodeticToECEF(observer.lat, observer.lon, observer.elevation) - const obsITRS: Vec3 = [obsECEF[0] / 1000, obsECEF[1] / 1000, obsECEF[2] / 1000] + const obsECEF = geodeticToECEF(observer.lat, observer.lon, observer.elevation); + const obsITRS: Vec3 = [obsECEF[0] / 1000, obsECEF[1] / 1000, obsECEF[2] / 1000]; // Convert to GCRS (inertial frame) — required for correct topocentric subtraction // GCRS body vectors (from SPK) and observer must be in the same frame before subtracting - const obsGCRS = itrsToGcrs(obsITRS, ts) + const obsGCRS = itrsToGcrs(obsITRS, ts); // Airless alt/az — required by Yallop/Odeh criteria - const moonAirless = computeAzAlt(moonGCRS, observer, ts, true) - const sunAirless = computeAzAlt(sunGCRS, observer, ts, true) + const moonAirless = computeAzAlt(moonGCRS, observer, ts, true); + const sunAirless = computeAzAlt(sunGCRS, observer, ts, true); // Apparent alt/az (with refraction) — for guidance text - const moonApparent = computeAzAlt(moonGCRS, observer, ts, false) + const moonApparent = computeAzAlt(moonGCRS, observer, ts, false); // Illumination and moon age - const illumData = computeIllumination(moonGCRS, sunGCRS) - const illumination = illumData.illumination * 100 - const prevNewMoonJD = nearestNewMoon(ts.jdTT - 15) - const moonAgeHours = (ts.jdTT - prevNewMoonJD) * 24 + const illumData = computeIllumination(moonGCRS, sunGCRS); + const illumination = illumData.illumination * 100; + const prevNewMoonJD = nearestNewMoon(ts.jdTT - 15); + const moonAgeHours = (ts.jdTT - prevNewMoonJD) * 24; // Topocentric vectors for crescent geometry (GCRS - observer GCRS) const moonTopo: Vec3 = [ moonGCRS[0] - obsGCRS[0], moonGCRS[1] - obsGCRS[1], moonGCRS[2] - obsGCRS[2], - ] - const sunTopo: Vec3 = [sunGCRS[0] - obsGCRS[0], sunGCRS[1] - obsGCRS[1], sunGCRS[2] - obsGCRS[2]] + ]; + const sunTopo: Vec3 = [sunGCRS[0] - obsGCRS[0], sunGCRS[1] - obsGCRS[1], sunGCRS[2] - obsGCRS[2]]; const geometry = computeCrescentGeometry( moonAirless, @@ -382,14 +383,14 @@ export async function getMoonSightingReport( sunTopo, sunsetUTC, moonsetUTC, - ) + ); - const { Wprime } = computeCrescentWidth(moonTopo, geometry.ARCL) - const yallop = computeYallop(geometry, Wprime) - const odeh = computeOdeh(geometry) + const { Wprime } = computeCrescentWidth(moonTopo, geometry.ARCL); + const yallop = computeYallop(geometry, Wprime); + const odeh = computeOdeh(geometry); - const moonAboveHorizon = moonAirless.altitude > 0 - const sightingPossible = moonAboveHorizon && lagMinutes > 0 + const moonAboveHorizon = moonAirless.altitude > 0; + const sightingPossible = moonAboveHorizon && lagMinutes > 0; const guidance = buildGuidanceText( yallop, @@ -398,7 +399,7 @@ export async function getMoonSightingReport( moonApparent.altitude, bestTimeUTC, lagMinutes, - ) + ); return { date, @@ -416,10 +417,10 @@ export async function getMoonSightingReport( yallop, odeh, guidance, - ephemerisSource: 'DE442S', + ephemerisSource: "DE442S", moonAboveHorizon, sightingPossible, - } + }; } /** Build a null report for cases where sighting geometry cannot be computed. */ @@ -427,7 +428,7 @@ function buildNullReport( date: Date, observer: Observer, events: SunMoonEvents, - source: 'DE442S' | 'approximate', + source: "DE442S" | "approximate", sightingPossible: boolean, ): MoonSightingReport { return { @@ -446,25 +447,25 @@ function buildNullReport( yallop: null, odeh: null, guidance: - 'Sighting not possible: sunset or moonset could not be determined for this date and location.', + "Sighting not possible: sunset or moonset could not be determined for this date and location.", ephemerisSource: source, moonAboveHorizon: null, sightingPossible, - } + }; } // ─── Phase display lookup ────────────────────────────────────────────────────── const PHASE_DISPLAY: Record = { - 'new-moon': { name: 'New Moon', symbol: '🌑' }, - 'waxing-crescent': { name: 'Waxing Crescent', symbol: '🌒' }, - 'first-quarter': { name: 'First Quarter', symbol: '🌓' }, - 'waxing-gibbous': { name: 'Waxing Gibbous', symbol: '🌔' }, - 'full-moon': { name: 'Full Moon', symbol: '🌕' }, - 'waning-gibbous': { name: 'Waning Gibbous', symbol: '🌖' }, - 'last-quarter': { name: 'Last Quarter', symbol: '🌗' }, - 'waning-crescent': { name: 'Waning Crescent', symbol: '🌘' }, -} + "new-moon": { name: "New Moon", symbol: "🌑" }, + "waxing-crescent": { name: "Waxing Crescent", symbol: "🌒" }, + "first-quarter": { name: "First Quarter", symbol: "🌓" }, + "waxing-gibbous": { name: "Waxing Gibbous", symbol: "🌔" }, + "full-moon": { name: "Full Moon", symbol: "🌕" }, + "waning-gibbous": { name: "Waning Gibbous", symbol: "🌖" }, + "last-quarter": { name: "Last Quarter", symbol: "🌗" }, + "waning-crescent": { name: "Waning Crescent", symbol: "🌘" }, +}; /** * Compute the Moon's current phase, illumination, and next phase times. @@ -485,23 +486,23 @@ const PHASE_DISPLAY: Record = { * ``` */ export function getMoonPhase(date = new Date()): MoonPhaseResult { - validateDate(date, 'getMoonPhase') - const ts = computeTimeScales(date) - const { moonGCRS, sunGCRS } = getMoonSunApproximate(ts.jdTT) + validateDate(date, "getMoonPhase"); + const ts = computeTimeScales(date); + const { moonGCRS, sunGCRS } = getMoonSunApproximate(ts.jdTT); - const { illumination, elongationDeg, isWaxing } = computeIllumination(moonGCRS, sunGCRS) - const illuminationPct = illumination * 100 + const { illumination, elongationDeg, isWaxing } = computeIllumination(moonGCRS, sunGCRS); + const illuminationPct = illumination * 100; // Age in hours since previous new moon // Search 15 days back/forward to ensure we clear the current lunation boundary - const prevNewMoonJD = nearestNewMoon(ts.jdTT - 15) - const age = (ts.jdTT - prevNewMoonJD) * 24 + const prevNewMoonJD = nearestNewMoon(ts.jdTT - 15); + const age = (ts.jdTT - prevNewMoonJD) * 24; - const phaseKey = elongationToPhase(elongationDeg, isWaxing) - const { name: phaseName, symbol: phaseSymbol } = PHASE_DISPLAY[phaseKey] + const phaseKey = elongationToPhase(elongationDeg, isWaxing); + const { name: phaseName, symbol: phaseSymbol } = PHASE_DISPLAY[phaseKey]; - const nextNewMoonJD = nearestNewMoon(ts.jdTT + 15) - const nextFullMoonJD = nearestFullMoon(ts.jdTT) + const nextNewMoonJD = nearestNewMoon(ts.jdTT + 15); + const nextFullMoonJD = nearestFullMoon(ts.jdTT); return { phase: phaseKey, @@ -514,7 +515,7 @@ export function getMoonPhase(date = new Date()): MoonPhaseResult { nextNewMoon: jdToDate(nextNewMoonJD), nextFullMoon: jdToDate(nextFullMoonJD), prevNewMoon: jdToDate(prevNewMoonJD), - } + }; } /** @@ -542,34 +543,34 @@ export function getMoonPosition( lon: number, elevation = 0, ): MoonPosition { - validateDate(date, 'getMoonPosition') - validateLatitude(lat, 'getMoonPosition') - validateLongitude(lon, 'getMoonPosition') - const ts = computeTimeScales(date) - const { moonGCRS } = getMoonSunApproximate(ts.jdTT) + validateDate(date, "getMoonPosition"); + validateLatitude(lat, "getMoonPosition"); + validateLongitude(lon, "getMoonPosition"); + const ts = computeTimeScales(date); + const { moonGCRS } = getMoonSunApproximate(ts.jdTT); // Apparent az/alt with Bennett refraction — uses existing observer pipeline - const observer: Observer = { lat, lon, elevation } - const azAlt = computeAzAlt(moonGCRS, observer, ts, false) + const observer: Observer = { lat, lon, elevation }; + const azAlt = computeAzAlt(moonGCRS, observer, ts, false); // Distance in km - const distance = Math.sqrt(moonGCRS[0] ** 2 + moonGCRS[1] ** 2 + moonGCRS[2] ** 2) + const distance = Math.sqrt(moonGCRS[0] ** 2 + moonGCRS[1] ** 2 + moonGCRS[2] ** 2); // Equatorial coordinates for parallactic angle - const RA_moon = Math.atan2(moonGCRS[1], moonGCRS[0]) - const dec_moon = Math.asin(Math.max(-1, Math.min(1, moonGCRS[2] / distance))) + const RA_moon = Math.atan2(moonGCRS[1], moonGCRS[0]); + const dec_moon = Math.asin(Math.max(-1, Math.min(1, moonGCRS[2] / distance))); // Hour angle: ERA(UT1) + longitude − right ascension - const era = computeERA(ts.jdUT1) - const HA = era + lon * DEG2RAD - RA_moon + const era = computeERA(ts.jdUT1); + const HA = era + lon * DEG2RAD - RA_moon; // Parallactic angle: signed angle between zenith and north pole as seen from the Moon const parallacticAngle = Math.atan2( Math.sin(HA), Math.cos(lat * DEG2RAD) * Math.tan(dec_moon) - Math.sin(lat * DEG2RAD) * Math.cos(HA), - ) + ); - return { azimuth: azAlt.azimuth, altitude: azAlt.altitude, distance, parallacticAngle } + return { azimuth: azAlt.azimuth, altitude: azAlt.altitude, distance, parallacticAngle }; } /** @@ -590,33 +591,33 @@ export function getMoonPosition( * ``` */ export function getMoonIllumination(date: Date = new Date()): MoonIlluminationResult { - validateDate(date, 'getMoonIllumination') - const ts = computeTimeScales(date) - const { moonGCRS, sunGCRS } = getMoonSunApproximate(ts.jdTT) + validateDate(date, "getMoonIllumination"); + const ts = computeTimeScales(date); + const { moonGCRS, sunGCRS } = getMoonSunApproximate(ts.jdTT); - const { illumination, elongationDeg, isWaxing } = computeIllumination(moonGCRS, sunGCRS) + const { illumination, elongationDeg, isWaxing } = computeIllumination(moonGCRS, sunGCRS); // Phase fraction: 0 = new moon, 0.25 = first quarter, 0.5 = full moon, 0.75 = last quarter - const phase = isWaxing ? elongationDeg / 360 : 1 - elongationDeg / 360 + const phase = isWaxing ? elongationDeg / 360 : 1 - elongationDeg / 360; // Position angle of the bright limb midpoint, measured eastward from north celestial pole. // PA = atan2(cos(dec_sun) * sin(RA_sun - RA_moon), // sin(dec_sun) * cos(dec_moon) - cos(dec_sun) * sin(dec_moon) * cos(RA_sun - RA_moon)) - const moonDist = Math.sqrt(moonGCRS[0] ** 2 + moonGCRS[1] ** 2 + moonGCRS[2] ** 2) - const sunDist = Math.sqrt(sunGCRS[0] ** 2 + sunGCRS[1] ** 2 + sunGCRS[2] ** 2) + const moonDist = Math.sqrt(moonGCRS[0] ** 2 + moonGCRS[1] ** 2 + moonGCRS[2] ** 2); + const sunDist = Math.sqrt(sunGCRS[0] ** 2 + sunGCRS[1] ** 2 + sunGCRS[2] ** 2); - const RA_moon = Math.atan2(moonGCRS[1], moonGCRS[0]) - const dec_moon = Math.asin(Math.max(-1, Math.min(1, moonGCRS[2] / moonDist))) - const RA_sun = Math.atan2(sunGCRS[1], sunGCRS[0]) - const dec_sun = Math.asin(Math.max(-1, Math.min(1, sunGCRS[2] / sunDist))) + const RA_moon = Math.atan2(moonGCRS[1], moonGCRS[0]); + const dec_moon = Math.asin(Math.max(-1, Math.min(1, moonGCRS[2] / moonDist))); + const RA_sun = Math.atan2(sunGCRS[1], sunGCRS[0]); + const dec_sun = Math.asin(Math.max(-1, Math.min(1, sunGCRS[2] / sunDist))); - const dRA = RA_sun - RA_moon + const dRA = RA_sun - RA_moon; const angle = Math.atan2( Math.cos(dec_sun) * Math.sin(dRA), Math.sin(dec_sun) * Math.cos(dec_moon) - Math.cos(dec_sun) * Math.sin(dec_moon) * Math.cos(dRA), - ) + ); - return { fraction: illumination, phase, angle, isWaxing } + return { fraction: illumination, phase, angle, isWaxing }; } /** @@ -650,54 +651,60 @@ export function getMoonVisibilityEstimate( lon: number, elevation = 0, ): MoonVisibilityEstimate { - validateDate(date, 'getMoonVisibilityEstimate') - validateLatitude(lat, 'getMoonVisibilityEstimate') - validateLongitude(lon, 'getMoonVisibilityEstimate') - const ts = computeTimeScales(date) - const { moonGCRS, sunGCRS } = getMoonSunApproximate(ts.jdTT) - const observer: Observer = { lat, lon, elevation } + validateDate(date, "getMoonVisibilityEstimate"); + validateLatitude(lat, "getMoonVisibilityEstimate"); + validateLongitude(lon, "getMoonVisibilityEstimate"); + const ts = computeTimeScales(date); + const { moonGCRS, sunGCRS } = getMoonSunApproximate(ts.jdTT); + const observer: Observer = { lat, lon, elevation }; // Airless positions — Odeh uses airless altitudes (no refraction) - const moonAirless = computeAzAlt(moonGCRS, observer, ts, true) - const sunAirless = computeAzAlt(sunGCRS, observer, ts, true) + const moonAirless = computeAzAlt(moonGCRS, observer, ts, true); + const sunAirless = computeAzAlt(sunGCRS, observer, ts, true); // ARCL = elongation (geocentric, degrees) - const { elongationDeg } = computeIllumination(moonGCRS, sunGCRS) - const ARCL = elongationDeg + const { elongationDeg } = computeIllumination(moonGCRS, sunGCRS); + const ARCL = elongationDeg; // ARCV = Moon airless altitude minus Sun airless altitude - const ARCV = moonAirless.altitude - sunAirless.altitude + const ARCV = moonAirless.altitude - sunAirless.altitude; // Topocentric Moon vector for crescent width - const obsECEF = geodeticToECEF(lat, lon, elevation) - const obsITRS: Vec3 = [obsECEF[0] / 1000, obsECEF[1] / 1000, obsECEF[2] / 1000] - const obsGCRS = itrsToGcrs(obsITRS, ts) + const obsECEF = geodeticToECEF(lat, lon, elevation); + const obsITRS: Vec3 = [obsECEF[0] / 1000, obsECEF[1] / 1000, obsECEF[2] / 1000]; + const obsGCRS = itrsToGcrs(obsITRS, ts); const moonTopo: Vec3 = [ moonGCRS[0] - obsGCRS[0], moonGCRS[1] - obsGCRS[1], moonGCRS[2] - obsGCRS[2], - ] + ]; - const { W } = computeCrescentWidth(moonTopo, ARCL) + const { W } = computeCrescentWidth(moonTopo, ARCL); // Odeh 2006: V = ARCV - arcv_minimum(W) - const V = ARCV - arcvMinimum(W) + const V = ARCV - arcvMinimum(W); const zone: OdehZone = - V >= ODEH_THRESHOLDS.A ? 'A' : V >= ODEH_THRESHOLDS.B ? 'B' : V >= ODEH_THRESHOLDS.C ? 'C' : 'D' + V >= ODEH_THRESHOLDS.A + ? "A" + : V >= ODEH_THRESHOLDS.B + ? "B" + : V >= ODEH_THRESHOLDS.C + ? "C" + : "D"; return { V, zone, description: ODEH_DESCRIPTIONS[zone], - isVisibleNakedEye: zone === 'A', - isVisibleWithOpticalAid: zone === 'A' || zone === 'B', + isVisibleNakedEye: zone === "A", + isVisibleWithOpticalAid: zone === "A" || zone === "B", ARCL, ARCV, W, moonAboveHorizon: moonAirless.altitude > 0, isApproximate: true, - } + }; } /** @@ -730,15 +737,15 @@ export function getMoon( lon: number, elevation = 0, ): MoonSnapshot { - validateDate(date, 'getMoon') - validateLatitude(lat, 'getMoon') - validateLongitude(lon, 'getMoon') + validateDate(date, "getMoon"); + validateLatitude(lat, "getMoon"); + validateLongitude(lon, "getMoon"); return { phase: getMoonPhase(date), position: getMoonPosition(date, lat, lon, elevation), illumination: getMoonIllumination(date), visibility: getMoonVisibilityEstimate(date, lat, lon, elevation), - } + }; } // ─── Internal helpers ───────────────────────────────────────────────────────── @@ -748,34 +755,34 @@ export function getMoon( * Full moon corrections differ from new moon; these are from Meeus Table 49.A. */ function nearestFullMoon(jdTT: number): number { - const Y = 2000 + (jdTT - J2000) / 365.25 - const kBase = Math.round((Y - 2000.0) * 12.3685) + const Y = 2000 + (jdTT - J2000) / 365.25; + const kBase = Math.round((Y - 2000.0) * 12.3685); // Check the full moons on either side of the nearest new moon (k ± 0.5) - const k1 = kBase - 0.5 - const k2 = kBase + 0.5 - const jde1 = fullMoonJDE(k1) - const jde2 = fullMoonJDE(k2) - const d1 = Math.abs(jde1 - jdTT) - const d2 = Math.abs(jde2 - jdTT) - return d1 < d2 ? jde1 : jde2 + const k1 = kBase - 0.5; + const k2 = kBase + 0.5; + const jde1 = fullMoonJDE(k1); + const jde2 = fullMoonJDE(k2); + const d1 = Math.abs(jde1 - jdTT); + const d2 = Math.abs(jde2 - jdTT); + return d1 < d2 ? jde1 : jde2; } /** Full moon JDE for a half-integer k (Meeus Ch. 49, Table 49.A). */ function fullMoonJDE(k: number): number { - const T = k / 1236.85 + const T = k / 1236.85; let JDE = 2451550.09766 + 29.530588861 * k + 0.00015437 * T * T - 0.00000015 * T * T * T + - 0.00000000073 * T * T * T * T + 0.00000000073 * T * T * T * T; - const M = (2.5534 + 29.1053567 * k - 0.0000014 * T * T) * DEG2RAD - const Mp = (201.5643 + 385.81693528 * k + 0.0107582 * T * T) * DEG2RAD - const Fc = (160.7108 + 390.67050284 * k - 0.0016118 * T * T) * DEG2RAD - const Om = (124.7746 - 1.56375588 * k + 0.0020672 * T * T) * DEG2RAD - const E = 1 - 0.002516 * T - 0.0000074 * T * T + const M = (2.5534 + 29.1053567 * k - 0.0000014 * T * T) * DEG2RAD; + const Mp = (201.5643 + 385.81693528 * k + 0.0107582 * T * T) * DEG2RAD; + const Fc = (160.7108 + 390.67050284 * k - 0.0016118 * T * T) * DEG2RAD; + const Om = (124.7746 - 1.56375588 * k + 0.0020672 * T * T) * DEG2RAD; + const E = 1 - 0.002516 * T - 0.0000074 * T * T; JDE += -0.40614 * Math.sin(Mp) + @@ -802,9 +809,9 @@ function fullMoonJDE(k: number): number { 0.00003 * Math.sin(Mp - M + 2 * Fc) - 0.00002 * Math.sin(Mp - M - 2 * Fc) - 0.00002 * Math.sin(3 * Mp + M) + - 0.00002 * Math.sin(4 * Mp) + 0.00002 * Math.sin(4 * Mp); - return JDE + return JDE; } /** @@ -812,12 +819,12 @@ function fullMoonJDE(k: number): number { * Boundaries: new (<5°), crescent (5-85°), quarter (85-95°), gibbous (95-175°), full (>175°). */ function elongationToPhase(elongationDeg: number, isWaxing: boolean): MoonPhaseName { - const e = elongationDeg - if (e < 5) return 'new-moon' - if (e > 175) return 'full-moon' - if (e < 85) return isWaxing ? 'waxing-crescent' : 'waning-crescent' - if (e < 95) return isWaxing ? 'first-quarter' : 'last-quarter' - return isWaxing ? 'waxing-gibbous' : 'waning-gibbous' + const e = elongationDeg; + if (e < 5) return "new-moon"; + if (e > 175) return "full-moon"; + if (e < 85) return isWaxing ? "waxing-crescent" : "waning-crescent"; + if (e < 95) return isWaxing ? "first-quarter" : "last-quarter"; + return isWaxing ? "waxing-gibbous" : "waning-gibbous"; } /** @@ -833,10 +840,10 @@ function elongationToPhase(elongationDeg: number, isWaxing: boolean): MoonPhaseN export async function getSunMoonEvents( date: Date, observer: Observer, - options?: Pick, + options?: Pick, ): Promise { - validateDate(date, 'getSunMoonEvents') - validateObserver(observer, 'getSunMoonEvents') - const kernel = await resolveKernel(options?.kernels) - return eventsGetSunMoonEvents(date, observer, kernel) + validateDate(date, "getSunMoonEvents"); + validateObserver(observer, "getSunMoonEvents"); + const kernel = await resolveKernel(options?.kernels); + return eventsGetSunMoonEvents(date, observer, kernel); } diff --git a/src/bodies/index.ts b/src/bodies/index.ts index 5dd9f22..f6a3f7b 100644 --- a/src/bodies/index.ts +++ b/src/bodies/index.ts @@ -16,21 +16,21 @@ * New Crescent Moon. NAO Technical Note 69. HM Nautical Almanac Office. */ -import type { StateVector, Vec3 } from '../types.js' -import type { SpkKernel } from '../spk/index.js' -import { NAIF_IDS } from '../spk/index.js' -import { J2000, DAYS_PER_JULIAN_CENTURY } from '../time/index.js' -import { DEG2RAD, vdot, vnorm } from '../math/index.js' +import type { StateVector, Vec3 } from "../types.js"; +import type { SpkKernel } from "../spk/index.js"; +import { NAIF_IDS } from "../spk/index.js"; +import { J2000, DAYS_PER_JULIAN_CENTURY } from "../time/index.js"; +import { DEG2RAD, vdot, vnorm } from "../math/index.js"; // ─── Constants ──────────────────────────────────────────────────────────────── -const AU_KM = 149597870.7 +const AU_KM = 149597870.7; /** Mean radius of the Moon in km (IAU 2015 nominal value) */ -const MOON_RADIUS_KM = 1737.4 +const MOON_RADIUS_KM = 1737.4; /** Mean radius of the Sun in km */ -const _SUN_RADIUS_KM = 696000.0 -void _SUN_RADIUS_KM // reserved for future solar semi-diameter calculations +const _SUN_RADIUS_KM = 696000.0; +void _SUN_RADIUS_KM; // reserved for future solar semi-diameter calculations // ─── Geocentric state ───────────────────────────────────────────────────────── @@ -45,7 +45,7 @@ void _SUN_RADIUS_KM // reserved for future solar semi-diameter calculations * @returns Moon state vector relative to Earth center, km and km/s, GCRS */ export function getMoonGeocentricState(kernel: SpkKernel, et: number): StateVector { - return kernel.getState(NAIF_IDS.MOON, NAIF_IDS.EARTH, et) + return kernel.getState(NAIF_IDS.MOON, NAIF_IDS.EARTH, et); } /** @@ -59,7 +59,7 @@ export function getMoonGeocentricState(kernel: SpkKernel, et: number): StateVect * @returns Sun state vector relative to Earth center, km and km/s, GCRS */ export function getSunGeocentricState(kernel: SpkKernel, et: number): StateVector { - return kernel.getState(NAIF_IDS.SUN, NAIF_IDS.EARTH, et) + return kernel.getState(NAIF_IDS.SUN, NAIF_IDS.EARTH, et); } // ─── Moon illumination ──────────────────────────────────────────────────────── @@ -85,12 +85,12 @@ export function computeIllumination( moonGCRS: Vec3, sunGCRS: Vec3, ): { illumination: number; phaseAngleDeg: number; elongationDeg: number; isWaxing: boolean } { - const rMoon = vnorm(moonGCRS) - const rSun = vnorm(sunGCRS) + const rMoon = vnorm(moonGCRS); + const rSun = vnorm(sunGCRS); // Elongation ψ: angle at Earth between Moon and Sun - const cosElong = vdot(moonGCRS, sunGCRS) / (rMoon * rSun) - const elongationDeg = Math.acos(Math.max(-1, Math.min(1, cosElong))) / DEG2RAD + const cosElong = vdot(moonGCRS, sunGCRS) / (rMoon * rSun); + const elongationDeg = Math.acos(Math.max(-1, Math.min(1, cosElong))) / DEG2RAD; // Phase angle i: angle at Moon between Earth and Sun // Vector from Moon to Earth: -moonGCRS @@ -99,21 +99,21 @@ export function computeIllumination( sunGCRS[0] - moonGCRS[0], sunGCRS[1] - moonGCRS[1], sunGCRS[2] - moonGCRS[2], - ] - const moonToEarth: Vec3 = [-moonGCRS[0], -moonGCRS[1], -moonGCRS[2]] - const rMoonToSun = vnorm(moonToSun) + ]; + const moonToEarth: Vec3 = [-moonGCRS[0], -moonGCRS[1], -moonGCRS[2]]; + const rMoonToSun = vnorm(moonToSun); - const cosPhase = vdot(moonToEarth, moonToSun) / (rMoon * rMoonToSun) - const phaseAngleDeg = Math.acos(Math.max(-1, Math.min(1, cosPhase))) / DEG2RAD + const cosPhase = vdot(moonToEarth, moonToSun) / (rMoon * rMoonToSun); + const phaseAngleDeg = Math.acos(Math.max(-1, Math.min(1, cosPhase))) / DEG2RAD; - const illumination = (1 + Math.cos(phaseAngleDeg * DEG2RAD)) / 2 + const illumination = (1 + Math.cos(phaseAngleDeg * DEG2RAD)) / 2; // Moon is waxing when it is east of the Sun (elongation increasing). // Cross product sunGCRS × moonGCRS z-component: positive when Moon is east of Sun. - const crossZ = sunGCRS[0] * moonGCRS[1] - sunGCRS[1] * moonGCRS[0] - const isWaxing = crossZ > 0 + const crossZ = sunGCRS[0] * moonGCRS[1] - sunGCRS[1] * moonGCRS[0]; + const isWaxing = crossZ > 0; - return { illumination, phaseAngleDeg, elongationDeg, isWaxing } + return { illumination, phaseAngleDeg, elongationDeg, isWaxing }; } /** @@ -140,17 +140,17 @@ export function computeCrescentWidth( moonTopoVec: Vec3, ARCL: number, ): { W: number; Wprime: number } { - const rMoon = Math.sqrt(moonTopoVec[0] ** 2 + moonTopoVec[1] ** 2 + moonTopoVec[2] ** 2) + const rMoon = Math.sqrt(moonTopoVec[0] ** 2 + moonTopoVec[1] ** 2 + moonTopoVec[2] ** 2); // Topocentric semi-diameter in arc minutes - const SDmoon_arcmin = (Math.atan(MOON_RADIUS_KM / rMoon) / DEG2RAD) * 60 + const SDmoon_arcmin = (Math.atan(MOON_RADIUS_KM / rMoon) / DEG2RAD) * 60; // Crescent width in arc minutes - const ARCL_rad = ARCL * DEG2RAD - const W = SDmoon_arcmin * (1 - Math.cos(ARCL_rad)) + const ARCL_rad = ARCL * DEG2RAD; + const W = SDmoon_arcmin * (1 - Math.cos(ARCL_rad)); // Wprime ≡ W for both Odeh and Yallop in this formulation - return { W, Wprime: W } + return { W, Wprime: W }; } // ─── Approximate positions (no kernel) ──────────────────────────────────────── @@ -169,48 +169,48 @@ export function computeCrescentWidth( * @returns Geocentric GCRS positions in km (approximate, light-time not corrected) */ export function getMoonSunApproximate(jdTT: number): { - moonGCRS: Vec3 - sunGCRS: Vec3 + moonGCRS: Vec3; + sunGCRS: Vec3; } { - const T = (jdTT - J2000) / DAYS_PER_JULIAN_CENTURY + const T = (jdTT - J2000) / DAYS_PER_JULIAN_CENTURY; // ── Sun (Meeus Ch. 25) ────────────────────────────────────────────────────── // Mean longitude L0 and mean anomaly M (degrees) - const L0 = 280.46646 + 36000.76983 * T + 0.0003032 * T * T - const M_sun = 357.52911 + 35999.05029 * T - 0.0001537 * T * T - const M_sun_rad = (M_sun % 360) * DEG2RAD + const L0 = 280.46646 + 36000.76983 * T + 0.0003032 * T * T; + const M_sun = 357.52911 + 35999.05029 * T - 0.0001537 * T * T; + const M_sun_rad = (M_sun % 360) * DEG2RAD; - const e_sun = 0.016708634 - 0.000042037 * T - 0.0000001267 * T * T + const e_sun = 0.016708634 - 0.000042037 * T - 0.0000001267 * T * T; // Equation of center (degrees) const C = (1.914602 - 0.004817 * T - 0.000014 * T * T) * Math.sin(M_sun_rad) + (0.019993 - 0.000101 * T) * Math.sin(2 * M_sun_rad) + - 0.000289 * Math.sin(3 * M_sun_rad) + 0.000289 * Math.sin(3 * M_sun_rad); // True longitude and anomaly - const sunLonDeg = L0 + C - const nu_rad = M_sun_rad + C * DEG2RAD + const sunLonDeg = L0 + C; + const nu_rad = M_sun_rad + C * DEG2RAD; // Geometric distance in AU - const R_AU = (1.000001018 * (1 - e_sun * e_sun)) / (1 + e_sun * Math.cos(nu_rad)) - const R_km = R_AU * AU_KM + const R_AU = (1.000001018 * (1 - e_sun * e_sun)) / (1 + e_sun * Math.cos(nu_rad)); + const R_km = R_AU * AU_KM; // Nutation correction for apparent longitude (simplified) - const omega = (125.04 - 1934.136 * T) * DEG2RAD - const sunLonApp = sunLonDeg - 0.00569 - 0.00478 * Math.sin(omega) - const sunLon_rad = sunLonApp * DEG2RAD + const omega = (125.04 - 1934.136 * T) * DEG2RAD; + const sunLonApp = sunLonDeg - 0.00569 - 0.00478 * Math.sin(omega); + const sunLon_rad = sunLonApp * DEG2RAD; // Mean obliquity of the ecliptic (IAU 1980 approximation, degrees) const eps = - (23.439291111 - 0.013004167 * T - 0.0000001638 * T * T + 0.0000005036 * T * T * T) * DEG2RAD + (23.439291111 - 0.013004167 * T - 0.0000001638 * T * T + 0.0000005036 * T * T * T) * DEG2RAD; const sunGCRS: Vec3 = [ R_km * Math.cos(sunLon_rad), R_km * Math.sin(sunLon_rad) * Math.cos(eps), R_km * Math.sin(sunLon_rad) * Math.sin(eps), - ] + ]; // ── Moon (Meeus Ch. 47) ───────────────────────────────────────────────────── @@ -220,40 +220,40 @@ export function getMoonSunApproximate(jdTT: number): { 481267.88123421 * T - 0.0015786 * T * T + (T * T * T) / 538841 - - (T * T * T * T) / 65194000 + (T * T * T * T) / 65194000; const D = 297.8501921 + 445267.1114034 * T - 0.0018819 * T * T + (T * T * T) / 545868 - - (T * T * T * T) / 113065000 - const M = 357.5291092 + 35999.0502909 * T - 0.0001536 * T * T + (T * T * T) / 24490000 + (T * T * T * T) / 113065000; + const M = 357.5291092 + 35999.0502909 * T - 0.0001536 * T * T + (T * T * T) / 24490000; const Mp = 134.9633964 + 477198.8675055 * T + 0.0087414 * T * T + (T * T * T) / 69699 - - (T * T * T * T) / 14712000 + (T * T * T * T) / 14712000; const F = 93.272095 + 483202.0175233 * T - 0.0036539 * T * T - (T * T * T) / 3526000 + - (T * T * T * T) / 863310000 + (T * T * T * T) / 863310000; // Additive terms for longitude/latitude - const A1 = (119.75 + 131.849 * T) * DEG2RAD - const A2 = (53.09 + 479264.29 * T) * DEG2RAD - const A3 = (313.45 + 481266.484 * T) * DEG2RAD + const A1 = (119.75 + 131.849 * T) * DEG2RAD; + const A2 = (53.09 + 479264.29 * T) * DEG2RAD; + const A3 = (313.45 + 481266.484 * T) * DEG2RAD; // Convert to radians for accumulation - const D_r = (D % 360) * DEG2RAD - const M_r = (M % 360) * DEG2RAD - const Mp_r = (Mp % 360) * DEG2RAD - const F_r = (F % 360) * DEG2RAD + const D_r = (D % 360) * DEG2RAD; + const M_r = (M % 360) * DEG2RAD; + const Mp_r = (Mp % 360) * DEG2RAD; + const F_r = (F % 360) * DEG2RAD; // Eccentricity correction for terms involving M (Earth's orbital eccentricity) - const E = 1 - 0.002516 * T - 0.0000074 * T * T + const E = 1 - 0.002516 * T - 0.0000074 * T * T; // Longitude and distance accumulation — 30 main terms from Meeus Table 47.A // [d, m, mp, f, Σl (0.000001°), Σr (0.001 km)] @@ -288,19 +288,19 @@ export function getMoonSunApproximate(jdTT: number): { [0, 1, -2, 0, -2689, -7003], [2, 0, -1, 2, -2602, 0], [2, -1, -2, 0, 2390, 10056], - ] + ]; let Sl = 0, - Sr = 0 + Sr = 0; for (const [d, m, mp, f, sl, sr] of LD) { - const arg = d * D_r + m * M_r + mp * Mp_r + f * F_r - const eCorr = Math.abs(m) === 2 ? E * E : Math.abs(m) === 1 ? E : 1 - Sl += sl * eCorr * Math.sin(arg) - Sr += sr * eCorr * Math.cos(arg) + const arg = d * D_r + m * M_r + mp * Mp_r + f * F_r; + const eCorr = Math.abs(m) === 2 ? E * E : Math.abs(m) === 1 ? E : 1; + Sl += sl * eCorr * Math.sin(arg); + Sr += sr * eCorr * Math.cos(arg); } // Additive longitude corrections (Meeus §47) - Sl += 3958 * Math.sin(A1) + 1962 * Math.sin((Lp - F) * DEG2RAD) + 318 * Math.sin(A2) + Sl += 3958 * Math.sin(A1) + 1962 * Math.sin((Lp - F) * DEG2RAD) + 318 * Math.sin(A2); // Latitude accumulation — 20 main terms from Meeus Table 47.B // [d, m, mp, f, Σb (0.000001°)] @@ -325,13 +325,13 @@ export function getMoonSunApproximate(jdTT: number): { [0, 1, -1, -1, -1870], [4, 0, -1, -1, 1828], [0, 1, 0, 1, -1794], - ] + ]; - let Sb = 0 + let Sb = 0; for (const [d, m, mp, f, sb] of FB) { - const arg = d * D_r + m * M_r + mp * Mp_r + f * F_r - const eCorr = Math.abs(m) === 2 ? E * E : Math.abs(m) === 1 ? E : 1 - Sb += sb * eCorr * Math.sin(arg) + const arg = d * D_r + m * M_r + mp * Mp_r + f * F_r; + const eCorr = Math.abs(m) === 2 ? E * E : Math.abs(m) === 1 ? E : 1; + Sb += sb * eCorr * Math.sin(arg); } // Additive latitude corrections @@ -341,15 +341,15 @@ export function getMoonSunApproximate(jdTT: number): { 175 * Math.sin(A1 - F_r) + 175 * Math.sin(A1 + F_r) + 127 * Math.sin((Lp - Mp) * DEG2RAD) - - 115 * Math.sin((Lp + Mp) * DEG2RAD) + 115 * Math.sin((Lp + Mp) * DEG2RAD); // Moon ecliptic coordinates - const moonLonDeg = Lp + Sl * 1e-6 - const moonLatDeg = Sb * 1e-6 - const moonDistKm = 385000.56 + Sr * 0.001 + const moonLonDeg = Lp + Sl * 1e-6; + const moonLatDeg = Sb * 1e-6; + const moonDistKm = 385000.56 + Sr * 0.001; - const moonLon_rad = moonLonDeg * DEG2RAD - const moonLat_rad = moonLatDeg * DEG2RAD + const moonLon_rad = moonLonDeg * DEG2RAD; + const moonLat_rad = moonLatDeg * DEG2RAD; // Ecliptic to equatorial (GCRS ≈ J2000 equatorial for this accuracy level) const moonGCRS: Vec3 = [ @@ -360,9 +360,9 @@ export function getMoonSunApproximate(jdTT: number): { moonDistKm * (Math.sin(eps) * Math.cos(moonLat_rad) * Math.sin(moonLon_rad) + Math.cos(eps) * Math.sin(moonLat_rad)), - ] + ]; - return { moonGCRS, sunGCRS } + return { moonGCRS, sunGCRS }; } /** @@ -374,11 +374,11 @@ export function getMoonSunApproximate(jdTT: number): { */ export function nearestNewMoon(jdTT: number): number { // Convert JD to approximate decimal year - const Y = 2000.0 + (jdTT - J2000) / 365.25 + const Y = 2000.0 + (jdTT - J2000) / 365.25; // k = approximate lunation number (0 = Jan 6, 2000 new moon) - const k = Math.round((Y - 2000.0) * 12.3685) - const T = k / 1236.85 + const k = Math.round((Y - 2000.0) * 12.3685); + const T = k / 1236.85; // JDE of mean new moon (Meeus Eq. 49.1) let JDE = @@ -386,16 +386,16 @@ export function nearestNewMoon(jdTT: number): number { 29.530588861 * k + 0.00015437 * T * T - 0.00000015 * T * T * T + - 0.00000000073 * T * T * T * T + 0.00000000073 * T * T * T * T; // Fundamental arguments for the corrections (degrees → radians) - const M = (2.5534 + 29.1053567 * k - 0.0000014 * T * T - 0.00000011 * T * T * T) * DEG2RAD - const Mp = (201.5643 + 385.81693528 * k + 0.0107582 * T * T + 0.00001238 * T * T * T) * DEG2RAD - const Fc = (160.7108 + 390.67050284 * k - 0.0016118 * T * T - 0.00000227 * T * T * T) * DEG2RAD - const Om = (124.7746 - 1.56375588 * k + 0.0020672 * T * T + 0.00000215 * T * T * T) * DEG2RAD + const M = (2.5534 + 29.1053567 * k - 0.0000014 * T * T - 0.00000011 * T * T * T) * DEG2RAD; + const Mp = (201.5643 + 385.81693528 * k + 0.0107582 * T * T + 0.00001238 * T * T * T) * DEG2RAD; + const Fc = (160.7108 + 390.67050284 * k - 0.0016118 * T * T - 0.00000227 * T * T * T) * DEG2RAD; + const Om = (124.7746 - 1.56375588 * k + 0.0020672 * T * T + 0.00000215 * T * T * T) * DEG2RAD; // Eccentricity of Earth's orbit - const E = 1 - 0.002516 * T - 0.0000074 * T * T + const E = 1 - 0.002516 * T - 0.0000074 * T * T; // Corrections from Meeus Table 49.A (new moon) JDE += @@ -423,7 +423,7 @@ export function nearestNewMoon(jdTT: number): number { 0.00003 * Math.sin(Mp - M + 2 * Fc) - 0.00002 * Math.sin(Mp - M - 2 * Fc) - 0.00002 * Math.sin(3 * Mp + M) + - 0.00002 * Math.sin(4 * Mp) + 0.00002 * Math.sin(4 * Mp); - return JDE + return JDE; } diff --git a/src/cli/index.ts b/src/cli/index.ts index d5f0cdf..6b797c2 100644 --- a/src/cli/index.ts +++ b/src/cli/index.ts @@ -15,31 +15,31 @@ import { verifyKernels, getMoonSightingReport, getMoonPhase, -} from '../api/index.js' +} from "../api/index.js"; -const args = process.argv.slice(2) -const command = args[0] +const args = process.argv.slice(2); +const command = args[0]; async function main() { switch (command) { - case 'download-kernels': - await cmdDownloadKernels() - break - case 'verify-kernels': - await cmdVerifyKernels() - break - case 'sighting': - await cmdSighting(args.slice(1)) - break - case 'phase': - cmdPhase(args[1]) - break - case 'benchmark': - await cmdBenchmark() - break + case "download-kernels": + await cmdDownloadKernels(); + break; + case "verify-kernels": + await cmdVerifyKernels(); + break; + case "sighting": + await cmdSighting(args.slice(1)); + break; + case "phase": + cmdPhase(args[1]); + break; + case "benchmark": + await cmdBenchmark(); + break; default: - printHelp() - process.exit(command ? 1 : 0) + printHelp(); + process.exit(command ? 1 : 0); } } @@ -57,130 +57,130 @@ Examples: moon-sighting download-kernels moon-sighting sighting 51.5 -0.1 2025-03-29 moon-sighting sighting 21.4 39.8 # Mecca - moon-sighting phase 2025-03-01`) + moon-sighting phase 2025-03-01`); } async function cmdDownloadKernels() { - await downloadKernels() - console.log('Kernels ready.') + await downloadKernels(); + console.log("Kernels ready."); } async function cmdVerifyKernels() { - const result = await verifyKernels() + const result = await verifyKernels(); if (result.ok) { - console.log('Kernels OK.') + console.log("Kernels OK."); } else { - for (const err of result.errors) console.error(err) - process.exit(1) + for (const err of result.errors) console.error(err); + process.exit(1); } } async function cmdSighting(cmdArgs: string[]) { - const lat = parseFloat(cmdArgs[0] ?? '') - const lon = parseFloat(cmdArgs[1] ?? '') - const dateStr = cmdArgs[2] ?? new Date().toISOString().slice(0, 10) + const lat = parseFloat(cmdArgs[0] ?? ""); + const lon = parseFloat(cmdArgs[1] ?? ""); + const dateStr = cmdArgs[2] ?? new Date().toISOString().slice(0, 10); if (isNaN(lat) || isNaN(lon)) { - console.error('Usage: moon-sighting sighting [YYYY-MM-DD]') - process.exit(1) + console.error("Usage: moon-sighting sighting [YYYY-MM-DD]"); + process.exit(1); } - const date = new Date(`${dateStr}T00:00:00Z`) + const date = new Date(`${dateStr}T00:00:00Z`); if (isNaN(date.getTime())) { - console.error(`Invalid date: ${dateStr}. Use YYYY-MM-DD format.`) - process.exit(1) + console.error(`Invalid date: ${dateStr}. Use YYYY-MM-DD format.`); + process.exit(1); } - console.log(`Computing sighting report for ${lat}°N ${lon}°E on ${dateStr}...`) - await initKernels() + console.log(`Computing sighting report for ${lat}°N ${lon}°E on ${dateStr}...`); + await initKernels(); - const report = await getMoonSightingReport(date, { lat, lon, elevation: 0 }) + const report = await getMoonSightingReport(date, { lat, lon, elevation: 0 }); - console.log('') - console.log(`Sunset: ${fmtDate(report.sunsetUTC)}`) - console.log(`Moonset: ${fmtDate(report.moonsetUTC)}`) - console.log(`Best time: ${fmtDate(report.bestTimeUTC)}`) + console.log(""); + console.log(`Sunset: ${fmtDate(report.sunsetUTC)}`); + console.log(`Moonset: ${fmtDate(report.moonsetUTC)}`); + console.log(`Best time: ${fmtDate(report.bestTimeUTC)}`); console.log( - `Lag: ${report.lagMinutes !== null ? Math.round(report.lagMinutes) + ' min' : 'N/A'}`, - ) - console.log('') + `Lag: ${report.lagMinutes !== null ? Math.round(report.lagMinutes) + " min" : "N/A"}`, + ); + console.log(""); if (report.geometry) { - const g = report.geometry - console.log(`ARCL: ${g.ARCL.toFixed(2)}°`) - console.log(`ARCV: ${g.ARCV.toFixed(2)}°`) - console.log(`DAZ: ${g.DAZ.toFixed(2)}°`) - console.log(`W: ${g.W.toFixed(3)} arcmin`) + const g = report.geometry; + console.log(`ARCL: ${g.ARCL.toFixed(2)}°`); + console.log(`ARCV: ${g.ARCV.toFixed(2)}°`); + console.log(`DAZ: ${g.DAZ.toFixed(2)}°`); + console.log(`W: ${g.W.toFixed(3)} arcmin`); } if (report.yallop && report.odeh) { - console.log('') - console.log(`Yallop: ${report.yallop.category} — ${report.yallop.description}`) - console.log(`Odeh: ${report.odeh.zone} — ${report.odeh.description}`) + console.log(""); + console.log(`Yallop: ${report.yallop.category} — ${report.yallop.description}`); + console.log(`Odeh: ${report.odeh.zone} — ${report.odeh.description}`); } - console.log('') - console.log(report.guidance) + console.log(""); + console.log(report.guidance); } function cmdPhase(dateStr?: string) { - const date = dateStr ? new Date(`${dateStr}T00:00:00Z`) : new Date() + const date = dateStr ? new Date(`${dateStr}T00:00:00Z`) : new Date(); if (isNaN(date.getTime())) { - console.error(`Invalid date: ${dateStr}. Use YYYY-MM-DD format.`) - process.exit(1) + console.error(`Invalid date: ${dateStr}. Use YYYY-MM-DD format.`); + process.exit(1); } - const phase = getMoonPhase(date) + const phase = getMoonPhase(date); - console.log(`Moon phase for ${date.toISOString().slice(0, 10)}:`) - console.log(` Phase: ${phase.phase}`) - console.log(` Illumination: ${phase.illumination.toFixed(1)}%`) - console.log(` Age: ${phase.age.toFixed(1)} hours`) - console.log(` Elongation: ${phase.elongationDeg.toFixed(1)}°`) - console.log(` Waxing: ${phase.isWaxing}`) - console.log(` Prev new: ${phase.prevNewMoon.toISOString().slice(0, 16)} UTC`) - console.log(` Next new: ${phase.nextNewMoon.toISOString().slice(0, 16)} UTC`) - console.log(` Next full: ${phase.nextFullMoon.toISOString().slice(0, 16)} UTC`) + console.log(`Moon phase for ${date.toISOString().slice(0, 10)}:`); + console.log(` Phase: ${phase.phase}`); + console.log(` Illumination: ${phase.illumination.toFixed(1)}%`); + console.log(` Age: ${phase.age.toFixed(1)} hours`); + console.log(` Elongation: ${phase.elongationDeg.toFixed(1)}°`); + console.log(` Waxing: ${phase.isWaxing}`); + console.log(` Prev new: ${phase.prevNewMoon.toISOString().slice(0, 16)} UTC`); + console.log(` Next new: ${phase.nextNewMoon.toISOString().slice(0, 16)} UTC`); + console.log(` Next full: ${phase.nextFullMoon.toISOString().slice(0, 16)} UTC`); } async function cmdBenchmark() { - console.log('moon-sighting benchmark\n') + console.log("moon-sighting benchmark\n"); // Benchmark 1: getMoonPhase (no kernel needed) - const N_PHASE = 10000 - const phaseStart = performance.now() + const N_PHASE = 10000; + const phaseStart = performance.now(); for (let i = 0; i < N_PHASE; i++) { - getMoonPhase(new Date(Date.UTC(2025, 2, 1 + (i % 28)))) + getMoonPhase(new Date(Date.UTC(2025, 2, 1 + (i % 28)))); } - const phaseMs = performance.now() - phaseStart + const phaseMs = performance.now() - phaseStart; console.log( `getMoonPhase × ${N_PHASE}: ${phaseMs.toFixed(1)} ms (${((phaseMs / N_PHASE) * 1000).toFixed(1)} µs/call)`, - ) + ); // Benchmark 2: kernel load - const loadStart = performance.now() - await initKernels() - const loadMs = performance.now() - loadStart - console.log(`initKernels (cold/cached): ${loadMs.toFixed(1)} ms`) + const loadStart = performance.now(); + await initKernels(); + const loadMs = performance.now() - loadStart; + console.log(`initKernels (cold/cached): ${loadMs.toFixed(1)} ms`); // Benchmark 3: single sighting report - const observer = { lat: 51.5074, lon: -0.1278, elevation: 10 } - const reportStart = performance.now() - await getMoonSightingReport(new Date('2025-03-29T00:00:00Z'), observer) - const reportMs = performance.now() - reportStart - console.log(`getMoonSightingReport (single): ${reportMs.toFixed(1)} ms`) + const observer = { lat: 51.5074, lon: -0.1278, elevation: 10 }; + const reportStart = performance.now(); + await getMoonSightingReport(new Date("2025-03-29T00:00:00Z"), observer); + const reportMs = performance.now() - reportStart; + console.log(`getMoonSightingReport (single): ${reportMs.toFixed(1)} ms`); } /** Format a nullable Date as a short UTC string. */ function fmtDate(d: Date | null): string { - if (!d) return 'N/A' + if (!d) return "N/A"; return d .toISOString() - .replace('T', ' ') - .replace(/\.\d+Z$/, ' UTC') + .replace("T", " ") + .replace(/\.\d+Z$/, " UTC"); } main().catch((err) => { - console.error(err instanceof Error ? err.message : String(err)) - process.exit(1) -}) + console.error(err instanceof Error ? err.message : String(err)); + process.exit(1); +}); diff --git a/src/events/index.ts b/src/events/index.ts index 8997998..b19ebfb 100644 --- a/src/events/index.ts +++ b/src/events/index.ts @@ -19,11 +19,11 @@ * Reference: Yallop (1997) NAO TN 69 §2.4; Odeh (2006) §3 */ -import type { Vec3, Observer, SunMoonEvents, TimeScales } from '../types.js' -import type { SpkKernel } from '../spk/index.js' -import { NAIF_IDS } from '../spk/index.js' -import { brentRoot, vdot, vnorm } from '../math/index.js' -import { arcvMinimum } from '../visibility/index.js' +import type { Vec3, Observer, SunMoonEvents, TimeScales } from "../types.js"; +import type { SpkKernel } from "../spk/index.js"; +import { NAIF_IDS } from "../spk/index.js"; +import { brentRoot, vdot, vnorm } from "../math/index.js"; +import { arcvMinimum } from "../visibility/index.js"; import { J2000, SECONDS_PER_DAY, @@ -33,14 +33,14 @@ import { jdTTtoET, getDeltaAT, TT_MINUS_TAI, -} from '../time/index.js' +} from "../time/index.js"; import { getMoonGeocentricState, getSunGeocentricState, computeCrescentWidth, -} from '../bodies/index.js' -import { geodeticToECEF, computeAzAlt } from '../observer/index.js' -import { itrsToGcrs } from '../frames/index.js' +} from "../bodies/index.js"; +import { geodeticToECEF, computeAzAlt } from "../observer/index.js"; +import { itrsToGcrs } from "../frames/index.js"; // ─── Altitude threshold constants ───────────────────────────────────────────── @@ -49,14 +49,14 @@ import { itrsToGcrs } from '../frames/index.js' * Accounts for: standard refraction at horizon (34') + solar semi-diameter (16') * Total: −50' = −0.8333° */ -export const SUN_ALTITUDE_THRESHOLD = -0.8333 +export const SUN_ALTITUDE_THRESHOLD = -0.8333; /** * Standard threshold altitude for moonset/moonrise. * Accounts for: standard refraction at horizon (34') + lunar semi-diameter (~16') * Note: Moon's SD varies with distance (14.7'–16.8'). Use 0.2725° as mean. */ -export const MOON_ALTITUDE_THRESHOLD = -0.8333 // refined per actual distance in implementation +export const MOON_ALTITUDE_THRESHOLD = -0.8333; // refined per actual distance in implementation // ─── Internal helpers ───────────────────────────────────────────────────────── @@ -66,14 +66,14 @@ export const MOON_ALTITUDE_THRESHOLD = -0.8333 // refined per actual distance in */ function etToTS(et: number): TimeScales { // ET ≈ (jdTT - J2000) × 86400, so jdTT ≈ J2000 + et / 86400 - const jdTT = J2000 + et / SECONDS_PER_DAY + const jdTT = J2000 + et / SECONDS_PER_DAY; // Approximate UTC: TT - UTC ≈ deltaAT + 32.184s (typically ~69s currently) // Use a rough correction; getDeltaAT needs a UTC JD, so iterate once - const jdUTC_est = jdTT - 70.0 / SECONDS_PER_DAY - const deltaAT = getDeltaAT(jdUTC_est) - const jdUTC = jdTT - (deltaAT + TT_MINUS_TAI) / SECONDS_PER_DAY - const utc = jdToDate(jdUTC) - return computeTimeScales(utc) + const jdUTC_est = jdTT - 70.0 / SECONDS_PER_DAY; + const deltaAT = getDeltaAT(jdUTC_est); + const jdUTC = jdTT - (deltaAT + TT_MINUS_TAI) / SECONDS_PER_DAY; + const utc = jdToDate(jdUTC); + return computeTimeScales(utc); } /** @@ -81,9 +81,9 @@ function etToTS(et: number): TimeScales { */ function bodyPositionAt(kernel: SpkKernel, naifId: number, et: number): Vec3 { if (naifId === NAIF_IDS.SUN) { - return getSunGeocentricState(kernel, et).position + return getSunGeocentricState(kernel, et).position; } - return getMoonGeocentricState(kernel, et).position + return getMoonGeocentricState(kernel, et).position; } /** @@ -97,10 +97,10 @@ function altitudeMinusThreshold( et: number, threshold: number, ): number { - const ts = etToTS(et) - const bodyGCRS = bodyPositionAt(kernel, naifId, et) - const azAlt = computeAzAlt(bodyGCRS, observer, ts, true) - return azAlt.altitude - threshold + const ts = etToTS(et); + const bodyGCRS = bodyPositionAt(kernel, naifId, et); + const azAlt = computeAzAlt(bodyGCRS, observer, ts, true); + return azAlt.altitude - threshold; } // ─── Event finding ──────────────────────────────────────────────────────────── @@ -131,36 +131,36 @@ export function findAltitudeCrossing( threshold: number, rising: boolean, ): Date | null { - void ts // ts not needed — etToTS computes from ET directly + void ts; // ts not needed — etToTS computes from ET directly - const f = (et: number) => altitudeMinusThreshold(kernel, naifId, observer, et, threshold) + const f = (et: number) => altitudeMinusThreshold(kernel, naifId, observer, et, threshold); - const STEP_S = 600 // 10-minute coarse sampling - const nSteps = Math.ceil((endET - startET) / STEP_S) + const STEP_S = 600; // 10-minute coarse sampling + const nSteps = Math.ceil((endET - startET) / STEP_S); - let prevET = startET - let prevF = f(prevET) + let prevET = startET; + let prevF = f(prevET); for (let i = 1; i <= nSteps; i++) { - const currET = Math.min(startET + i * STEP_S, endET) - const currF = f(currET) + const currET = Math.min(startET + i * STEP_S, endET); + const currF = f(currET); - const isRisingCross = rising && prevF < 0 && currF >= 0 - const isSettingCross = !rising && prevF >= 0 && currF < 0 + const isRisingCross = rising && prevF < 0 && currF >= 0; + const isSettingCross = !rising && prevF >= 0 && currF < 0; if (isRisingCross || isSettingCross) { - const etCross = brentRoot(f, prevET, currET, 0.5) // 0.5s precision + const etCross = brentRoot(f, prevET, currET, 0.5); // 0.5s precision if (etCross !== null) { - const tsCross = etToTS(etCross) - return tsCross.utc + const tsCross = etToTS(etCross); + return tsCross.utc; } } - prevET = currET - prevF = currF + prevET = currET; + prevF = currF; } - return null + return null; } /** @@ -176,13 +176,13 @@ export function findAltitudeCrossing( */ export function getSunMoonEvents(date: Date, observer: Observer, kernel: SpkKernel): SunMoonEvents { // Anchor search at UTC midnight of the input date - const midnight = new Date(Date.UTC(date.getUTCFullYear(), date.getUTCMonth(), date.getUTCDate())) - const jdMidnight = dateToJD(midnight) + const midnight = new Date(Date.UTC(date.getUTCFullYear(), date.getUTCMonth(), date.getUTCDate())); + const jdMidnight = dateToJD(midnight); // Approximate ET at midnight - const etStart = jdTTtoET(jdMidnight + 70.0 / SECONDS_PER_DAY) // rough TT≈UTC+70s - const etEnd = etStart + 28 * 3600 // 28-hour window + const etStart = jdTTtoET(jdMidnight + 70.0 / SECONDS_PER_DAY); // rough TT≈UTC+70s + const etEnd = etStart + 28 * 3600; // 28-hour window - const ts0 = computeTimeScales(midnight) + const ts0 = computeTimeScales(midnight); // Sun events const sunriseUTC = findAltitudeCrossing( @@ -194,7 +194,7 @@ export function getSunMoonEvents(date: Date, observer: Observer, kernel: SpkKern etEnd, SUN_ALTITUDE_THRESHOLD, true, - ) + ); const sunsetUTC = findAltitudeCrossing( kernel, NAIF_IDS.SUN, @@ -204,7 +204,7 @@ export function getSunMoonEvents(date: Date, observer: Observer, kernel: SpkKern etEnd, SUN_ALTITUDE_THRESHOLD, false, - ) + ); // Twilight events (Sun setting below -6°, -12°, -18°) const civilTwilightEndUTC = findAltitudeCrossing( @@ -216,7 +216,7 @@ export function getSunMoonEvents(date: Date, observer: Observer, kernel: SpkKern etEnd, -6, false, - ) + ); const nauticalTwilightEndUTC = findAltitudeCrossing( kernel, NAIF_IDS.SUN, @@ -226,7 +226,7 @@ export function getSunMoonEvents(date: Date, observer: Observer, kernel: SpkKern etEnd, -12, false, - ) + ); const astronomicalTwilightEndUTC = findAltitudeCrossing( kernel, NAIF_IDS.SUN, @@ -236,7 +236,7 @@ export function getSunMoonEvents(date: Date, observer: Observer, kernel: SpkKern etEnd, -18, false, - ) + ); // Moon events const moonriseUTC = findAltitudeCrossing( @@ -248,7 +248,7 @@ export function getSunMoonEvents(date: Date, observer: Observer, kernel: SpkKern etEnd, MOON_ALTITUDE_THRESHOLD, true, - ) + ); const moonsetUTC = findAltitudeCrossing( kernel, NAIF_IDS.MOON, @@ -258,7 +258,7 @@ export function getSunMoonEvents(date: Date, observer: Observer, kernel: SpkKern etEnd, MOON_ALTITUDE_THRESHOLD, false, - ) + ); return { sunriseUTC, @@ -268,7 +268,7 @@ export function getSunMoonEvents(date: Date, observer: Observer, kernel: SpkKern civilTwilightEndUTC, nauticalTwilightEndUTC, astronomicalTwilightEndUTC, - } + }; } // ─── Best-time computation ──────────────────────────────────────────────────── @@ -290,16 +290,16 @@ export function bestTimeHeuristic( sunsetUTC: Date, moonsetUTC: Date, ): { bestTimeUTC: Date; lagMinutes: number } | null { - const lagMs = moonsetUTC.getTime() - sunsetUTC.getTime() - if (lagMs <= 0) return null // Moon sets before Sun — no sighting possible + const lagMs = moonsetUTC.getTime() - sunsetUTC.getTime(); + if (lagMs <= 0) return null; // Moon sets before Sun — no sighting possible - const lagMinutes = lagMs / 60000 - const bestTimeMs = sunsetUTC.getTime() + (4 / 9) * lagMs + const lagMinutes = lagMs / 60000; + const bestTimeMs = sunsetUTC.getTime() + (4 / 9) * lagMs; return { bestTimeUTC: new Date(bestTimeMs), lagMinutes, - } + }; } /** @@ -323,59 +323,59 @@ export function bestTimeOptimized( observer: Observer, steps = 90, ): { bestTimeUTC: Date; lagMinutes: number; maxV: number } | null { - const lagMs = moonsetUTC.getTime() - sunsetUTC.getTime() - if (lagMs <= 0) return null + const lagMs = moonsetUTC.getTime() - sunsetUTC.getTime(); + if (lagMs <= 0) return null; - const lagMinutes = lagMs / 60000 + const lagMinutes = lagMs / 60000; // Observer ITRS position (km) — fixed on Earth, computed once outside the loop - const obsECEF = geodeticToECEF(observer.lat, observer.lon, observer.elevation) - const obsITRS: Vec3 = [obsECEF[0] / 1000, obsECEF[1] / 1000, obsECEF[2] / 1000] + const obsECEF = geodeticToECEF(observer.lat, observer.lon, observer.elevation); + const obsITRS: Vec3 = [obsECEF[0] / 1000, obsECEF[1] / 1000, obsECEF[2] / 1000]; - let bestTimeUTC = sunsetUTC - let maxV = -Infinity + let bestTimeUTC = sunsetUTC; + let maxV = -Infinity; for (let i = 0; i <= steps; i++) { - const t = new Date(sunsetUTC.getTime() + (lagMs * i) / steps) - const ts = computeTimeScales(t) - const et = jdTTtoET(ts.jdTT) + const t = new Date(sunsetUTC.getTime() + (lagMs * i) / steps); + const ts = computeTimeScales(t); + const et = jdTTtoET(ts.jdTT); - const moonGCRS = getMoonGeocentricState(kernel, et).position - const sunGCRS = getSunGeocentricState(kernel, et).position + const moonGCRS = getMoonGeocentricState(kernel, et).position; + const sunGCRS = getSunGeocentricState(kernel, et).position; // Convert observer ITRS → GCRS at this timestep (Earth rotation changes per step) - const obsGCRS = itrsToGcrs(obsITRS, ts) + const obsGCRS = itrsToGcrs(obsITRS, ts); // Airless altitudes via the full pipeline - const moonAzAlt = computeAzAlt(moonGCRS, observer, ts, true) - const sunAzAlt = computeAzAlt(sunGCRS, observer, ts, true) + const moonAzAlt = computeAzAlt(moonGCRS, observer, ts, true); + const sunAzAlt = computeAzAlt(sunGCRS, observer, ts, true); - const ARCV = moonAzAlt.altitude - sunAzAlt.altitude + const ARCV = moonAzAlt.altitude - sunAzAlt.altitude; // Topocentric ARCL (Sun-Moon angular separation — all vectors in GCRS) const moonTopo: Vec3 = [ moonGCRS[0] - obsGCRS[0], moonGCRS[1] - obsGCRS[1], moonGCRS[2] - obsGCRS[2], - ] + ]; const sunTopo: Vec3 = [ sunGCRS[0] - obsGCRS[0], sunGCRS[1] - obsGCRS[1], sunGCRS[2] - obsGCRS[2], - ] - const cosARCL = vdot(moonTopo, sunTopo) / (vnorm(moonTopo) * vnorm(sunTopo)) - const ARCL = Math.acos(Math.max(-1, Math.min(1, cosARCL))) * (180 / Math.PI) + ]; + const cosARCL = vdot(moonTopo, sunTopo) / (vnorm(moonTopo) * vnorm(sunTopo)); + const ARCL = Math.acos(Math.max(-1, Math.min(1, cosARCL))) * (180 / Math.PI); - const { W } = computeCrescentWidth(moonTopo, ARCL) - const V = ARCV - arcvMinimum(W) + const { W } = computeCrescentWidth(moonTopo, ARCL); + const V = ARCV - arcvMinimum(W); if (V > maxV) { - maxV = V - bestTimeUTC = t + maxV = V; + bestTimeUTC = t; } } - return { bestTimeUTC, lagMinutes, maxV } + return { bestTimeUTC, lagMinutes, maxV }; } /** @@ -387,6 +387,6 @@ export function bestTimeOptimized( * @returns [start, end] UTC Date pair */ export function computeObservationWindow(bestTimeUTC: Date, windowMinutes = 20): [Date, Date] { - const windowMs = windowMinutes * 60000 - return [new Date(bestTimeUTC.getTime() - windowMs), new Date(bestTimeUTC.getTime() + windowMs)] + const windowMs = windowMinutes * 60000; + return [new Date(bestTimeUTC.getTime() - windowMs), new Date(bestTimeUTC.getTime() + windowMs)]; } diff --git a/src/frames/index.ts b/src/frames/index.ts index 1bc09e8..7daf326 100644 --- a/src/frames/index.ts +++ b/src/frames/index.ts @@ -21,18 +21,18 @@ * Capitaine et al. (2003), Astronomy & Astrophysics 412, 567-586 */ -import type { Vec3, TimeScales } from '../types.js' -import type { Mat3 } from '../math/index.js' -import { mvmul, mmmul, mtranspose, rotX, rotY, rotZ } from '../math/index.js' -import { J2000, DAYS_PER_JULIAN_CENTURY } from '../time/index.js' +import type { Vec3, TimeScales } from "../types.js"; +import type { Mat3 } from "../math/index.js"; +import { mvmul, mmmul, mtranspose, rotX, rotY, rotZ } from "../math/index.js"; +import { J2000, DAYS_PER_JULIAN_CENTURY } from "../time/index.js"; // ─── Constants ──────────────────────────────────────────────────────────────── /** Arcseconds to radians */ -const ARCSEC_RAD = Math.PI / (180 * 3600) +const ARCSEC_RAD = Math.PI / (180 * 3600); /** 0.1 microarcseconds to arcseconds (units of nutation coefficients) */ -const UAS01_TO_ARCSEC = 1e-7 +const UAS01_TO_ARCSEC = 1e-7; // ─── IAU 2000B Nutation Series ──────────────────────────────────────────────── // @@ -207,7 +207,7 @@ const NUT_2000B: ReadonlyArray< [0, 0, 4, -2, 4, 1341.0, 0.0, 0.0, -577.0, 0.0, 0.0], // 77 [0, 0, 2, -2, 4, -1316.0, 0.0, 0.0, 567.0, 0.0, 0.0], -] +]; // ─── Fundamental arguments (Delaunay) ──────────────────────────────────────── // Source: SOFA iauFal03, iauFalp03, iauFaf03, iauFad03, iauFaom03 @@ -215,43 +215,43 @@ const NUT_2000B: ReadonlyArray< /** Reduce arcseconds to [0, 2π) radians */ function arcsecToRad(arcsec: number): number { - const r = (arcsec * ARCSEC_RAD) % (2 * Math.PI) - return r >= 0 ? r : r + 2 * Math.PI + const r = (arcsec * ARCSEC_RAD) % (2 * Math.PI); + return r >= 0 ? r : r + 2 * Math.PI; } /** Mean anomaly of the Moon l (IAU 2003) */ function fundamentalL(T: number): number { return arcsecToRad( 485868.249036 + T * (1717915923.2178 + T * (31.8792 + T * (0.051635 + T * -0.0002447))), - ) + ); } /** Mean anomaly of the Sun l' (IAU 2003) */ function fundamentalLp(T: number): number { return arcsecToRad( 1287104.793048 + T * (129596581.0481 + T * (-0.5532 + T * (0.000136 + T * -0.00001149))), - ) + ); } /** Moon's argument of latitude F = L - Ω (IAU 2003) */ function fundamentalF(T: number): number { return arcsecToRad( 335779.526232 + T * (1739527262.8478 + T * (-12.7512 + T * (-0.001037 + T * 0.00000417))), - ) + ); } /** Mean elongation of the Moon D (IAU 2003) */ function fundamentalD(T: number): number { return arcsecToRad( 1072260.703692 + T * (1602961601.209 + T * (-6.3706 + T * (0.006593 + T * -0.00003169))), - ) + ); } /** Longitude of Moon's ascending node Ω (IAU 2003) */ function fundamentalOm(T: number): number { return arcsecToRad( 450160.398036 + T * (-6962890.5431 + T * (7.4722 + T * (0.007702 + T * -0.00005939))), - ) + ); } // ─── CIP coordinates ───────────────────────────────────────────────────────── @@ -273,29 +273,29 @@ function fundamentalOm(T: number): number { * @returns { X, Y, s } in radians */ export function computeCIPXYs(jdTT: number): { X: number; Y: number; s: number } { - const T = (jdTT - J2000) / DAYS_PER_JULIAN_CENTURY + const T = (jdTT - J2000) / DAYS_PER_JULIAN_CENTURY; // Delaunay fundamental arguments - const l = fundamentalL(T) - const lp = fundamentalLp(T) - const F = fundamentalF(T) - const D = fundamentalD(T) - const Om = fundamentalOm(T) + const l = fundamentalL(T); + const lp = fundamentalLp(T); + const F = fundamentalF(T); + const D = fundamentalD(T); + const Om = fundamentalOm(T); // Accumulate nutation in longitude (dpsi) and obliquity (deps) — units: 0.1 uas - let dpsi = 0.0 - let deps = 0.0 + let dpsi = 0.0; + let deps = 0.0; for (const [nl, nlp, nF, nD, nOm, ps, pst, pc, ec, ect, es] of NUT_2000B) { - const arg = nl * l + nlp * lp + nF * F + nD * D + nOm * Om - const sinA = Math.sin(arg) - const cosA = Math.cos(arg) - dpsi += (ps + pst * T) * sinA + pc * cosA - deps += (ec + ect * T) * cosA + es * sinA + const arg = nl * l + nlp * lp + nF * F + nD * D + nOm * Om; + const sinA = Math.sin(arg); + const cosA = Math.cos(arg); + dpsi += (ps + pst * T) * sinA + pc * cosA; + deps += (ec + ect * T) * cosA + es * sinA; } // Convert 0.1 uas → arcseconds → radians - const dpsiRad = dpsi * UAS01_TO_ARCSEC * ARCSEC_RAD - const depsRad = deps * UAS01_TO_ARCSEC * ARCSEC_RAD + const dpsiRad = dpsi * UAS01_TO_ARCSEC * ARCSEC_RAD; + const depsRad = deps * UAS01_TO_ARCSEC * ARCSEC_RAD; // Mean obliquity eps0 (IAU 2006, arcseconds → radians) // Reference: IERS Conventions (2010) Table 5.1 @@ -304,30 +304,30 @@ export function computeCIPXYs(jdTT: number): { X: number; Y: number; s: number } T * (-46.836769 + T * (-0.0001831 + T * (0.0020034 + T * (-0.000000576 + T * -0.0000000434))))) * - ARCSEC_RAD + ARCSEC_RAD; // IAU 2006 precession polynomial for X (arcseconds) // Reference: IERS Conventions (2010) Table 5.2a, polynomial s_X const Xarcsec = -0.016617 + - T * (2004.191898 + T * (-0.4297829 + T * (-0.19861834 + T * (0.000007578 + T * 0.0000059285)))) + T * (2004.191898 + T * (-0.4297829 + T * (-0.19861834 + T * (0.000007578 + T * 0.0000059285)))); // IAU 2006 precession polynomial for Y (arcseconds) // Reference: IERS Conventions (2010) Table 5.2a, polynomial s_Y const Yarcsec = -0.006951 + - T * (-0.025896 + T * (-22.4072747 + T * (0.00190059 + T * (0.001112526 + T * 0.0000001358)))) + T * (-0.025896 + T * (-22.4072747 + T * (0.00190059 + T * (0.001112526 + T * 0.0000001358)))); // CIP X, Y: precession polynomial + first-order nutation correction - const X = Xarcsec * ARCSEC_RAD + dpsiRad * Math.sin(eps0) - const Y = Yarcsec * ARCSEC_RAD - depsRad + const X = Xarcsec * ARCSEC_RAD + dpsiRad * Math.sin(eps0); + const Y = Yarcsec * ARCSEC_RAD - depsRad; // CIO locator s ≈ -X·Y/2 + small polynomial (IERS Conventions 2010 Eq. 5.9) // Polynomial term: s_poly ≈ -0.041775"·T (arcseconds) - const sPoly = -0.041775 * T * ARCSEC_RAD - const s = (-X * Y) / 2 + sPoly + const sPoly = -0.041775 * T * ARCSEC_RAD; + const s = (-X * Y) / 2 + sPoly; - return { X, Y, s } + return { X, Y, s }; } // ─── Earth Rotation Angle ──────────────────────────────────────────────────── @@ -345,9 +345,9 @@ export function computeCIPXYs(jdTT: number): { X: number; Y: number; s: number } * Reference: IAU 2000 Resolution B1.8; IERS Conventions (2010) §5.4.4 */ export function computeERA(jdUT1: number): number { - const Du = jdUT1 - 2451545.0 - const era = 2 * Math.PI * (0.779057273264 + 1.0027378119113546 * Du) - return ((era % (2 * Math.PI)) + 2 * Math.PI) % (2 * Math.PI) + const Du = jdUT1 - 2451545.0; + const era = 2 * Math.PI * (0.779057273264 + 1.0027378119113546 * Du); + return ((era % (2 * Math.PI)) + 2 * Math.PI) % (2 * Math.PI); } // ─── Frame rotation matrices ────────────────────────────────────────────────── @@ -363,10 +363,10 @@ export function computeERA(jdUT1: number): number { * Reference: SOFA iauC2ixys; IERS Conventions (2010) Eq. 5.7 */ export function celestialMotionMatrix(X: number, Y: number, s: number): Mat3 { - const r2 = X * X + Y * Y - const e = r2 > 0 ? Math.atan2(Y, X) : 0 - const d = Math.asin(Math.sqrt(r2)) - return mmmul(rotZ(-(e + s)), mmmul(rotY(d), rotZ(e))) + const r2 = X * X + Y * Y; + const e = r2 > 0 ? Math.atan2(Y, X) : 0; + const d = Math.asin(Math.sqrt(r2)); + return mmmul(rotZ(-(e + s)), mmmul(rotY(d), rotZ(e))); } /** @@ -374,7 +374,7 @@ export function celestialMotionMatrix(X: number, Y: number, s: number): Mat3 { * Simple rotation about the CIP pole (z-axis) by ERA. */ export function earthRotationMatrix(era: number): Mat3 { - return rotZ(era) + return rotZ(era); } /** @@ -391,7 +391,7 @@ export function earthRotationMatrix(era: number): Mat3 { * @param yp - Pole y-offset in radians (from IERS Bulletin A) */ export function polarMotionMatrix(xp: number, yp: number): Mat3 { - return mmmul(rotY(xp), rotX(-yp)) + return mmmul(rotY(xp), rotX(-yp)); } // ─── Full transformation ────────────────────────────────────────────────────── @@ -408,14 +408,14 @@ export function polarMotionMatrix(xp: number, yp: number): Mat3 { * @returns Vector in ITRS frame (km) */ export function gcrsToItrs(gcrsVec: Vec3, ts: TimeScales, xp = 0, yp = 0): Vec3 { - const { X, Y, s } = computeCIPXYs(ts.jdTT) - const Q = celestialMotionMatrix(X, Y, s) - const era = computeERA(ts.jdUT1) - const R = earthRotationMatrix(era) - const W = polarMotionMatrix(xp, yp) + const { X, Y, s } = computeCIPXYs(ts.jdTT); + const Q = celestialMotionMatrix(X, Y, s); + const era = computeERA(ts.jdUT1); + const R = earthRotationMatrix(era); + const W = polarMotionMatrix(xp, yp); // Apply Q first (GCRS→CIRS), then R (CIRS→TIRS), then W (TIRS→ITRS) - const combined = mmmul(W, mmmul(R, Q)) - return mvmul(combined, gcrsVec) + const combined = mmmul(W, mmmul(R, Q)); + return mvmul(combined, gcrsVec); } /** @@ -430,12 +430,12 @@ export function gcrsToItrs(gcrsVec: Vec3, ts: TimeScales, xp = 0, yp = 0): Vec3 * @returns Vector in GCRS frame (km) */ export function itrsToGcrs(itrsVec: Vec3, ts: TimeScales, xp = 0, yp = 0): Vec3 { - const { X, Y, s } = computeCIPXYs(ts.jdTT) - const Q = celestialMotionMatrix(X, Y, s) - const era = computeERA(ts.jdUT1) - const R = earthRotationMatrix(era) - const W = polarMotionMatrix(xp, yp) + const { X, Y, s } = computeCIPXYs(ts.jdTT); + const Q = celestialMotionMatrix(X, Y, s); + const era = computeERA(ts.jdUT1); + const R = earthRotationMatrix(era); + const W = polarMotionMatrix(xp, yp); // [GCRS] = Qᵀ · Rᵀ · Wᵀ · [ITRS] - const combined = mmmul(mtranspose(Q), mmmul(mtranspose(R), mtranspose(W))) - return mvmul(combined, itrsVec) + const combined = mmmul(mtranspose(Q), mmmul(mtranspose(R), mtranspose(W))); + return mvmul(combined, itrsVec); } diff --git a/src/index.ts b/src/index.ts index 7d09b93..a6366b1 100644 --- a/src/index.ts +++ b/src/index.ts @@ -27,7 +27,7 @@ export { initKernels, downloadKernels, verifyKernels, -} from './api/index.js' +} from "./api/index.js"; // ─── Types ──────────────────────────────────────────────────────────────────── @@ -61,7 +61,7 @@ export type { // Ephemeris internals (for advanced use) SpkSegment, ChebRecord, -} from './types.js' +} from "./types.js"; // ─── Constants ──────────────────────────────────────────────────────────────── @@ -71,4 +71,4 @@ export { ODEH_THRESHOLDS, ODEH_DESCRIPTIONS, WGS84, -} from './types.js' +} from "./types.js"; diff --git a/src/math/index.ts b/src/math/index.ts index a5a2264..3db601f 100644 --- a/src/math/index.ts +++ b/src/math/index.ts @@ -5,57 +5,57 @@ * Uses Float64Array for coefficient storage to match JS engine optimization paths. */ -import type { Vec3 } from '../types.js' +import type { Vec3 } from "../types.js"; // ─── Vector operations ──────────────────────────────────────────────────────── /** Add two 3-vectors */ export function vadd(a: Vec3, b: Vec3): Vec3 { - return [a[0] + b[0], a[1] + b[1], a[2] + b[2]] + return [a[0] + b[0], a[1] + b[1], a[2] + b[2]]; } /** Subtract b from a */ export function vsub(a: Vec3, b: Vec3): Vec3 { - return [a[0] - b[0], a[1] - b[1], a[2] - b[2]] + return [a[0] - b[0], a[1] - b[1], a[2] - b[2]]; } /** Scale a 3-vector */ export function vscale(a: Vec3, s: number): Vec3 { - return [a[0] * s, a[1] * s, a[2] * s] + return [a[0] * s, a[1] * s, a[2] * s]; } /** Dot product */ export function vdot(a: Vec3, b: Vec3): number { - return a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; } /** Euclidean norm */ export function vnorm(a: Vec3): number { - return Math.sqrt(vdot(a, a)) + return Math.sqrt(vdot(a, a)); } /** Cross product */ export function vcross(a: Vec3, b: Vec3): Vec3 { - return [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]] + return [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]]; } /** Unit vector (normalized) */ export function vunit(a: Vec3): Vec3 { - const n = vnorm(a) - if (n === 0) throw new RangeError('Cannot normalize a zero vector') - return vscale(a, 1 / n) + const n = vnorm(a); + if (n === 0) throw new RangeError("Cannot normalize a zero vector"); + return vscale(a, 1 / n); } /** Angular separation between two direction vectors in radians */ export function angularSep(a: Vec3, b: Vec3): number { - const cosAngle = Math.max(-1, Math.min(1, vdot(vunit(a), vunit(b)))) - return Math.acos(cosAngle) + const cosAngle = Math.max(-1, Math.min(1, vdot(vunit(a), vunit(b)))); + return Math.acos(cosAngle); } // ─── 3×3 matrix operations ──────────────────────────────────────────────────── /** 3×3 matrix stored row-major as a 9-element tuple */ -export type Mat3 = [number, number, number, number, number, number, number, number, number] +export type Mat3 = [number, number, number, number, number, number, number, number, number]; /** Multiply 3×3 matrix by 3-vector */ export function mvmul(m: Mat3, v: Vec3): Vec3 { @@ -63,7 +63,7 @@ export function mvmul(m: Mat3, v: Vec3): Vec3 { m[0] * v[0] + m[1] * v[1] + m[2] * v[2], m[3] * v[0] + m[4] * v[1] + m[5] * v[2], m[6] * v[0] + m[7] * v[1] + m[8] * v[2], - ] + ]; } /** Multiply two 3×3 matrices */ @@ -78,12 +78,12 @@ export function mmmul(a: Mat3, b: Mat3): Mat3 { a[6] * b[0] + a[7] * b[3] + a[8] * b[6], a[6] * b[1] + a[7] * b[4] + a[8] * b[7], a[6] * b[2] + a[7] * b[5] + a[8] * b[8], - ] + ]; } /** Transpose a 3×3 matrix */ export function mtranspose(m: Mat3): Mat3 { - return [m[0], m[3], m[6], m[1], m[4], m[7], m[2], m[5], m[8]] + return [m[0], m[3], m[6], m[1], m[4], m[7], m[2], m[5], m[8]]; } /** @@ -91,27 +91,27 @@ export function mtranspose(m: Mat3): Mat3 { * Follows right-hand rule. */ export function rotX(theta: number): Mat3 { - const c = Math.cos(theta) - const s = Math.sin(theta) - return [1, 0, 0, 0, c, s, 0, -s, c] + const c = Math.cos(theta); + const s = Math.sin(theta); + return [1, 0, 0, 0, c, s, 0, -s, c]; } /** * Rotation matrix around the Y axis by angle θ (radians). */ export function rotY(theta: number): Mat3 { - const c = Math.cos(theta) - const s = Math.sin(theta) - return [c, 0, -s, 0, 1, 0, s, 0, c] + const c = Math.cos(theta); + const s = Math.sin(theta); + return [c, 0, -s, 0, 1, 0, s, 0, c]; } /** * Rotation matrix around the Z axis by angle θ (radians). */ export function rotZ(theta: number): Mat3 { - const c = Math.cos(theta) - const s = Math.sin(theta) - return [c, s, 0, -s, c, 0, 0, 0, 1] + const c = Math.cos(theta); + const s = Math.sin(theta); + return [c, s, 0, -s, c, 0, 0, 0, 1]; } // ─── Chebyshev polynomial evaluation ───────────────────────────────────────── @@ -128,22 +128,22 @@ export function rotZ(theta: number): Mat3 { * @returns Polynomial value at x */ export function chebyshevEval(coeffs: Float64Array, x: number): number { - const n = coeffs.length - if (n === 0) return 0 - if (n === 1) return coeffs[0]! // Float64Array length checked above; index 0 is valid + const n = coeffs.length; + if (n === 0) return 0; + if (n === 1) return coeffs[0]!; // Float64Array length checked above; index 0 is valid // Double-x for Clenshaw efficiency - const x2 = 2 * x - let b2 = 0 - let b1 = 0 + const x2 = 2 * x; + let b2 = 0; + let b1 = 0; for (let k = n - 1; k >= 1; k--) { - const b0 = coeffs[k]! + x2 * b1 - b2 // k in [1, n-1]; all valid indices - b2 = b1 - b1 = b0 + const b0 = coeffs[k]! + x2 * b1 - b2; // k in [1, n-1]; all valid indices + b2 = b1; + b1 = b0; } - return coeffs[0]! + x * b1 - b2 // index 0 valid; n >= 2 at this point + return coeffs[0]! + x * b1 - b2; // index 0 valid; n >= 2 at this point } /** @@ -160,29 +160,29 @@ export function chebyshevEvalWithDerivative( x: number, radius: number, ): [number, number] { - const n = coeffs.length - if (n === 0) return [0, 0] - if (n === 1) return [coeffs[0]!, 0] // length checked; index 0 valid + const n = coeffs.length; + if (n === 0) return [0, 0]; + if (n === 1) return [coeffs[0]!, 0]; // length checked; index 0 valid - const x2 = 2 * x - let b2 = 0 - let b1 = 0 - let db2 = 0 - let db1 = 0 + const x2 = 2 * x; + let b2 = 0; + let b1 = 0; + let db2 = 0; + let db1 = 0; for (let k = n - 1; k >= 1; k--) { - const b0 = coeffs[k]! + x2 * b1 - b2 // k in [1, n-1]; all valid indices - const db0 = 2 * b1 + x2 * db1 - db2 - b2 = b1 - b1 = b0 - db2 = db1 - db1 = db0 + const b0 = coeffs[k]! + x2 * b1 - b2; // k in [1, n-1]; all valid indices + const db0 = 2 * b1 + x2 * db1 - db2; + b2 = b1; + b1 = b0; + db2 = db1; + db1 = db0; } - const value = coeffs[0]! + x * b1 - b2 // index 0 valid; n >= 2 at this point - const dvalue = b1 + x * db1 - db2 + const value = coeffs[0]! + x * b1 - b2; // index 0 valid; n >= 2 at this point + const dvalue = b1 + x * db1 - db2; // Scale derivative from normalized domain back to seconds - return [value, dvalue / radius] + return [value, dvalue / radius]; } // ─── Root finding ───────────────────────────────────────────────────────────── @@ -208,71 +208,71 @@ export function brentRoot( tol = 1e-9, maxIter = 64, ): number | null { - let fa = f(a) - let fb = f(b) + let fa = f(a); + let fb = f(b); // No sign change in bracket - if (fa * fb > 0) return null + if (fa * fb > 0) return null; // Swap so |f(b)| <= |f(a)| if (Math.abs(fa) < Math.abs(fb)) { - ;[a, b] = [b, a] - ;[fa, fb] = [fb, fa] + [a, b] = [b, a]; + [fa, fb] = [fb, fa]; } - let c = a - let fc = fa - let mflag = true - let s: number - let d = 0 + let c = a; + let fc = fa; + let mflag = true; + let s: number; + let d = 0; for (let i = 0; i < maxIter; i++) { - if (Math.abs(b - a) < tol) return b + if (Math.abs(b - a) < tol) return b; if (fa !== fc && fb !== fc) { // Inverse quadratic interpolation s = (a * fb * fc) / ((fa - fb) * (fa - fc)) + (b * fa * fc) / ((fb - fa) * (fb - fc)) + - (c * fa * fb) / ((fc - fa) * (fc - fb)) + (c * fa * fb) / ((fc - fa) * (fc - fb)); } else { // Secant method - s = b - fb * ((b - a) / (fb - fa)) + s = b - fb * ((b - a) / (fb - fa)); } - const cond1 = s < (3 * a + b) / 4 || s > b - const cond2 = mflag && Math.abs(s - b) >= Math.abs(b - c) / 2 - const cond3 = !mflag && Math.abs(s - b) >= Math.abs(c - d) / 2 - const cond4 = mflag && Math.abs(b - c) < tol - const cond5 = !mflag && Math.abs(c - d) < tol + const cond1 = s < (3 * a + b) / 4 || s > b; + const cond2 = mflag && Math.abs(s - b) >= Math.abs(b - c) / 2; + const cond3 = !mflag && Math.abs(s - b) >= Math.abs(c - d) / 2; + const cond4 = mflag && Math.abs(b - c) < tol; + const cond5 = !mflag && Math.abs(c - d) < tol; if (cond1 || cond2 || cond3 || cond4 || cond5) { - s = (a + b) / 2 - mflag = true + s = (a + b) / 2; + mflag = true; } else { - mflag = false + mflag = false; } - const fs = f(s) - d = c - c = b - fc = fb + const fs = f(s); + d = c; + c = b; + fc = fb; if (fa * fs < 0) { - b = s - fb = fs + b = s; + fb = fs; } else { - a = s - fa = fs + a = s; + fa = fs; } if (Math.abs(fa) < Math.abs(fb)) { - ;[a, b] = [b, a] - ;[fa, fb] = [fb, fa] + [a, b] = [b, a]; + [fa, fb] = [fb, fa]; } } - return b + return b; } /** @@ -286,54 +286,55 @@ export function brentRoot( * @returns Array of root locations */ export function findRoots(f: (t: number) => number, a: number, b: number, steps = 48): number[] { - const dt = (b - a) / steps - const roots: number[] = [] - let tPrev = a - let fPrev = f(a) + const dt = (b - a) / steps; + const roots: number[] = []; + let tPrev = a; + let fPrev = f(a); for (let i = 1; i <= steps; i++) { - const t = a + i * dt - const ft = f(t) + const t = a + i * dt; + const ft = f(t); if (fPrev * ft <= 0) { // Sign change in [tPrev, t] — apply Brent's method - const root = brentRoot(f, tPrev, t, 1e-9, 64) + const root = brentRoot(f, tPrev, t, 1e-9, 64); if (root !== null) { // Deduplicate roots that are too close together - if (roots.length === 0 || Math.abs(root - roots[roots.length - 1]!) > 1e-6) { // length > 0 checked via || - roots.push(root) + if (roots.length === 0 || Math.abs(root - roots[roots.length - 1]!) > 1e-6) { + // length > 0 checked via || + roots.push(root); } } } - tPrev = t - fPrev = ft + tPrev = t; + fPrev = ft; } - return roots + return roots; } // ─── Angle utilities ───────────────────────────────────────────────────────── /** Convert degrees to radians */ -export const DEG2RAD = Math.PI / 180 +export const DEG2RAD = Math.PI / 180; /** Convert radians to degrees */ -export const RAD2DEG = 180 / Math.PI +export const RAD2DEG = 180 / Math.PI; /** Normalize an angle to [0, 2π) */ export function mod2pi(angle: number): number { - const twoPi = 2 * Math.PI - return ((angle % twoPi) + twoPi) % twoPi + const twoPi = 2 * Math.PI; + return ((angle % twoPi) + twoPi) % twoPi; } /** Normalize an angle in degrees to [0, 360) */ export function mod360(deg: number): number { - return ((deg % 360) + 360) % 360 + return ((deg % 360) + 360) % 360; } /** Normalize an angle in degrees to [-180, 180) */ export function normalizeDeg180(deg: number): number { - deg = mod360(deg) - return deg >= 180 ? deg - 360 : deg + deg = mod360(deg); + return deg >= 180 ? deg - 360 : deg; } diff --git a/src/observer/index.ts b/src/observer/index.ts index 4dbdf79..b22cf0f 100644 --- a/src/observer/index.ts +++ b/src/observer/index.ts @@ -18,10 +18,10 @@ * Saemundsson (1986), Sky & Telescope 72, 70 */ -import type { Vec3, Observer, AzAlt, TimeScales } from '../types.js' -import { WGS84 } from '../types.js' -import { gcrsToItrs } from '../frames/index.js' -import { vdot } from '../math/index.js' +import type { Vec3, Observer, AzAlt, TimeScales } from "../types.js"; +import { WGS84 } from "../types.js"; +import { gcrsToItrs } from "../frames/index.js"; +import { vdot } from "../math/index.js"; // ─── Geodetic ↔ ECEF ───────────────────────────────────────────────────────── @@ -35,17 +35,17 @@ import { vdot } from '../math/index.js' * @returns ECEF position vector in meters */ export function geodeticToECEF(lat: number, lon: number, elev: number): Vec3 { - const phi = (lat * Math.PI) / 180 - const lam = (lon * Math.PI) / 180 - const sinPhi = Math.sin(phi) - const cosPhi = Math.cos(phi) + const phi = (lat * Math.PI) / 180; + const lam = (lon * Math.PI) / 180; + const sinPhi = Math.sin(phi); + const cosPhi = Math.cos(phi); // Prime vertical radius of curvature - const N = WGS84.a / Math.sqrt(1 - WGS84.e2 * sinPhi * sinPhi) + const N = WGS84.a / Math.sqrt(1 - WGS84.e2 * sinPhi * sinPhi); return [ (N + elev) * cosPhi * Math.cos(lam), (N + elev) * cosPhi * Math.sin(lam), (N * (1 - WGS84.e2) + elev) * sinPhi, - ] + ]; } /** @@ -55,27 +55,27 @@ export function geodeticToECEF(lat: number, lon: number, elev: number): Vec3 { * @returns { lat, lon, h } — latitude/longitude in degrees, height in meters */ export function ecefToGeodetic(ecef: Vec3): { lat: number; lon: number; h: number } { - const [X, Y, Z] = ecef - const p = Math.sqrt(X * X + Y * Y) - const lon = Math.atan2(Y, X) + const [X, Y, Z] = ecef; + const p = Math.sqrt(X * X + Y * Y); + const lon = Math.atan2(Y, X); // Bowring iteration for geodetic latitude - let lat = Math.atan2(Z, p * (1 - WGS84.e2)) + let lat = Math.atan2(Z, p * (1 - WGS84.e2)); for (let i = 0; i < 4; i++) { - const sinLat = Math.sin(lat) - const N = WGS84.a / Math.sqrt(1 - WGS84.e2 * sinLat * sinLat) - lat = Math.atan2(Z + WGS84.e2 * N * sinLat, p) + const sinLat = Math.sin(lat); + const N = WGS84.a / Math.sqrt(1 - WGS84.e2 * sinLat * sinLat); + lat = Math.atan2(Z + WGS84.e2 * N * sinLat, p); } - const sinLat = Math.sin(lat) - const N = WGS84.a / Math.sqrt(1 - WGS84.e2 * sinLat * sinLat) - const h = p / Math.cos(lat) - N + const sinLat = Math.sin(lat); + const N = WGS84.a / Math.sqrt(1 - WGS84.e2 * sinLat * sinLat); + const h = p / Math.cos(lat) - N; return { lat: (lat * 180) / Math.PI, lon: (lon * 180) / Math.PI, h, - } + }; } // ─── Topocentric ENU ───────────────────────────────────────────────────────── @@ -92,18 +92,18 @@ export function ecefToGeodetic(ecef: Vec3): { lat: number; lon: number; h: numbe * @param lon - Longitude in degrees */ export function computeENUBasis(lat: number, lon: number): { east: Vec3; north: Vec3; up: Vec3 } { - const phi = (lat * Math.PI) / 180 - const lam = (lon * Math.PI) / 180 + const phi = (lat * Math.PI) / 180; + const lam = (lon * Math.PI) / 180; const sinPhi = Math.sin(phi), - cosPhi = Math.cos(phi) + cosPhi = Math.cos(phi); const sinLam = Math.sin(lam), - cosLam = Math.cos(lam) + cosLam = Math.cos(lam); - const east: Vec3 = [-sinLam, cosLam, 0] - const north: Vec3 = [-sinPhi * cosLam, -sinPhi * sinLam, cosPhi] - const up: Vec3 = [cosPhi * cosLam, cosPhi * sinLam, sinPhi] + const east: Vec3 = [-sinLam, cosLam, 0]; + const north: Vec3 = [-sinPhi * cosLam, -sinPhi * sinLam, cosPhi]; + const up: Vec3 = [cosPhi * cosLam, cosPhi * sinLam, sinPhi]; - return { east, north, up } + return { east, north, up }; } /** @@ -116,8 +116,8 @@ export function computeENUBasis(lat: number, lon: number): { east: Vec3; north: * @returns ENU vector [east, north, up] in the same units as input */ export function ecefToENU(ecefDelta: Vec3, lat: number, lon: number): Vec3 { - const { east, north, up } = computeENUBasis(lat, lon) - return [vdot(ecefDelta, east), vdot(ecefDelta, north), vdot(ecefDelta, up)] + const { east, north, up } = computeENUBasis(lat, lon); + return [vdot(ecefDelta, east), vdot(ecefDelta, north), vdot(ecefDelta, up)]; } /** @@ -130,13 +130,13 @@ export function ecefToENU(ecefDelta: Vec3, lat: number, lon: number): Vec3 { * @returns Azimuth (degrees, [0, 360)) and altitude (degrees) */ export function enuToAzAlt(enu: Vec3): AzAlt { - const [e, n, u] = enu - const horiz = Math.sqrt(e * e + n * n) - const altitude = (Math.atan2(u, horiz) * 180) / Math.PI + const [e, n, u] = enu; + const horiz = Math.sqrt(e * e + n * n); + const altitude = (Math.atan2(u, horiz) * 180) / Math.PI; // atan2(east, north) gives bearing from North; convert to [0, 360) - let azimuth = (Math.atan2(e, n) * 180) / Math.PI - if (azimuth < 0) azimuth += 360 - return { azimuth, altitude } + let azimuth = (Math.atan2(e, n) * 180) / Math.PI; + if (azimuth < 0) azimuth += 360; + return { azimuth, altitude }; } // ─── Topocentric parallax ───────────────────────────────────────────────────── @@ -158,7 +158,7 @@ export function topocentricPosition(bodyGCRS: Vec3, observerGCRS: Vec3): Vec3 { bodyGCRS[0] - observerGCRS[0], bodyGCRS[1] - observerGCRS[1], bodyGCRS[2] - observerGCRS[2], - ] + ]; } // ─── Full pipeline: GCRS → az/alt ──────────────────────────────────────────── @@ -188,27 +188,31 @@ export function computeAzAlt( airless: boolean, ): AzAlt { // 1. Convert body position from GCRS to ITRS (km) - const bodyITRS = gcrsToItrs(bodyGCRS, ts) + const bodyITRS = gcrsToItrs(bodyGCRS, ts); // 2. Observer position in ITRS: geodeticToECEF returns meters → convert to km - const obsECEF = geodeticToECEF(observer.lat, observer.lon, observer.elevation) - const obsITRS: Vec3 = [obsECEF[0] / 1000, obsECEF[1] / 1000, obsECEF[2] / 1000] + const obsECEF = geodeticToECEF(observer.lat, observer.lon, observer.elevation); + const obsITRS: Vec3 = [obsECEF[0] / 1000, obsECEF[1] / 1000, obsECEF[2] / 1000]; // 3. Displacement vector from observer to body in ITRS (km — magnitude doesn't matter) - const delta: Vec3 = [bodyITRS[0] - obsITRS[0], bodyITRS[1] - obsITRS[1], bodyITRS[2] - obsITRS[2]] + const delta: Vec3 = [ + bodyITRS[0] - obsITRS[0], + bodyITRS[1] - obsITRS[1], + bodyITRS[2] - obsITRS[2], + ]; // 4. Project onto local ENU basis at the observer's location - const enu = ecefToENU(delta, observer.lat, observer.lon) + const enu = ecefToENU(delta, observer.lat, observer.lon); // 5. Convert ENU to azimuth + altitude - const azAlt = enuToAzAlt(enu) + const azAlt = enuToAzAlt(enu); // 6. Refraction correction if (!airless) { - azAlt.altitude = applyRefraction(azAlt.altitude, observer.pressure, observer.temperature) + azAlt.altitude = applyRefraction(azAlt.altitude, observer.pressure, observer.temperature); } - return azAlt + return azAlt; } // ─── Atmospheric refraction ─────────────────────────────────────────────────── @@ -237,17 +241,17 @@ export function bennettRefraction( temperature = 15, ): number { // No refraction below the geometric horizon (Bennett formula diverges below ~−1°) - if (altitudeDeg < -1) return 0 + if (altitudeDeg < -1) return 0; // Convert altitude argument to radians for the cot computation - const h = altitudeDeg - const argDeg = h + 7.31 / (h + 4.4) - const argRad = (argDeg * Math.PI) / 180 - const R = 1 / (Math.tan(argRad) * 60) // degrees + const h = altitudeDeg; + const argDeg = h + 7.31 / (h + 4.4); + const argRad = (argDeg * Math.PI) / 180; + const R = 1 / (Math.tan(argRad) * 60); // degrees // Pressure and temperature correction - const corrected = R * (pressure / 1010) * (283 / (273 + temperature)) - return Math.max(0, corrected) + const corrected = R * (pressure / 1010) * (283 / (273 + temperature)); + return Math.max(0, corrected); } /** @@ -255,7 +259,7 @@ export function bennettRefraction( * Returns the apparent (observed) altitude. */ export function applyRefraction(airlessAlt: number, pressure = 1013.25, temperature = 15): number { - return airlessAlt + bennettRefraction(airlessAlt, pressure, temperature) + return airlessAlt + bennettRefraction(airlessAlt, pressure, temperature); } /** @@ -268,9 +272,9 @@ export function removeRefraction( temperature = 15, ): number { // Start from the apparent altitude and iterate - let airless = apparentAlt + let airless = apparentAlt; for (let i = 0; i < 4; i++) { - airless = apparentAlt - bennettRefraction(airless, pressure, temperature) + airless = apparentAlt - bennettRefraction(airless, pressure, temperature); } - return airless + return airless; } diff --git a/src/spk/index.ts b/src/spk/index.ts index 2bcc67e..e584b8a 100644 --- a/src/spk/index.ts +++ b/src/spk/index.ts @@ -24,8 +24,8 @@ * SPICE source: SPKE02.f, DAFRA.f */ -import type { SpkSegment, StateVector } from '../types.js' -import { chebyshevEvalWithDerivative } from '../math/index.js' +import type { SpkSegment, StateVector } from "../types.js"; +import { chebyshevEvalWithDerivative } from "../math/index.js"; // ─── NAIF body IDs ──────────────────────────────────────────────────────────── @@ -44,14 +44,14 @@ export const NAIF_IDS = { SUN: 10, MOON: 301, EARTH: 399, -} as const +} as const; /** Frame code for ICRF/J2000 (the inertial reference frame used by DE442S) */ -export const FRAME_J2000 = 1 +export const FRAME_J2000 = 1; /** DAF record size in bytes */ -const DAF_RECORD_SIZE = 1024 -const BYTES_PER_DOUBLE = 8 +const DAF_RECORD_SIZE = 1024; +const BYTES_PER_DOUBLE = 8; // ─── SPK Kernel ─────────────────────────────────────────────────────────────── @@ -59,39 +59,39 @@ const BYTES_PER_DOUBLE = 8 * A loaded SPK kernel with segment index. */ export class SpkKernel { - private readonly buffer: ArrayBuffer - private readonly segments: SpkSegment[] - private readonly index: Map - private readonly le: boolean + private readonly buffer: ArrayBuffer; + private readonly segments: SpkSegment[]; + private readonly index: Map; + private readonly le: boolean; private constructor(buffer: ArrayBuffer, segments: SpkSegment[], le: boolean) { - this.buffer = buffer - this.segments = segments - this.le = le - this.index = new Map() + this.buffer = buffer; + this.segments = segments; + this.le = le; + this.index = new Map(); for (const seg of segments) { - const key = `${seg.target}:${seg.center}` - const list = this.index.get(key) ?? [] - list.push(seg) - this.index.set(key, list) + const key = `${seg.target}:${seg.center}`; + const list = this.index.get(key) ?? []; + list.push(seg); + this.index.set(key, list); } } /** Load a kernel from a binary ArrayBuffer. */ static fromBuffer(buffer: ArrayBuffer): SpkKernel { - const { nd, ni, fward, le } = parseDafFileRecord(buffer) - const segments = parseSummaryRecords(buffer, fward, nd, ni, le) - return new SpkKernel(buffer, segments, le) + const { nd, ni, fward, le } = parseDafFileRecord(buffer); + const segments = parseSummaryRecords(buffer, fward, nd, ni, le); + return new SpkKernel(buffer, segments, le); } /** Load a kernel from a file path (Node.js only). */ static fromFile(path: string): SpkKernel { // eslint-disable-next-line @typescript-eslint/no-require-imports - const fs = require('fs') as typeof import('fs') - const buf = fs.readFileSync(path) - const ab = buf.buffer.slice(buf.byteOffset, buf.byteOffset + buf.byteLength) - return SpkKernel.fromBuffer(ab as ArrayBuffer) + const fs = require("fs") as typeof import("fs"); + const buf = fs.readFileSync(path); + const ab = buf.buffer.slice(buf.byteOffset, buf.byteOffset + buf.byteLength); + return SpkKernel.fromBuffer(ab as ArrayBuffer); } /** @@ -99,104 +99,104 @@ export class SpkKernel { * Uses segment chaining when no direct segment exists. */ getState(target: number, center: number, et: number): StateVector { - const direct = this.findSeg(target, center, et) - if (direct) return evaluateSegment(this.buffer, direct, et, this.le) - return this.getChained(target, center, et) + const direct = this.findSeg(target, center, et); + if (direct) return evaluateSegment(this.buffer, direct, et, this.le); + return this.getChained(target, center, et); } private findSeg(target: number, center: number, et: number): SpkSegment | null { - const candidates = this.index.get(`${target}:${center}`) - if (!candidates) return null - return candidates.find((s) => et >= s.startET && et <= s.endET) ?? null + const candidates = this.index.get(`${target}:${center}`); + if (!candidates) return null; + return candidates.find((s) => et >= s.startET && et <= s.endET) ?? null; } private getChained(target: number, center: number, et: number): StateVector { - const ssb = NAIF_IDS.SSB - const emb = NAIF_IDS.EMB + const ssb = NAIF_IDS.SSB; + const emb = NAIF_IDS.EMB; // Moon relative to Earth: Moon-EMB minus Earth-EMB if (target === NAIF_IDS.MOON && center === NAIF_IDS.EARTH) { - const s1 = this.findSeg(NAIF_IDS.MOON, emb, et) - const s2 = this.findSeg(NAIF_IDS.EARTH, emb, et) + const s1 = this.findSeg(NAIF_IDS.MOON, emb, et); + const s2 = this.findSeg(NAIF_IDS.EARTH, emb, et); if (s1 && s2) { return subtractSV( evaluateSegment(this.buffer, s1, et, this.le), evaluateSegment(this.buffer, s2, et, this.le), - ) + ); } } // Earth relative to Moon (inverse) if (target === NAIF_IDS.EARTH && center === NAIF_IDS.MOON) { - const s1 = this.findSeg(NAIF_IDS.EARTH, emb, et) - const s2 = this.findSeg(NAIF_IDS.MOON, emb, et) + const s1 = this.findSeg(NAIF_IDS.EARTH, emb, et); + const s2 = this.findSeg(NAIF_IDS.MOON, emb, et); if (s1 && s2) { return subtractSV( evaluateSegment(this.buffer, s1, et, this.le), evaluateSegment(this.buffer, s2, et, this.le), - ) + ); } } // Sun relative to Earth if (target === NAIF_IDS.SUN && center === NAIF_IDS.EARTH) { - const sSunSsb = this.findSeg(NAIF_IDS.SUN, ssb, et) - const sEmbSsb = this.findSeg(emb, ssb, et) - const sEarthEmb = this.findSeg(NAIF_IDS.EARTH, emb, et) + const sSunSsb = this.findSeg(NAIF_IDS.SUN, ssb, et); + const sEmbSsb = this.findSeg(emb, ssb, et); + const sEarthEmb = this.findSeg(NAIF_IDS.EARTH, emb, et); if (sSunSsb && sEmbSsb && sEarthEmb) { - const svSunSsb = evaluateSegment(this.buffer, sSunSsb, et, this.le) - const svEmbSsb = evaluateSegment(this.buffer, sEmbSsb, et, this.le) - const svEarthEmb = evaluateSegment(this.buffer, sEarthEmb, et, this.le) + const svSunSsb = evaluateSegment(this.buffer, sSunSsb, et, this.le); + const svEmbSsb = evaluateSegment(this.buffer, sEmbSsb, et, this.le); + const svEarthEmb = evaluateSegment(this.buffer, sEarthEmb, et, this.le); // Earth/SSB = EMB/SSB - Earth/EMB - const earthSsb = subtractSV(svEmbSsb, svEarthEmb) - return subtractSV(svSunSsb, earthSsb) + const earthSsb = subtractSV(svEmbSsb, svEarthEmb); + return subtractSV(svSunSsb, earthSsb); } } // Generic two-hop via SSB - const sTargetSsb = this.findSeg(target, ssb, et) - const sCenterSsb = this.findSeg(center, ssb, et) + const sTargetSsb = this.findSeg(target, ssb, et); + const sCenterSsb = this.findSeg(center, ssb, et); if (sTargetSsb && sCenterSsb) { return subtractSV( evaluateSegment(this.buffer, sTargetSsb, et, this.le), evaluateSegment(this.buffer, sCenterSsb, et, this.le), - ) + ); } - throw new Error(`SpkKernel: no path for target=${target} center=${center} et=${et}`) + throw new Error(`SpkKernel: no path for target=${target} center=${center} et=${et}`); } getSegments(): ReadonlyArray { - return this.segments + return this.segments; } } // ─── DAF parsing ────────────────────────────────────────────────────────────── function parseDafFileRecord(buffer: ArrayBuffer): { - nd: number - ni: number - fward: number - bward: number - free: number - le: boolean + nd: number; + ni: number; + fward: number; + bward: number; + free: number; + le: boolean; } { - const dv = new DataView(buffer) + const dv = new DataView(buffer); // Detect endianness by reading ND (should be 2 for DE442S SPK) - let le = true - let nd = dv.getInt32(8, true) + let le = true; + let nd = dv.getInt32(8, true); if (nd < 1 || nd > 100) { - nd = dv.getInt32(8, false) - le = false + nd = dv.getInt32(8, false); + le = false; } - const ni = dv.getInt32(12, le) - const fward = dv.getInt32(256, le) - const bward = dv.getInt32(260, le) - const free = dv.getInt32(264, le) + const ni = dv.getInt32(12, le); + const fward = dv.getInt32(256, le); + const bward = dv.getInt32(260, le); + const free = dv.getInt32(264, le); - return { nd, ni, fward, bward, free, le } + return { nd, ni, fward, bward, free, le }; } function parseSummaryRecords( @@ -206,46 +206,46 @@ function parseSummaryRecords( ni: number, le: boolean, ): SpkSegment[] { - const dv = new DataView(buffer) - const segments: SpkSegment[] = [] - const summaryBytes = nd * BYTES_PER_DOUBLE + ni * 4 + const dv = new DataView(buffer); + const segments: SpkSegment[] = []; + const summaryBytes = nd * BYTES_PER_DOUBLE + ni * 4; - let recordNum = fward + let recordNum = fward; while (recordNum !== 0) { - const recOffset = (recordNum - 1) * DAF_RECORD_SIZE + const recOffset = (recordNum - 1) * DAF_RECORD_SIZE; // Control area: 3 doubles at start of record - const nextRecord = dv.getFloat64(recOffset, le) - const nSummaries = Math.round(dv.getFloat64(recOffset + 16, le)) + const nextRecord = dv.getFloat64(recOffset, le); + const nSummaries = Math.round(dv.getFloat64(recOffset + 16, le)); - let offset = recOffset + 24 // skip 3 control doubles (24 bytes) + let offset = recOffset + 24; // skip 3 control doubles (24 bytes) for (let i = 0; i < nSummaries; i++) { - if (offset + summaryBytes > buffer.byteLength) break + if (offset + summaryBytes > buffer.byteLength) break; - const startET = dv.getFloat64(offset, le) - const endET = dv.getFloat64(offset + 8, le) - offset += nd * BYTES_PER_DOUBLE + const startET = dv.getFloat64(offset, le); + const endET = dv.getFloat64(offset + 8, le); + offset += nd * BYTES_PER_DOUBLE; - const target = dv.getInt32(offset, le) - const center = dv.getInt32(offset + 4, le) - const frame = dv.getInt32(offset + 8, le) - const dataType = dv.getInt32(offset + 12, le) as 2 | 3 - const beginAddr = dv.getInt32(offset + 16, le) - const endAddr = dv.getInt32(offset + 20, le) - offset += ni * 4 + const target = dv.getInt32(offset, le); + const center = dv.getInt32(offset + 4, le); + const frame = dv.getInt32(offset + 8, le); + const dataType = dv.getInt32(offset + 12, le) as 2 | 3; + const beginAddr = dv.getInt32(offset + 16, le); + const endAddr = dv.getInt32(offset + 20, le); + offset += ni * 4; - const dataOffset = (beginAddr - 1) * BYTES_PER_DOUBLE - const dataSize = endAddr - beginAddr + 1 + const dataOffset = (beginAddr - 1) * BYTES_PER_DOUBLE; + const dataSize = endAddr - beginAddr + 1; - segments.push({ target, center, frame, dataType, startET, endET, dataOffset, dataSize }) + segments.push({ target, center, frame, dataType, startET, endET, dataOffset, dataSize }); } - recordNum = Math.round(nextRecord) + recordNum = Math.round(nextRecord); } - return segments + return segments; } // ─── Segment evaluation ─────────────────────────────────────────────────────── @@ -256,9 +256,9 @@ function evaluateSegment( et: number, le: boolean, ): StateVector { - if (seg.dataType === 2) return evaluateType2(buffer, seg, et, le) - if (seg.dataType === 3) return evaluateType3(buffer, seg, et, le) - throw new Error(`Unsupported SPK segment type: ${seg.dataType}`) + if (seg.dataType === 2) return evaluateType2(buffer, seg, et, le); + if (seg.dataType === 3) return evaluateType3(buffer, seg, et, le); + throw new Error(`Unsupported SPK segment type: ${seg.dataType}`); } /** @@ -270,36 +270,36 @@ export function evaluateType2( et: number, le = true, ): StateVector { - const dv = new DataView(buffer) - const endOffset = seg.dataOffset + seg.dataSize * BYTES_PER_DOUBLE + const dv = new DataView(buffer); + const endOffset = seg.dataOffset + seg.dataSize * BYTES_PER_DOUBLE; // Directory at end of data (4 doubles before the data end) - const N = dv.getFloat64(endOffset - BYTES_PER_DOUBLE, le) - const rsize = dv.getFloat64(endOffset - 2 * BYTES_PER_DOUBLE, le) - const intlen = dv.getFloat64(endOffset - 3 * BYTES_PER_DOUBLE, le) - const init = dv.getFloat64(endOffset - 4 * BYTES_PER_DOUBLE, le) + const N = dv.getFloat64(endOffset - BYTES_PER_DOUBLE, le); + const rsize = dv.getFloat64(endOffset - 2 * BYTES_PER_DOUBLE, le); + const intlen = dv.getFloat64(endOffset - 3 * BYTES_PER_DOUBLE, le); + const init = dv.getFloat64(endOffset - 4 * BYTES_PER_DOUBLE, le); // degree = (rsize - 2) / 3 (Type 2 stores 3 components) - const degree = Math.round((rsize - 2) / 3) - const nCoeffs = degree + 1 + const degree = Math.round((rsize - 2) / 3); + const nCoeffs = degree + 1; - let recIdx = Math.floor((et - init) / intlen) - recIdx = Math.max(0, Math.min(Math.round(N) - 1, recIdx)) + let recIdx = Math.floor((et - init) / intlen); + recIdx = Math.max(0, Math.min(Math.round(N) - 1, recIdx)); - const recOffset = seg.dataOffset + recIdx * rsize * BYTES_PER_DOUBLE - const mid = dv.getFloat64(recOffset, le) - const radius = dv.getFloat64(recOffset + BYTES_PER_DOUBLE, le) - const x = (et - mid) / radius + const recOffset = seg.dataOffset + recIdx * rsize * BYTES_PER_DOUBLE; + const mid = dv.getFloat64(recOffset, le); + const radius = dv.getFloat64(recOffset + BYTES_PER_DOUBLE, le); + const x = (et - mid) / radius; - const xC = readCoeffs(dv, recOffset + 2 * BYTES_PER_DOUBLE, nCoeffs, le) - const yC = readCoeffs(dv, recOffset + (2 + nCoeffs) * BYTES_PER_DOUBLE, nCoeffs, le) - const zC = readCoeffs(dv, recOffset + (2 + 2 * nCoeffs) * BYTES_PER_DOUBLE, nCoeffs, le) + const xC = readCoeffs(dv, recOffset + 2 * BYTES_PER_DOUBLE, nCoeffs, le); + const yC = readCoeffs(dv, recOffset + (2 + nCoeffs) * BYTES_PER_DOUBLE, nCoeffs, le); + const zC = readCoeffs(dv, recOffset + (2 + 2 * nCoeffs) * BYTES_PER_DOUBLE, nCoeffs, le); - const [px, vx] = chebyshevEvalWithDerivative(xC, x, radius) - const [py, vy] = chebyshevEvalWithDerivative(yC, x, radius) - const [pz, vz] = chebyshevEvalWithDerivative(zC, x, radius) + const [px, vx] = chebyshevEvalWithDerivative(xC, x, radius); + const [py, vy] = chebyshevEvalWithDerivative(yC, x, radius); + const [pz, vz] = chebyshevEvalWithDerivative(zC, x, radius); - return { position: [px, py, pz], velocity: [vx, vy, vz] } + return { position: [px, py, pz], velocity: [vx, vy, vz] }; } /** @@ -312,51 +312,51 @@ export function evaluateType3( et: number, le = true, ): StateVector { - const dv = new DataView(buffer) - const endOffset = seg.dataOffset + seg.dataSize * BYTES_PER_DOUBLE + const dv = new DataView(buffer); + const endOffset = seg.dataOffset + seg.dataSize * BYTES_PER_DOUBLE; - const N = dv.getFloat64(endOffset - BYTES_PER_DOUBLE, le) - const rsize = dv.getFloat64(endOffset - 2 * BYTES_PER_DOUBLE, le) - const intlen = dv.getFloat64(endOffset - 3 * BYTES_PER_DOUBLE, le) - const init = dv.getFloat64(endOffset - 4 * BYTES_PER_DOUBLE, le) + const N = dv.getFloat64(endOffset - BYTES_PER_DOUBLE, le); + const rsize = dv.getFloat64(endOffset - 2 * BYTES_PER_DOUBLE, le); + const intlen = dv.getFloat64(endOffset - 3 * BYTES_PER_DOUBLE, le); + const init = dv.getFloat64(endOffset - 4 * BYTES_PER_DOUBLE, le); - const degree = Math.round((rsize - 2) / 6) - const nCoeffs = degree + 1 + const degree = Math.round((rsize - 2) / 6); + const nCoeffs = degree + 1; - let recIdx = Math.floor((et - init) / intlen) - recIdx = Math.max(0, Math.min(Math.round(N) - 1, recIdx)) + let recIdx = Math.floor((et - init) / intlen); + recIdx = Math.max(0, Math.min(Math.round(N) - 1, recIdx)); - const recOffset = seg.dataOffset + recIdx * rsize * BYTES_PER_DOUBLE - const mid = dv.getFloat64(recOffset, le) - const radius = dv.getFloat64(recOffset + BYTES_PER_DOUBLE, le) - const x = (et - mid) / radius + const recOffset = seg.dataOffset + recIdx * rsize * BYTES_PER_DOUBLE; + const mid = dv.getFloat64(recOffset, le); + const radius = dv.getFloat64(recOffset + BYTES_PER_DOUBLE, le); + const x = (et - mid) / radius; - const xPC = readCoeffs(dv, recOffset + 2 * BYTES_PER_DOUBLE, nCoeffs, le) - const yPC = readCoeffs(dv, recOffset + (2 + nCoeffs) * BYTES_PER_DOUBLE, nCoeffs, le) - const zPC = readCoeffs(dv, recOffset + (2 + 2 * nCoeffs) * BYTES_PER_DOUBLE, nCoeffs, le) - const xVC = readCoeffs(dv, recOffset + (2 + 3 * nCoeffs) * BYTES_PER_DOUBLE, nCoeffs, le) - const yVC = readCoeffs(dv, recOffset + (2 + 4 * nCoeffs) * BYTES_PER_DOUBLE, nCoeffs, le) - const zVC = readCoeffs(dv, recOffset + (2 + 5 * nCoeffs) * BYTES_PER_DOUBLE, nCoeffs, le) + const xPC = readCoeffs(dv, recOffset + 2 * BYTES_PER_DOUBLE, nCoeffs, le); + const yPC = readCoeffs(dv, recOffset + (2 + nCoeffs) * BYTES_PER_DOUBLE, nCoeffs, le); + const zPC = readCoeffs(dv, recOffset + (2 + 2 * nCoeffs) * BYTES_PER_DOUBLE, nCoeffs, le); + const xVC = readCoeffs(dv, recOffset + (2 + 3 * nCoeffs) * BYTES_PER_DOUBLE, nCoeffs, le); + const yVC = readCoeffs(dv, recOffset + (2 + 4 * nCoeffs) * BYTES_PER_DOUBLE, nCoeffs, le); + const zVC = readCoeffs(dv, recOffset + (2 + 5 * nCoeffs) * BYTES_PER_DOUBLE, nCoeffs, le); - const px = chebyshevEvalWithDerivative(xPC, x, radius)[0] - const py = chebyshevEvalWithDerivative(yPC, x, radius)[0] - const pz = chebyshevEvalWithDerivative(zPC, x, radius)[0] + const px = chebyshevEvalWithDerivative(xPC, x, radius)[0]; + const py = chebyshevEvalWithDerivative(yPC, x, radius)[0]; + const pz = chebyshevEvalWithDerivative(zPC, x, radius)[0]; // Type 3: velocity polynomial evaluated at x gives km/s directly - const vx = chebyshevEvalWithDerivative(xVC, x, radius)[0] - const vy = chebyshevEvalWithDerivative(yVC, x, radius)[0] - const vz = chebyshevEvalWithDerivative(zVC, x, radius)[0] + const vx = chebyshevEvalWithDerivative(xVC, x, radius)[0]; + const vy = chebyshevEvalWithDerivative(yVC, x, radius)[0]; + const vz = chebyshevEvalWithDerivative(zVC, x, radius)[0]; - return { position: [px, py, pz], velocity: [vx, vy, vz] } + return { position: [px, py, pz], velocity: [vx, vy, vz] }; } // ─── Helpers ────────────────────────────────────────────────────────────────── function readCoeffs(dv: DataView, offset: number, n: number, le: boolean): Float64Array { - const arr = new Float64Array(n) + const arr = new Float64Array(n); for (let k = 0; k < n; k++) { - arr[k] = dv.getFloat64(offset + k * BYTES_PER_DOUBLE, le) + arr[k] = dv.getFloat64(offset + k * BYTES_PER_DOUBLE, le); } - return arr + return arr; } function subtractSV(a: StateVector, b: StateVector): StateVector { @@ -371,7 +371,7 @@ function subtractSV(a: StateVector, b: StateVector): StateVector { a.velocity[1] - b.velocity[1], a.velocity[2] - b.velocity[2], ], - } + }; } // ─── Leap-second kernel ─────────────────────────────────────────────────────── @@ -381,12 +381,12 @@ function subtractSV(a: StateVector, b: StateVector): StateVector { * Extracts DELTET/DELTA_AT pairs and converts to (JD_UTC, deltaAT) pairs. */ export function parseLsk(text: string): ReadonlyArray { - const results: [number, number][] = [] + const results: [number, number][] = []; - const match = text.match(/DELTET\/DELTA_AT\s*=\s*\(\s*([\s\S]*?)\)/m) - if (!match) return results + const match = text.match(/DELTET\/DELTA_AT\s*=\s*\(\s*([\s\S]*?)\)/m); + if (!match) return results; - const block = match[1]! // regex has one capture group; match[1] is present when match succeeds + const block = match[1]!; // regex has one capture group; match[1] is present when match succeeds const months: Record = { JAN: 1, FEB: 2, @@ -400,22 +400,22 @@ export function parseLsk(text: string): ReadonlyArray OCT: 10, NOV: 11, DEC: 12, - } + }; - const pairRe = /(-?\d+(?:\.\d+)?)\s*,\s*@(\d{4})-([A-Z]{3})-(\d{1,2})/g - let m: RegExpExecArray | null + const pairRe = /(-?\d+(?:\.\d+)?)\s*,\s*@(\d{4})-([A-Z]{3})-(\d{1,2})/g; + let m: RegExpExecArray | null; while ((m = pairRe.exec(block)) !== null) { // m[1]..m[4] are capture groups; regex guarantees they are present when exec matches - const deltaAT = parseFloat(m[1]!) - const year = parseInt(m[2]!) - const month = months[m[3]!] ?? 1 - const day = parseInt(m[4]!) + const deltaAT = parseFloat(m[1]!); + const year = parseInt(m[2]!); + const month = months[m[3]!] ?? 1; + const day = parseInt(m[4]!); // Gregorian to JD (noon = integer JD) - const a = Math.floor((14 - month) / 12) - const y = year + 4800 - a - const mo = month + 12 * a - 3 + const a = Math.floor((14 - month) / 12); + const y = year + 4800 - a; + const mo = month + 12 * a - 3; const jdNoon = day + Math.floor((153 * mo + 2) / 5) + @@ -423,10 +423,10 @@ export function parseLsk(text: string): ReadonlyArray Math.floor(y / 4) - Math.floor(y / 100) + Math.floor(y / 400) - - 32045 + 32045; // Midnight = JD - 0.5 - results.push([jdNoon - 0.5, deltaAT]) + results.push([jdNoon - 0.5, deltaAT]); } - return results + return results; } diff --git a/src/time/index.ts b/src/time/index.ts index a31b484..f4e29bf 100644 --- a/src/time/index.ts +++ b/src/time/index.ts @@ -19,21 +19,21 @@ * Espenak & Meeus — ΔT polynomial expressions */ -import type { TimeScales } from '../types.js' +import type { TimeScales } from "../types.js"; // ─── Constants ──────────────────────────────────────────────────────────────── /** Julian Date of J2000.0 epoch (2000 Jan 1, 12:00 TT) */ -export const J2000 = 2451545.0 +export const J2000 = 2451545.0; /** TT - TAI offset in seconds (exact, by definition) */ -export const TT_MINUS_TAI = 32.184 +export const TT_MINUS_TAI = 32.184; /** Seconds per day */ -export const SECONDS_PER_DAY = 86400.0 +export const SECONDS_PER_DAY = 86400.0; /** Days per Julian century */ -export const DAYS_PER_JULIAN_CENTURY = 36525.0 +export const DAYS_PER_JULIAN_CENTURY = 36525.0; // ─── Leap-second table ──────────────────────────────────────────────────────── @@ -73,19 +73,19 @@ export const LEAP_SECOND_TABLE: ReadonlyArray = [ [2456109.5, 35], // 2012 Jul 1 [2457204.5, 36], // 2015 Jul 1 [2457754.5, 37], // 2017 Jan 1 -] as const +] as const; /** * Get the current leap second count (TAI - UTC) for a given JD in UTC. * Returns 10 for dates before 1972 (the first leap second era). */ export function getDeltaAT(jdUTC: number): number { - let deltaAT = 10 + let deltaAT = 10; for (const [jd, dat] of LEAP_SECOND_TABLE) { - if (jdUTC >= jd) deltaAT = dat - else break + if (jdUTC >= jd) deltaAT = dat; + else break; } - return deltaAT + return deltaAT; } // ─── Julian Date ───────────────────────────────────────────────────────────── @@ -95,14 +95,14 @@ export function getDeltaAT(jdUTC: number): number { * Uses the standard formula; valid for dates after the Gregorian reform. */ export function dateToJD(date: Date): number { - return date.getTime() / 86400000 + 2440587.5 + return date.getTime() / 86400000 + 2440587.5; } /** * Convert a Julian Date in UTC to a JavaScript Date. */ export function jdToDate(jd: number): Date { - return new Date((jd - 2440587.5) * 86400000) + return new Date((jd - 2440587.5) * 86400000); } /** @@ -110,7 +110,7 @@ export function jdToDate(jd: number): Date { * Used as the standard argument for precession and nutation polynomials. */ export function jdTTtoT(jdTT: number): number { - return (jdTT - J2000) / DAYS_PER_JULIAN_CENTURY + return (jdTT - J2000) / DAYS_PER_JULIAN_CENTURY; } // ─── Time scale conversions ─────────────────────────────────────────────────── @@ -129,33 +129,33 @@ export function computeTimeScales( ut1utcOverride?: number, deltaTOverride?: number, ): TimeScales { - const jdUTC = dateToJD(utc) - const deltaAT = getDeltaAT(jdUTC) + const jdUTC = dateToJD(utc); + const deltaAT = getDeltaAT(jdUTC); // UTC → TAI → TT - const jdTAI = jdUTC + deltaAT / SECONDS_PER_DAY - const jdTT = jdTAI + TT_MINUS_TAI / SECONDS_PER_DAY + const jdTAI = jdUTC + deltaAT / SECONDS_PER_DAY; + const jdTT = jdTAI + TT_MINUS_TAI / SECONDS_PER_DAY; // TT → TDB (periodic correction, sub-millisecond) - const tdbCorrection = tdbMinusTT(jdTT) / SECONDS_PER_DAY - const jdTDB = jdTT + tdbCorrection + const tdbCorrection = tdbMinusTT(jdTT) / SECONDS_PER_DAY; + const jdTDB = jdTT + tdbCorrection; // UT1 - let jdUT1: number - let deltaT: number + let jdUT1: number; + let deltaT: number; if (ut1utcOverride !== undefined) { - jdUT1 = jdUTC + ut1utcOverride / SECONDS_PER_DAY - deltaT = (jdTT - jdUT1) * SECONDS_PER_DAY + jdUT1 = jdUTC + ut1utcOverride / SECONDS_PER_DAY; + deltaT = (jdTT - jdUT1) * SECONDS_PER_DAY; } else if (deltaTOverride !== undefined) { - deltaT = deltaTOverride - jdUT1 = jdTT - deltaT / SECONDS_PER_DAY + deltaT = deltaTOverride; + jdUT1 = jdTT - deltaT / SECONDS_PER_DAY; } else { - deltaT = deltaTPolynomial(jdTT) - jdUT1 = jdTT - deltaT / SECONDS_PER_DAY + deltaT = deltaTPolynomial(jdTT); + jdUT1 = jdTT - deltaT / SECONDS_PER_DAY; } - return { utc, jdUTC, jdTT, jdTDB, jdUT1, deltaT, deltaAT } + return { utc, jdUTC, jdTT, jdTDB, jdUT1, deltaT, deltaAT }; } /** @@ -165,8 +165,8 @@ export function computeTimeScales( * gives the proper SPICE-compatible ET value. */ export function jdTTtoET(jdTT: number): number { - const tdbCorr = tdbMinusTT(jdTT) - return (jdTT - J2000) * SECONDS_PER_DAY + tdbCorr + const tdbCorr = tdbMinusTT(jdTT); + return (jdTT - J2000) * SECONDS_PER_DAY + tdbCorr; } /** @@ -179,11 +179,11 @@ export function jdTTtoET(jdTT: number): number { * Maximum error: ~30 microseconds (acceptable for crescent work). */ export function tdbMinusTT(jdTT: number): number { - const d = jdTT - J2000 + const d = jdTT - J2000; // Mean anomaly of the Sun (degrees) - const gDeg = 357.53 + 0.9856003 * d - const g = (gDeg * Math.PI) / 180 - return 0.001658 * Math.sin(g) + 0.000014 * Math.sin(2 * g) + const gDeg = 357.53 + 0.9856003 * d; + const g = (gDeg * Math.PI) / 180; + return 0.001658 * Math.sin(g) + 0.000014 * Math.sin(2 * g); } /** @@ -194,13 +194,13 @@ export function tdbMinusTT(jdTT: number): number { */ export function deltaTPolynomial(jdTT: number): number { // Convert JD to decimal year - const y = 2000 + (jdTT - J2000) / 365.25 + const y = 2000 + (jdTT - J2000) / 365.25; if (y < -500) { - const u = (y - 1820) / 100 - return -20 + 32 * u * u + const u = (y - 1820) / 100; + return -20 + 32 * u * u; } else if (y < 500) { - const u = y / 100 + const u = y / 100; return ( 10583.6 - 1014.41 * u + @@ -209,9 +209,9 @@ export function deltaTPolynomial(jdTT: number): number { 0.1798452 * u ** 4 + 0.022174192 * u ** 5 + 0.0090316521 * u ** 6 - ) + ); } else if (y < 1600) { - const u = (y - 1000) / 100 + const u = (y - 1000) / 100; return ( 1574.2 - 556.01 * u + @@ -220,15 +220,15 @@ export function deltaTPolynomial(jdTT: number): number { 0.8503463 * u ** 4 - 0.005050998 * u ** 5 + 0.0083572073 * u ** 6 - ) + ); } else if (y < 1700) { - const t = y - 1600 - return 120 - 0.9808 * t - 0.01532 * t * t + t ** 3 / 7129 + const t = y - 1600; + return 120 - 0.9808 * t - 0.01532 * t * t + t ** 3 / 7129; } else if (y < 1800) { - const t = y - 1700 - return 8.83 + 0.1603 * t - 0.0059285 * t * t + 0.00013336 * t ** 3 - t ** 4 / 1174000 + const t = y - 1700; + return 8.83 + 0.1603 * t - 0.0059285 * t * t + 0.00013336 * t ** 3 - t ** 4 / 1174000; } else if (y < 1860) { - const t = y - 1800 + const t = y - 1800; return ( 13.72 - 0.332447 * t + @@ -238,9 +238,9 @@ export function deltaTPolynomial(jdTT: number): number { 0.0000121272 * t ** 5 - 0.0000001699 * t ** 6 + 0.000000000875 * t ** 7 - ) + ); } else if (y < 1900) { - const t = y - 1860 + const t = y - 1860; return ( 7.62 + 0.5737 * t - @@ -248,21 +248,21 @@ export function deltaTPolynomial(jdTT: number): number { 0.01680668 * t ** 3 - 0.0004473624 * t ** 4 + t ** 5 / 233174 - ) + ); } else if (y < 1920) { - const t = y - 1900 - return -2.79 + 1.494119 * t - 0.0598939 * t * t + 0.0061966 * t ** 3 - 0.000197 * t ** 4 + const t = y - 1900; + return -2.79 + 1.494119 * t - 0.0598939 * t * t + 0.0061966 * t ** 3 - 0.000197 * t ** 4; } else if (y < 1941) { - const t = y - 1920 - return 21.2 + 0.84493 * t - 0.0761 * t * t + 0.0020936 * t ** 3 + const t = y - 1920; + return 21.2 + 0.84493 * t - 0.0761 * t * t + 0.0020936 * t ** 3; } else if (y < 1961) { - const t = y - 1950 - return 29.07 + 0.407 * t - (t * t) / 233 + t ** 3 / 2547 + const t = y - 1950; + return 29.07 + 0.407 * t - (t * t) / 233 + t ** 3 / 2547; } else if (y < 1986) { - const t = y - 1975 - return 45.45 + 1.067 * t - (t * t) / 260 - t ** 3 / 718 + const t = y - 1975; + return 45.45 + 1.067 * t - (t * t) / 260 - t ** 3 / 718; } else if (y < 2005) { - const t = y - 2000 + const t = y - 2000; return ( 63.86 + 0.3345 * t - @@ -270,14 +270,14 @@ export function deltaTPolynomial(jdTT: number): number { 0.0017275 * t ** 3 + 0.000651814 * t ** 4 + 0.00002373599 * t ** 5 - ) + ); } else if (y < 2050) { - const t = y - 2000 - return 62.92 + 0.32217 * t + 0.005589 * t * t + const t = y - 2000; + return 62.92 + 0.32217 * t + 0.005589 * t * t; } else if (y < 2150) { - return -20 + 32 * ((y - 1820) / 100) ** 2 - 0.5628 * (2150 - y) + return -20 + 32 * ((y - 1820) / 100) ** 2 - 0.5628 * (2150 - y); } else { - const u = (y - 1820) / 100 - return -20 + 32 * u * u + const u = (y - 1820) / 100; + return -20 + 32 * u * u; } } diff --git a/src/types.ts b/src/types.ts index 5669891..34ae588 100644 --- a/src/types.ts +++ b/src/types.ts @@ -1,20 +1,20 @@ // ─── Primitive geometry ────────────────────────────────────────────────────── /** 3-element position or velocity vector in km or km/s */ -export type Vec3 = [number, number, number] +export type Vec3 = [number, number, number]; /** Position + velocity state vector from the ephemeris */ export interface StateVector { - position: Vec3 // km, in the frame determined by context - velocity: Vec3 // km/s + position: Vec3; // km, in the frame determined by context + velocity: Vec3; // km/s } /** Azimuth + altitude in degrees */ export interface AzAlt { /** Degrees from North, measured clockwise (0 = N, 90 = E, 180 = S, 270 = W) */ - azimuth: number + azimuth: number; /** Degrees above the horizon (negative = below) */ - altitude: number + altitude: number; } // ─── Kernel-free moon results ───────────────────────────────────────────────── @@ -26,17 +26,17 @@ export interface AzAlt { */ export interface MoonPosition { /** Azimuth in degrees from North, measured clockwise (0 = N, 90 = E, 180 = S, 270 = W) */ - azimuth: number + azimuth: number; /** Apparent altitude in degrees above the horizon (atmospheric refraction applied) */ - altitude: number + altitude: number; /** Distance from Earth center to Moon center, km */ - distance: number + distance: number; /** * Parallactic angle in radians. * The angle between the great circle through the Moon and zenith, and the great circle * through the Moon and the north celestial pole. Positive east of the meridian. */ - parallacticAngle: number + parallacticAngle: number; } /** @@ -46,38 +46,38 @@ export interface MoonPosition { */ export interface MoonIlluminationResult { /** Illuminated fraction of the Moon disk, 0 (new moon) to 1 (full moon) */ - fraction: number + fraction: number; /** * Phase cycle fraction in [0, 1): * 0 = new moon, 0.25 = first quarter, 0.5 = full moon, 0.75 = last quarter */ - phase: number + phase: number; /** * Position angle of the midpoint of the bright limb, measured eastward from * the north celestial pole, in radians. Matches the suncalc convention. */ - angle: number + angle: number; /** True while elongation is increasing (new moon toward full moon) */ - isWaxing: boolean + isWaxing: boolean; } // ─── Time ──────────────────────────────────────────────────────────────────── /** All relevant time scale values for a single moment */ export interface TimeScales { - utc: Date + utc: Date; /** Julian Date in UTC */ - jdUTC: number + jdUTC: number; /** Julian Date in Terrestrial Time (TT = TAI + 32.184s) */ - jdTT: number + jdTT: number; /** Julian Date in Barycentric Dynamical Time (used by JPL ephemerides) */ - jdTDB: number + jdTDB: number; /** Julian Date in UT1 (Earth rotation time) */ - jdUT1: number + jdUT1: number; /** TT - UT1 in seconds (delta-T) */ - deltaT: number + deltaT: number; /** TAI - UTC in seconds (leap seconds count) */ - deltaAT: number + deltaAT: number; } // ─── Observer ──────────────────────────────────────────────────────────────── @@ -85,28 +85,28 @@ export interface TimeScales { /** Observer location and environmental parameters */ export interface Observer { /** Geodetic latitude in degrees (north positive) */ - lat: number + lat: number; /** Longitude in degrees (east positive) */ - lon: number + lon: number; /** Height above WGS84 ellipsoid in meters */ - elevation: number + elevation: number; /** Optional label for the location */ - name?: string + name?: string; /** * Override TT - UT1 in seconds. * When provided, used directly. Otherwise the built-in polynomial is used. * For maximum accuracy, supply the current IERS value (typically within ±0.9s). */ - deltaT?: number + deltaT?: number; /** * Override UT1 - UTC in seconds (from IERS Bulletin A). * Takes precedence over deltaT when both are provided. */ - ut1utc?: number + ut1utc?: number; /** Atmospheric pressure in millibars (default 1013.25) */ - pressure?: number + pressure?: number; /** Ambient temperature in Celsius (default 15) */ - temperature?: number + temperature?: number; } // ─── Crescent geometry ─────────────────────────────────────────────────────── @@ -117,31 +117,31 @@ export interface Observer { */ export interface CrescentGeometry { /** Arc of light: topocentric Sun-Moon angular separation (elongation), degrees */ - ARCL: number + ARCL: number; /** * Arc of vision: Moon airless altitude minus Sun airless altitude, degrees. * Used as the primary visibility discriminant in both Yallop and Odeh. */ - ARCV: number + ARCV: number; /** * Relative azimuth: Sun azimuth minus Moon azimuth, normalized to [-180, 180], degrees. * Positive = Moon north of Sun. */ - DAZ: number + DAZ: number; /** * Topocentric crescent width in arc minutes. * Used directly in Odeh's polynomial V expression. */ - W: number + W: number; /** Moonset minus sunset in minutes. Negative = Moon sets before Sun (no sighting possible). */ - lag: number + lag: number; } // ─── Yallop q-test ─────────────────────────────────────────────────────────── /** Yallop q-test visibility category (NAO Technical Note 69) */ /** Yallop visibility category (A = easily visible, F = below Danjon limit). */ -export type YallopCategory = 'A' | 'B' | 'C' | 'D' | 'E' | 'F' +export type YallopCategory = "A" | "B" | "C" | "D" | "E" | "F"; /** * Published q thresholds (Yallop 1997, NAO TN 69): @@ -158,43 +158,43 @@ export const YALLOP_THRESHOLDS = { C: -0.16, D: -0.232, E: -0.293, -} as const +} as const; /** * Human-readable descriptions for each Yallop visibility category (A–F). * Sourced from Yallop (NAO TN 69, 1997). */ export const YALLOP_DESCRIPTIONS: Record = { - A: 'Easily visible to the naked eye', - B: 'Visible under perfect conditions', - C: 'May need optical aid to find; naked eye possible', - D: 'Optical aid needed; naked eye not possible', - E: 'Not visible even with telescope under good conditions', - F: 'Below Danjon limit — crescent cannot form', -} + A: "Easily visible to the naked eye", + B: "Visible under perfect conditions", + C: "May need optical aid to find; naked eye possible", + D: "Optical aid needed; naked eye not possible", + E: "Not visible even with telescope under good conditions", + F: "Below Danjon limit — crescent cannot form", +}; export interface YallopResult { /** The continuous q parameter (higher = more visible) */ - q: number + q: number; /** Visibility category A through F */ - category: YallopCategory + category: YallopCategory; /** Human-readable interpretation */ - description: string + description: string; /** True for categories A and B */ - isVisibleNakedEye: boolean + isVisibleNakedEye: boolean; /** True for categories C and D */ - requiresOpticalAid: boolean + requiresOpticalAid: boolean; /** True for category F */ - isBelowDanjonLimit: boolean + isBelowDanjonLimit: boolean; /** Topocentric crescent width W' used in the q formula, arc minutes */ - Wprime: number + Wprime: number; } // ─── Odeh criterion ────────────────────────────────────────────────────────── /** Odeh visibility zone (Experimental Astronomy 2006) */ /** Odeh visibility zone (A = naked eye visible, D = not visible with any aid). */ -export type OdehZone = 'A' | 'B' | 'C' | 'D' +export type OdehZone = "A" | "B" | "C" | "D"; /** * Published V thresholds (Odeh 2006): @@ -207,33 +207,33 @@ export const ODEH_THRESHOLDS = { A: 5.65, B: 2.0, C: -0.96, -} as const +} as const; /** * Human-readable descriptions for each Odeh visibility zone (A–D). * Sourced from Odeh (Experimental Astronomy, 2006). */ export const ODEH_DESCRIPTIONS: Record = { - A: 'Visible with naked eye', - B: 'Visible with optical aid; may be seen with naked eye under excellent conditions', - C: 'Visible with optical aid only', - D: 'Not visible even with optical aid', -} + A: "Visible with naked eye", + B: "Visible with optical aid; may be seen with naked eye under excellent conditions", + C: "Visible with optical aid only", + D: "Not visible even with optical aid", +}; export interface OdehResult { /** * Continuous visibility parameter V = ARCV - (arcv_minimum(W)). * Positive = crescent exceeds minimum visibility threshold. */ - V: number + V: number; /** Visibility zone A through D */ - zone: OdehZone + zone: OdehZone; /** Human-readable interpretation */ - description: string + description: string; /** True for zone A */ - isVisibleNakedEye: boolean + isVisibleNakedEye: boolean; /** True for zones A and B */ - isVisibleWithOpticalAid: boolean + isVisibleWithOpticalAid: boolean; } // ─── Kernel-free visibility estimate ───────────────────────────────────────── @@ -248,25 +248,25 @@ export interface MoonVisibilityEstimate { * Odeh V parameter: V = ARCV − f(W). * Positive = crescent exceeds minimum visibility threshold. */ - V: number + V: number; /** Visibility zone A through D */ - zone: OdehZone + zone: OdehZone; /** Human-readable zone description */ - description: string + description: string; /** True for zone A */ - isVisibleNakedEye: boolean + isVisibleNakedEye: boolean; /** True for zones A and B */ - isVisibleWithOpticalAid: boolean + isVisibleWithOpticalAid: boolean; /** Arc of light (Sun-Moon elongation) in degrees */ - ARCL: number + ARCL: number; /** Arc of vision (Moon airless altitude minus Sun airless altitude) in degrees */ - ARCV: number + ARCV: number; /** Topocentric crescent width in arc minutes */ - W: number + W: number; /** True when Moon is above the horizon at the given time */ - moonAboveHorizon: boolean + moonAboveHorizon: boolean; /** Always true: computed via Meeus approximation, not DE442S */ - isApproximate: true + isApproximate: true; } /** @@ -276,118 +276,118 @@ export interface MoonVisibilityEstimate { */ export interface MoonSnapshot { /** Phase name, illumination, age, and next events */ - phase: MoonPhaseResult + phase: MoonPhaseResult; /** Topocentric az/alt, distance, parallactic angle */ - position: MoonPosition + position: MoonPosition; /** Illumination fraction, phase cycle, bright limb angle, waxing/waning */ - illumination: MoonIlluminationResult + illumination: MoonIlluminationResult; /** Quick Odeh-based crescent visibility estimate */ - visibility: MoonVisibilityEstimate + visibility: MoonVisibilityEstimate; } // ─── Moon phase ────────────────────────────────────────────────────────────── export type MoonPhaseName = - | 'new-moon' - | 'waxing-crescent' - | 'first-quarter' - | 'waxing-gibbous' - | 'full-moon' - | 'waning-gibbous' - | 'last-quarter' - | 'waning-crescent' + | "new-moon" + | "waxing-crescent" + | "first-quarter" + | "waxing-gibbous" + | "full-moon" + | "waning-gibbous" + | "last-quarter" + | "waning-crescent"; export interface MoonPhaseResult { /** Named phase based on illumination and waxing/waning state */ - phase: MoonPhaseName + phase: MoonPhaseName; /** Human-readable phase name, e.g. "Waxing Crescent" */ - phaseName: string + phaseName: string; /** Moon phase emoji symbol, e.g. "🌒" */ - phaseSymbol: string + phaseSymbol: string; /** Illuminated fraction 0-100 (percent) */ - illumination: number + illumination: number; /** Hours since last new moon */ - age: number + age: number; /** Ecliptic longitude of the Moon minus the Sun, degrees [0, 360) */ - elongationDeg: number + elongationDeg: number; /** True when Moon is moving away from the Sun (illumination increasing) */ - isWaxing: boolean + isWaxing: boolean; /** UTC date of the next new moon */ - nextNewMoon: Date + nextNewMoon: Date; /** UTC date of the next full moon */ - nextFullMoon: Date + nextFullMoon: Date; /** UTC date of the previous new moon */ - prevNewMoon: Date + prevNewMoon: Date; } // ─── Event times ───────────────────────────────────────────────────────────── export interface SunMoonEvents { /** UTC time of sunset for the given date at the observer's location */ - sunsetUTC: Date | null + sunsetUTC: Date | null; /** UTC time of moonset for the given date at the observer's location */ - moonsetUTC: Date | null + moonsetUTC: Date | null; /** UTC time of sunrise */ - sunriseUTC: Date | null + sunriseUTC: Date | null; /** UTC time of moonrise */ - moonriseUTC: Date | null + moonriseUTC: Date | null; /** UTC time when civil twilight ends (Sun at -6°) */ - civilTwilightEndUTC: Date | null + civilTwilightEndUTC: Date | null; /** UTC time when nautical twilight ends (Sun at -12°) */ - nauticalTwilightEndUTC: Date | null + nauticalTwilightEndUTC: Date | null; /** UTC time when astronomical twilight ends (Sun at -18°) */ - astronomicalTwilightEndUTC: Date | null + astronomicalTwilightEndUTC: Date | null; } // ─── Full moon sighting report ──────────────────────────────────────────────── export interface MoonSightingReport { /** Date for which the sighting report was computed */ - date: Date + date: Date; /** Observer location used */ - observer: Observer + observer: Observer; // Event times - sunsetUTC: Date | null - moonsetUTC: Date | null + sunsetUTC: Date | null; + moonsetUTC: Date | null; /** Moonset minus sunset, in minutes. Null if either event is null. */ - lagMinutes: number | null + lagMinutes: number | null; /** Best observation time (Odeh/Yallop: T_s + 4/9 * Lag) */ - bestTimeUTC: Date | null + bestTimeUTC: Date | null; /** Conservative observation window [bestTime - 20min, bestTime + 20min] */ - bestTimeWindowUTC: [Date, Date] | null + bestTimeWindowUTC: [Date, Date] | null; // At best time /** Topocentric Moon position at best time */ - moonPosition: AzAlt | null + moonPosition: AzAlt | null; /** Topocentric Sun position at best time */ - sunPosition: AzAlt | null + sunPosition: AzAlt | null; /** Moon illumination percent at best time */ - illumination: number | null + illumination: number | null; /** Hours since conjunction (new moon) */ - moonAge: number | null + moonAge: number | null; // Crescent geometry at best time - geometry: CrescentGeometry | null + geometry: CrescentGeometry | null; // Visibility criteria results - yallop: YallopResult | null - odeh: OdehResult | null + yallop: YallopResult | null; + odeh: OdehResult | null; // Sighting guidance /** * Plain-language direction for observers. * Includes where to look (azimuth, altitude), when (best time), and what to expect. */ - guidance: string + guidance: string; // Metadata /** Source ephemeris used for this calculation */ - ephemerisSource: 'DE442S' | 'approximate' + ephemerisSource: "DE442S" | "approximate"; /** Whether the Moon is even above the horizon at best time */ - moonAboveHorizon: boolean | null + moonAboveHorizon: boolean | null; /** Whether sighting is geometrically possible (lag > 0, Moon above horizon at best time) */ - sightingPossible: boolean + sightingPossible: boolean; } // ─── Kernel configuration ───────────────────────────────────────────────────── @@ -397,40 +397,40 @@ export interface MoonSightingReport { * Used for both the planetary SPK (de442s.bsp) and leap-second kernel (naif0012.tls). */ export type KernelSource = - | { type: 'file'; path: string } - | { type: 'buffer'; data: ArrayBuffer; name: string } - | { type: 'url'; url: string } - | { type: 'auto' } // auto-download from NAIF, cache in ~/.cache/moon-sighting + | { type: "file"; path: string } + | { type: "buffer"; data: ArrayBuffer; name: string } + | { type: "url"; url: string } + | { type: "auto" }; // auto-download from NAIF, cache in ~/.cache/moon-sighting export interface KernelConfig { /** Planetary SPK kernel — defaults to de442s.bsp via auto-download */ - planetary?: KernelSource + planetary?: KernelSource; /** Leap-second kernel — defaults to naif0012.tls via auto-download */ - leapSeconds?: KernelSource + leapSeconds?: KernelSource; /** * Directory for the download cache. * Defaults to ~/.cache/moon-sighting on POSIX, %LOCALAPPDATA%\moon-sighting on Windows. */ - cacheDir?: string + cacheDir?: string; /** * SHA-256 checksum of de442s.bsp for download verification. * Bundled default matches the NAIF distribution as of 2024. */ - checksumOverride?: string + checksumOverride?: string; } // ─── Top-level options ──────────────────────────────────────────────────────── export interface SightingOptions { /** Kernel acquisition configuration. Defaults to auto-download. */ - kernels?: KernelConfig + kernels?: KernelConfig; /** * Best-time computation method. * 'heuristic' — T_b = T_sunset + 4/9 * Lag (Odeh/Yallop approximation, fast) * 'optimized' — scan sunset-to-moonset interval, maximize Odeh V parameter * Default: 'heuristic' */ - bestTimeMethod?: 'heuristic' | 'optimized' + bestTimeMethod?: "heuristic" | "optimized"; } // ─── WGS84 constants ───────────────────────────────────────────────────────── @@ -445,44 +445,44 @@ export const WGS84 = { f: 1 / 298.257223563, /** Semi-minor axis in meters */ get b() { - return this.a * (1 - this.f) + return this.a * (1 - this.f); }, /** First eccentricity squared */ get e2() { - return 2 * this.f - this.f * this.f + return 2 * this.f - this.f * this.f; }, -} as const +} as const; // ─── Internal ephemeris types ───────────────────────────────────────────────── /** A segment in a JPL SPK (DAF) kernel file */ export interface SpkSegment { /** NAIF body ID of the target body */ - target: number + target: number; /** NAIF body ID of the center body */ - center: number + center: number; /** Reference frame code */ - frame: number + frame: number; /** SPK data type (2 = Chebyshev position only, 3 = Chebyshev position + velocity) */ - dataType: 2 | 3 + dataType: 2 | 3; /** Segment start time in ET seconds past J2000 */ - startET: number + startET: number; /** Segment end time in ET seconds past J2000 */ - endET: number + endET: number; /** Byte offset of the data array in the file */ - dataOffset: number + dataOffset: number; /** Number of double-precision values in the data array */ - dataSize: number + dataSize: number; } /** A decoded Chebyshev record from a Type 2 or Type 3 SPK segment */ export interface ChebRecord { /** Midpoint of the record interval in ET seconds past J2000 */ - mid: number + mid: number; /** Half-width of the record interval in seconds */ - radius: number + radius: number; /** Chebyshev coefficients for X, Y, Z [3][degree+1] */ - coeffs: Float64Array[] + coeffs: Float64Array[]; /** Degree of the polynomial */ - degree: number + degree: number; } diff --git a/src/visibility/index.ts b/src/visibility/index.ts index 88c64f3..0321448 100644 --- a/src/visibility/index.ts +++ b/src/visibility/index.ts @@ -30,15 +30,15 @@ import type { YallopCategory, OdehResult, OdehZone, -} from '../types.js' +} from "../types.js"; import { YALLOP_THRESHOLDS, YALLOP_DESCRIPTIONS, ODEH_THRESHOLDS, ODEH_DESCRIPTIONS, -} from '../types.js' -import { angularSep } from '../math/index.js' -import { computeCrescentWidth } from '../bodies/index.js' +} from "../types.js"; +import { angularSep } from "../math/index.js"; +import { computeCrescentWidth } from "../bodies/index.js"; // ─── Shared polynomial ──────────────────────────────────────────────────────── @@ -56,7 +56,7 @@ import { computeCrescentWidth } from '../bodies/index.js' * @returns Minimum ARCV required for detection, in degrees */ export function arcvMinimum(W: number): number { - return 11.8371 - 6.3226 * W + 0.7319 * W * W - 0.1018 * W * W * W + return 11.8371 - 6.3226 * W + 0.7319 * W * W - 0.1018 * W * W * W; } // ─── Yallop q-test ──────────────────────────────────────────────────────────── @@ -74,7 +74,7 @@ export function arcvMinimum(W: number): number { * @returns q parameter (continuous) */ export function computeYallopQ(ARCV: number, Wprime: number): number { - return (ARCV - arcvMinimum(Wprime)) / 10 + return (ARCV - arcvMinimum(Wprime)) / 10; } /** @@ -89,12 +89,12 @@ export function computeYallopQ(ARCV: number, Wprime: number): number { * F: q <= -0.293 */ export function yallopCategory(q: number): YallopCategory { - if (q > YALLOP_THRESHOLDS.A) return 'A' - if (q > YALLOP_THRESHOLDS.B) return 'B' - if (q > YALLOP_THRESHOLDS.C) return 'C' - if (q > YALLOP_THRESHOLDS.D) return 'D' - if (q > YALLOP_THRESHOLDS.E) return 'E' - return 'F' + if (q > YALLOP_THRESHOLDS.A) return "A"; + if (q > YALLOP_THRESHOLDS.B) return "B"; + if (q > YALLOP_THRESHOLDS.C) return "C"; + if (q > YALLOP_THRESHOLDS.D) return "D"; + if (q > YALLOP_THRESHOLDS.E) return "E"; + return "F"; } /** @@ -104,18 +104,18 @@ export function yallopCategory(q: number): YallopCategory { * @param Wprime - Topocentric crescent width in arc minutes (may differ from geometry.W) */ export function computeYallop(geometry: CrescentGeometry, Wprime: number): YallopResult { - const q = computeYallopQ(geometry.ARCV, Wprime) - const category = yallopCategory(q) + const q = computeYallopQ(geometry.ARCV, Wprime); + const category = yallopCategory(q); return { q, category, description: YALLOP_DESCRIPTIONS[category], - isVisibleNakedEye: category === 'A' || category === 'B', - requiresOpticalAid: category === 'C' || category === 'D', - isBelowDanjonLimit: category === 'F', + isVisibleNakedEye: category === "A" || category === "B", + requiresOpticalAid: category === "C" || category === "D", + isBelowDanjonLimit: category === "F", Wprime, - } + }; } // ─── Odeh criterion ─────────────────────────────────────────────────────────── @@ -135,7 +135,7 @@ export function computeYallop(geometry: CrescentGeometry, Wprime: number): Yallo * @param W - Topocentric crescent width in arc minutes */ export function computeOdehV(ARCV: number, W: number): number { - return ARCV - arcvMinimum(W) + return ARCV - arcvMinimum(W); } /** @@ -148,10 +148,10 @@ export function computeOdehV(ARCV: number, W: number): number { * D: V < -0.96 — Not visible even with optical aid */ export function odehZone(V: number): OdehZone { - if (V >= ODEH_THRESHOLDS.A) return 'A' - if (V >= ODEH_THRESHOLDS.B) return 'B' - if (V >= ODEH_THRESHOLDS.C) return 'C' - return 'D' + if (V >= ODEH_THRESHOLDS.A) return "A"; + if (V >= ODEH_THRESHOLDS.B) return "B"; + if (V >= ODEH_THRESHOLDS.C) return "C"; + return "D"; } /** @@ -159,16 +159,16 @@ export function odehZone(V: number): OdehZone { * Uses geometry.W directly as the Odeh topocentric crescent width. */ export function computeOdeh(geometry: CrescentGeometry): OdehResult { - const V = computeOdehV(geometry.ARCV, geometry.W) - const zone = odehZone(V) + const V = computeOdehV(geometry.ARCV, geometry.W); + const zone = odehZone(V); return { V, zone, description: ODEH_DESCRIPTIONS[zone], - isVisibleNakedEye: zone === 'A', - isVisibleWithOpticalAid: zone === 'A' || zone === 'B', - } + isVisibleNakedEye: zone === "A", + isVisibleWithOpticalAid: zone === "A" || zone === "B", + }; } // ─── Geometry computation ───────────────────────────────────────────────────── @@ -190,30 +190,30 @@ export function computeOdeh(geometry: CrescentGeometry): OdehResult { export function computeCrescentGeometry( moonAirlessAzAlt: { azimuth: number; altitude: number }, sunAirlessAzAlt: { azimuth: number; altitude: number }, - moonGCRS: import('../types.js').Vec3, - sunGCRS: import('../types.js').Vec3, + moonGCRS: import("../types.js").Vec3, + sunGCRS: import("../types.js").Vec3, sunsetUTC: Date, moonsetUTC: Date, ): CrescentGeometry { // ARCV: airless arc of vision (Moon altitude minus Sun altitude) - const ARCV = moonAirlessAzAlt.altitude - sunAirlessAzAlt.altitude + const ARCV = moonAirlessAzAlt.altitude - sunAirlessAzAlt.altitude; // DAZ: Sun azimuth minus Moon azimuth, normalized to (−180, 180] - let DAZ = sunAirlessAzAlt.azimuth - moonAirlessAzAlt.azimuth - if (DAZ > 180) DAZ -= 360 - if (DAZ < -180) DAZ += 360 + let DAZ = sunAirlessAzAlt.azimuth - moonAirlessAzAlt.azimuth; + if (DAZ > 180) DAZ -= 360; + if (DAZ < -180) DAZ += 360; // ARCL: topocentric Sun-Moon angular separation in degrees // angularSep returns radians; both vectors must be topocentric for accurate ARCL - const ARCL = angularSep(moonGCRS, sunGCRS) * (180 / Math.PI) + const ARCL = angularSep(moonGCRS, sunGCRS) * (180 / Math.PI); // W: topocentric crescent width in arc minutes - const { W } = computeCrescentWidth(moonGCRS, ARCL) + const { W } = computeCrescentWidth(moonGCRS, ARCL); // lag: moonset minus sunset in minutes (negative = Moon sets before Sun) - const lag = (moonsetUTC.getTime() - sunsetUTC.getTime()) / 60000 + const lag = (moonsetUTC.getTime() - sunsetUTC.getTime()) / 60000; - return { ARCL, ARCV, DAZ, W, lag } + return { ARCL, ARCV, DAZ, W, lag }; } // ─── Guidance text ──────────────────────────────────────────────────────────── @@ -237,22 +237,22 @@ export function buildGuidanceText( bestTimeUTC: Date, lagMinutes: number, ): string { - const direction = azimuthToCardinal(moonAz) + const direction = azimuthToCardinal(moonAz); const timeStr = bestTimeUTC .toISOString() - .replace('T', ' ') - .replace(/\.\d+Z$/, ' UTC') - const lagStr = `${Math.round(lagMinutes)} min after sunset` + .replace("T", " ") + .replace(/\.\d+Z$/, " UTC"); + const lagStr = `${Math.round(lagMinutes)} min after sunset`; - let visibility: string + let visibility: string; if (yallop.isVisibleNakedEye && odeh.isVisibleNakedEye) { - visibility = 'should be visible to the naked eye' + visibility = "should be visible to the naked eye"; } else if (odeh.isVisibleWithOpticalAid) { - visibility = 'may require binoculars or a telescope to spot' + visibility = "may require binoculars or a telescope to spot"; } else if (yallop.isBelowDanjonLimit) { - visibility = 'is too close to the Sun to form a visible crescent (below Danjon limit)' + visibility = "is too close to the Sun to form a visible crescent (below Danjon limit)"; } else { - visibility = 'is not expected to be visible even with optical aid' + visibility = "is not expected to be visible even with optical aid"; } return ( @@ -261,29 +261,29 @@ export function buildGuidanceText( `The crescent ${visibility}. ` + `Yallop: ${yallop.category} (${yallop.description}). ` + `Odeh: ${odeh.zone} (${odeh.description}).` - ) + ); } /** Convert azimuth degrees to a cardinal/intercardinal direction label */ function azimuthToCardinal(az: number): string { const dirs = [ - 'North', - 'NNE', - 'NE', - 'ENE', - 'East', - 'ESE', - 'SE', - 'SSE', - 'South', - 'SSW', - 'SW', - 'WSW', - 'West', - 'WNW', - 'NW', - 'NNW', - ] - const idx = Math.round(az / 22.5) % 16 - return dirs[(idx + 16) % 16]! // dirs has 16 elements; (idx+16)%16 is always 0..15 + "North", + "NNE", + "NE", + "ENE", + "East", + "ESE", + "SE", + "SSE", + "South", + "SSW", + "SW", + "WSW", + "West", + "WNW", + "NW", + "NNW", + ]; + const idx = Math.round(az / 22.5) % 16; + return dirs[(idx + 16) % 16]!; // dirs has 16 elements; (idx+16)%16 is always 0..15 }