aric/moon-sighting
1
0
Fork 0
mirror of https://github.com/acamarata/moon-sighting.git synced 2026-06-30 19:04:24 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Initial release: moon-calc v1.0.0

Browse source 
High-accuracy lunar crescent visibility using JPL DE442S ephemerides.
Implements Yallop q-test and Odeh V-parameter with all five crescent
geometry quantities (ARCL, ARCV, DAZ, W, lag). Full time scale chain
(UTC → TDB), IERS Q·R·W frame transforms, WGS84 observer model,
Bennett refraction, and kernel-free moon phase via Meeus approximation.
This commit is contained in:
Aric Camarata 2026-02-25 15:45:41 -05:00
commit b46ba2a74c
38 changed files with 7205 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

18
.editorconfig Normal file
View file

@ -0,0 +1,18 @@
root = true
[*]
end_of_line = lf
charset = utf-8
trim_trailing_whitespace = true
insert_final_newline = true
[*.{js,mjs,cjs,ts,mts,cts,json,yaml,yml,md}]
indent_style = space
indent_size = 2
[*.{c,h}]
indent_style = space
indent_size = 4
[Makefile]
indent_style = tab

70
.github/workflows/ci.yml vendored Normal file
View file

@ -0,0 +1,70 @@
name: CI
on:
push:
branches: [main]
pull_request:
branches: [main]
jobs:
test:
name: Test (Node ${{ matrix.node }})
runs-on: ubuntu-latest
strategy:
matrix:
node: [20, 22, 24]
steps:
- uses: actions/checkout@v4
- uses: pnpm/action-setup@v4
with:
version: 10
- uses: actions/setup-node@v4
with:
node-version: ${{ matrix.node }}
cache: pnpm
- run: pnpm install
- run: pnpm build
- run: node test.mjs
- run: node test-cjs.cjs
typecheck:
name: TypeScript
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4
- uses: pnpm/action-setup@v4
with:
version: 10
- uses: actions/setup-node@v4
with:
node-version: 24
cache: pnpm
- run: pnpm install
- run: pnpm run typecheck
pack-check:
name: Pack check
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4
- uses: pnpm/action-setup@v4
with:
version: 10
- uses: actions/setup-node@v4
with:
node-version: 24
cache: pnpm
- run: pnpm install
- run: pnpm build
- name: Verify pack contents
run: |
npm pack --dry-run 2>&1 | tee pack-output.txt
# Verify expected files are present
grep -q "dist/index.cjs" pack-output.txt || (echo "Missing dist/index.cjs" && exit 1)
grep -q "dist/index.mjs" pack-output.txt || (echo "Missing dist/index.mjs" && exit 1)
grep -q "dist/index.d.ts" pack-output.txt || (echo "Missing dist/index.d.ts" && exit 1)
grep -q "README.md" pack-output.txt || (echo "Missing README.md" && exit 1)
# Verify no test files or .gitignore items are included
! grep -q "test.mjs" pack-output.txt || (echo "test.mjs should not be in pack" && exit 1)
! grep -q "node_modules" pack-output.txt || (echo "node_modules should not be in pack" && exit 1)
echo "Pack contents verified."

23
.github/workflows/wiki-sync.yml vendored Normal file
View file

@ -0,0 +1,23 @@
name: Sync Wiki
on:
push:
branches: [main]
paths:
- '.wiki/**'
jobs:
sync:
name: Sync .wiki/ to GitHub Wiki
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4
- name: Sync wiki pages
uses: newrelic/wiki-sync-action@v1.0.1
with:
source: .wiki
destination: wiki
token: ${{ secrets.GITHUB_TOKEN }}
gitAuthorName: github-actions[bot]
gitAuthorEmail: 41898282+github-actions[bot]@users.noreply.github.com

18
.gitignore vendored Normal file
View file

@ -0,0 +1,18 @@
node_modules/
dist/
*.tgz
*.log
.DS_Store
.env
.env.*
# Ephemeris kernel cache (large binary files)
.kernel-cache/
*.bsp
*.tls
# AI agent directories
.claude/
.cursor/
.aider/
.continue/

0
.npmrc Normal file
View file

1
.nvmrc Normal file
View file

@ -0,0 +1 @@
24

249
.wiki/API-Reference.md Normal file
View file

@ -0,0 +1,249 @@
# API Reference
## Functions
### `initKernels(config?)`
Load the JPL DE442S ephemeris kernel. Must be called before `getMoonSightingReport()` or `getSunMoonEvents()`. Can be called multiple times; subsequent calls replace the loaded kernel.
```ts
async function initKernels(config?: KernelConfig): Promise<void>
```
**KernelConfig:**
```ts
interface KernelConfig {
planetary?: KernelSource // DE442S source. Default: auto-download
leapSeconds?: KernelSource // LSK source. Default: auto-download
cacheDir?: string // Cache dir. Default: ~/.cache/moon-calc
checksumOverride?: string // SHA-256 override for de442s.bsp
}
type KernelSource =
| { type: 'auto' }
| { type: 'file'; path: string }
| { type: 'buffer'; data: ArrayBuffer; name: string }
| { type: 'url'; url: string }
```
---
### `getMoonSightingReport(date, observer, options?)`
Compute a complete moon sighting report. Requires kernel.
```ts
async function getMoonSightingReport(
date: Date,
observer: Observer,
options?: SightingOptions,
): Promise<MoonSightingReport>
```
**Observer:**
```ts
interface Observer {
lat: number // Geodetic latitude, degrees (north positive)
lon: number // Longitude, degrees (east positive)
elevation: number // Height above WGS84 ellipsoid, meters
name?: string // Optional label
deltaT?: number // Override TT - UT1, seconds
ut1utc?: number // Override UT1 - UTC, seconds (takes precedence over deltaT)
pressure?: number // Atmospheric pressure, mbar (default 1013.25)
temperature?: number // Temperature, Celsius (default 15)
}
```
**SightingOptions:**
```ts
interface SightingOptions {
kernels?: KernelConfig
bestTimeMethod?: 'heuristic' | 'optimized' // default: 'heuristic'
}
```
**MoonSightingReport:**
```ts
interface MoonSightingReport {
date: Date
observer: Observer
// Event times
sunsetUTC: Date | null
moonsetUTC: Date | null
lagMinutes: number | null // moonset - sunset, minutes
bestTimeUTC: Date | null // T_sunset + 4/9 × lag
bestTimeWindowUTC: [Date, Date] | null // ±20 min around best time
// Body positions at best time
moonPosition: AzAlt | null // { azimuth, altitude }
sunPosition: AzAlt | null
illumination: number | null // percent, 0–100
moonAge: number | null // hours since conjunction
// Crescent geometry at best time
geometry: CrescentGeometry | null
// Visibility criteria
yallop: YallopResult | null
odeh: OdehResult | null
// Guidance
guidance: string
// Metadata
ephemerisSource: 'DE442S' | 'approximate'
moonAboveHorizon: boolean | null
sightingPossible: boolean
}
```
---
### `getMoonPhase(date?)`
Compute moon phase data. Works without a kernel.
```ts
function getMoonPhase(date?: Date): MoonPhaseResult
```
**MoonPhaseResult:**
```ts
interface MoonPhaseResult {
phase: MoonPhaseName // 'new-moon' | 'waxing-crescent' | ... | 'waning-crescent'
illumination: number // 0–100 percent
age: number // hours since last new moon
elongationDeg: number // Moon - Sun ecliptic longitude, [0, 360)
isWaxing: boolean
nextNewMoon: Date
nextFullMoon: Date
prevNewMoon: Date
}
```
---
### `getSunMoonEvents(date, observer, options?)`
Rise, set, and twilight times. Requires kernel.
```ts
async function getSunMoonEvents(
date: Date,
observer: Observer,
options?: Pick<SightingOptions, 'kernels'>,
): Promise<SunMoonEvents>
```
**SunMoonEvents:**
```ts
interface SunMoonEvents {
sunsetUTC: Date | null
moonsetUTC: Date | null
sunriseUTC: Date | null
moonriseUTC: Date | null
civilTwilightEndUTC: Date | null // Sun at -6°
nauticalTwilightEndUTC: Date | null // Sun at -12°
astronomicalTwilightEndUTC: Date | null // Sun at -18°
}
```
---
### `downloadKernels(config?)`
Download DE442S and naif0012.tls to the local cache (Node.js only).
```ts
async function downloadKernels(config?: KernelConfig): Promise<{
planetaryPath: string
leapSecondsPath: string
}>
```
---
### `verifyKernels(config?)`
Verify locally cached kernels by SHA-256 checksum.
```ts
async function verifyKernels(config?: KernelConfig): Promise<{
ok: boolean
errors: string[]
}>
```
---
## Types
### `CrescentGeometry`
```ts
interface CrescentGeometry {
ARCL: number // Elongation (Sun-Moon angular separation), degrees
ARCV: number // Moon altitude - Sun altitude (airless), degrees
DAZ: number // Sun azimuth - Moon azimuth, [-180, 180], degrees
W: number // Topocentric crescent width, arc minutes
lag: number // Moonset - sunset, minutes
}
```
### `YallopResult`
```ts
interface YallopResult {
q: number // Continuous q parameter
category: YallopCategory // 'A' | 'B' | 'C' | 'D' | 'E' | 'F'
description: string
isVisibleNakedEye: boolean // A or B
requiresOpticalAid: boolean // C or D
isBelowDanjonLimit: boolean // F
Wprime: number // W' used in q formula, arc minutes
}
```
### `OdehResult`
```ts
interface OdehResult {
V: number // Continuous V parameter
zone: OdehZone // 'A' | 'B' | 'C' | 'D'
description: string
isVisibleNakedEye: boolean // A
isVisibleWithOpticalAid: boolean // A or B
}
```
### `AzAlt`
```ts
interface AzAlt {
azimuth: number // Degrees from North, clockwise (0–360)
altitude: number // Degrees above horizon (negative = below)
}
```
---
## Constants
```ts
YALLOP_THRESHOLDS // { A: 0.216, B: -0.014, C: -0.160, D: -0.232, E: -0.293 }
ODEH_THRESHOLDS // { A: 5.65, B: 2.00, C: -0.96 }
WGS84 // { a: 6378137.0, invF: 298.257223563, f, b, e2 }
YALLOP_DESCRIPTIONS // Record<YallopCategory, string>
ODEH_DESCRIPTIONS // Record<OdehZone, string>
```
---
*Previous: [Home](Home) | Next: [Architecture](Architecture)*

139
.wiki/Architecture.md Normal file
View file

@ -0,0 +1,139 @@
# Architecture
## Module layout
```text
src/
math/ Numerical utilities
time/ Time scale conversions
spk/ JPL kernel reader
frames/ Earth orientation transforms
observer/ Observer model and refraction
bodies/ Moon/Sun state computation
events/ Rise/set and event finding
visibility/ Crescent visibility criteria
api/ User-facing API and kernel management
cli/ Command-line interface
```
Each module has zero application-level side effects and no circular dependencies. The layering is strict: lower layers never import from higher ones.
```text
math <── time <── spk <── frames <── observer <── bodies <── events <── visibility <── api
└── cli
```
## Data flow for a sighting report
```text
User calls: getMoonSightingReport(date, observer)
│
├── computeTimeScales(date)
│ UTC → JD(UTC) → ΔAT → JD(TAI) → JD(TT) → JD(TDB) → JD(UT1)
│
├── getSunMoonEvents(date, observer, kernel)
│ For 48 time samples across the civil day:
│ computeAzAlt(Moon, observer, ts)
│ computeAzAlt(Sun, observer, ts)
│ Brent root-finding on altitude crossings
│ → sunsetUTC, moonsetUTC, lag, twilight times
│
├── bestTimeHeuristic(sunsetUTC, moonsetUTC)
│ T_b = T_sunset + (4/9) × Lag
│
├── At best time T_b:
│ computeTimeScales(T_b)
│ getMoonGeocentricState(kernel, ET) ← DE442S SPK evaluation
│ getSunGeocentricState(kernel, ET) ← DE442S SPK evaluation
│ gcrsToItrs(moonGCRS, ts) ← IERS Q·R·W chain
│ geodeticToECEF(observer) ← WGS84
│ topocentricPosition(moon, observer) ← parallax correction
│ enuToAzAlt(moonENU) ← local horizon coords
│
├── computeCrescentGeometry(moonAzAlt, sunAzAlt, moonVec, sunVec)
│ → { ARCL, ARCV, DAZ, W, lag }
│
├── computeYallop(geometry, W') ← q parameter, category A–F
├── computeOdeh(geometry) ← V parameter, zone A–D
│
└── buildGuidanceText(...)
→ "Best time to look: ..."
```
## Ephemeris evaluation
The DE442S kernel is the most important external dependency. It is a JPL SPK file containing Chebyshev polynomial fits to the positions of the Sun, Moon, and planets from 1849 to 2150. Reading it requires:
1. **DAF parser:** reads the binary Double Precision Array File header, iterates summary records, and builds a segment index keyed by (target, center) NAIF ID pairs.
2. **Segment chaining:** DE442S does not store Moon-relative-to-Earth directly. Instead it stores Moon-relative-to-EMB and Earth-relative-to-EMB; the code adds these vectors to get Moon-relative-to-Earth. Similarly for Sun-relative-to-Earth via SSB segments.
3. **Type 2 Chebyshev evaluation:** each time interval has a fixed set of Chebyshev coefficients. Given an ET, the code locates the correct record, computes the normalized time argument x ∈ [−1, 1], and evaluates the degree-n polynomial via the Clenshaw recurrence for each of X, Y, Z.
See [Ephemeris](Ephemeris) for the full technical description.
## Frame transformation
JPL ephemerides are expressed in the ICRF (International Celestial Reference Frame), essentially the J2000 inertial frame. Topocentric alt/az requires Earth-fixed (ITRS) coordinates. The transformation chain (IERS Conventions 2010) is:
```text
[ITRS] = W(t) · R(t) · Q(t) · [GCRS]
```
- **Q(t):** IAU 2006 precession + IAU 2000A nutation (celestial motion), parameterized by CIP coordinates X, Y and CIO locator s.
- **R(t):** Earth rotation angle (ERA) from UT1, a simple rotation about the pole.
- **W(t):** polar motion (xp, yp), typically < 0.5 arcsec; defaults to zero.
This is the largest code in the project (large nutation series tables), but far smaller than the ephemeris data.
See [Reference Frames](Reference-Frames) for implementation details.
## Observer model
Observer position is computed in three stages:
1. **Geodetic → ECEF:** WGS84 ellipsoid, exact formula using N(φ) (prime vertical radius of curvature). Output in meters.
2. **ECEF → GCRS:** Apply the inverse of the Q·R·W transform (since ECEF = ITRS, and GCRS is the inertial frame at that epoch).
3. **Topocentric parallax:** Subtract observer GCRS vector from body GCRS vector to get the topocentric direction.
4. **ENU → az/alt:** Project topocentric vector onto local East-North-Up basis, then compute azimuth and altitude.
Atmospheric refraction (Bennett 1982) is applied as a post-processing step on the computed altitude. For the Yallop/Odeh criteria, airless (refraction-free) altitudes are used. For rise/set times and practical "where to look" output, refracted altitudes are used.
See [Observer Model](Observer-Model) for details.
## Two operating modes
**Full mode** (kernel loaded): all features available, DE442S accuracy.
**Lite mode** (no kernel): `getMoonPhase()` only. Uses Meeus Ch. 25 (Sun) and Ch. 47 (Moon) low-accuracy positions. Error is < 1° in ecliptic longitude. Not suitable for crescent sighting reports.
The API is designed so `getMoonPhase()` never throws for missing kernel; it always works.
## Performance design
- Segment index is built once at kernel load time: O(1) lookup by (target, center, ET).
- Chebyshev records are cached per segment (last-used-record cache). Repeated evaluations in the same time interval cost only the polynomial evaluation, not a binary search.
- `Float64Array` for coefficient storage enables V8/SpiderMonkey typed array optimization paths.
- Clenshaw recurrence avoids the instability of naive power series and is numerically identical to the SPICE SPKE02 implementation.
- The rise/set solver samples at 30-minute intervals before applying Brent, meaning each event finder costs ~48 ephemeris evaluations for bracketing plus ~10 iterations of Brent per root. Total per event: ~60–100 evaluations, each taking tens of microseconds.
Target: a full sighting report (sunset + moonset + best-time geometry + Yallop + Odeh) in 5–15 ms on Node.js.
## Error budget
A crescent sighting report's accuracy is limited by the worst source in the chain:
| Source | Contribution |
| ------ | ------------ |
| DE442S position error | < 1 km (~0.001 arcsec at Moon distance) |
| IERS Q·R·W transform (with user-supplied EOP) | < 1 mas |
| IERS Q·R·W transform (polynomial ΔT approximation) | < 5 arcsec |
| WGS84 observer position | < 1 m (negligible in angle) |
| Bennett refraction (standard atmosphere) | < 1 arcmin for alt > 5° |
| Bennett refraction (non-standard conditions) | up to 15 arcmin near horizon |
In practice, refraction uncertainty dominates all other error sources for crescent sighting near the horizon.
---
*Previous: [API Reference](API-Reference) | Next: [Crescent Visibility](Crescent-Visibility)*

149
.wiki/Crescent-Visibility.md Normal file
View file

@ -0,0 +1,149 @@
# Crescent Visibility Criteria
The new crescent moon is visible when it has moved far enough from the Sun, climbed high enough above the horizon at sunset, and formed a wide enough arc to overcome sky brightness and atmospheric extinction. No criterion can guarantee a sighting. Real-world results depend on atmospheric clarity, observer acuity, and optical equipment.
moon-calc implements two complementary published criteria and outputs both simultaneously, so applications can explain the prediction in terms of either model.
## The five geometric quantities
All major crescent criteria reduce to five observational quantities computed at a canonical "best time" shortly after sunset.
### ARCL: Arc of light (elongation)
The angular separation between the Sun and Moon, measured topocentrically (from the observer, not from Earth's center). This is the true elongation. It determines whether a crescent has begun to form at all: below the Danjon limit (roughly 7°), the Moon is too close to the Sun for the illuminated sliver to survive atmospheric and physiological scattering. ARCL drives the crescent width W.
### ARCV: Arc of vision
Moon altitude minus Sun altitude, both computed without refraction ("airless"). When the Sun is below the horizon, this is approximately the Moon's altitude above the horizon. It controls how dark the sky background is at the moment the crescent is observed. A higher ARCV means a darker sky and a more favorable contrast ratio. This is the primary discriminant in both Yallop and Odeh.
### DAZ: Relative azimuth
Sun azimuth minus Moon azimuth. Positive when the Moon is north of the Sun; negative when south. This affects the geometry of the crescent's orientation and how much sky brightness surrounds it, but plays a secondary role in both criteria.
### W: Topocentric crescent width
The linear width of the illuminated crescent in arc minutes. It serves as a proxy for intrinsic crescent brightness: a thicker crescent is easier to see than a thin one at the same ARCV. Both Yallop and Odeh express the minimum detectable ARCV as a polynomial in W. W is geometrically related to ARCL; as elongation grows, the crescent fattens.
### Lag
The time difference in minutes between moonset and sunset. A positive lag means the Moon sets after the Sun, creating a window for observation. The Yallop/Odeh heuristic computes best time as `T_sunset + 4/9 × Lag`.
## Yallop q-test (NAO TN 69, 1997)
Bernard Yallop's criterion was published in 1997 as Royal Greenwich Observatory Technical Note 69. It defined a standard set of computational procedures and produced the q-parameter, a continuous score that maps to six letter categories.
### Formula
```text
q = (ARCV - arcv_min(W')) / 10
```
Where `W'` is the topocentric crescent width in arc minutes and `arcv_min` is the empirically derived polynomial:
```text
arcv_min(W) = 11.8371 - 6.3226·W + 0.7319·W² - 0.1018·W³
```
This polynomial represents the minimum ARCV observed in historical crescent sightings as a function of crescent width. When `q > 0`, the observed ARCV exceeds the historical minimum, meaning the crescent is potentially visible. The division by 10 is a scaling convention; the thresholds are chosen so that the categories span a meaningful range.
### Categories
| Category | q range | Meaning |
| -------- | ------- | ------- |
| A | q > +0.216 | Easily visible to the naked eye |
| B | q > −0.014 | Visible under perfect conditions |
| C | q > −0.160 | May need optical aid to locate; naked eye possible |
| D | q > −0.232 | Optical aid necessary; naked eye not possible |
| E | q > −0.293 | Not visible even with telescope |
| F | q ≤ −0.293 | Below Danjon limit; crescent cannot form |
Category F corresponds to ARCL below roughly 7° (the Danjon limit), where the Moon is geometrically too close to the Sun for the crescent arc to sustain itself.
### W' vs W
Yallop defines two variants of crescent width. The geocentric W uses the semi-diameter of the Moon and Sun at their geocentric distances. The topocentric W' (W-prime) applies a parallax-based correction because the Moon's apparent diameter changes by ~1% between the geocenter and a surface observer at typical latitudes. moon-calc computes W' directly from the topocentric state vector.
### Implementation notes
- ARCV must be the airless (refraction-free) altitude difference.
- W' must be in arc minutes (the polynomial constants assume this).
- The q thresholds are stated to three decimal places in Yallop's Table 1.
## Odeh criterion (Experimental Astronomy, 2006)
Mohammad Odeh published an updated criterion in 2006 using a larger observational database, including ICOP (Islamic Crescent Observation Project) records collected over multiple decades.
### Odeh formula
```text
V = ARCV - arcv_min(W)
```
The polynomial `arcv_min(W)` is identical to Yallop's formulation. The difference is in the zone boundaries (finer calibration) and in which form of W to use.
```text
V = ARCV - (11.8371 - 6.3226·W + 0.7319·W² - 0.1018·W³)
```
### Zones
| Zone | V range | Meaning |
| ---- | ------- | ------- |
| A | V ≥ 5.65 | Visible with naked eye |
| B | V ≥ 2.00 | Visible with optical aid; may be naked eye under excellent conditions |
| C | V ≥ −0.96 | Visible with optical aid only |
| D | V < −0.96 | Not visible even with optical aid |
### Key differences from Yallop
Odeh provides:
1. **Finer zone boundaries** based on a larger and more systematic observational dataset.
2. **Best-time formula:** the T_b = T_s + 4/9 × Lag expression appears explicitly in Odeh's derivation, though it originates with Yallop.
3. **No F category:** Odeh's Zone D subsumes the Danjon-limit cases.
### Which criterion to trust?
Both criteria agree in the clear cases: large ARCV and wide crescent = visible; very small ARCV or negative W' = not visible. They diverge in the boundary region (Zones B/C, Categories B/C/D), which is inherently uncertain regardless of the model.
In practice, use both: if Yallop says A or B and Odeh says A or B, the crescent is very likely visible. If they disagree near the boundary, treat the result as uncertain and specify weather and optical aid conditions explicitly.
## Best-time computation
### Heuristic (default)
```text
T_b = T_sunset + (4/9) × Lag
```
For a typical 90-minute lag this gives approximately 40 minutes after sunset, which is when the sky has darkened enough for crescent contrast while the Moon is still well above the horizon.
### Optimized (optional)
moon-calc can scan the interval from sunset to moonset, computing the Odeh V parameter at each step and finding the time that maximizes it. This handles high-latitude cases where the Moon's altitude changes quickly and the heuristic can be significantly off.
## Observation window
The library returns a `bestTimeWindowUTC` of ±20 minutes around the best time, giving observers a practical range to watch for the crescent rather than a single instant.
## Limitations
- Refraction uncertainty dominates below 5° altitude. The Bennett formula assumes a standard atmosphere; real conditions can differ by 5–15 arcminutes near the horizon.
- Atmospheric extinction (aerosols, dust) is not modeled. The Schaefer physics-based approach exists in the literature but is not included in this library.
- Urban light pollution reduces effective visibility but is not accounted for by any criterion.
- Observer acuity varies: the criteria are calibrated for average human vision without optical correction.
## Physics-based model (Schaefer, not included)
The Schaefer (1988, 1991) and Doggett-Schaefer contrast threshold model provides a more rigorous treatment of crescent visibility by computing:
- Sky background luminance at the Moon's position
- Atmospheric extinction as a function of aerosol optical depth
- Crescent surface brightness
- Contrast threshold relative to sky background
This approach requires additional atmospheric inputs (aerosol optical depth, humidity proxy) and produces probabilistic visibility estimates rather than deterministic zones. It is not part of this library, but the Yallop and Odeh criteria remain the standard for practical Islamic calendar work.
---
*Previous: [Architecture](Architecture) | Next: [Ephemeris](Ephemeris)*

155
.wiki/Ephemeris.md Normal file
View file

@ -0,0 +1,155 @@
# Ephemeris
## What is DE442S?
DE442 is the current JPL planetary ephemeris, integrated in May 2024. It incorporates Uranus occultation data and additional spacecraft ranging from Mars orbiters and Juno (Jupiter). DE442S is the "short" variant, covering 1849-12-26 to 2150-01-22 (TDB) and weighing 31 MB, practical for download and distribution.
The full DE442 file covers a much longer time span (~114 MB) and is not needed for moon sighting applications. For dates before 1849 or after 2150, you would need the full kernel.
The kernel is distributed by NASA's Navigation and Ancillary Information Facility (NAIF) at:
`https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de442s.bsp`
## SPK file format
SPK (Spacecraft and Planet Kernel) files use the DAF (Double Precision Array File) binary format. DAF is a flexible array-of-records format designed by NAIF for efficient random access.
### File structure
```
├── File record (1024 bytes)
│ Magic: 'NAIF/DAF'
│ ND: number of double-precision summary components (2 for SPK)
│ NI: number of integer summary components (6 for SPK)
│ FWARD: record number of first summary record
│ BWARD: record number of last summary record
│
├── Summary record 1 (1024 bytes)
│ Next/Prev record pointers
│ N: count of summaries in this record
│ [Summary 1: 2 doubles + 6 ints per segment]
│ doubles: start_ET, end_ET
│ ints: target, center, frame, type, begin_addr, end_addr
│
├── Summary record 2 ...
│
└── Data records (variable, Chebyshev coefficients)
```
### Endianness
NAIF kernels are created on the platform that generated them (typically IEEE little-endian for modern kernels). The file record contains an endianness flag; moon-calc reads this and byte-swaps doubles as needed.
### Summary record navigation
Summary records form a doubly-linked list. FWARD and BWARD in the file record give the first and last summary record numbers. Each summary record links to the next via its header. Segment index construction reads all summary records at load time.
## SPK segment types
DE442S uses Type 2 (Chebyshev polynomial position) for all planetary segments.
### Type 2: Chebyshev position-only
Each record covers a fixed time interval and stores coefficients for X, Y, Z:
```
[0] = MID (interval midpoint, ET seconds past J2000)
[1] = RADIUS (interval half-width, seconds)
[2..n+1] = X Chebyshev coefficients C_0..C_n
[n+2..2n+1] = Y Chebyshev coefficients C_0..C_n
[2n+2..3n+1] = Z Chebyshev coefficients C_0..C_n
```
The polynomial degree n is derived from RSIZE (record size in doubles):
```
n = (RSIZE - 2) / 3 - 1
```
Velocity is computed by differentiating the Chebyshev polynomial analytically (using the recurrence relation for Chebyshev derivatives), not by finite differencing.
### Type 3: Chebyshev position + velocity
Type 3 stores separate Chebyshev fits for position and velocity. The structure has 6 coefficient arrays (X, Y, Z for position; X, Y, Z for velocity). Used for some satellite ephemerides; moon-calc implements it for forward compatibility.
## Chebyshev evaluation
The Clenshaw recurrence evaluates a degree-n Chebyshev series in O(n) operations with good numerical stability:
```
Evaluate T_k(x) for x ∈ [-1, 1]:
b_{n+1} = 0
b_{n+2} = 0
for k = n downto 0:
b_k = c_k + 2x·b_{k+1} - b_{k+2}
result = (b_0 - b_2) / 2 [for standard Clenshaw]
```
For moon-calc's implementation (where the constant term c_0 appears differently), we use the SPICE convention:
```
result = c_0 + x·b_1 - b_2
```
This produces the position. Velocity requires the derivative d(result)/dt, computed via the Chebyshev derivative recurrence, not by finite differencing, which would lose accuracy.
Transforming from normalized domain back to physical time:
```
x = (et - MID) / RADIUS
dx/dt = 1/RADIUS
d(result)/dt = d(result)/dx · dx/dt
```
## Segment chaining
DE442S does not provide Moon-relative-to-Earth directly. It provides:
```
Moon (301) relative to Earth-Moon Barycenter (3)
Earth (399) relative to Earth-Moon Barycenter (3)
Earth-Moon Barycenter (3) relative to SSB (0)
Sun (10) relative to SSB (0)
```
The required vectors are assembled by vector addition:
```
Moon relative to Earth = [Moon - EMB] - [Earth - EMB]
Sun relative to Earth = [Sun - SSB] - [EMB - SSB] - [Earth - EMB]
(Earth relative to SSB)
```
The SPK kernel class handles this automatically via a `getState(target, center, et)` interface that recursively chains segments.
## Kernel acquisition
The library supports four acquisition modes:
**auto** (default): Downloads from NAIF on first use, caches locally, verifies by SHA-256. The checksum is bundled in the library and corresponds to the NAIF distribution as of the library release.
**file**: Provide a local path. Useful when the kernel is already downloaded or for offline environments.
**buffer**: Provide an `ArrayBuffer`. Works in browsers and any runtime. The user is responsible for fetching the kernel.
**url**: Provide a custom URL. The library streams the download and verifies the checksum.
## Performance
On first load, the segment index is built in one pass over the summary records. This is fast (~5 ms for DE442S). The data records are not read at load time; they are accessed on demand.
Subsequent calls cache the last-used Chebyshev record per segment, so repeated evaluations in the same time interval cost only the polynomial evaluation (~microseconds).
For batch evaluation (e.g., scanning 1000 times across a year for calendar generation), the segment cache means most evaluations hit the same record and incur no I/O. The total cost scales with the number of distinct time intervals, not the number of evaluations.
## Verification
moon-calc's ephemeris evaluations can be verified against:
1. **SPICE:** the reference implementation from NAIF. Using the same kernel and the same (target, center, frame, ET) arguments, SPICE should produce identical results (to floating-point precision) because both use the same binary data and the same Chebyshev algorithm.
2. **JPL Horizons:** NASA's online ephemeris service, which uses the same JPL ephemerides. Provides tabular output for arbitrary times and locations.
See [Validation](Validation) for the test methodology.
---
*Previous: [Crescent Visibility](Crescent-Visibility) | Next: [Time Scales](Time-Scales)*

183
.wiki/Getting-Started.md Normal file
View file

@ -0,0 +1,183 @@
# Getting Started
## Requirements
- Node.js 20 or later (for full mode with kernel)
- ~35 MB disk space for the kernel cache
- Internet access on first run (for kernel download)
Moon phase queries work in browsers and any runtime without a kernel or network access.
## Install
```bash
npm install moon-calc
# or
pnpm add moon-calc
```
## Download the kernel (one time)
moon-calc uses the JPL DE442S planetary ephemeris. The binary kernel file is 31 MB and is not bundled with the npm package.
```bash
npx moon-calc download-kernels
```
This downloads two files:
- `de442s.bsp` (31 MB): planetary ephemeris, covering 1849–2150
- `naif0012.tls` (4 KB): leap-second table
Default cache location:
- Linux/macOS: `~/.cache/moon-calc/`
- Windows: `%LOCALAPPDATA%\moon-calc\`
To use a custom cache directory:
```ts
await initKernels({ cacheDir: '/my/data/dir' })
```
To verify the download afterward:
```bash
npx moon-calc verify-kernels
```
## First sighting report
```ts
import { initKernels, getMoonSightingReport } from 'moon-calc'
// Load the kernel once per process
await initKernels()
const observer = {
lat: 51.5074, // London
lon: -0.1278,
elevation: 10, // meters above WGS84 ellipsoid
name: 'London, UK',
}
const report = await getMoonSightingReport(new Date('2025-03-29'), observer)
// Summary
console.log(report.yallop.category) // 'A' through 'F'
console.log(report.odeh.zone) // 'A' through 'D'
console.log(report.guidance)
// Event times
console.log(report.sunsetUTC)
console.log(report.moonsetUTC)
console.log(report.lagMinutes)
console.log(report.bestTimeUTC)
// Crescent geometry
console.log(report.geometry)
// { ARCL: 12.3, ARCV: 8.1, DAZ: -2.4, W: 0.21, lag: 67 }
// Where to look
console.log(report.moonPosition)
// { azimuth: 258.3, altitude: 7.9 }
```
## Moon phase (no kernel needed)
```ts
import { getMoonPhase } from 'moon-calc'
const phase = getMoonPhase()
console.log(phase.phase) // 'waxing-crescent'
console.log(phase.illumination) // 23.4
console.log(phase.age) // 4.2 (hours since last new moon)
console.log(phase.nextFullMoon) // Date
// For a specific date
const past = getMoonPhase(new Date('2024-01-01'))
```
## Rise and set times
```ts
import { initKernels, getSunMoonEvents } from 'moon-calc'
await initKernels()
const events = await getSunMoonEvents(new Date('2025-03-29'), {
lat: 21.4225, lon: 39.8262, elevation: 300, name: 'Mecca'
})
console.log(events.sunsetUTC)
console.log(events.moonsetUTC)
console.log(events.civilTwilightEndUTC)
```
## Supplying the kernel manually
### File path (Node.js)
```ts
await initKernels({
planetary: { type: 'file', path: '/data/kernels/de442s.bsp' },
leapSeconds: { type: 'file', path: '/data/kernels/naif0012.tls' },
})
```
### ArrayBuffer (browser)
```ts
const response = await fetch('/kernels/de442s.bsp')
const data = await response.arrayBuffer()
await initKernels({
planetary: { type: 'buffer', data, name: 'de442s.bsp' },
})
```
### URL (streaming load)
```ts
await initKernels({
planetary: { type: 'url', url: 'https://example.com/kernels/de442s.bsp' },
})
```
## Accuracy tips
**Supply the current delta-T from IERS for maximum accuracy.** The built-in polynomial can be off by up to 5 seconds near the present date, introducing a few arcseconds of error in azimuth and altitude.
```ts
// IERS Bulletin A value for UT1-UTC (current, as of 2025-03)
await getMoonSightingReport(date, {
...observer,
ut1utc: 0.0341, // seconds, from IERS Bulletin A
})
```
**The biggest uncertainty in practice is atmospheric refraction**, not the ephemeris. At 5° altitude, refraction uncertainty at non-standard conditions can exceed 10 arcminutes, far larger than any ephemeris error.
## CLI
All features are accessible from the command line:
```bash
# Download kernels
npx moon-calc download-kernels
# Sighting report
npx moon-calc sighting 51.5 -0.1 2025-03-29
# Moon phase
npx moon-calc phase 2025-03-01
# Verify kernel integrity
npx moon-calc verify-kernels
# Performance benchmark
npx moon-calc benchmark
```
---
*Previous: [Home](Home) | Next: [API Reference](API-Reference)*

52
.wiki/Home.md Normal file
View file

@ -0,0 +1,52 @@
# moon-calc
High-accuracy lunar crescent visibility and moon sighting calculations using JPL DE442S ephemerides.
## What this library does
moon-calc computes whether the new crescent moon will be visible at a specific location on a specific date. It produces the geometric quantities used by astronomers and Islamic lunar calendar authorities worldwide, including the Yallop q-test and Odeh visibility zones.
It also computes moon phase data, rise/set times, and twilight periods for any location.
## Key design decisions
**JPL DE442S, not VSOP87 or Meeus.** The library uses the same planetary ephemeris as professional observatories. It reads an SPK binary kernel from NASA's NAIF, covering 1849–2150 with sub-arcsecond accuracy. This makes the position calculations verifiable against SPICE and JPL Horizons.
**Two authoritative visibility criteria.** The Yallop q-test (NAO TN 69, 1997) and Odeh V-parameter (Experimental Astronomy, 2006) are both implemented exactly as published, with all five required geometric variables (ARCL, ARCV, DAZ, W, Lag).
**Kernel not bundled.** The DE442S kernel is 31 MB. It is downloaded once to a local cache, verified by checksum, and reused. This keeps the npm package small.
**Lite mode without kernel.** Moon phase, illumination, and next new/full moon work immediately, no kernel needed. These use Meeus approximations (accurate to ~1°).
## Pages
- [Getting Started](Getting-Started): Installation, kernel setup, first sighting report
- [API Reference](API-Reference): Complete function and type documentation
- [Architecture](Architecture): Module structure, data flow, design rationale
- [Crescent Visibility](Crescent-Visibility): Yallop and Odeh criteria explained in depth
- [Ephemeris](Ephemeris): DE442S kernel format (DAF/SPK), segment chaining, evaluation
- [Time Scales](Time-Scales): UTC, TAI, TT, TDB, UT1, delta-T, leap seconds
- [Reference Frames](Reference-Frames): IERS Q·R·W chain, precession, nutation, ERA
- [Observer Model](Observer-Model): WGS84, topocentric transforms, atmospheric refraction
- [Validation](Validation): Accuracy goals, test methodology, SPICE/Horizons comparison
## Quick example
```ts
import { initKernels, getMoonSightingReport } from 'moon-calc'
await initKernels()
const report = await getMoonSightingReport(new Date('2025-03-29'), {
lat: 51.5074, lon: -0.1278, elevation: 10
})
console.log(report.yallop.category) // 'A'
console.log(report.guidance)
```
## Related packages
- [nrel-spa](https://github.com/acamarata/nrel-spa): Solar position algorithm (zero deps)
- [pray-calc](https://github.com/acamarata/pray-calc): Islamic prayer times
- [luxon-hijri](https://github.com/acamarata/luxon-hijri): Hijri/Gregorian conversion

129
.wiki/Observer-Model.md Normal file
View file

@ -0,0 +1,129 @@
# Observer Model and Refraction
## WGS84 ellipsoid
The library uses the WGS84 (World Geodetic System 1984) reference ellipsoid, which is the standard for GPS coordinates, Google Maps, and most modern mapping systems.
Key constants:
```
a = 6378137.0 m (semi-major axis, equatorial radius)
1/f = 298.257223563 (inverse flattening)
b = a × (1 − f) (semi-minor axis, polar radius ≈ 6356752 m)
e² = 2f − f² (first eccentricity squared ≈ 0.00669438)
```
The shape is an oblate spheroid: slightly wider at the equator than at the poles (~21 km difference).
## Geodetic to ECEF
Given geodetic latitude φ (degrees, north positive), longitude λ (degrees, east positive), and height h (meters above ellipsoid):
```
N(φ) = a / sqrt(1 − e² sin²φ) (prime vertical radius of curvature)
X = (N + h) cos φ cos λ
Y = (N + h) cos φ sin λ
Z = (N(1 − e²) + h) sin φ
```
The result is in meters in the ECEF (Earth-Centered, Earth-Fixed) frame, which is identical to the ITRS (International Terrestrial Reference System) at a given epoch.
Note: geodetic latitude φ is NOT the same as geocentric latitude. For a point on the equator (φ = 0°), both agree. At φ = 45°, the geocentric latitude is about 0.2° less. The distinction matters because the Moon's parallax correction (topocentric vs geocentric position) depends on the observer's actual 3D position.
## ECEF to geodetic
The inverse transformation (ECEF → geodetic) uses Bowring's iterative method, which converges in 2–3 iterations to full double-precision accuracy. This is needed when mapping ECEF coordinates back to lat/lon.
## Local ENU basis
The local East-North-Up (ENU) frame is the observer's natural coordinate system. Its basis vectors in ECEF are:
```
East = (−sin λ, cos λ, 0)
North = (−sin φ cos λ, −sin φ sin λ, cos φ)
Up = (cos φ cos λ, cos φ sin λ, sin φ)
```
These are unit vectors. To convert a topocentric ECEF displacement Δ (in meters) to ENU:
```
e = East · Δ
n = North · Δ
u = Up · Δ
```
## ENU to azimuth/altitude
```
azimuth = atan2(e, n) [radians, convert to degrees, normalize to [0°, 360°)]
altitude = atan2(u, sqrt(e² + n²)) [degrees]
```
Azimuth is measured from North, clockwise: 0° = North, 90° = East, 180° = South, 270° = West.
Altitude is the angle above the horizontal plane: 0° = horizon, 90° = zenith, negative = below horizon.
## Topocentric parallax
The Moon's geocentric position (from the ephemeris) differs from its topocentric position (as seen by a surface observer) because the Moon is close enough that the baseline between Earth's center and the observer's surface position is significant. This is the diurnal parallax.
The correction is simply:
```
topocentric_direction = moon_ITRS − observer_ITRS
```
The Moon's horizontal parallax is approximately 57 arcminutes, meaning the topocentric and geocentric azimuths can differ by up to ~57' (about 1°). For the Sun, the parallax is ~8.7 arcseconds, small but not negligible.
## Atmospheric refraction
The atmosphere bends light from celestial objects upward as they approach the horizon, making them appear higher than their geometric position. The effect is:
- ~34 arcminutes at the geometric horizon (altitude = 0°)
- ~10 arcminutes at altitude 5°
- ~1 arcminute at altitude 20°
- ~0.1 arcminute above 45°
### Bennett (1982) formula
The Bennett formula is the standard practical approximation, accepted by the IAU and widely used in software:
```
R = cot(h + 7.31 / (h + 4.4)) / 60 [degrees]
```
Where h is the geometric (airless) altitude in degrees, and R is the refraction correction to add to h. The constants were derived by fitting to observational data.
With pressure P (millibars) and temperature T (Celsius) corrections:
```
R_actual = R × (P / 1010) × (283 / (273 + T))
```
Standard conditions: P = 1013.25 mbar, T = 15°C. The correction factors adjust for non-standard density.
### Accuracy limits
The Bennett formula is accurate to:
- ~0.1 arcminute for h > 5°
- ~0.5 arcminute for h = 2°–5°
- ~1–2 arcminutes for h < 2°
- Fails below h ≈ −0.5° (below the geometric horizon, refraction becomes strongly non-linear)
In practice, the dominant uncertainty near the horizon is atmospheric variability: temperature inversions, humidity gradients, and dust can shift refraction by 5–15 arcminutes from the standard formula. No formula based solely on pressure and temperature can capture this.
This is why crescent sighting criteria use "airless" (refraction-free) altitudes for ARCV: the criteria were calibrated from historical observations without correcting individual atmospheric conditions, so applying a standard refraction would introduce a systematic bias.
### When to apply refraction
| Use case | Mode |
|----------|------|
| Yallop ARCV input | Airless |
| Odeh ARCV input | Airless |
| Sunset/moonset threshold | Standard refraction |
| "Where to look" altitude output | Standard refraction |
| Civil/nautical/astronomical twilight | Standard refraction |
moon-calc computes both airless and apparent altitudes for each body position and uses the appropriate one for each purpose.
---
*Previous: [Reference Frames](Reference-Frames) | Next: [Validation](Validation)*

108
.wiki/Reference-Frames.md Normal file
View file

@ -0,0 +1,108 @@
# Reference Frames
## The problem
The JPL ephemeris gives Moon and Sun positions in the ICRF (International Celestial Reference Frame), an inertial frame aligned to distant quasars. An observer on Earth's surface needs positions in the local horizon frame (azimuth and altitude). Getting from one to the other requires knowing exactly how Earth is oriented in inertial space at the moment of observation.
## The IERS Q·R·W chain
The IERS Conventions (2010) define the standard transformation:
```
[ITRS] = W(t) · R(t) · Q(t) · [GCRS]
```
Where:
- **GCRS** = Geocentric Celestial Reference System (essentially the inertial J2000 frame at Earth's center)
- **ITRS** = International Terrestrial Reference System (Earth-fixed frame, rotates with the solid Earth)
- **Q(t)** = celestial motion matrix (precession + nutation)
- **R(t)** = Earth rotation matrix
- **W(t)** = polar motion matrix
The matrices are applied right-to-left: first Q (inertial → intermediate frame), then R (apply Earth's rotation), then W (correct for pole wobble).
## Q(t): Celestial motion
Q(t) captures the slow drift of Earth's rotation axis due to gravitational torques from the Moon and Sun (precession, ~26,000-year period) and higher-frequency oscillations (nutation, dominant period ~18.6 years).
The IAU 2006 precession model and IAU 2000A nutation model together parameterize Q(t) via three quantities:
- **X, Y:** celestial intermediate pole (CIP) coordinates in radians
- **s:** CIO locator, a small angle that ensures continuity of the CIO position
The CIP X,Y series has:
- A polynomial part (degree 5 in T = Julian centuries from J2000.0)
- 1,306 luni-solar nutation terms
- 687 planetary nutation terms
This is the largest tabular computation in the library. The IAU 2000B model reduces to 77 terms with < 1 milliarcsecond error, an option for lower-accuracy "lite" builds.
The Q matrix from X, Y, s (from IERS Conventions eq. 5.7):
```
Q = Q₁ · Rz(s)
```
Where Q₁ is built from X and Y directly using the exact (non-small-angle) formula, and Rz(s) rotates by the small CIO locator angle s (< 1 arcsec throughout the century).
## R(t): Earth rotation
R(t) is a simple rotation about the CIP (z-axis direction) by the Earth Rotation Angle:
```
ERA(UT1) = 2π × (0.7790572732640 + 1.00273781191135448 × Du)
Du = JD(UT1) − 2451545.0
```
ERA replaced Greenwich Mean Sidereal Time (GMST) in the IAU 2000 model. It is defined directly from UT1 (Earth's rotation) rather than being derived from precession and nutation models, which is conceptually cleaner.
The rate multiplier 1.00273781191135448 (slightly more than 1 revolution per solar day) accounts for Earth's orbital motion around the Sun.
## W(t): Polar motion
The Earth's rotation pole (CIP) does not coincide exactly with the conventional terrestrial pole. The wobble (Chandler wobble + annual wobble + longer-period terms) is described by two angles:
- **xp:** pole x-offset in arcseconds (toward 0° longitude)
- **yp:** pole y-offset in arcseconds (toward 270° longitude)
Typical magnitudes: ≤ 0.3 arcseconds. At the Earth's surface, this introduces errors of < 30 meters if ignored, negligible for crescent sighting work (Moon angular diameter ~0.5°).
moon-calc defaults to xp = yp = 0. Supply current values from IERS Bulletin A for maximum accuracy.
## IAU 2000A vs 2000B
| Feature | 2000A | 2000B |
|---------|-------|-------|
| Luni-solar terms | 1,306 | 77 |
| Planetary terms | 687 | 0 |
| Max error | < 0.1 mas | < 1 mas |
| Computation | ~2× slower | fast |
| Suitable for | moon sighting | approximate work |
For crescent sighting at the horizon where refraction dominates, 2000B is more than sufficient. moon-calc defaults to 2000A for correctness; 2000B will be available as a compile-time option for size-sensitive builds.
## From GCRS to local alt/az
After applying Q·R·W:
1. **GCRS → ITRS**: Apply Q, R, W to get Earth-fixed Cartesian coordinates.
2. **Observer ECEF position**: Computed from WGS84 geodetic coordinates.
3. **Topocentric vector**: Subtract observer ECEF from body ECEF (in ITRS). The result is the geocentric vector reduced to the observer's position, the topocentric parallax correction.
4. **ECEF → ENU**: Project the topocentric vector onto the observer's local East-North-Up basis.
5. **ENU → az/alt**: Azimuth = atan2(east, north); altitude = atan2(up, horizontal_distance).
See [Observer Model](Observer-Model) for the WGS84 and ENU computation details.
## Accuracy
With user-supplied EOP (Earth orientation parameters from IERS Bulletin A):
- Azimuth/altitude accuracy: < 0.1 arcsecond (dominated by nutation model error)
With polynomial ΔT approximation (no user EOP):
- Azimuth/altitude accuracy: typically < 5 arcseconds, occasionally up to ~30 arcseconds in pathological ΔT errors
For comparison, the Moon's angular diameter is ~1800 arcseconds, and refraction uncertainty near the horizon is 600–900 arcseconds. The frame transform is not the limiting factor for crescent sighting.
---
*Previous: [Time Scales](Time-Scales) | Next: [Observer Model](Observer-Model)*

109
.wiki/Time-Scales.md Normal file
View file

@ -0,0 +1,109 @@
# Time Scales
Getting time right is the most error-prone part of an astronomy library. moon-calc implements a complete conversion chain from the user's familiar UTC input to the internal TDB time argument required by the JPL ephemeris.
## The time scale chain
```text
User input (Date / ISO string)
│
▼
UTC — Coordinated Universal Time
│ + ΔAT (leap seconds, integer steps)
▼
TAI — International Atomic Time
│ + 32.184 s (exact, by definition)
▼
TT — Terrestrial Time (formerly ET)
│ + periodic correction (~1.7 ms max)
▼
TDB — Barycentric Dynamical Time ← used by JPL DE442S
Also needed:
UTC + (UT1 - UTC) ← from IERS Bulletin A
▼
UT1 — Universal Time 1 ← used for Earth Rotation Angle
```
## UTC: Coordinated Universal Time
UTC is the international time standard, synchronized to TAI via occasional leap second insertions. JavaScript's `Date` object stores Unix time (seconds since 1970-01-01 00:00:00 UTC), which is equivalent to UTC ignoring leap seconds.
## TAI: International Atomic Time
TAI = UTC + ΔAT, where ΔAT is the cumulative leap second count. ΔAT grows in 1-second steps whenever the IERS determines that Earth's rotation has slowed enough. As of 2024, ΔAT = 37 seconds.
moon-calc ships the complete leap-second table (through 2017) and parses the NAIF LSK kernel (`naif0012.tls`) to stay current when the user downloads it.
## TT: Terrestrial Time
TT = TAI + 32.184 seconds exactly. This exact offset is the definition. TT replaced the old Ephemeris Time (ET) designation in IAU 1991; for practical purposes TT and the old ET are identical.
TT is the primary argument for Earth-based positional calculations: precession-nutation series, mean orbital elements, and topocentric corrections all use TT (expressed as Julian centuries from J2000.0).
```text
T = (JD_TT − 2451545.0) / 36525.0 (Julian centuries from J2000.0)
```
## TDB: Barycentric Dynamical Time
TDB is the time coordinate for the JPL ephemeris. It differs from TT by a periodic term (relativistic clock correction for the eccentricity of Earth's orbit) that never exceeds 1.7 milliseconds:
```text
TDB − TT ≈ 0.001658 s × sin(g) + 0.000014 s × sin(2g)
g = 357.53° + 0.9856003° × (JD_TT − 2451545.0) (Sun's mean anomaly)
```
For crescent sighting purposes this ~1 ms difference is negligible (Moon moves ~30 arcsec/second, so 1 ms → ~0.03 arcsec position error). moon-calc applies the correction anyway to be consistent with SPICE.
In the SPICE system, the internal time argument is called ET (Ephemeris Time) and is expressed as seconds past J2000.0 TDB:
```text
ET = (JD_TDB − 2451545.0) × 86400.0
```
## UT1: Universal Time 1
UT1 is the measure of Earth's actual rotation angle. It tracks mean solar time at the Greenwich meridian and is needed to compute the Earth Rotation Angle (ERA), which determines the relationship between the inertial celestial frame and the Earth-fixed terrestrial frame.
UT1 differs from UTC by at most ±0.9 seconds (the IERS inserts leap seconds before this threshold is reached). The difference UT1 − UTC is published in IERS Bulletin A, updated weekly.
For applications requiring maximum accuracy, supply the current `ut1utc` value from IERS Bulletin A to `Observer.ut1utc`. This eliminates the primary source of error in the Earth orientation model.
## ΔT: the historical problem
ΔT = TT − UT1 combines the accumulated leap seconds plus the sub-second UT1 − UTC offset. It grows roughly parabolically over historical timescales as Earth's rotation gradually slows.
In 2024, ΔT ≈ 69.2 seconds. The built-in polynomial (Espenak-Meeus) approximates ΔT without IERS data:
- Error within ±5 seconds for 1950–2050
- Error grows to ±30 seconds at 1800 and ±100 seconds at 1600
- Forward extrapolation beyond 2050 is unreliable (ΔT growth rate varies)
For historical crescent analysis (e.g., verifying ancient sightings), the ΔT uncertainty is often the dominant error source, not the ephemeris.
## Julian Day
All internal time arithmetic uses Julian Day (JD), a continuous count of days since noon on January 1, 4713 BC (Julian calendar). This avoids calendar ambiguity and makes time differences trivial.
J2000.0 epoch:
```text
JD = 2451545.0 (2000 Jan 1, 12:00:00 TT)
```
Converting JavaScript Date to JD:
```text
JD(UTC) = Date.getTime() / 86400000 + 2440587.5
```
The constant 2440587.5 is JD for 1970-01-01 00:00:00 UTC.
## Leap-second kernel
The NAIF LSK (`naif0012.tls`) is a plain-text file in NAIF text kernel format. It contains the ΔAT table and the constants needed to compute TDB from TT. moon-calc parses this file when the user downloads kernels, ensuring the library always reflects the latest leap seconds. The hardcoded table in `time/index.ts` serves as a fallback when no LSK is provided.
---
*Previous: [Ephemeris](Ephemeris) | Next: [Reference Frames](Reference-Frames)*

101
.wiki/Validation.md Normal file
View file

@ -0,0 +1,101 @@
# Validation
moon-calc is validated at multiple levels: kernel parsing, ephemeris evaluation, frame transforms, and full sighting reports. The goal is to make each component independently verifiable against authoritative references.
## Validation philosophy
Crescent visibility criteria depend on a chain of computations, each of which can introduce errors. Validating the end result (Yallop category) alone is insufficient. A wrong category can result from an error in any step. The validation strategy separates:
1. **Ephemeris evaluation:** does the kernel parsing produce the right position?
2. **Frame transforms:** does the IERS chain produce the right ITRS position?
3. **Topocentric position:** does the observer model produce the right az/alt?
4. **Crescent geometry:** do ARCL, ARCV, DAZ, W match reference implementations?
5. **Visibility criteria:** do the q and V parameters match published examples?
## Reference implementations
### SPICE (CSPICE / SpiceyPy)
NASA NAIF's SPICE toolkit is the authoritative reference for reading JPL ephemerides. It is written in C (CSPICE) with Python bindings (SpiceyPy). Using the same kernel (`de442s.bsp`) and the same time/frame arguments, SPICE should produce positions that are bit-identical to moon-calc's output (to double-precision floating-point).
Any deviation in the SPK Chebyshev evaluation from SPICE indicates a parsing or algorithm error in moon-calc.
**How to compare:**
```python
import spiceypy as spice
spice.furnsh('de442s.bsp')
spice.furnsh('naif0012.tls')
# Moon relative to Earth center in J2000 at ET
et = spice.str2et('2025-03-29 20:00:00 UTC')
state, lt = spice.spkezr('301', et, 'J2000', 'NONE', '399')
print(state[:3]) # position in km
```
The moon-calc equivalent:
```ts
const kernel = SpkKernel.fromFile('de442s.bsp')
const ts = computeTimeScales(new Date('2025-03-29T20:00:00Z'))
const state = kernel.getState(NAIF_IDS.MOON, NAIF_IDS.EARTH, jdTTtoET(ts.jdTT))
```
Expected agreement: < 1 meter (floating-point evaluation precision).
### JPL Horizons
JPL Horizons is the online solar system ephemeris service. It uses the same JPL ephemerides and provides tabular output for:
- Apparent RA/Dec and az/alt for any observer and time
- Observer-centered quantities (elongation, illumination, phase angle)
- Rise/transit/set times
Horizons uses SPICE internally, so it represents an independent end-to-end validation of the full pipeline including frame transforms and refraction.
**How to use for validation:**
Go to https://ssd.jpl.nasa.gov/horizons/, select:
- Target body: Moon (or Sun)
- Observer location: user-defined geodetic lat/lon/elevation
- Time span: the date of interest
- Output quantities: Observer az/alt, Illuminated fraction, Elongation
Compare Horizons' output with moon-calc's topocentric az/alt. Differences of < 30 arcseconds indicate the frame transforms are correct.
## Acceptance thresholds
| Quantity | Expected error vs SPICE | Notes |
|----------|------------------------|-------|
| Geocentric position | < 1 m (< 0.001 arcsec) | SPK parsing precision |
| Topocentric az/alt (with EOP) | < 0.1 arcsec | Frame transform precision |
| Topocentric az/alt (polynomial ΔT) | < 30 arcsec | ΔT polynomial error |
| ARCL | < 1 arcsec | Derived from positions |
| ARCV | < 30 arcsec | Dominated by ΔT uncertainty |
| Yallop q | < 0.005 | q is dimensionless; <0.005 difference = same category in most cases |
| Sunset/moonset | < 10 seconds | Root-finding convergence |
## Validation suite
The test harness will:
1. Load `de442s.bsp` and evaluate Moon/Sun states at 1000 randomly sampled times within 1849–2150.
2. Compare each result against SPICE output for the same kernel, time, and frames.
3. For 50 geographically diverse locations × 12 months, compute full sighting reports and compare ARCL, ARCV, DAZ, W, q, and V against reference implementations.
4. For the same 50 locations, compare Horizons tabular output with moon-calc's az/alt at the same times.
The suite is run in CI on major releases. Raw comparison data (SPICE reference outputs) is stored as CSV files in the test fixtures directory.
## Spot-check cases
Historical crescent sightings from the ICOP database can serve as sanity checks. A positive sighting report from ICOP should correspond to Yallop A or B, and Odeh A or B, at the recording location and date. A failed sighting should correspond to Yallop C–F / Odeh C–D.
Be cautious: ICOP records include weather and observer acuity information that the criteria cannot account for. A clear Zone A (easily visible) sighting that was missed due to clouds or poor horizon is not a failure of the criterion.
## Known limitations
**Near-horizon refraction.** The Bennett formula assumes standard atmospheric conditions. At altitudes below 3°, real refraction can vary by ±20% from the standard value. The validation harness cannot test this without actual atmospheric measurements.
**Historical delta-T.** For dates before 1800, the Espenak-Meeus polynomial error can exceed 1 minute. Sighting analysis for ancient dates requires user-supplied ΔT from Morrison & Stephenson or similar.
**Polar and extreme latitudes.** In Arctic regions, the Sun can be below the horizon for days or weeks. The event finder handles circumpolar conditions (no rise or set during the search window) by returning null for the missing events.
---
*Previous: [Observer Model](Observer-Model) | Next: [API Reference](API-Reference)*

35
CHANGELOG.md Normal file
View file

@ -0,0 +1,35 @@
# Changelog
All notable changes to moon-calc are documented here.
Format follows [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
## [1.0.0] - 2026-02-25
### Added
- Core type system: `Observer`, `MoonSightingReport`, `CrescentGeometry`, `YallopResult`, `OdehResult`, `MoonPhaseResult`, `SunMoonEvents`, `KernelConfig`, `SightingOptions`
- Module architecture: `math/`, `time/`, `spk/`, `frames/`, `observer/`, `bodies/`, `events/`, `visibility/`, `api/`, `cli/`
- DAF/SPK Type 2 parser and Chebyshev evaluator for DE442S
- Kernel auto-download and SHA-256 checksum verification (`initKernels`, `downloadKernels`, `verifyKernels`)
- Full time scale chain: UTC → TAI → TT → TDB with delta-T polynomial, leap-second table, ERA computation
- IERS Q·R·W frame transforms: IAU 2006/2000A precession, nutation, polar motion
- WGS84 geodetic ↔ ECEF conversion, topocentric ENU projection
- Bennett (1982) atmospheric refraction formula with pressure/temperature correction
- Topocentric Moon/Sun state computation via DE442S segment chaining
- Meeus approximate positions for kernel-free `getMoonPhase()`
- Rise/set event solver using Brent's method over the altitude function
- Twilight computation: civil (−6°), nautical (−12°), astronomical (−18°)
- Full crescent geometry: ARCL, ARCV, DAZ, W (arc minutes), lag
- Yallop q-test constants and category thresholds (NAO TN 69)
- Odeh zone constants, V-parameter thresholds, and best-time optimizer (Experimental Astronomy 2006)
- WGS84 ellipsoid constants and Clenshaw Chebyshev evaluation
- Vector/matrix math utilities, Brent root-finding
- Best-time heuristic: T_b = T_sunset + (4/9) × Lag
- Observation window computation (±20 min around best time)
- Odeh-based and Yallop-based guidance text generation
- `getMoonSightingReport()` full pipeline
- `getMoonPhase()` (Meeus approximation, kernel-free)
- `getSunMoonEvents()` full pipeline
- CLI commands: `download-kernels`, `verify-kernels`, `sighting`, `phase`, `benchmark`
- CI: Node matrix (20/22/24), typecheck, pack-check
- Wiki: Architecture, API Reference, Crescent Visibility, Ephemeris, Time Scales, Reference Frames, Observer Model, Validation, Getting Started

21
LICENSE Normal file
View file

@ -0,0 +1,21 @@
MIT License
Copyright (c) 2026 Aric Camarata
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

264
README.md Normal file
View file

@ -0,0 +1,264 @@
# moon-calc
[![npm version](https://img.shields.io/npm/v/moon-calc.svg)](https://www.npmjs.com/package/moon-calc)
[![CI](https://github.com/acamarata/moon-calc/actions/workflows/ci.yml/badge.svg)](https://github.com/acamarata/moon-calc/actions/workflows/ci.yml)
[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](LICENSE)
High-accuracy lunar crescent visibility and moon sighting library for Node.js and browsers. Uses JPL DE442S ephemerides with full IERS Earth orientation for sub-arcsecond topocentric Moon and Sun positions.
Implements the Yallop (NAO TN 69) and Odeh (Experimental Astronomy 2006) crescent visibility criteria, the two most widely used models in Islamic crescent sighting workflows.
## Installation
```bash
npm install moon-calc
```
After installing, download the JPL ephemeris kernel (31 MB, one-time setup):
```bash
npx moon-calc download-kernels
```
This fetches `de442s.bsp` and `naif0012.tls` from NASA's NAIF server and caches them locally. The download is verified by SHA-256 checksum.
## Quick start
```ts
import { initKernels, getMoonSightingReport, getMoonPhase } from 'moon-calc'
// One-time setup: load the ephemeris kernel
await initKernels()
// Full sighting report for a date and location
const report = await getMoonSightingReport(new Date('2025-03-29'), {
lat: 51.5074,
lon: -0.1278,
elevation: 10,
name: 'London, UK'
})
console.log(report.yallop.category) // 'A' (easily visible to the naked eye)
console.log(report.odeh.zone) // 'A' (visible with naked eye)
console.log(report.guidance)
// "Best time to look: 2025-03-29 20:14 UTC (73 min after sunset).
// Look West at 8° above the horizon. The crescent should be visible
// to the naked eye. Yallop: A (Easily visible to the naked eye).
// Odeh: A (Visible with naked eye)."
// Moon phase works without a kernel
const phase = getMoonPhase()
console.log(phase.phase) // 'waxing-crescent'
console.log(phase.illumination) // 14.3
console.log(phase.nextFullMoon) // Date
```
## API
### `initKernels(config?)`
Load the ephemeris kernel. Required before `getMoonSightingReport()`.
```ts
// Auto-download and cache (default)
await initKernels()
// User-supplied file path (Node.js)
await initKernels({ planetary: { type: 'file', path: '/data/de442s.bsp' } })
// ArrayBuffer (browser)
await initKernels({ planetary: { type: 'buffer', data: buf, name: 'de442s.bsp' } })
```
### `getMoonSightingReport(date, observer, options?)`
Returns a complete moon sighting report.
| Parameter | Type | Description |
| --------- | ---- | ----------- |
| `date` | `Date` | Civil date to check (UTC) |
| `observer.lat` | `number` | Geodetic latitude, degrees (north positive) |
| `observer.lon` | `number` | Longitude, degrees (east positive) |
| `observer.elevation` | `number` | Height above WGS84 ellipsoid, meters |
| `observer.deltaT` | `number?` | Override TT - UT1 in seconds (IERS value) |
| `observer.pressure` | `number?` | Atmospheric pressure, mbar (default 1013.25) |
| `observer.temperature` | `number?` | Temperature, Celsius (default 15) |
**Returns** `MoonSightingReport`:
| Field | Type | Description |
| ----- | ---- | ----------- |
| `sunsetUTC` | `Date` | Sunset time |
| `moonsetUTC` | `Date` | Moonset time |
| `lagMinutes` | `number` | Moonset - sunset, minutes |
| `bestTimeUTC` | `Date` | Optimal observation time |
| `geometry.ARCL` | `number` | Arc of light (elongation), degrees |
| `geometry.ARCV` | `number` | Arc of vision (airless), degrees |
| `geometry.DAZ` | `number` | Relative azimuth, degrees |
| `geometry.W` | `number` | Crescent width, arc minutes |
| `yallop.category` | `'A'`–`'F'` | Yallop visibility class |
| `yallop.q` | `number` | Continuous q parameter |
| `odeh.zone` | `'A'`–`'D'` | Odeh visibility zone |
| `odeh.V` | `number` | Continuous V parameter |
| `moonPosition` | `AzAlt` | Moon azimuth/altitude at best time |
| `guidance` | `string` | Plain-language sighting instructions |
### `getMoonPhase(date?)`
Compute the Moon's current phase. Works without a kernel.
| Field | Type | Description |
| ----- | ---- | ----------- |
| `phase` | `string` | Phase name (e.g., `'waxing-crescent'`) |
| `illumination` | `number` | Illuminated fraction, 0–100 |
| `age` | `number` | Hours since last new moon |
| `isWaxing` | `boolean` | True when illumination is increasing |
| `prevNewMoon` | `Date` | Time of previous new moon |
| `nextNewMoon` | `Date` | Time of next new moon |
| `nextFullMoon` | `Date` | Time of next full moon |
### `getSunMoonEvents(date, observer)`
Get rise, set, and twilight times. Requires kernel.
| Field | Description |
| ----- | ----------- |
| `sunsetUTC` | Sunset |
| `moonsetUTC` | Moonset |
| `sunriseUTC` | Sunrise |
| `moonriseUTC` | Moonrise |
| `civilTwilightEndUTC` | Civil twilight (Sun at -6°) |
| `nauticalTwilightEndUTC` | Nautical twilight (Sun at -12°) |
| `astronomicalTwilightEndUTC` | Astronomical twilight (Sun at -18°) |
### `downloadKernels(config?)`
Download DE442S and naif0012.tls to the local cache (Node.js).
### `verifyKernels(config?)`
Verify cached kernels by SHA-256 checksum.
## Visibility criteria
### Yallop (A–F)
| Category | q range | Interpretation |
| -------- | ------- | -------------- |
| A | q > +0.216 | Easily visible to the naked eye |
| B | q > -0.014 | Visible under perfect conditions |
| C | q > -0.160 | May need optical aid to find |
| D | q > -0.232 | Needs optical aid; not visible naked eye |
| E | q > -0.293 | Not visible even with telescope |
| F | q ≤ -0.293 | Below Danjon limit |
### Odeh (A–D)
| Zone | V range | Interpretation |
| ---- | ------- | -------------- |
| A | V ≥ 5.65 | Visible with naked eye |
| B | V ≥ 2.00 | Visible with optical aid; may be naked eye |
| C | V ≥ -0.96 | Visible with optical aid only |
| D | V < -0.96 | Not visible even with optical aid |
## Architecture
```text
src/
math/ Vector/matrix, Chebyshev evaluation, root-finding
time/ UTC → TAI → TT → TDB conversions, leap seconds, Julian Day
spk/ JPL DAF/SPK kernel parser, Chebyshev segment evaluator
frames/ IERS Q·R·W chain: precession + nutation + ERA + polar motion
observer/ WGS84 geodetic → ECEF, topocentric ENU, Bennett refraction
bodies/ Moon/Sun state computation, illumination, crescent width
events/ Rise/set solver, twilight, best-time computation
visibility/ Yallop q-test, Odeh zones, crescent geometry
api/ User-facing functions, kernel management
cli/ Command-line interface
```
See the [Architecture wiki page](https://github.com/acamarata/moon-calc/wiki/Architecture) for a full technical description.
## CLI
```bash
# Setup (one-time)
npx moon-calc download-kernels
# Sighting report
npx moon-calc sighting 51.5 -0.1 2025-03-29
npx moon-calc sighting 21.4 39.8 # Mecca
# Moon phase
npx moon-calc phase 2025-03-01
# Verify downloaded kernels
npx moon-calc verify-kernels
# Benchmark
npx moon-calc benchmark
```
## Compatibility
| Environment | Support |
| ----------- | ------- |
| Node.js 20+ | Full (all features) |
| Node.js 22, 24 | Full |
| Browser | Partial (no auto-download; supply kernel buffer) |
| ESM | `import` from `moon-calc` |
| CommonJS | `require('moon-calc')` |
| TypeScript | Full type definitions included |
## TypeScript
```ts
import type {
Observer,
MoonSightingReport,
YallopCategory,
OdehZone,
KernelConfig,
} from 'moon-calc'
```
## Documentation
Full documentation is on the [GitHub Wiki](https://github.com/acamarata/moon-calc/wiki):
- [Getting Started](https://github.com/acamarata/moon-calc/wiki/Getting-Started)
- [API Reference](https://github.com/acamarata/moon-calc/wiki/API-Reference)
- [Architecture](https://github.com/acamarata/moon-calc/wiki/Architecture)
- [Crescent Visibility Criteria](https://github.com/acamarata/moon-calc/wiki/Crescent-Visibility)
- [Ephemeris and Kernel Setup](https://github.com/acamarata/moon-calc/wiki/Ephemeris)
- [Time Scales](https://github.com/acamarata/moon-calc/wiki/Time-Scales)
- [Reference Frames](https://github.com/acamarata/moon-calc/wiki/Reference-Frames)
- [Observer Model and Refraction](https://github.com/acamarata/moon-calc/wiki/Observer-Model)
- [Validation](https://github.com/acamarata/moon-calc/wiki/Validation)
## Related
- [nrel-spa](https://github.com/acamarata/nrel-spa): Pure JS solar position algorithm (zero deps)
- [pray-calc](https://github.com/acamarata/pray-calc): Islamic prayer times with dynamic angle algorithm
- [luxon-hijri](https://github.com/acamarata/luxon-hijri): Hijri/Gregorian calendar conversion
## Acknowledgments
Crescent visibility criteria implemented from:
- B.D. Yallop, "A Method for Predicting the First Sighting of the New Crescent Moon," NAO Technical Note No. 69, Royal Greenwich Observatory, 1997.
- M.Sh. Odeh, "New Criterion for Lunar Crescent Visibility," Experimental Astronomy 18(1), 39–64, 2006.
Planetary ephemeris data from:
- JPL DE442S. Jet Propulsion Laboratory, NASA. Ryan S. Park et al. (2021). "The JPL Planetary and Lunar Ephemerides DE440 and DE441." Astronomical Journal 161(3), 105. [doi:10.3847/1538-3881/abd414](https://doi.org/10.3847/1538-3881/abd414)
NAIF SPICE toolkit concepts:
- Navigation and Ancillary Information Facility (NAIF), Jet Propulsion Laboratory.
## License
MIT. See [LICENSE](LICENSE).
The DE442S kernel data is provided by NASA/JPL and is not redistributed with this package. It is downloaded separately from the NAIF public server.

74
package.json Normal file
View file

@ -0,0 +1,74 @@
{
"name": "moon-calc",
"version": "1.0.0",
"description": "High-accuracy lunar crescent visibility and moon sighting calculations using JPL DE442S ephemerides. Implements Yallop and Odeh criteria for Islamic crescent sighting workflows.",
"author": "Aric Camarata",
"license": "MIT",
"main": "./dist/index.cjs",
"module": "./dist/index.mjs",
"types": "./dist/index.d.ts",
"exports": {
".": {
"types": "./dist/index.d.ts",
"import": "./dist/index.mjs",
"require": "./dist/index.cjs"
}
},
"bin": {
"moon-calc": "./dist/cli/index.cjs"
},
"sideEffects": false,
"files": [
"dist",
"README.md",
"CHANGELOG.md",
"LICENSE"
],
"engines": {
"node": ">=20"
},
"scripts": {
"build": "tsup",
"typecheck": "tsc --noEmit",
"pretest": "tsup",
"test": "node test.mjs && node test-cjs.cjs",
"prepublishOnly": "tsup",
"cli": "node dist/cli/index.cjs"
},
"devDependencies": {
"@types/node": "latest",
"tsup": "latest",
"typescript": "latest"
},
"publishConfig": {
"access": "public",
"registry": "https://registry.npmjs.org/"
},
"repository": {
"type": "git",
"url": "git+https://github.com/acamarata/moon-calc.git"
},
"homepage": "https://github.com/acamarata/moon-calc#readme",
"bugs": {
"url": "https://github.com/acamarata/moon-calc/issues"
},
"keywords": [
"moon",
"lunar",
"crescent",
"moon-sighting",
"hilal",
"yallop",
"odeh",
"islamic",
"ephemeris",
"jpl",
"de442s",
"astronomy",
"moon-phase",
"visibility",
"sighting",
"ramadan",
"prayer-times"
]
}

925
pnpm-lock.yaml Normal file
View file

@ -0,0 +1,925 @@
lockfileVersion: '9.0'
settings:
autoInstallPeers: true
excludeLinksFromLockfile: false
importers:
.:
devDependencies:
'@types/node':
specifier: latest
version: 25.3.0
tsup:
specifier: latest
version: 8.5.1(typescript@5.9.3)
typescript:
specifier: latest
version: 5.9.3
packages:
'@esbuild/aix-ppc64@0.27.3':
resolution: {integrity: sha512-9fJMTNFTWZMh5qwrBItuziu834eOCUcEqymSH7pY+zoMVEZg3gcPuBNxH1EvfVYe9h0x/Ptw8KBzv7qxb7l8dg==}
engines: {node: '>=18'}
cpu: [ppc64]
os: [aix]
'@esbuild/android-arm64@0.27.3':
resolution: {integrity: sha512-YdghPYUmj/FX2SYKJ0OZxf+iaKgMsKHVPF1MAq/P8WirnSpCStzKJFjOjzsW0QQ7oIAiccHdcqjbHmJxRb/dmg==}
engines: {node: '>=18'}
cpu: [arm64]
os: [android]
'@esbuild/android-arm@0.27.3':
resolution: {integrity: sha512-i5D1hPY7GIQmXlXhs2w8AWHhenb00+GxjxRncS2ZM7YNVGNfaMxgzSGuO8o8SJzRc/oZwU2bcScvVERk03QhzA==}
engines: {node: '>=18'}
cpu: [arm]
os: [android]
'@esbuild/android-x64@0.27.3':
resolution: {integrity: sha512-IN/0BNTkHtk8lkOM8JWAYFg4ORxBkZQf9zXiEOfERX/CzxW3Vg1ewAhU7QSWQpVIzTW+b8Xy+lGzdYXV6UZObQ==}
engines: {node: '>=18'}
cpu: [x64]
os: [android]
'@esbuild/darwin-arm64@0.27.3':
resolution: {integrity: sha512-Re491k7ByTVRy0t3EKWajdLIr0gz2kKKfzafkth4Q8A5n1xTHrkqZgLLjFEHVD+AXdUGgQMq+Godfq45mGpCKg==}
engines: {node: '>=18'}
cpu: [arm64]
os: [darwin]
'@esbuild/darwin-x64@0.27.3':
resolution: {integrity: sha512-vHk/hA7/1AckjGzRqi6wbo+jaShzRowYip6rt6q7VYEDX4LEy1pZfDpdxCBnGtl+A5zq8iXDcyuxwtv3hNtHFg==}
engines: {node: '>=18'}
cpu: [x64]
os: [darwin]
'@esbuild/freebsd-arm64@0.27.3':
resolution: {integrity: sha512-ipTYM2fjt3kQAYOvo6vcxJx3nBYAzPjgTCk7QEgZG8AUO3ydUhvelmhrbOheMnGOlaSFUoHXB6un+A7q4ygY9w==}
engines: {node: '>=18'}
cpu: [arm64]
os: [freebsd]
'@esbuild/freebsd-x64@0.27.3':
resolution: {integrity: sha512-dDk0X87T7mI6U3K9VjWtHOXqwAMJBNN2r7bejDsc+j03SEjtD9HrOl8gVFByeM0aJksoUuUVU9TBaZa2rgj0oA==}
engines: {node: '>=18'}
cpu: [x64]
os: [freebsd]
'@esbuild/linux-arm64@0.27.3':
resolution: {integrity: sha512-sZOuFz/xWnZ4KH3YfFrKCf1WyPZHakVzTiqji3WDc0BCl2kBwiJLCXpzLzUBLgmp4veFZdvN5ChW4Eq/8Fc2Fg==}
engines: {node: '>=18'}
cpu: [arm64]
os: [linux]
'@esbuild/linux-arm@0.27.3':
resolution: {integrity: sha512-s6nPv2QkSupJwLYyfS+gwdirm0ukyTFNl3KTgZEAiJDd+iHZcbTPPcWCcRYH+WlNbwChgH2QkE9NSlNrMT8Gfw==}
engines: {node: '>=18'}
cpu: [arm]
os: [linux]
'@esbuild/linux-ia32@0.27.3':
resolution: {integrity: sha512-yGlQYjdxtLdh0a3jHjuwOrxQjOZYD/C9PfdbgJJF3TIZWnm/tMd/RcNiLngiu4iwcBAOezdnSLAwQDPqTmtTYg==}
engines: {node: '>=18'}
cpu: [ia32]
os: [linux]
'@esbuild/linux-loong64@0.27.3':
resolution: {integrity: sha512-WO60Sn8ly3gtzhyjATDgieJNet/KqsDlX5nRC5Y3oTFcS1l0KWba+SEa9Ja1GfDqSF1z6hif/SkpQJbL63cgOA==}
engines: {node: '>=18'}
cpu: [loong64]
os: [linux]
'@esbuild/linux-mips64el@0.27.3':
resolution: {integrity: sha512-APsymYA6sGcZ4pD6k+UxbDjOFSvPWyZhjaiPyl/f79xKxwTnrn5QUnXR5prvetuaSMsb4jgeHewIDCIWljrSxw==}
engines: {node: '>=18'}
cpu: [mips64el]
os: [linux]
'@esbuild/linux-ppc64@0.27.3':
resolution: {integrity: sha512-eizBnTeBefojtDb9nSh4vvVQ3V9Qf9Df01PfawPcRzJH4gFSgrObw+LveUyDoKU3kxi5+9RJTCWlj4FjYXVPEA==}
engines: {node: '>=18'}
cpu: [ppc64]
os: [linux]
'@esbuild/linux-riscv64@0.27.3':
resolution: {integrity: sha512-3Emwh0r5wmfm3ssTWRQSyVhbOHvqegUDRd0WhmXKX2mkHJe1SFCMJhagUleMq+Uci34wLSipf8Lagt4LlpRFWQ==}
engines: {node: '>=18'}
cpu: [riscv64]
os: [linux]
'@esbuild/linux-s390x@0.27.3':
resolution: {integrity: sha512-pBHUx9LzXWBc7MFIEEL0yD/ZVtNgLytvx60gES28GcWMqil8ElCYR4kvbV2BDqsHOvVDRrOxGySBM9Fcv744hw==}
engines: {node: '>=18'}
cpu: [s390x]
os: [linux]
'@esbuild/linux-x64@0.27.3':
resolution: {integrity: sha512-Czi8yzXUWIQYAtL/2y6vogER8pvcsOsk5cpwL4Gk5nJqH5UZiVByIY8Eorm5R13gq+DQKYg0+JyQoytLQas4dA==}
engines: {node: '>=18'}
cpu: [x64]
os: [linux]
'@esbuild/netbsd-arm64@0.27.3':
resolution: {integrity: sha512-sDpk0RgmTCR/5HguIZa9n9u+HVKf40fbEUt+iTzSnCaGvY9kFP0YKBWZtJaraonFnqef5SlJ8/TiPAxzyS+UoA==}
engines: {node: '>=18'}
cpu: [arm64]
os: [netbsd]
'@esbuild/netbsd-x64@0.27.3':
resolution: {integrity: sha512-P14lFKJl/DdaE00LItAukUdZO5iqNH7+PjoBm+fLQjtxfcfFE20Xf5CrLsmZdq5LFFZzb5JMZ9grUwvtVYzjiA==}
engines: {node: '>=18'}
cpu: [x64]
os: [netbsd]
'@esbuild/openbsd-arm64@0.27.3':
resolution: {integrity: sha512-AIcMP77AvirGbRl/UZFTq5hjXK+2wC7qFRGoHSDrZ5v5b8DK/GYpXW3CPRL53NkvDqb9D+alBiC/dV0Fb7eJcw==}
engines: {node: '>=18'}
cpu: [arm64]
os: [openbsd]
'@esbuild/openbsd-x64@0.27.3':
resolution: {integrity: sha512-DnW2sRrBzA+YnE70LKqnM3P+z8vehfJWHXECbwBmH/CU51z6FiqTQTHFenPlHmo3a8UgpLyH3PT+87OViOh1AQ==}
engines: {node: '>=18'}
cpu: [x64]
os: [openbsd]
'@esbuild/openharmony-arm64@0.27.3':
resolution: {integrity: sha512-NinAEgr/etERPTsZJ7aEZQvvg/A6IsZG/LgZy+81wON2huV7SrK3e63dU0XhyZP4RKGyTm7aOgmQk0bGp0fy2g==}
engines: {node: '>=18'}
cpu: [arm64]
os: [openharmony]
'@esbuild/sunos-x64@0.27.3':
resolution: {integrity: sha512-PanZ+nEz+eWoBJ8/f8HKxTTD172SKwdXebZ0ndd953gt1HRBbhMsaNqjTyYLGLPdoWHy4zLU7bDVJztF5f3BHA==}
engines: {node: '>=18'}
cpu: [x64]
os: [sunos]
'@esbuild/win32-arm64@0.27.3':
resolution: {integrity: sha512-B2t59lWWYrbRDw/tjiWOuzSsFh1Y/E95ofKz7rIVYSQkUYBjfSgf6oeYPNWHToFRr2zx52JKApIcAS/D5TUBnA==}
engines: {node: '>=18'}
cpu: [arm64]
os: [win32]
'@esbuild/win32-ia32@0.27.3':
resolution: {integrity: sha512-QLKSFeXNS8+tHW7tZpMtjlNb7HKau0QDpwm49u0vUp9y1WOF+PEzkU84y9GqYaAVW8aH8f3GcBck26jh54cX4Q==}
engines: {node: '>=18'}
cpu: [ia32]
os: [win32]
'@esbuild/win32-x64@0.27.3':
resolution: {integrity: sha512-4uJGhsxuptu3OcpVAzli+/gWusVGwZZHTlS63hh++ehExkVT8SgiEf7/uC/PclrPPkLhZqGgCTjd0VWLo6xMqA==}
engines: {node: '>=18'}
cpu: [x64]
os: [win32]
'@jridgewell/gen-mapping@0.3.13':
resolution: {integrity: sha512-2kkt/7niJ6MgEPxF0bYdQ6etZaA+fQvDcLKckhy1yIQOzaoKjBBjSj63/aLVjYE3qhRt5dvM+uUyfCg6UKCBbA==}
'@jridgewell/resolve-uri@3.1.2':
resolution: {integrity: sha512-bRISgCIjP20/tbWSPWMEi54QVPRZExkuD9lJL+UIxUKtwVJA8wW1Trb1jMs1RFXo1CBTNZ/5hpC9QvmKWdopKw==}
engines: {node: '>=6.0.0'}
'@jridgewell/sourcemap-codec@1.5.5':
resolution: {integrity: sha512-cYQ9310grqxueWbl+WuIUIaiUaDcj7WOq5fVhEljNVgRfOUhY9fy2zTvfoqWsnebh8Sl70VScFbICvJnLKB0Og==}
'@jridgewell/trace-mapping@0.3.31':
resolution: {integrity: sha512-zzNR+SdQSDJzc8joaeP8QQoCQr8NuYx2dIIytl1QeBEZHJ9uW6hebsrYgbz8hJwUQao3TWCMtmfV8Nu1twOLAw==}
'@rollup/rollup-android-arm-eabi@4.59.0':
resolution: {integrity: sha512-upnNBkA6ZH2VKGcBj9Fyl9IGNPULcjXRlg0LLeaioQWueH30p6IXtJEbKAgvyv+mJaMxSm1l6xwDXYjpEMiLMg==}
cpu: [arm]
os: [android]
'@rollup/rollup-android-arm64@4.59.0':
resolution: {integrity: sha512-hZ+Zxj3SySm4A/DylsDKZAeVg0mvi++0PYVceVyX7hemkw7OreKdCvW2oQ3T1FMZvCaQXqOTHb8qmBShoqk69Q==}
cpu: [arm64]
os: [android]
'@rollup/rollup-darwin-arm64@4.59.0':
resolution: {integrity: sha512-W2Psnbh1J8ZJw0xKAd8zdNgF9HRLkdWwwdWqubSVk0pUuQkoHnv7rx4GiF9rT4t5DIZGAsConRE3AxCdJ4m8rg==}
cpu: [arm64]
os: [darwin]
'@rollup/rollup-darwin-x64@4.59.0':
resolution: {integrity: sha512-ZW2KkwlS4lwTv7ZVsYDiARfFCnSGhzYPdiOU4IM2fDbL+QGlyAbjgSFuqNRbSthybLbIJ915UtZBtmuLrQAT/w==}
cpu: [x64]
os: [darwin]
'@rollup/rollup-freebsd-arm64@4.59.0':
resolution: {integrity: sha512-EsKaJ5ytAu9jI3lonzn3BgG8iRBjV4LxZexygcQbpiU0wU0ATxhNVEpXKfUa0pS05gTcSDMKpn3Sx+QB9RlTTA==}
cpu: [arm64]
os: [freebsd]
'@rollup/rollup-freebsd-x64@4.59.0':
resolution: {integrity: sha512-d3DuZi2KzTMjImrxoHIAODUZYoUUMsuUiY4SRRcJy6NJoZ6iIqWnJu9IScV9jXysyGMVuW+KNzZvBLOcpdl3Vg==}
cpu: [x64]
os: [freebsd]
'@rollup/rollup-linux-arm-gnueabihf@4.59.0':
resolution: {integrity: sha512-t4ONHboXi/3E0rT6OZl1pKbl2Vgxf9vJfWgmUoCEVQVxhW6Cw/c8I6hbbu7DAvgp82RKiH7TpLwxnJeKv2pbsw==}
cpu: [arm]
os: [linux]
libc: [glibc]
'@rollup/rollup-linux-arm-musleabihf@4.59.0':
resolution: {integrity: sha512-CikFT7aYPA2ufMD086cVORBYGHffBo4K8MQ4uPS/ZnY54GKj36i196u8U+aDVT2LX4eSMbyHtyOh7D7Zvk2VvA==}
cpu: [arm]
os: [linux]
libc: [musl]
'@rollup/rollup-linux-arm64-gnu@4.59.0':
resolution: {integrity: sha512-jYgUGk5aLd1nUb1CtQ8E+t5JhLc9x5WdBKew9ZgAXg7DBk0ZHErLHdXM24rfX+bKrFe+Xp5YuJo54I5HFjGDAA==}
cpu: [arm64]
os: [linux]
libc: [glibc]
'@rollup/rollup-linux-arm64-musl@4.59.0':
resolution: {integrity: sha512-peZRVEdnFWZ5Bh2KeumKG9ty7aCXzzEsHShOZEFiCQlDEepP1dpUl/SrUNXNg13UmZl+gzVDPsiCwnV1uI0RUA==}
cpu: [arm64]
os: [linux]
libc: [musl]
'@rollup/rollup-linux-loong64-gnu@4.59.0':
resolution: {integrity: sha512-gbUSW/97f7+r4gHy3Jlup8zDG190AuodsWnNiXErp9mT90iCy9NKKU0Xwx5k8VlRAIV2uU9CsMnEFg/xXaOfXg==}
cpu: [loong64]
os: [linux]
libc: [glibc]
'@rollup/rollup-linux-loong64-musl@4.59.0':
resolution: {integrity: sha512-yTRONe79E+o0FWFijasoTjtzG9EBedFXJMl888NBEDCDV9I2wGbFFfJQQe63OijbFCUZqxpHz1GzpbtSFikJ4Q==}
cpu: [loong64]
os: [linux]
libc: [musl]
'@rollup/rollup-linux-ppc64-gnu@4.59.0':
resolution: {integrity: sha512-sw1o3tfyk12k3OEpRddF68a1unZ5VCN7zoTNtSn2KndUE+ea3m3ROOKRCZxEpmT9nsGnogpFP9x6mnLTCaoLkA==}
cpu: [ppc64]
os: [linux]
libc: [glibc]
'@rollup/rollup-linux-ppc64-musl@4.59.0':
resolution: {integrity: sha512-+2kLtQ4xT3AiIxkzFVFXfsmlZiG5FXYW7ZyIIvGA7Bdeuh9Z0aN4hVyXS/G1E9bTP/vqszNIN/pUKCk/BTHsKA==}
cpu: [ppc64]
os: [linux]
libc: [musl]
'@rollup/rollup-linux-riscv64-gnu@4.59.0':
resolution: {integrity: sha512-NDYMpsXYJJaj+I7UdwIuHHNxXZ/b/N2hR15NyH3m2qAtb/hHPA4g4SuuvrdxetTdndfj9b1WOmy73kcPRoERUg==}
cpu: [riscv64]
os: [linux]
libc: [glibc]
'@rollup/rollup-linux-riscv64-musl@4.59.0':
resolution: {integrity: sha512-nLckB8WOqHIf1bhymk+oHxvM9D3tyPndZH8i8+35p/1YiVoVswPid2yLzgX7ZJP0KQvnkhM4H6QZ5m0LzbyIAg==}
cpu: [riscv64]
os: [linux]
libc: [musl]
'@rollup/rollup-linux-s390x-gnu@4.59.0':
resolution: {integrity: sha512-oF87Ie3uAIvORFBpwnCvUzdeYUqi2wY6jRFWJAy1qus/udHFYIkplYRW+wo+GRUP4sKzYdmE1Y3+rY5Gc4ZO+w==}
cpu: [s390x]
os: [linux]
libc: [glibc]
'@rollup/rollup-linux-x64-gnu@4.59.0':
resolution: {integrity: sha512-3AHmtQq/ppNuUspKAlvA8HtLybkDflkMuLK4DPo77DfthRb71V84/c4MlWJXixZz4uruIH4uaa07IqoAkG64fg==}
cpu: [x64]
os: [linux]
libc: [glibc]
'@rollup/rollup-linux-x64-musl@4.59.0':
resolution: {integrity: sha512-2UdiwS/9cTAx7qIUZB/fWtToJwvt0Vbo0zmnYt7ED35KPg13Q0ym1g442THLC7VyI6JfYTP4PiSOWyoMdV2/xg==}
cpu: [x64]
os: [linux]
libc: [musl]
'@rollup/rollup-openbsd-x64@4.59.0':
resolution: {integrity: sha512-M3bLRAVk6GOwFlPTIxVBSYKUaqfLrn8l0psKinkCFxl4lQvOSz8ZrKDz2gxcBwHFpci0B6rttydI4IpS4IS/jQ==}
cpu: [x64]
os: [openbsd]
'@rollup/rollup-openharmony-arm64@4.59.0':
resolution: {integrity: sha512-tt9KBJqaqp5i5HUZzoafHZX8b5Q2Fe7UjYERADll83O4fGqJ49O1FsL6LpdzVFQcpwvnyd0i+K/VSwu/o/nWlA==}
cpu: [arm64]
os: [openharmony]
'@rollup/rollup-win32-arm64-msvc@4.59.0':
resolution: {integrity: sha512-V5B6mG7OrGTwnxaNUzZTDTjDS7F75PO1ae6MJYdiMu60sq0CqN5CVeVsbhPxalupvTX8gXVSU9gq+Rx1/hvu6A==}
cpu: [arm64]
os: [win32]
'@rollup/rollup-win32-ia32-msvc@4.59.0':
resolution: {integrity: sha512-UKFMHPuM9R0iBegwzKF4y0C4J9u8C6MEJgFuXTBerMk7EJ92GFVFYBfOZaSGLu6COf7FxpQNqhNS4c4icUPqxA==}
cpu: [ia32]
os: [win32]
'@rollup/rollup-win32-x64-gnu@4.59.0':
resolution: {integrity: sha512-laBkYlSS1n2L8fSo1thDNGrCTQMmxjYY5G0WFWjFFYZkKPjsMBsgJfGf4TLxXrF6RyhI60L8TMOjBMvXiTcxeA==}
cpu: [x64]
os: [win32]
'@rollup/rollup-win32-x64-msvc@4.59.0':
resolution: {integrity: sha512-2HRCml6OztYXyJXAvdDXPKcawukWY2GpR5/nxKp4iBgiO3wcoEGkAaqctIbZcNB6KlUQBIqt8VYkNSj2397EfA==}
cpu: [x64]
os: [win32]
'@types/estree@1.0.8':
resolution: {integrity: sha512-dWHzHa2WqEXI/O1E9OjrocMTKJl2mSrEolh1Iomrv6U+JuNwaHXsXx9bLu5gG7BUWFIN0skIQJQ/L1rIex4X6w==}
'@types/node@25.3.0':
resolution: {integrity: sha512-4K3bqJpXpqfg2XKGK9bpDTc6xO/xoUP/RBWS7AtRMug6zZFaRekiLzjVtAoZMquxoAbzBvy5nxQ7veS5eYzf8A==}
acorn@8.16.0:
resolution: {integrity: sha512-UVJyE9MttOsBQIDKw1skb9nAwQuR5wuGD3+82K6JgJlm/Y+KI92oNsMNGZCYdDsVtRHSak0pcV5Dno5+4jh9sw==}
engines: {node: '>=0.4.0'}
hasBin: true
any-promise@1.3.0:
resolution: {integrity: sha512-7UvmKalWRt1wgjL1RrGxoSJW/0QZFIegpeGvZG9kjp8vrRu55XTHbwnqq2GpXm9uLbcuhxm3IqX9OB4MZR1b2A==}
bundle-require@5.1.0:
resolution: {integrity: sha512-3WrrOuZiyaaZPWiEt4G3+IffISVC9HYlWueJEBWED4ZH4aIAC2PnkdnuRrR94M+w6yGWn4AglWtJtBI8YqvgoA==}
engines: {node: ^12.20.0 || ^14.13.1 || >=16.0.0}
peerDependencies:
esbuild: '>=0.18'
cac@6.7.14:
resolution: {integrity: sha512-b6Ilus+c3RrdDk+JhLKUAQfzzgLEPy6wcXqS7f/xe1EETvsDP6GORG7SFuOs6cID5YkqchW/LXZbX5bc8j7ZcQ==}
engines: {node: '>=8'}
chokidar@4.0.3:
resolution: {integrity: sha512-Qgzu8kfBvo+cA4962jnP1KkS6Dop5NS6g7R5LFYJr4b8Ub94PPQXUksCw9PvXoeXPRRddRNC5C1JQUR2SMGtnA==}
engines: {node: '>= 14.16.0'}
commander@4.1.1:
resolution: {integrity: sha512-NOKm8xhkzAjzFx8B2v5OAHT+u5pRQc2UCa2Vq9jYL/31o2wi9mxBA7LIFs3sV5VSC49z6pEhfbMULvShKj26WA==}
engines: {node: '>= 6'}
confbox@0.1.8:
resolution: {integrity: sha512-RMtmw0iFkeR4YV+fUOSucriAQNb9g8zFR52MWCtl+cCZOFRNL6zeB395vPzFhEjjn4fMxXudmELnl/KF/WrK6w==}
consola@3.4.2:
resolution: {integrity: sha512-5IKcdX0nnYavi6G7TtOhwkYzyjfJlatbjMjuLSfE2kYT5pMDOilZ4OvMhi637CcDICTmz3wARPoyhqyX1Y+XvA==}
engines: {node: ^14.18.0 || >=16.10.0}
debug@4.4.3:
resolution: {integrity: sha512-RGwwWnwQvkVfavKVt22FGLw+xYSdzARwm0ru6DhTVA3umU5hZc28V3kO4stgYryrTlLpuvgI9GiijltAjNbcqA==}
engines: {node: '>=6.0'}
peerDependencies:
supports-color: '*'
peerDependenciesMeta:
supports-color:
optional: true
esbuild@0.27.3:
resolution: {integrity: sha512-8VwMnyGCONIs6cWue2IdpHxHnAjzxnw2Zr7MkVxB2vjmQ2ivqGFb4LEG3SMnv0Gb2F/G/2yA8zUaiL1gywDCCg==}
engines: {node: '>=18'}
hasBin: true
fdir@6.5.0:
resolution: {integrity: sha512-tIbYtZbucOs0BRGqPJkshJUYdL+SDH7dVM8gjy+ERp3WAUjLEFJE+02kanyHtwjWOnwrKYBiwAmM0p4kLJAnXg==}
engines: {node: '>=12.0.0'}
peerDependencies:
picomatch: ^3 || ^4
peerDependenciesMeta:
picomatch:
optional: true
fix-dts-default-cjs-exports@1.0.1:
resolution: {integrity: sha512-pVIECanWFC61Hzl2+oOCtoJ3F17kglZC/6N94eRWycFgBH35hHx0Li604ZIzhseh97mf2p0cv7vVrOZGoqhlEg==}
fsevents@2.3.3:
resolution: {integrity: sha512-5xoDfX+fL7faATnagmWPpbFtwh/R77WmMMqqHGS65C3vvB0YHrgF+B1YmZ3441tMj5n63k0212XNoJwzlhffQw==}
engines: {node: ^8.16.0 || ^10.6.0 || >=11.0.0}
os: [darwin]
joycon@3.1.1:
resolution: {integrity: sha512-34wB/Y7MW7bzjKRjUKTa46I2Z7eV62Rkhva+KkopW7Qvv/OSWBqvkSY7vusOPrNuZcUG3tApvdVgNB8POj3SPw==}
engines: {node: '>=10'}
lilconfig@3.1.3:
resolution: {integrity: sha512-/vlFKAoH5Cgt3Ie+JLhRbwOsCQePABiU3tJ1egGvyQ+33R/vcwM2Zl2QR/LzjsBeItPt3oSVXapn+m4nQDvpzw==}
engines: {node: '>=14'}
lines-and-columns@1.2.4:
resolution: {integrity: sha512-7ylylesZQ/PV29jhEDl3Ufjo6ZX7gCqJr5F7PKrqc93v7fzSymt1BpwEU8nAUXs8qzzvqhbjhK5QZg6Mt/HkBg==}
load-tsconfig@0.2.5:
resolution: {integrity: sha512-IXO6OCs9yg8tMKzfPZ1YmheJbZCiEsnBdcB03l0OcfK9prKnJb96siuHCr5Fl37/yo9DnKU+TLpxzTUspw9shg==}
engines: {node: ^12.20.0 || ^14.13.1 || >=16.0.0}
magic-string@0.30.21:
resolution: {integrity: sha512-vd2F4YUyEXKGcLHoq+TEyCjxueSeHnFxyyjNp80yg0XV4vUhnDer/lvvlqM/arB5bXQN5K2/3oinyCRyx8T2CQ==}
mlly@1.8.0:
resolution: {integrity: sha512-l8D9ODSRWLe2KHJSifWGwBqpTZXIXTeo8mlKjY+E2HAakaTeNpqAyBZ8GSqLzHgw4XmHmC8whvpjJNMbFZN7/g==}
ms@2.1.3:
resolution: {integrity: sha512-6FlzubTLZG3J2a/NVCAleEhjzq5oxgHyaCU9yYXvcLsvoVaHJq/s5xXI6/XXP6tz7R9xAOtHnSO/tXtF3WRTlA==}
mz@2.7.0:
resolution: {integrity: sha512-z81GNO7nnYMEhrGh9LeymoE4+Yr0Wn5McHIZMK5cfQCl+NDX08sCZgUc9/6MHni9IWuFLm1Z3HTCXu2z9fN62Q==}
object-assign@4.1.1:
resolution: {integrity: sha512-rJgTQnkUnH1sFw8yT6VSU3zD3sWmu6sZhIseY8VX+GRu3P6F7Fu+JNDoXfklElbLJSnc3FUQHVe4cU5hj+BcUg==}
engines: {node: '>=0.10.0'}
pathe@2.0.3:
resolution: {integrity: sha512-WUjGcAqP1gQacoQe+OBJsFA7Ld4DyXuUIjZ5cc75cLHvJ7dtNsTugphxIADwspS+AraAUePCKrSVtPLFj/F88w==}
picocolors@1.1.1:
resolution: {integrity: sha512-xceH2snhtb5M9liqDsmEw56le376mTZkEX/jEb/RxNFyegNul7eNslCXP9FDj/Lcu0X8KEyMceP2ntpaHrDEVA==}
picomatch@4.0.3:
resolution: {integrity: sha512-5gTmgEY/sqK6gFXLIsQNH19lWb4ebPDLA4SdLP7dsWkIXHWlG66oPuVvXSGFPppYZz8ZDZq0dYYrbHfBCVUb1Q==}
engines: {node: '>=12'}
pirates@4.0.7:
resolution: {integrity: sha512-TfySrs/5nm8fQJDcBDuUng3VOUKsd7S+zqvbOTiGXHfxX4wK31ard+hoNuvkicM/2YFzlpDgABOevKSsB4G/FA==}
engines: {node: '>= 6'}
pkg-types@1.3.1:
resolution: {integrity: sha512-/Jm5M4RvtBFVkKWRu2BLUTNP8/M2a+UwuAX+ae4770q1qVGtfjG+WTCupoZixokjmHiry8uI+dlY8KXYV5HVVQ==}
postcss-load-config@6.0.1:
resolution: {integrity: sha512-oPtTM4oerL+UXmx+93ytZVN82RrlY/wPUV8IeDxFrzIjXOLF1pN+EmKPLbubvKHT2HC20xXsCAH2Z+CKV6Oz/g==}
engines: {node: '>= 18'}
peerDependencies:
jiti: '>=1.21.0'
postcss: '>=8.0.9'
tsx: ^4.8.1
yaml: ^2.4.2
peerDependenciesMeta:
jiti:
optional: true
postcss:
optional: true
tsx:
optional: true
yaml:
optional: true
readdirp@4.1.2:
resolution: {integrity: sha512-GDhwkLfywWL2s6vEjyhri+eXmfH6j1L7JE27WhqLeYzoh/A3DBaYGEj2H/HFZCn/kMfim73FXxEJTw06WtxQwg==}
engines: {node: '>= 14.18.0'}
resolve-from@5.0.0:
resolution: {integrity: sha512-qYg9KP24dD5qka9J47d0aVky0N+b4fTU89LN9iDnjB5waksiC49rvMB0PrUJQGoTmH50XPiqOvAjDfaijGxYZw==}
engines: {node: '>=8'}
rollup@4.59.0:
resolution: {integrity: sha512-2oMpl67a3zCH9H79LeMcbDhXW/UmWG/y2zuqnF2jQq5uq9TbM9TVyXvA4+t+ne2IIkBdrLpAaRQAvo7YI/Yyeg==}
engines: {node: '>=18.0.0', npm: '>=8.0.0'}
hasBin: true
source-map@0.7.6:
resolution: {integrity: sha512-i5uvt8C3ikiWeNZSVZNWcfZPItFQOsYTUAOkcUPGd8DqDy1uOUikjt5dG+uRlwyvR108Fb9DOd4GvXfT0N2/uQ==}
engines: {node: '>= 12'}
sucrase@3.35.1:
resolution: {integrity: sha512-DhuTmvZWux4H1UOnWMB3sk0sbaCVOoQZjv8u1rDoTV0HTdGem9hkAZtl4JZy8P2z4Bg0nT+YMeOFyVr4zcG5Tw==}
engines: {node: '>=16 || 14 >=14.17'}
hasBin: true
thenify-all@1.6.0:
resolution: {integrity: sha512-RNxQH/qI8/t3thXJDwcstUO4zeqo64+Uy/+sNVRBx4Xn2OX+OZ9oP+iJnNFqplFra2ZUVeKCSa2oVWi3T4uVmA==}
engines: {node: '>=0.8'}
thenify@3.3.1:
resolution: {integrity: sha512-RVZSIV5IG10Hk3enotrhvz0T9em6cyHBLkH/YAZuKqd8hRkKhSfCGIcP2KUY0EPxndzANBmNllzWPwak+bheSw==}
tinyexec@0.3.2:
resolution: {integrity: sha512-KQQR9yN7R5+OSwaK0XQoj22pwHoTlgYqmUscPYoknOoWCWfj/5/ABTMRi69FrKU5ffPVh5QcFikpWJI/P1ocHA==}
tinyglobby@0.2.15:
resolution: {integrity: sha512-j2Zq4NyQYG5XMST4cbs02Ak8iJUdxRM0XI5QyxXuZOzKOINmWurp3smXu3y5wDcJrptwpSjgXHzIQxR0omXljQ==}
engines: {node: '>=12.0.0'}
tree-kill@1.2.2:
resolution: {integrity: sha512-L0Orpi8qGpRG//Nd+H90vFB+3iHnue1zSSGmNOOCh1GLJ7rUKVwV2HvijphGQS2UmhUZewS9VgvxYIdgr+fG1A==}
hasBin: true
ts-interface-checker@0.1.13:
resolution: {integrity: sha512-Y/arvbn+rrz3JCKl9C4kVNfTfSm2/mEp5FSz5EsZSANGPSlQrpRI5M4PKF+mJnE52jOO90PnPSc3Ur3bTQw0gA==}
tsup@8.5.1:
resolution: {integrity: sha512-xtgkqwdhpKWr3tKPmCkvYmS9xnQK3m3XgxZHwSUjvfTjp7YfXe5tT3GgWi0F2N+ZSMsOeWeZFh7ZZFg5iPhing==}
engines: {node: '>=18'}
hasBin: true
peerDependencies:
'@microsoft/api-extractor': ^7.36.0
'@swc/core': ^1
postcss: ^8.4.12
typescript: '>=4.5.0'
peerDependenciesMeta:
'@microsoft/api-extractor':
optional: true
'@swc/core':
optional: true
postcss:
optional: true
typescript:
optional: true
typescript@5.9.3:
resolution: {integrity: sha512-jl1vZzPDinLr9eUt3J/t7V6FgNEw9QjvBPdysz9KfQDD41fQrC2Y4vKQdiaUpFT4bXlb1RHhLpp8wtm6M5TgSw==}
engines: {node: '>=14.17'}
hasBin: true
ufo@1.6.3:
resolution: {integrity: sha512-yDJTmhydvl5lJzBmy/hyOAA0d+aqCBuwl818haVdYCRrWV84o7YyeVm4QlVHStqNrrJSTb6jKuFAVqAFsr+K3Q==}
undici-types@7.18.2:
resolution: {integrity: sha512-AsuCzffGHJybSaRrmr5eHr81mwJU3kjw6M+uprWvCXiNeN9SOGwQ3Jn8jb8m3Z6izVgknn1R0FTCEAP2QrLY/w==}
snapshots:
'@esbuild/aix-ppc64@0.27.3':
optional: true
'@esbuild/android-arm64@0.27.3':
optional: true
'@esbuild/android-arm@0.27.3':
optional: true
'@esbuild/android-x64@0.27.3':
optional: true
'@esbuild/darwin-arm64@0.27.3':
optional: true
'@esbuild/darwin-x64@0.27.3':
optional: true
'@esbuild/freebsd-arm64@0.27.3':
optional: true
'@esbuild/freebsd-x64@0.27.3':
optional: true
'@esbuild/linux-arm64@0.27.3':
optional: true
'@esbuild/linux-arm@0.27.3':
optional: true
'@esbuild/linux-ia32@0.27.3':
optional: true
'@esbuild/linux-loong64@0.27.3':
optional: true
'@esbuild/linux-mips64el@0.27.3':
optional: true
'@esbuild/linux-ppc64@0.27.3':
optional: true
'@esbuild/linux-riscv64@0.27.3':
optional: true
'@esbuild/linux-s390x@0.27.3':
optional: true
'@esbuild/linux-x64@0.27.3':
optional: true
'@esbuild/netbsd-arm64@0.27.3':
optional: true
'@esbuild/netbsd-x64@0.27.3':
optional: true
'@esbuild/openbsd-arm64@0.27.3':
optional: true
'@esbuild/openbsd-x64@0.27.3':
optional: true
'@esbuild/openharmony-arm64@0.27.3':
optional: true
'@esbuild/sunos-x64@0.27.3':
optional: true
'@esbuild/win32-arm64@0.27.3':
optional: true
'@esbuild/win32-ia32@0.27.3':
optional: true
'@esbuild/win32-x64@0.27.3':
optional: true
'@jridgewell/gen-mapping@0.3.13':
dependencies:
'@jridgewell/sourcemap-codec': 1.5.5
'@jridgewell/trace-mapping': 0.3.31
'@jridgewell/resolve-uri@3.1.2': {}
'@jridgewell/sourcemap-codec@1.5.5': {}
'@jridgewell/trace-mapping@0.3.31':
dependencies:
'@jridgewell/resolve-uri': 3.1.2
'@jridgewell/sourcemap-codec': 1.5.5
'@rollup/rollup-android-arm-eabi@4.59.0':
optional: true
'@rollup/rollup-android-arm64@4.59.0':
optional: true
'@rollup/rollup-darwin-arm64@4.59.0':
optional: true
'@rollup/rollup-darwin-x64@4.59.0':
optional: true
'@rollup/rollup-freebsd-arm64@4.59.0':
optional: true
'@rollup/rollup-freebsd-x64@4.59.0':
optional: true
'@rollup/rollup-linux-arm-gnueabihf@4.59.0':
optional: true
'@rollup/rollup-linux-arm-musleabihf@4.59.0':
optional: true
'@rollup/rollup-linux-arm64-gnu@4.59.0':
optional: true
'@rollup/rollup-linux-arm64-musl@4.59.0':
optional: true
'@rollup/rollup-linux-loong64-gnu@4.59.0':
optional: true
'@rollup/rollup-linux-loong64-musl@4.59.0':
optional: true
'@rollup/rollup-linux-ppc64-gnu@4.59.0':
optional: true
'@rollup/rollup-linux-ppc64-musl@4.59.0':
optional: true
'@rollup/rollup-linux-riscv64-gnu@4.59.0':
optional: true
'@rollup/rollup-linux-riscv64-musl@4.59.0':
optional: true
'@rollup/rollup-linux-s390x-gnu@4.59.0':
optional: true
'@rollup/rollup-linux-x64-gnu@4.59.0':
optional: true
'@rollup/rollup-linux-x64-musl@4.59.0':
optional: true
'@rollup/rollup-openbsd-x64@4.59.0':
optional: true
'@rollup/rollup-openharmony-arm64@4.59.0':
optional: true
'@rollup/rollup-win32-arm64-msvc@4.59.0':
optional: true
'@rollup/rollup-win32-ia32-msvc@4.59.0':
optional: true
'@rollup/rollup-win32-x64-gnu@4.59.0':
optional: true
'@rollup/rollup-win32-x64-msvc@4.59.0':
optional: true
'@types/estree@1.0.8': {}
'@types/node@25.3.0':
dependencies:
undici-types: 7.18.2
acorn@8.16.0: {}
any-promise@1.3.0: {}
bundle-require@5.1.0(esbuild@0.27.3):
dependencies:
esbuild: 0.27.3
load-tsconfig: 0.2.5
cac@6.7.14: {}
chokidar@4.0.3:
dependencies:
readdirp: 4.1.2
commander@4.1.1: {}
confbox@0.1.8: {}
consola@3.4.2: {}
debug@4.4.3:
dependencies:
ms: 2.1.3
esbuild@0.27.3:
optionalDependencies:
'@esbuild/aix-ppc64': 0.27.3
'@esbuild/android-arm': 0.27.3
'@esbuild/android-arm64': 0.27.3
'@esbuild/android-x64': 0.27.3
'@esbuild/darwin-arm64': 0.27.3
'@esbuild/darwin-x64': 0.27.3
'@esbuild/freebsd-arm64': 0.27.3
'@esbuild/freebsd-x64': 0.27.3
'@esbuild/linux-arm': 0.27.3
'@esbuild/linux-arm64': 0.27.3
'@esbuild/linux-ia32': 0.27.3
'@esbuild/linux-loong64': 0.27.3
'@esbuild/linux-mips64el': 0.27.3
'@esbuild/linux-ppc64': 0.27.3
'@esbuild/linux-riscv64': 0.27.3
'@esbuild/linux-s390x': 0.27.3
'@esbuild/linux-x64': 0.27.3
'@esbuild/netbsd-arm64': 0.27.3
'@esbuild/netbsd-x64': 0.27.3
'@esbuild/openbsd-arm64': 0.27.3
'@esbuild/openbsd-x64': 0.27.3
'@esbuild/openharmony-arm64': 0.27.3
'@esbuild/sunos-x64': 0.27.3
'@esbuild/win32-arm64': 0.27.3
'@esbuild/win32-ia32': 0.27.3
'@esbuild/win32-x64': 0.27.3
fdir@6.5.0(picomatch@4.0.3):
optionalDependencies:
picomatch: 4.0.3
fix-dts-default-cjs-exports@1.0.1:
dependencies:
magic-string: 0.30.21
mlly: 1.8.0
rollup: 4.59.0
fsevents@2.3.3:
optional: true
joycon@3.1.1: {}
lilconfig@3.1.3: {}
lines-and-columns@1.2.4: {}
load-tsconfig@0.2.5: {}
magic-string@0.30.21:
dependencies:
'@jridgewell/sourcemap-codec': 1.5.5
mlly@1.8.0:
dependencies:
acorn: 8.16.0
pathe: 2.0.3
pkg-types: 1.3.1
ufo: 1.6.3
ms@2.1.3: {}
mz@2.7.0:
dependencies:
any-promise: 1.3.0
object-assign: 4.1.1
thenify-all: 1.6.0
object-assign@4.1.1: {}
pathe@2.0.3: {}
picocolors@1.1.1: {}
picomatch@4.0.3: {}
pirates@4.0.7: {}
pkg-types@1.3.1:
dependencies:
confbox: 0.1.8
mlly: 1.8.0
pathe: 2.0.3
postcss-load-config@6.0.1:
dependencies:
lilconfig: 3.1.3
readdirp@4.1.2: {}
resolve-from@5.0.0: {}
rollup@4.59.0:
dependencies:
'@types/estree': 1.0.8
optionalDependencies:
'@rollup/rollup-android-arm-eabi': 4.59.0
'@rollup/rollup-android-arm64': 4.59.0
'@rollup/rollup-darwin-arm64': 4.59.0
'@rollup/rollup-darwin-x64': 4.59.0
'@rollup/rollup-freebsd-arm64': 4.59.0
'@rollup/rollup-freebsd-x64': 4.59.0
'@rollup/rollup-linux-arm-gnueabihf': 4.59.0
'@rollup/rollup-linux-arm-musleabihf': 4.59.0
'@rollup/rollup-linux-arm64-gnu': 4.59.0
'@rollup/rollup-linux-arm64-musl': 4.59.0
'@rollup/rollup-linux-loong64-gnu': 4.59.0
'@rollup/rollup-linux-loong64-musl': 4.59.0
'@rollup/rollup-linux-ppc64-gnu': 4.59.0
'@rollup/rollup-linux-ppc64-musl': 4.59.0
'@rollup/rollup-linux-riscv64-gnu': 4.59.0
'@rollup/rollup-linux-riscv64-musl': 4.59.0
'@rollup/rollup-linux-s390x-gnu': 4.59.0
'@rollup/rollup-linux-x64-gnu': 4.59.0
'@rollup/rollup-linux-x64-musl': 4.59.0
'@rollup/rollup-openbsd-x64': 4.59.0
'@rollup/rollup-openharmony-arm64': 4.59.0
'@rollup/rollup-win32-arm64-msvc': 4.59.0
'@rollup/rollup-win32-ia32-msvc': 4.59.0
'@rollup/rollup-win32-x64-gnu': 4.59.0
'@rollup/rollup-win32-x64-msvc': 4.59.0
fsevents: 2.3.3
source-map@0.7.6: {}
sucrase@3.35.1:
dependencies:
'@jridgewell/gen-mapping': 0.3.13
commander: 4.1.1
lines-and-columns: 1.2.4
mz: 2.7.0
pirates: 4.0.7
tinyglobby: 0.2.15
ts-interface-checker: 0.1.13
thenify-all@1.6.0:
dependencies:
thenify: 3.3.1
thenify@3.3.1:
dependencies:
any-promise: 1.3.0
tinyexec@0.3.2: {}
tinyglobby@0.2.15:
dependencies:
fdir: 6.5.0(picomatch@4.0.3)
picomatch: 4.0.3
tree-kill@1.2.2: {}
ts-interface-checker@0.1.13: {}
tsup@8.5.1(typescript@5.9.3):
dependencies:
bundle-require: 5.1.0(esbuild@0.27.3)
cac: 6.7.14
chokidar: 4.0.3
consola: 3.4.2
debug: 4.4.3
esbuild: 0.27.3
fix-dts-default-cjs-exports: 1.0.1
joycon: 3.1.1
picocolors: 1.1.1
postcss-load-config: 6.0.1
resolve-from: 5.0.0
rollup: 4.59.0
source-map: 0.7.6
sucrase: 3.35.1
tinyexec: 0.3.2
tinyglobby: 0.2.15
tree-kill: 1.2.2
optionalDependencies:
typescript: 5.9.3
transitivePeerDependencies:
- jiti
- supports-color
- tsx
- yaml
typescript@5.9.3: {}
ufo@1.6.3: {}
undici-types@7.18.2: {}

2
pnpm-workspace.yaml Normal file
View file

@ -0,0 +1,2 @@
onlyBuiltDependencies:
- esbuild

575
src/api/index.ts Normal file
View file

@ -0,0 +1,575 @@
/**
* api — User-facing functions and kernel management.
*
* This is the only module users need to import directly.
* Everything else in src/ is internal plumbing.
*
* Two operating modes:
*
* 1. Lite mode (no kernel): getMoonPhase() works without DE442S.
* Uses approximate Meeus algorithms. Accurate to ~1° for phase purposes.
* Not suitable for crescent sighting reports.
*
* 2. Full mode (kernel loaded): getMoonSightingReport() and getSunMoonEvents()
* require DE442S. Call initKernels() (or downloadKernels()) first.
*/
import type {
Observer,
SightingOptions,
MoonSightingReport,
MoonPhaseResult,
MoonPhaseName,
SunMoonEvents,
KernelConfig,
Vec3,
} from '../types.js'
import { SpkKernel } from '../spk/index.js'
import {
computeTimeScales,
jdTTtoET,
J2000,
} from '../time/index.js'
import {
getMoonGeocentricState,
getSunGeocentricState,
computeIllumination,
computeCrescentWidth,
getMoonSunApproximate,
nearestNewMoon,
} from '../bodies/index.js'
import {
geodeticToECEF,
computeAzAlt,
} from '../observer/index.js'
import { itrsToGcrs } from '../frames/index.js'
import {
getSunMoonEvents as eventsGetSunMoonEvents,
bestTimeHeuristic,
bestTimeOptimized,
computeObservationWindow,
} from '../events/index.js'
import {
computeCrescentGeometry,
computeYallop,
computeOdeh,
buildGuidanceText,
} from '../visibility/index.js'
// ─── Module-level kernel singleton ─────────────────────────────────────────────
let activeKernel: SpkKernel | null = null
// ─── Cache directory resolution ────────────────────────────────────────────────
/**
* 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-calc`
}
return `${env['HOME'] ?? '/tmp'}/.cache/moon-calc`
}
// ─── 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'
// ─── Kernel lifecycle ─────────────────────────────────────────────────────────
/**
* Initialize the kernel engine from an already-downloaded kernel.
* Must be called before getMoonSightingReport() or getSunMoonEvents().
*
* Supports three source modes:
* - File path (Node.js): initKernels({ planetary: { type: 'file', path: '/path/to/de442s.bsp' } })
* - ArrayBuffer (browser): initKernels({ planetary: { type: 'buffer', data: buf, name: 'de442s.bsp' } })
* - Auto (Node.js): initKernels() — downloads and caches automatically
*
* @param config - Kernel source configuration. Defaults to auto-download.
*/
export async function initKernels(config?: KernelConfig): Promise<void> {
const source = config?.planetary ?? { type: 'auto' as const }
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 (!res.ok) 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)
}
activeKernel = SpkKernel.fromBuffer(buffer)
}
/** Read a file into an ArrayBuffer (Node.js only). */
async function readFileAsBuffer(filePath: string): Promise<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
}
/**
* Download the DE442S planetary kernel and naif0012.tls leap-second kernel
* to the local cache directory. Verifies the download by SHA-256 checksum
* when a checksum is supplied via config.checksumOverride.
*
* @param config - Optional kernel config (to customize cache directory or checksum)
* @returns Paths where kernels were saved
*/
export async function downloadKernels(config?: KernelConfig): Promise<{
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')
await mkdir(cacheDir, { recursive: true })
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)`)
if (config?.checksumOverride) {
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.')
}
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.')
} else {
console.log('naif0012.tls already cached.')
}
return { planetaryPath, leapSecondsPath }
}
/** Compute the SHA-256 hex digest of a local file. */
async function sha256File(filePath: string): Promise<string> {
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')
}
/**
* Verify that locally cached kernels exist (and match checksums if supplied).
*
* @param config - Optional kernel config (to customize cache directory or checksum)
* @returns { ok, errors[] } — ok is true when all checks pass
*/
export async function verifyKernels(config?: KernelConfig): Promise<{
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 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.`)
} else if (config?.checksumOverride) {
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.`)
}
return { ok: errors.length === 0, errors }
}
// ─── Kernel resolution ─────────────────────────────────────────────────────────
/**
* Resolve a SpkKernel from the given config or module singleton.
* If neither is available, triggers auto-download.
*/
async function resolveKernel(config?: KernelConfig): Promise<SpkKernel> {
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())
}
}
if (activeKernel) return activeKernel
// auto-init as last resort
await initKernels(config)
if (!activeKernel) throw new Error('Kernel failed to initialize. Call initKernels() before computing.')
return activeKernel
}
// ─── Primary API ──────────────────────────────────────────────────────────────
/**
* Compute a complete moon sighting report for a date and location.
*
* Returns all quantities needed for Islamic crescent sighting determination:
* sunset/moonset times, best observation time, crescent geometry (ARCL, ARCV,
* DAZ, W, Lag), Yallop category, Odeh zone, and plain-language guidance.
*
* Requires initKernels() to have been called first (or pass kernels in options).
*
* @param date - Date to check (UTC midnight or any time on that civil day)
* @param observer - Observer location and environmental parameters
* @param options - Optional configuration for refraction model, best-time method, etc.
* @returns Complete MoonSightingReport
*
* @example
* ```ts
* await initKernels()
* const report = await getMoonSightingReport(new Date('2025-03-01'), {
* lat: 51.5074, lon: -0.1278, elevation: 10, name: 'London'
* })
* console.log(report.yallop.category) // 'A' through 'F'
* console.log(report.guidance) // "Best time to look: ..."
* ```
*/
export async function getMoonSightingReport(
date: Date,
observer: Observer,
options?: SightingOptions,
): Promise<MoonSightingReport> {
const kernel = await resolveKernel(options?.kernels)
// Event times (sunset, moonset, twilight, rise)
const events = eventsGetSunMoonEvents(date, observer, kernel)
const { sunsetUTC, moonsetUTC } = events
if (!sunsetUTC || !moonsetUTC) {
return buildNullReport(date, observer, events, 'DE442S', false)
}
// Best observation time
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 (!bestTimeResult) {
bestTimeResult = bestTimeHeuristic(sunsetUTC, moonsetUTC)
}
if (!bestTimeResult) {
return buildNullReport(date, observer, events, 'DE442S', false)
}
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)
// Body positions in GCRS (geocentric)
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]
// 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)
// Airless alt/az — required by Yallop/Odeh criteria
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)
// 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
// 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 geometry = computeCrescentGeometry(
moonAirless,
sunAirless,
moonTopo,
sunTopo,
sunsetUTC,
moonsetUTC,
)
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 guidance = buildGuidanceText(
yallop,
odeh,
moonApparent.azimuth,
moonApparent.altitude,
bestTimeUTC,
lagMinutes,
)
return {
date,
observer,
sunsetUTC,
moonsetUTC,
lagMinutes,
bestTimeUTC,
bestTimeWindowUTC,
moonPosition: moonApparent,
sunPosition: sunAirless,
illumination,
moonAge: moonAgeHours,
geometry,
yallop,
odeh,
guidance,
ephemerisSource: 'DE442S',
moonAboveHorizon,
sightingPossible,
}
}
/** Build a null report for cases where sighting geometry cannot be computed. */
function buildNullReport(
date: Date,
observer: Observer,
events: SunMoonEvents,
source: 'DE442S' | 'approximate',
sightingPossible: boolean,
): MoonSightingReport {
return {
date,
observer,
sunsetUTC: events.sunsetUTC,
moonsetUTC: events.moonsetUTC,
lagMinutes: null,
bestTimeUTC: null,
bestTimeWindowUTC: null,
moonPosition: null,
sunPosition: null,
illumination: null,
moonAge: null,
geometry: null,
yallop: null,
odeh: null,
guidance: 'Sighting not possible: sunset or moonset could not be determined for this date and location.',
ephemerisSource: source,
moonAboveHorizon: null,
sightingPossible,
}
}
/**
* Compute the Moon's current phase, illumination, and next phase times.
*
* Works WITHOUT a kernel (uses Meeus approximation).
*
* @param date - Date to compute phase for (default: now)
* @returns MoonPhaseResult with illumination, phase name, age, and next events
*
* @example
* ```ts
* const phase = getMoonPhase(new Date())
* console.log(phase.phase) // 'waxing-crescent'
* console.log(phase.illumination) // 14.3 (percent)
* console.log(phase.nextFullMoon) // Date object
* ```
*/
export function getMoonPhase(date = new Date()): MoonPhaseResult {
const ts = computeTimeScales(date)
const { moonGCRS, sunGCRS } = getMoonSunApproximate(ts.jdTT)
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 phaseName = elongationToPhase(elongationDeg, isWaxing)
const nextNewMoonJD = nearestNewMoon(ts.jdTT + 15)
const nextFullMoonJD = nearestFullMoon(ts.jdTT)
return {
phase: phaseName,
illumination: illuminationPct,
age,
elongationDeg,
isWaxing,
nextNewMoon: jdToJSDate(nextNewMoonJD),
nextFullMoon: jdToJSDate(nextFullMoonJD),
prevNewMoon: jdToJSDate(prevNewMoonJD),
}
}
// ─── Internal helpers ─────────────────────────────────────────────────────────
/** Convert JD to a UTC Date. */
function jdToJSDate(jd: number): Date {
return new Date((jd - 2440587.5) * 86400000)
}
/**
* Approximate the nearest full moon JD using Meeus Ch. 49 (full moon k = n + 0.5).
* 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)
// 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
}
/** 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 DEG = Math.PI / 180
let JDE = 2451550.09766
+ 29.530588861 * k
+ 0.00015437 * T * T
- 0.000000150 * T * T * T
+ 0.00000000073 * T * T * T * T
const M = (2.5534 + 29.10535670 * k - 0.0000014 * T * T) * DEG
const Mp = (201.5643 + 385.81693528 * k + 0.0107582 * T * T) * DEG
const Fc = (160.7108 + 390.67050284 * k - 0.0016118 * T * T) * DEG
const Om = (124.7746 - 1.56375588 * k + 0.0020672 * T * T) * DEG
const E = 1 - 0.002516 * T - 0.0000074 * T * T
JDE +=
-0.40614 * Math.sin(Mp)
+ 0.17302 * E * Math.sin(M)
+ 0.01614 * Math.sin(2 * Mp)
+ 0.01043 * Math.sin(2 * Fc)
+ 0.00734 * E * Math.sin(Mp - M)
- 0.00515 * E * Math.sin(Mp + M)
+ 0.00209 * E * E * Math.sin(2 * M)
- 0.00111 * Math.sin(Mp - 2 * Fc)
- 0.00057 * Math.sin(Mp + 2 * Fc)
+ 0.00056 * E * Math.sin(2 * Mp + M)
- 0.00042 * Math.sin(3 * Mp)
+ 0.00042 * E * Math.sin(M + 2 * Fc)
+ 0.00038 * E * Math.sin(M - 2 * Fc)
- 0.00024 * E * Math.sin(2 * Mp - M)
- 0.00017 * Math.sin(Om)
- 0.00007 * Math.sin(Mp + 2 * M)
+ 0.00004 * Math.sin(2 * Mp - 2 * Fc)
+ 0.00004 * Math.sin(3 * M)
+ 0.00003 * Math.sin(Mp + M - 2 * Fc)
+ 0.00003 * Math.sin(2 * Mp + 2 * Fc)
- 0.00003 * Math.sin(Mp + M + 2 * Fc)
+ 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)
return JDE
}
/**
* Map elongation and waxing direction to a named phase.
* 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'
}
/**
* Get rise, set, and twilight times for the Sun and Moon on a given date.
*
* Requires initKernels() for accurate results.
*
* @param date - Date to compute events for
* @param observer - Observer location
* @param options - Optional kernel configuration
* @returns SunMoonEvents with all times in UTC
*/
export async function getSunMoonEvents(
date: Date,
observer: Observer,
options?: Pick<SightingOptions, 'kernels'>,
): Promise<SunMoonEvents> {
const kernel = await resolveKernel(options?.kernels)
return eventsGetSunMoonEvents(date, observer, kernel)
}

402
src/bodies/index.ts Normal file
View file

@ -0,0 +1,402 @@
/**
* bodies — Moon and Sun state computation and illumination geometry.
*
* This module assembles the full pipeline from ephemeris evaluation to
* geocentric and topocentric body positions. It uses the SpkKernel's
* segment-chaining getState() to retrieve Moon/Sun positions relative to Earth.
*
* When the DE442S kernel is not available, getMoonSunApproximate() provides
* low-accuracy positions using Meeus Ch. 25 (Sun) and Ch. 47 (Moon).
*
* References:
* Meeus, J. (1998). Astronomical Algorithms, 2nd ed. Willmann-Bell.
* Odeh, M. (2006). New Criterion for Lunar Crescent Visibility.
* Experimental Astronomy, 17(1-3), 117-138.
* Yallop, B.D. (1997). A Method for Predicting the First Sighting of the
* 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'
// ─── Constants ────────────────────────────────────────────────────────────────
const DEG = Math.PI / 180
const AU_KM = 149597870.7
/** Mean radius of the Moon in km (IAU 2015 nominal value) */
const MOON_RADIUS_KM = 1737.4
/** Mean radius of the Sun in km */
const SUN_RADIUS_KM = 696000.0
// ─── Geocentric state ─────────────────────────────────────────────────────────
/**
* Compute the geocentric state of the Moon in GCRS at the given ET.
*
* The SpkKernel handles DE442S segment chaining automatically:
* Moon→Earth = Moon→EMB − Earth→EMB
*
* @param kernel - Loaded DE442S kernel
* @param et - Ephemeris time, seconds past J2000 TDB
* @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)
}
/**
* Compute the geocentric state of the Sun in GCRS at the given ET.
*
* The SpkKernel handles DE442S segment chaining automatically:
* Sun→Earth = Sun→SSB − (EMB→SSB − Earth→EMB)
*
* @param kernel - Loaded DE442S kernel
* @param et - Ephemeris time, seconds past J2000 TDB
* @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)
}
// ─── Moon illumination ────────────────────────────────────────────────────────
/**
* Compute the Moon's illumination fraction and related phase quantities.
*
* Phase angle i = angle at Moon between Earth and Sun directions.
* Illumination k = (1 + cos(i)) / 2 — 0 at new moon, 1 at full moon.
* Elongation ψ = angle at Earth between Moon and Sun.
*
* Since the Sun is ~400× farther from Earth than the Moon:
* i ≈ π − ψ (they are approximately supplementary)
* cos(i) ≈ −cos(ψ) → k = (1 − cos(ψ)) / 2
*
* Reference: Meeus §48.
*
* @param moonGCRS - Moon geocentric position (km)
* @param sunGCRS - Sun geocentric position (km)
* @returns illumination [0-1], phaseAngleDeg, elongationDeg, isWaxing
*/
export function computeIllumination(
moonGCRS: Vec3,
sunGCRS: Vec3,
): { illumination: number; phaseAngleDeg: number; elongationDeg: number; isWaxing: boolean } {
const dot = (a: Vec3, b: Vec3) => a[0]*b[0] + a[1]*b[1] + a[2]*b[2]
const norm = (v: Vec3) => Math.sqrt(dot(v, v))
const rMoon = norm(moonGCRS)
const rSun = norm(sunGCRS)
// Elongation ψ: angle at Earth between Moon and Sun
const cosElong = dot(moonGCRS, sunGCRS) / (rMoon * rSun)
const elongationDeg = Math.acos(Math.max(-1, Math.min(1, cosElong))) / DEG
// Phase angle i: angle at Moon between Earth and Sun
// Vector from Moon to Earth: -moonGCRS
// Vector from Moon to Sun: sunGCRS - moonGCRS
const moonToSun: Vec3 = [
sunGCRS[0] - moonGCRS[0],
sunGCRS[1] - moonGCRS[1],
sunGCRS[2] - moonGCRS[2],
]
const moonToEarth: Vec3 = [-moonGCRS[0], -moonGCRS[1], -moonGCRS[2]]
const rMoonToSun = norm(moonToSun)
const cosPhase = dot(moonToEarth, moonToSun) / (rMoon * rMoonToSun)
const phaseAngleDeg = Math.acos(Math.max(-1, Math.min(1, cosPhase))) / DEG
const illumination = (1 + Math.cos(phaseAngleDeg * DEG)) / 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
return { illumination, phaseAngleDeg, elongationDeg, isWaxing }
}
/**
* Compute the topocentric crescent width W (Odeh) and W' (Yallop) in arc minutes.
*
* The crescent width is the angular thickness of the illuminated limb at its widest,
* measured in the direction perpendicular to the cusps axis.
*
* For a sphere of topocentric semi-diameter SD and topocentric elongation ARCL:
* W = SD × (1 − cos ARCL) [SD and W in the same angular units]
*
* This is exact for a spherical Moon model and gives:
* W = 0 at new moon (ARCL = 0°) — correct, no crescent
* W = SD at ARCL = 90°
* W = 2·SD at full moon (ARCL = 180°)
*
* Both Odeh W and Yallop W' use this formula with topocentric ARCL and SD.
*
* @param moonTopoVec - Topocentric Moon position vector (km)
* @param ARCL - Topocentric Sun-Moon angular separation (degrees)
* @returns W and Wprime: topocentric crescent width in arc minutes
*/
export function computeCrescentWidth(
moonTopoVec: Vec3,
ARCL: number,
): { W: number; Wprime: number } {
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) / DEG * 60
// Crescent width in arc minutes
const ARCL_rad = ARCL * DEG
const W = SDmoon_arcmin * (1 - Math.cos(ARCL_rad))
// Wprime ≡ W for both Odeh and Yallop in this formulation
return { W, Wprime: W }
}
// ─── Approximate positions (no kernel) ────────────────────────────────────────
/**
* Low-accuracy Sun and Moon positions using Meeus Ch. 25 (Sun) and Ch. 47 (Moon).
*
* Error budget:
* Sun: < 0.01° in ecliptic longitude (main terms only)
* Moon: < 0.3° in ecliptic longitude, < 0.2° in latitude
*
* Not suitable for crescent sighting reports. Intended for moon phase displays
* and for bootstrapping event-time searches.
*
* @param jdTT - Julian Date in TT
* @returns Geocentric GCRS positions in km (approximate, light-time not corrected)
*/
export function getMoonSunApproximate(jdTT: number): {
moonGCRS: Vec3
sunGCRS: Vec3
} {
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) * DEG
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)
// True longitude and anomaly
const sunLonDeg = L0 + C
const nu_rad = M_sun_rad + C * DEG
// 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
// Nutation correction for apparent longitude (simplified)
const omega = (125.04 - 1934.136 * T) * DEG
const sunLonApp = sunLonDeg - 0.00569 - 0.00478 * Math.sin(omega)
const sunLon_rad = sunLonApp * DEG
// 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) * DEG
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) ─────────────────────────────────────────────────────
// Fundamental arguments (degrees)
const Lp = 218.3164477 + 481267.88123421 * T - 0.0015786 * T * T + T * T * T / 538841 - 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
const Mp = 134.9633964 + 477198.8675055 * T + 0.0087414 * T * T + T * T * T / 69699 - T * T * T * T / 14712000
const F = 93.2720950 + 483202.0175233 * T - 0.0036539 * T * T - T * T * T / 3526000 + T * T * T * T / 863310000
// Additive terms for longitude/latitude
const A1 = (119.75 + 131.849 * T) * DEG
const A2 = ( 53.09 + 479264.290 * T) * DEG
const A3 = (313.45 + 481266.484 * T) * DEG
// Convert to radians for accumulation
const D_r = (D % 360) * DEG
const M_r = (M % 360) * DEG
const Mp_r = (Mp % 360) * DEG
const F_r = (F % 360) * DEG
// Eccentricity correction for terms involving M (Earth's orbital eccentricity)
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)]
const LD: ReadonlyArray<readonly [number,number,number,number,number,number]> = [
[ 0, 0, 1, 0, 6288774, -20905355],
[ 2, 0,-1, 0, 1274027, -3699111],
[ 2, 0, 0, 0, 658314, -2955968],
[ 0, 0, 2, 0, 213618, -569925],
[ 0, 1, 0, 0, -185116, 48888],
[ 0, 0, 0, 2, -114332, -3149],
[ 2, 0,-2, 0, 58793, 246158],
[ 2,-1,-1, 0, 57066, -152138],
[ 2, 0, 1, 0, 53322, -170733],
[ 2,-1, 0, 0, 45758, -204586],
[ 0, 1,-1, 0, -40923, -129620],
[ 1, 0, 0, 0, -34720, 108743],
[ 0, 1, 1, 0, -30383, 104755],
[ 2, 0, 0,-2, 15327, 10321],
[ 0, 0, 1, 2, -12528, 0],
[ 0, 0, 1,-2, 10980, 79661],
[ 4, 0,-1, 0, 10675, -34782],
[ 0, 0, 3, 0, 10034, -23210],
[ 4, 0,-2, 0, 8548, -21636],
[ 2, 1,-1, 0, -7888, 24208],
[ 2, 1, 0, 0, -6766, 30824],
[ 1, 0,-1, 0, -5163, -8379],
[ 1, 1, 0, 0, 4987, -16675],
[ 2,-1, 1, 0, 4036, -12831],
[ 2, 0, 2, 0, 3994, -10445],
[ 4, 0, 0, 0, 3861, -11650],
[ 2, 0,-3, 0, 3665, 14403],
[ 0, 1,-2, 0, -2689, -7003],
[ 2, 0,-1, 2, -2602, 0],
[ 2,-1,-2, 0, 2390, 10056],
]
let Sl = 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)
}
// Additive longitude corrections (Meeus §47)
Sl += 3958 * Math.sin(A1) + 1962 * Math.sin((Lp - F) * DEG) + 318 * Math.sin(A2)
// Latitude accumulation — 20 main terms from Meeus Table 47.B
// [d, m, mp, f, Σb (0.000001°)]
const FB: ReadonlyArray<readonly [number,number,number,number,number]> = [
[ 0, 0, 0, 1, 5128122],
[ 0, 0, 1, 1, 280602],
[ 0, 0, 1,-1, 277693],
[ 2, 0, 0,-1, 173237],
[ 2, 0,-1, 1, 55413],
[ 2, 0,-1,-1, 46271],
[ 2, 0, 0, 1, 32573],
[ 0, 0, 2, 1, 17198],
[ 2, 0, 1,-1, 9266],
[ 0, 0, 2,-1, 8822],
[ 2,-1, 0,-1, 8216],
[ 2, 0,-2,-1, 4324],
[ 2, 0, 1, 1, 4200],
[ 2, 1, 0,-1, -3359],
[ 2,-1,-1, 1, 2463],
[ 2,-1, 0, 1, 2211],
[ 2,-1,-1,-1, 2065],
[ 0, 1,-1,-1, -1870],
[ 4, 0,-1,-1, 1828],
[ 0, 1, 0, 1, -1794],
]
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)
}
// Additive latitude corrections
Sb += -2235 * Math.sin(Lp * DEG) + 382 * Math.sin(A3) + 175 * Math.sin(A1 - F_r)
+ 175 * Math.sin(A1 + F_r) + 127 * Math.sin((Lp - Mp) * DEG) - 115 * Math.sin((Lp + Mp) * DEG)
// Moon ecliptic coordinates
const moonLonDeg = Lp + Sl * 1e-6
const moonLatDeg = Sb * 1e-6
const moonDistKm = 385000.56 + Sr * 0.001
const moonLon_rad = moonLonDeg * DEG
const moonLat_rad = moonLatDeg * DEG
// Ecliptic to equatorial (GCRS ≈ J2000 equatorial for this accuracy level)
const moonGCRS: Vec3 = [
moonDistKm * Math.cos(moonLat_rad) * Math.cos(moonLon_rad),
moonDistKm * (Math.cos(eps) * Math.cos(moonLat_rad) * Math.sin(moonLon_rad) - Math.sin(eps) * Math.sin(moonLat_rad)),
moonDistKm * (Math.sin(eps) * Math.cos(moonLat_rad) * Math.sin(moonLon_rad) + Math.cos(eps) * Math.sin(moonLat_rad)),
]
return { moonGCRS, sunGCRS }
}
/**
* Estimate the time of the nearest new moon using Meeus Ch. 49.
* Accurate to within ~2 hours; sufficient for phase age calculations.
*
* @param jdTT - Julian Date in TT near the desired new moon
* @returns Julian Date in TT of the nearest new moon
*/
export function nearestNewMoon(jdTT: number): number {
// Convert JD to approximate decimal year
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
// JDE of mean new moon (Meeus Eq. 49.1)
let JDE = 2451550.09766
+ 29.530588861 * k
+ 0.00015437 * T * T
- 0.000000150 * T * T * T
+ 0.00000000073 * T * T * T * T
// Fundamental arguments for the corrections (degrees → radians)
const M = (2.5534 + 29.10535670 * k - 0.0000014 * T * T - 0.00000011 * T * T * T) * DEG
const Mp = (201.5643 + 385.81693528 * k + 0.0107582 * T * T + 0.00001238 * T * T * T) * DEG
const Fc = (160.7108 + 390.67050284 * k - 0.0016118 * T * T - 0.00000227 * T * T * T) * DEG
const Om = (124.7746 - 1.56375588 * k + 0.0020672 * T * T + 0.00000215 * T * T * T) * DEG
// Eccentricity of Earth's orbit
const E = 1 - 0.002516 * T - 0.0000074 * T * T
// Corrections from Meeus Table 49.A (new moon)
JDE +=
-0.40720 * Math.sin(Mp)
+ 0.17241 * E * Math.sin(M)
+ 0.01608 * Math.sin(2 * Mp)
+ 0.01039 * Math.sin(2 * Fc)
+ 0.00739 * E * Math.sin(Mp - M)
- 0.00514 * E * Math.sin(Mp + M)
+ 0.00208 * E * E * Math.sin(2 * M)
- 0.00111 * Math.sin(Mp - 2 * Fc)
- 0.00057 * Math.sin(Mp + 2 * Fc)
+ 0.00056 * E * Math.sin(2 * Mp + M)
- 0.00042 * Math.sin(3 * Mp)
+ 0.00042 * E * Math.sin(M + 2 * Fc)
+ 0.00038 * E * Math.sin(M - 2 * Fc)
- 0.00024 * E * Math.sin(2 * Mp - M)
- 0.00017 * Math.sin(Om)
- 0.00007 * Math.sin(Mp + 2 * M)
+ 0.00004 * Math.sin(2 * Mp - 2 * Fc)
+ 0.00004 * Math.sin(3 * M)
+ 0.00003 * Math.sin(Mp + M - 2 * Fc)
+ 0.00003 * Math.sin(2 * Mp + 2 * Fc)
- 0.00003 * Math.sin(Mp + M + 2 * Fc)
+ 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)
return JDE
}

179
src/cli/index.ts Normal file
View file

@ -0,0 +1,179 @@
/**
* moon-calc CLI
*
* Commands:
* moon-calc download-kernels Download DE442S and naif0012.tls to cache
* moon-calc verify-kernels Verify cached kernels by SHA-256 checksum
* moon-calc sighting <lat> <lon> [date] Print a sighting report
* moon-calc phase [date] Print current moon phase
* moon-calc benchmark Run performance benchmark
*/
import {
initKernels,
downloadKernels,
verifyKernels,
getMoonSightingReport,
getMoonPhase,
} from '../api/index.js'
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
default:
printHelp()
process.exit(command ? 1 : 0)
}
}
function printHelp() {
console.log(`moon-calc — Lunar crescent visibility calculator
Commands:
download-kernels Download DE442S and naif0012.tls to cache
verify-kernels Verify cached kernels by SHA-256
sighting <lat> <lon> [date] Print sighting report (date: YYYY-MM-DD, default today)
phase [date] Print moon phase (date: YYYY-MM-DD, default today)
benchmark Run performance benchmark
Examples:
moon-calc download-kernels
moon-calc sighting 51.5 -0.1 2025-03-29
moon-calc sighting 21.4 39.8 # Mecca
moon-calc phase 2025-03-01`)
}
async function cmdDownloadKernels() {
await downloadKernels()
console.log('Kernels ready.')
}
async function cmdVerifyKernels() {
const result = await verifyKernels()
if (result.ok) {
console.log('Kernels OK.')
} else {
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)
if (isNaN(lat) || isNaN(lon)) {
console.error('Usage: moon-calc sighting <lat> <lon> [YYYY-MM-DD]')
process.exit(1)
}
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.log(`Computing sighting report for ${lat}°N ${lon}°E on ${dateStr}...`)
await initKernels()
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(`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`)
}
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(report.guidance)
}
function cmdPhase(dateStr?: string) {
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)
}
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`)
}
async function cmdBenchmark() {
console.log('moon-calc benchmark\n')
// Benchmark 1: getMoonPhase (no kernel needed)
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))))
}
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`)
// 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`)
}
/** Format a nullable Date as a short UTC string. */
function fmtDate(d: Date | null): string {
if (!d) return 'N/A'
return d.toISOString().replace('T', ' ').replace(/\.\d+Z$/, ' UTC')
}
main().catch(err => {
console.error(err instanceof Error ? err.message : String(err))
process.exit(1)
})

359
src/events/index.ts Normal file
View file

@ -0,0 +1,359 @@
/**
* events — Rise/set, twilight, and best-observation-time computation.
*
* Finding when the Sun and Moon cross the horizon is a root-finding problem:
*
* f(t) = apparent_altitude(t) − h0 = 0
*
* Where h0 is the threshold altitude (accounting for refraction, semi-diameter,
* and parallax). For sunset, the standard threshold is −0.8333° (34' refraction
* + 16' solar semi-diameter ≈ −50' = −0.8333°).
*
* The solver brackets the crossing by sampling altitude at coarse time steps,
* then applies Brent's method to each sign-change bracket for precision.
*
* Best-time computation:
* Yallop/Odeh approximation: T_b = T_sunset + (4/9) × Lag
* where Lag = T_moonset − T_sunset in minutes
*
* 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 } from '../math/index.js'
import {
J2000,
SECONDS_PER_DAY,
computeTimeScales,
jdToDate,
dateToJD,
jdTTtoET,
getDeltaAT,
TT_MINUS_TAI,
} 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'
// ─── Altitude threshold constants ─────────────────────────────────────────────
/**
* Standard threshold altitude for sunset/sunrise.
* Accounts for: standard refraction at horizon (34') + solar semi-diameter (16')
* Total: −50' = −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
// ─── Internal helpers ─────────────────────────────────────────────────────────
/**
* Convert ET (seconds past J2000 TDB) to approximate TimeScales.
* Accurate to within ~1 second for UTC — sufficient for event finding.
*/
function etToTS(et: number): TimeScales {
// ET ≈ (jdTT - J2000) × 86400, so jdTT ≈ J2000 + et / 86400
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)
}
/**
* Get the geocentric GCRS position of a body at the given ET.
*/
function bodyPositionAt(kernel: SpkKernel, naifId: number, et: number): Vec3 {
if (naifId === NAIF_IDS.SUN) {
return getSunGeocentricState(kernel, et).position
}
return getMoonGeocentricState(kernel, et).position
}
/**
* Compute the airless altitude of a body at a given ET.
* Returns altitude minus threshold so the zero crossing equals the event time.
*/
function altitudeMinusThreshold(
kernel: SpkKernel,
naifId: number,
observer: Observer,
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
}
// ─── Event finding ────────────────────────────────────────────────────────────
/**
* Find the time when a body (Sun or Moon) crosses a given altitude threshold.
* Returns the first crossing in the requested direction within the search window.
*
* Algorithm: coarse sample (10-min steps) to bracket sign changes, then Brent.
*
* @param kernel - DE442S kernel
* @param naifId - NAIF body ID (10=Sun, 301=Moon)
* @param observer - Observer location
* @param ts - Time scales at the search start (used only for context; not interpolated)
* @param startET - Search start (ET seconds past J2000)
* @param endET - Search end (ET seconds past J2000)
* @param threshold - Altitude threshold in degrees
* @param rising - True to find a rising event, false for setting
* @returns Event time as Date (UTC), or null if no crossing found
*/
export function findAltitudeCrossing(
kernel: SpkKernel,
naifId: number,
observer: Observer,
ts: TimeScales,
startET: number,
endET: number,
threshold: number,
rising: boolean,
): Date | null {
void ts // ts not needed — etToTS computes from ET directly
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)
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 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
if (etCross !== null) {
const tsCross = etToTS(etCross)
return tsCross.utc
}
}
prevET = currET
prevF = currF
}
return null
}
/**
* Find all Sun and Moon rise/set and twilight times for a given date and observer.
*
* Searches a 28-hour window from the start of the UTC date to cover events
* at all longitudes for the corresponding civil date.
*
* @param date - Civil date (any time on the desired UTC day)
* @param observer - Observer location
* @param kernel - DE442S kernel
* @returns SunMoonEvents with all event times in UTC, or null if event doesn't occur
*/
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)
// 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 ts0 = computeTimeScales(midnight)
// Sun events
const sunriseUTC = findAltitudeCrossing(
kernel, NAIF_IDS.SUN, observer, ts0, etStart, etEnd, SUN_ALTITUDE_THRESHOLD, true,
)
const sunsetUTC = findAltitudeCrossing(
kernel, NAIF_IDS.SUN, observer, ts0, etStart, etEnd, SUN_ALTITUDE_THRESHOLD, false,
)
// Twilight events (Sun setting below -6°, -12°, -18°)
const civilTwilightEndUTC = findAltitudeCrossing(
kernel, NAIF_IDS.SUN, observer, ts0, etStart, etEnd, -6, false,
)
const nauticalTwilightEndUTC = findAltitudeCrossing(
kernel, NAIF_IDS.SUN, observer, ts0, etStart, etEnd, -12, false,
)
const astronomicalTwilightEndUTC = findAltitudeCrossing(
kernel, NAIF_IDS.SUN, observer, ts0, etStart, etEnd, -18, false,
)
// Moon events
const moonriseUTC = findAltitudeCrossing(
kernel, NAIF_IDS.MOON, observer, ts0, etStart, etEnd, MOON_ALTITUDE_THRESHOLD, true,
)
const moonsetUTC = findAltitudeCrossing(
kernel, NAIF_IDS.MOON, observer, ts0, etStart, etEnd, MOON_ALTITUDE_THRESHOLD, false,
)
return {
sunriseUTC,
sunsetUTC,
moonriseUTC,
moonsetUTC,
civilTwilightEndUTC,
nauticalTwilightEndUTC,
astronomicalTwilightEndUTC,
}
}
// ─── Best-time computation ────────────────────────────────────────────────────
/**
* Compute the best observation time using the Yallop/Odeh heuristic.
*
* T_b = T_sunset + (4/9) × Lag
*
* This gives approximately 40 minutes after sunset for a typical 90-minute lag,
* which is the empirically optimal time for crescent visibility (Yallop 1997 §2.4,
* confirmed by Odeh's ICOP observation database analysis).
*
* @param sunsetUTC - UTC time of sunset
* @param moonsetUTC - UTC time of moonset
* @returns Best observation time (UTC) and lag in minutes, or null if lag ≤ 0
*/
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 lagMinutes = lagMs / 60000
const bestTimeMs = sunsetUTC.getTime() + (4 / 9) * lagMs
return {
bestTimeUTC: new Date(bestTimeMs),
lagMinutes,
}
}
/**
* Odeh arcv minimum polynomial (Odeh 2006, Eq. 1).
* Returns the minimum ARCV needed for visibility at crescent width W (arc minutes).
*/
function odehArcvMin(W: number): number {
return 11.8371 - 6.3226 * W + 0.7319 * W * W - 0.1018 * W * W * W
}
/**
* Find the optimal observation time by maximizing the Odeh V parameter
* over the interval [sunset, moonset].
*
* More accurate than the heuristic at high latitudes or unusual geometries
* where the (4/9)×Lag formula can be misleading.
*
* @param sunsetUTC - UTC time of sunset
* @param moonsetUTC - UTC time of moonset
* @param kernel - DE442S kernel
* @param observer - Observer location
* @param steps - Number of scan steps (default 90, giving ~1-min resolution for typical lag)
* @returns Optimal time, its V parameter, and lag in minutes
*/
export function bestTimeOptimized(
sunsetUTC: Date,
moonsetUTC: Date,
kernel: SpkKernel,
observer: Observer,
steps = 90,
): { bestTimeUTC: Date; lagMinutes: number; maxV: number } | null {
const lagMs = moonsetUTC.getTime() - sunsetUTC.getTime()
if (lagMs <= 0) return null
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 dot = (a: Vec3, b: Vec3) => a[0]*b[0] + a[1]*b[1] + a[2]*b[2]
const norm = (v: Vec3) => Math.sqrt(dot(v, v))
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 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)
// Airless altitudes via the full pipeline
const moonAzAlt = computeAzAlt(moonGCRS, observer, ts, true)
const sunAzAlt = computeAzAlt(sunGCRS, observer, ts, true)
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 = dot(moonTopo, sunTopo) / (norm(moonTopo) * norm(sunTopo))
const ARCL = Math.acos(Math.max(-1, Math.min(1, cosARCL))) * (180 / Math.PI)
const { W } = computeCrescentWidth(moonTopo, ARCL)
const V = ARCV - odehArcvMin(W)
if (V > maxV) {
maxV = V
bestTimeUTC = t
}
}
return { bestTimeUTC, lagMinutes, maxV }
}
/**
* Compute a time window around the best time where the crescent may be visible.
* Default window: ±20 minutes around best time (practical approximation).
*
* @param bestTimeUTC - Best observation time
* @param windowMinutes - Half-width of window in minutes (default 20)
* @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),
]
}

466
src/frames/index.ts Normal file
View file

@ -0,0 +1,466 @@
/**
* frames — GCRS to ITRS coordinate transformation (IERS Q·R·W chain).
*
* To compute topocentric alt/az from inertial (GCRS/ICRF) coordinates, we apply
* the IAU/IERS reference frame transformation:
*
* [ITRS] = W(t) · R(t) · Q(t) · [GCRS]
*
* Where:
* Q(t) — Celestial motion: precession + nutation (IAU 2006/2000B)
* R(t) — Earth rotation: Earth Rotation Angle (ERA) from UT1
* W(t) — Polar motion: small offsets xp, yp from the CIP (typically < 1 arcsec)
*
* The CIP (Celestial Intermediate Pole) coordinates X, Y and the CIO locator s
* parameterize Q(t). These are computed from IAU 2006 precession polynomials plus
* the IAU 2000B 77-term nutation series (< 1 mas error, far below our 30" target).
*
* References:
* IERS Conventions (2010), Chapter 5
* IAU SOFA release 2023-10-11 (iauNut00b, iauEra00, iauPom00, iauC2ixys)
* 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'
// ─── Constants ────────────────────────────────────────────────────────────────
/** Arcseconds to radians */
const ARCSEC_RAD = Math.PI / (180 * 3600)
/** 0.1 microarcseconds to arcseconds (units of nutation coefficients) */
const UAS01_TO_ARCSEC = 1e-7
// ─── IAU 2000B Nutation Series ────────────────────────────────────────────────
//
// 77-term luni-solar nutation series from SOFA iauNut00b.c.
// Source: Mathews et al. (2002), J. Geophys. Res. 107(B4); SOFA release 2023-10-11.
//
// Format per row: [nl, nlp, nF, nD, nOm, ps, pst, pc, ec, ect, es]
// nl..nOm — integer multipliers for Delaunay arguments l, l', F, D, Ω
// ps, pst — sin(arg) coefficient for dpsi, and its T rate [0.1 uas, 0.1 uas/cy]
// pc — cos(arg) coefficient for dpsi [0.1 uas]
// ec, ect — cos(arg) coefficient for deps, and its T rate [0.1 uas, 0.1 uas/cy]
// es — sin(arg) coefficient for deps [0.1 uas]
//
// Accumulation formulas:
// dpsi += (ps + pst*T)*sin(arg) + pc*cos(arg)
// deps += (ec + ect*T)*cos(arg) + es*sin(arg)
const NUT_2000B: ReadonlyArray<readonly [
number,number,number,number,number,
number,number,number,
number,number,number
]> = [
// 1
[ 0, 0, 0, 0, 1, -172064161.0, -174666.0, 33386.0, 92052331.0, 9086.0, 15377.0],
// 2
[ 0, 0, 2,-2, 2, -13170906.0, -13696.0, -13238.0, 5730336.0, -3015.0, -4587.0],
// 3
[ 0, 0, 2, 0, 2, -2276413.0, -2353.0, 2796.0, 978459.0, -619.0, 645.0],
// 4
[ 0, 0, 0, 0, 2, 2074554.0, 2352.0, -2635.0, -897492.0, 307.0, -187.0],
// 5
[ 0, 1, 0, 0, 0, 1475877.0, -11817.0, 11817.0, 73871.0, -184.0, -1924.0],
// 6
[ 0, 1, 2,-2, 2, -516821.0, 1226.0, -524.0, 224386.0, -677.0, -174.0],
// 7
[ 1, 0, 0, 0, 0, 711159.0, 73.0, -872.0, -6750.0, 0.0, 358.0],
// 8
[ 0, 0, 2, 0, 1, -387298.0, -367.0, 380.0, 200728.0, 18.0, 318.0],
// 9
[ 1, 0, 2, 0, 2, -301461.0, -36.0, 816.0, 129025.0, -63.0, 367.0],
// 10
[ 0,-1, 2,-2, 2, 215829.0, -494.0, 111.0, -95929.0, 299.0, 132.0],
// 11
[ 0, 0, 2,-2, 1, 128227.0, 137.0, 181.0, -68982.0, -9.0, 39.0],
// 12
[-1, 0, 2, 0, 2, 123457.0, 11.0, 19.0, -53311.0, 32.0, -4.0],
// 13
[-1, 0, 0, 2, 0, 156994.0, 10.0, -168.0, -1235.0, 0.0, 82.0],
// 14
[ 1, 0, 0, 0, 1, 63110.0, 63.0, 27.0, -33228.0, 0.0, -9.0],
// 15
[-1, 0, 0, 0, 1, -57976.0, -63.0, -189.0, 31429.0, 0.0, -75.0],
// 16
[-1, 0, 2, 2, 2, -59641.0, -11.0, 149.0, 25543.0, -11.0, 66.0],
// 17
[ 1, 0, 2, 0, 1, -51613.0, -42.0, 129.0, 26366.0, 0.0, 78.0],
// 18
[-2, 0, 2, 0, 1, 45893.0, 50.0, 31.0, -24236.0, -10.0, 20.0],
// 19
[ 0, 0, 0, 2, 0, 63384.0, 11.0, -150.0, -1220.0, 0.0, 29.0],
// 20
[ 0, 0, 2, 2, 2, -38571.0, -1.0, 158.0, 16452.0, -11.0, 68.0],
// 21
[ 0,-2, 2,-2, 2, 32481.0, 0.0, 0.0, -13870.0, 0.0, 0.0],
// 22
[-2, 0, 0, 2, 0, -47722.0, 0.0, -18.0, 477.0, 0.0, -25.0],
// 23
[ 2, 0, 2, 0, 2, -31046.0, -1.0, 131.0, 13238.0, -11.0, 59.0],
// 24
[ 1, 0, 2,-2, 2, 28593.0, 0.0, -1.0, -12338.0, 10.0, -3.0],
// 25
[-1, 0, 2, 0, 1, 20441.0, 21.0, 10.0, -10758.0, 0.0, -3.0],
// 26
[ 2, 0, 0, 0, 0, 29243.0, 0.0, -74.0, -609.0, 0.0, 13.0],
// 27
[ 0, 0, 2, 0, 0, 25887.0, 0.0, -66.0, -550.0, 0.0, 11.0],
// 28
[ 0, 1, 0, 0, 1, -14053.0, -25.0, 79.0, 8551.0, -2.0, -45.0],
// 29
[-1, 0, 0, 2, 1, 15164.0, 10.0, 11.0, -8001.0, 0.0, -1.0],
// 30
[ 0, 2, 2,-2, 2, -15794.0, 72.0, -16.0, 6850.0, -42.0, -5.0],
// 31
[ 0, 0,-2, 2, 0, 21783.0, 0.0, 13.0, -167.0, 0.0, 13.0],
// 32
[ 1, 0, 0,-2, 1, -12873.0, -10.0, -37.0, 6953.0, 0.0, -14.0],
// 33
[ 0,-1, 0, 0, 1, -12654.0, 11.0, 63.0, 6415.0, 0.0, 26.0],
// 34
[-1, 0, 2, 2, 1, -10204.0, 0.0, 25.0, 5222.0, 0.0, 15.0],
// 35
[ 0, 2, 0, 0, 0, 16707.0, -85.0, -10.0, 168.0, -1.0, 10.0],
// 36
[ 1, 0, 2, 2, 2, -7691.0, 0.0, 44.0, 3268.0, 0.0, 19.0],
// 37
[-2, 0, 2, 0, 0, -11024.0, 0.0, -14.0, 104.0, 0.0, 2.0],
// 38
[ 0, 1, 2, 0, 2, 7566.0, -21.0, -11.0, -3250.0, 0.0, -5.0],
// 39
[ 0, 0, 2, 2, 1, -6637.0, -11.0, 25.0, 3353.0, 0.0, 14.0],
// 40
[ 0,-1, 2, 0, 2, -7141.0, 21.0, 8.0, 3070.0, 0.0, 4.0],
// 41
[ 0, 0, 0, 2, 1, -6302.0, -11.0, 2.0, 3272.0, 0.0, 4.0],
// 42
[ 1, 0, 2,-2, 1, 5800.0, 10.0, 2.0, -3045.0, 0.0, -1.0],
// 43
[ 2, 0, 2,-2, 2, 6443.0, 0.0, -7.0, -2768.0, 0.0, -4.0],
// 44
[-2, 0, 0, 2, 1, -5774.0, -11.0, -15.0, 3041.0, 0.0, -5.0],
// 45
[ 2, 0, 2, 0, 1, -5350.0, 0.0, 21.0, 2695.0, 0.0, 12.0],
// 46
[ 0,-1, 2,-2, 1, -4752.0, -11.0, -3.0, 2719.0, 0.0, -3.0],
// 47
[ 0, 0, 0,-2, 1, -4940.0, -11.0, -21.0, 2720.0, 0.0, -9.0],
// 48
[-1,-1, 0, 2, 0, 7350.0, 0.0, -8.0, -51.0, 0.0, 4.0],
// 49
[ 2, 0, 0,-2, 1, 4065.0, 0.0, 6.0, -2206.0, 0.0, 1.0],
// 50
[ 1, 0, 0, 2, 0, 6579.0, 0.0, -24.0, -199.0, 0.0, 2.0],
// 51
[ 0, 1, 2,-2, 1, 3579.0, 0.0, 5.0, -1900.0, 0.0, 1.0],
// 52
[ 1,-1, 0, 0, 0, 4725.0, 0.0, -6.0, -41.0, 0.0, 3.0],
// 53
[-2, 0, 2, 0, 2, -3075.0, 0.0, -2.0, 1313.0, 0.0, -1.0],
// 54
[ 3, 0, 2, 0, 2, -2904.0, 0.0, 15.0, 1233.0, 0.0, 7.0],
// 55
[ 0,-1, 0, 2, 0, 4348.0, 0.0, -10.0, -81.0, 0.0, 2.0],
// 56
[ 1,-1, 2, 0, 2, -2878.0, 0.0, 8.0, 1232.0, 0.0, 4.0],
// 57
[ 0, 0, 0, 1, 0, -4230.0, 0.0, 5.0, -20.0, 0.0, -2.0],
// 58
[-1,-1, 2, 2, 2, -2819.0, 0.0, 7.0, 1207.0, 0.0, 3.0],
// 59
[-1, 0, 2, 0, 0, -4056.0, 0.0, 5.0, 40.0, 0.0, -2.0],
// 60
[ 0,-1, 2, 2, 2, -2647.0, 0.0, 11.0, 1129.0, 0.0, 5.0],
// 61
[-2, 0, 0, 0, 1, -2294.0, 0.0, -10.0, 1266.0, 0.0, -4.0],
// 62
[ 1, 1, 2, 0, 2, 2481.0, 0.0, -7.0, -1062.0, 0.0, -3.0],
// 63
[ 2, 0, 0, 0, 1, 2179.0, 0.0, -2.0, -1129.0, 0.0, -2.0],
// 64
[-1, 1, 0, 1, 0, 3276.0, 0.0, 1.0, -9.0, 0.0, 0.0],
// 65
[ 1, 1, 0, 0, 0, -3389.0, 0.0, 5.0, 35.0, 0.0, -2.0],
// 66
[ 1, 0, 2, 0, 0, 3339.0, 0.0, -13.0, -107.0, 0.0, 1.0],
// 67
[-1, 0, 2,-2, 1, -1987.0, 0.0, -6.0, 1073.0, 0.0, -2.0],
// 68
[ 1, 0, 0, 0, 2, -1981.0, 0.0, 0.0, 854.0, 0.0, 0.0],
// 69
[-1, 0, 0, 1, 0, 4026.0, 0.0, -353.0, -553.0, 0.0, -139.0],
// 70
[ 0, 0, 2, 1, 2, 1660.0, 0.0, -5.0, -710.0, 0.0, -2.0],
// 71
[-1, 0, 2, 4, 2, -1521.0, 0.0, 9.0, 647.0, 0.0, 4.0],
// 72
[-1, 1, 0, 1, 1, 1464.0, 0.0, -11.0, -527.0, 0.0, -1.0],
// 73
[ 0,-2, 2,-2, 1, -1389.0, 0.0, 3.0, 656.0, 0.0, 1.0],
// 74
[ 1,-1, 2, 2, 2, -1377.0, 0.0, 8.0, 594.0, 0.0, 4.0],
// 75
[ 3, 0, 2,-2, 2, 1371.0, 0.0, -2.0, -588.0, 0.0, -1.0],
// 76
[ 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
// T = Julian centuries from J2000.0 (TT)
/** 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
}
/** 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.00024470))))
)
}
/** 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.2090 + 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 ─────────────────────────────────────────────────────────
/**
* Compute CIP X, Y and CIO locator s using IAU 2006 precession + IAU 2000B nutation.
*
* Method:
* 1. Compute dpsi, deps from the 77-term IAU 2000B luni-solar nutation series.
* 2. Add to IAU 2006 precession polynomials for X and Y (IERS Conventions 2010 Eq. 5.16).
* 3. Compute s ≈ -X·Y/2 + polynomial term.
*
* Accuracy: < 1 mas over ±50 years from J2000 — well within the 30" library target.
*
* Implements the IAU 2000B standard (77-term series, truncated from IAU 2000A's 1306 terms).
* Accuracy: ~1 mas in X and Y, sufficient for all practical crescent visibility work.
*
* @param jdTT - Julian Date in TT
* @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
// Delaunay fundamental arguments
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
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
}
// Convert 0.1 uas → arcseconds → radians
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
const eps0 = (
84381.406
+ T * (-46.836769
+ T * (-0.0001831
+ T * ( 0.00200340
+ T * (-0.000000576
+ T * (-0.0000000434)))))
) * 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))))
// 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))))
// CIP X, Y: precession polynomial + first-order nutation correction
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
return { X, Y, s }
}
// ─── Earth Rotation Angle ────────────────────────────────────────────────────
/**
* Compute the Earth Rotation Angle (ERA) in radians.
*
* ERA is the angle between the Celestial Intermediate Origin (CIO) and the
* Terrestrial Intermediate Origin (TIO), measured in the equatorial plane.
* It replaces GMST in the IAU 2000+ Earth rotation model.
*
* ERA(UT1) = 2π(0.7790572732640 + 1.00273781191135448 · Du)
* where Du = JD(UT1) − 2451545.0
*
* 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.7790572732640 + 1.00273781191135448 * Du)
return ((era % (2 * Math.PI)) + 2 * Math.PI) % (2 * Math.PI)
}
// ─── Frame rotation matrices ──────────────────────────────────────────────────
/**
* Celestial motion matrix Q(t) from CIP X, Y, s.
* Converts GCRS to the Celestial Intermediate Reference System (CIRS).
*
* Uses the exact SOFA iauC2ixys formula:
* Q = Rz(-(e+s)) · Ry(d) · Rz(e)
* where e = atan2(Y, X) and d = asin(sqrt(X²+Y²)).
*
* 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)))
}
/**
* Earth rotation matrix R(t) from the ERA.
* Simple rotation about the CIP pole (z-axis) by ERA.
*/
export function earthRotationMatrix(era: number): Mat3 {
return rotZ(era)
}
/**
* Polar motion matrix W(t) from CIP offsets xp, yp.
* Small correction (< 0.5 arcsec) for the wobble of the Earth's pole.
*
* W = Ry(xp) · Rx(−yp)
* (TIO locator sp is omitted — < 0.001" effect over one year)
*
* For moon sighting purposes, default xp = yp = 0 introduces < 30m error
* in the observer position — negligible compared to refraction uncertainty.
*
* @param xp - Pole x-offset in radians (from IERS Bulletin A)
* @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))
}
// ─── Full transformation ──────────────────────────────────────────────────────
/**
* Transform a vector from GCRS (inertial) to ITRS (Earth-fixed).
*
* Full chain: [ITRS] = W · R · Q · [GCRS]
*
* @param gcrsVec - 3-vector in GCRS frame (km)
* @param ts - Time scales for the epoch
* @param xp - Polar motion x (radians, default 0)
* @param yp - Polar motion y (radians, default 0)
* @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)
// Apply Q first (GCRS→CIRS), then R (CIRS→TIRS), then W (TIRS→ITRS)
const combined = mmmul(W, mmmul(R, Q))
return mvmul(combined, gcrsVec)
}
/**
* Transform a vector from ITRS (Earth-fixed) to GCRS (inertial).
* Inverse of gcrsToItrs — the combined rotation matrix is orthogonal, so its
* inverse equals its transpose.
*
* @param itrsVec - 3-vector in ITRS frame (km)
* @param ts - Time scales for the epoch
* @param xp - Polar motion x (radians, default 0)
* @param yp - Polar motion y (radians, default 0)
* @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)
// [GCRS] = Qᵀ · Rᵀ · Wᵀ · [ITRS]
const combined = mmmul(mtranspose(Q), mmmul(mtranspose(R), mtranspose(W)))
return mvmul(combined, itrsVec)
}

66
src/index.ts Normal file
View file

@ -0,0 +1,66 @@
/**
* moon-calc — High-accuracy lunar crescent visibility and moon sighting calculations.
*
* Uses JPL DE442S ephemerides with full IERS Earth orientation for precise
* topocentric Sun/Moon positions. Implements the Yallop (NAO TN 69) and
* Odeh (Experimental Astronomy 2006) crescent visibility criteria.
*
* Quick start:
* import { getMoonSightingReport } from 'moon-calc'
*
* const report = await getMoonSightingReport(new Date('2025-03-01'), {
* lat: 51.5, lon: -0.1, elevation: 20, name: 'London'
* })
* console.log(report.yallop.category, report.guidance)
*/
// ─── Primary API ──────────────────────────────────────────────────────────────
export {
getMoonSightingReport,
getMoonPhase,
getSunMoonEvents,
initKernels,
downloadKernels,
verifyKernels,
} from './api/index.js'
// ─── Types ────────────────────────────────────────────────────────────────────
export type {
// Observer
Observer,
// Results
MoonSightingReport,
MoonPhaseResult,
MoonPhaseName,
SunMoonEvents,
CrescentGeometry,
YallopResult,
YallopCategory,
OdehResult,
OdehZone,
// Configuration
KernelConfig,
KernelSource,
SightingOptions,
// Time
TimeScales,
// Geometry primitives
AzAlt,
Vec3,
StateVector,
// Ephemeris internals (for advanced use)
SpkSegment,
ChebRecord,
} from './types.js'
// ─── Constants ────────────────────────────────────────────────────────────────
export {
YALLOP_THRESHOLDS,
YALLOP_DESCRIPTIONS,
ODEH_THRESHOLDS,
ODEH_DESCRIPTIONS,
WGS84,
} from './types.js'

348
src/math/index.ts Normal file
View file

@ -0,0 +1,348 @@
/**
* math — Core numerical utilities.
*
* All computation in this module is pure (no I/O, no state).
* Uses Float64Array for coefficient storage to match JS engine optimization paths.
*/
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]]
}
/** 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]]
}
/** Scale a 3-vector */
export function vscale(a: Vec3, s: number): Vec3 {
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]
}
/** Euclidean norm */
export function vnorm(a: Vec3): number {
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],
]
}
/** 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)
}
/** 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)
}
// ─── 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,
]
/** Multiply 3×3 matrix by 3-vector */
export function mvmul(m: Mat3, v: Vec3): Vec3 {
return [
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 */
export function mmmul(a: Mat3, b: Mat3): Mat3 {
return [
a[0]*b[0] + a[1]*b[3] + a[2]*b[6],
a[0]*b[1] + a[1]*b[4] + a[2]*b[7],
a[0]*b[2] + a[1]*b[5] + a[2]*b[8],
a[3]*b[0] + a[4]*b[3] + a[5]*b[6],
a[3]*b[1] + a[4]*b[4] + a[5]*b[7],
a[3]*b[2] + a[4]*b[5] + a[5]*b[8],
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]]
}
/**
* Rotation matrix around the X axis by angle θ (radians).
* 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]
}
/**
* 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]
}
/**
* 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]
}
// ─── Chebyshev polynomial evaluation ─────────────────────────────────────────
/**
* Evaluate a Chebyshev polynomial at a normalized point x ∈ [-1, 1]
* using the Clenshaw recurrence algorithm.
*
* The Clenshaw algorithm is numerically superior to Horner's method for
* Chebyshev series and is the standard approach in SPICE's SPKE02.
*
* @param coeffs - Chebyshev coefficients c[0..n] (degree n polynomial)
* @param x - Evaluation point, must be in [-1, 1]
* @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]
// Double-x for Clenshaw efficiency
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
b2 = b1
b1 = b0
}
return coeffs[0] + x * b1 - b2
}
/**
* Evaluate a Chebyshev polynomial and its derivative simultaneously.
* Uses the extended Clenshaw recurrence (more efficient than two separate evaluations).
*
* @returns [value, derivative] where derivative is with respect to the original time variable
* @param coeffs - Chebyshev coefficients
* @param x - Evaluation point in [-1, 1]
* @param radius - Half-interval width in seconds (for scaling derivative to original time)
*/
export function chebyshevEvalWithDerivative(
coeffs: Float64Array,
x: number,
radius: number,
): [number, number] {
const n = coeffs.length
if (n === 0) return [0, 0]
if (n === 1) return [coeffs[0], 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
const db0 = 2 * b1 + x2 * db1 - db2
b2 = b1; b1 = b0
db2 = db1; db1 = db0
}
const value = coeffs[0] + x * b1 - b2
const dvalue = b1 + x * db1 - db2
// Scale derivative from normalized domain back to seconds
return [value, dvalue / radius]
}
// ─── Root finding ─────────────────────────────────────────────────────────────
/**
* Find a root of f(t) in [a, b] using Brent's method.
* Requires f(a) and f(b) to have opposite signs.
*
* Brent's method combines bisection, secant, and inverse quadratic interpolation
* for guaranteed convergence with superlinear speed in practice.
*
* @param f - Function to find root of
* @param a - Left bracket
* @param b - Right bracket
* @param tol - Tolerance (default 1e-9 seconds for astronomical work)
* @param maxIter - Maximum iterations (default 64)
* @returns Root location, or null if bracket does not contain a sign change
*/
export function brentRoot(
f: (t: number) => number,
a: number,
b: number,
tol = 1e-9,
maxIter = 64,
): number | null {
let fa = f(a)
let fb = f(b)
// No sign change in bracket
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]
}
let c = a
let fc = fa
let mflag = true
let s = 0
let d = 0
for (let i = 0; i < maxIter; i++) {
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))
} else {
// Secant method
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
if (cond1 || cond2 || cond3 || cond4 || cond5) {
s = (a + b) / 2
mflag = true
} else {
mflag = false
}
const fs = f(s)
d = c
c = b
fc = fb
if (fa * fs < 0) {
b = s
fb = fs
} else {
a = s
fa = fs
}
if (Math.abs(fa) < Math.abs(fb)) {
;[a, b] = [b, a]
;[fa, fb] = [fb, fa]
}
}
return b
}
/**
* Find all roots of f(t) in [a, b] by adaptive subdivision then Brent.
* Used for finding rise/set events where multiple crossings may occur in a day.
*
* @param f - Function to find roots of
* @param a - Start of search interval
* @param b - End of search interval
* @param steps - Number of initial subdivision steps (default 48 for 30-min resolution over a day)
* @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)
for (let i = 1; i <= steps; i++) {
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)
if (root !== null) {
// Deduplicate roots that are too close together
if (roots.length === 0 || Math.abs(root - roots[roots.length - 1]) > 1e-6) {
roots.push(root)
}
}
}
tPrev = t
fPrev = ft
}
return roots
}
// ─── Angle utilities ─────────────────────────────────────────────────────────
/** Convert degrees to radians */
export const DEG2RAD = Math.PI / 180
/** Convert radians to degrees */
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
}
/** Normalize an angle in degrees to [0, 360) */
export function mod360(deg: number): number {
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
}

287
src/observer/index.ts Normal file
View file

@ -0,0 +1,287 @@
/**
* observer — WGS84 geodetic to ECEF, topocentric transforms, and refraction.
*
* Observer position chain:
* geodetic (lat, lon, elev) → ECEF (ITRS, meters) → km → GCRS (via frames/) → topocentric ENU → az/alt
*
* The WGS84 ellipsoid defines the reference for geodetic coordinates:
* a = 6378137.0 m (semi-major axis)
* 1/f = 298.257223563 (inverse flattening)
*
* Atmospheric refraction is significant near the horizon (up to ~34' at altitude 0°).
* The Bennett (1982) formula is the standard practical approximation; it accepts
* pressure and temperature for improved accuracy.
*
* References:
* WGS84 definition: NIMA TR8350.2, 3rd edition
* Bennett (1982), Astronomical Refraction for Upper and Lower Limbs
* 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'
// ─── Geodetic ↔ ECEF ─────────────────────────────────────────────────────────
/**
* Convert geodetic coordinates to ECEF (Earth-Centered, Earth-Fixed) position.
* Result is in meters, consistent with the WGS84 ellipsoid.
*
* @param lat - Geodetic latitude in degrees (north positive)
* @param lon - Longitude in degrees (east positive)
* @param elev - Height above ellipsoid in meters
* @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)
// Prime vertical radius of curvature
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,
]
}
/**
* Convert ECEF position (meters) to geodetic coordinates.
* Uses Bowring's iterative method (converges in 3 iterations for < 1 mm error).
*
* @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)
// Bowring iteration for geodetic latitude
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)
const h = p / Math.cos(lat) - N
return {
lat: (lat * 180) / Math.PI,
lon: (lon * 180) / Math.PI,
h,
}
}
// ─── Topocentric ENU ─────────────────────────────────────────────────────────
/**
* Compute the East-North-Up (ENU) basis vectors at a geodetic location.
* These unit vectors define the local horizon coordinate frame in ECEF/ITRS.
*
* East = (−sin λ, cos λ, 0)
* North = (−sin φ·cos λ, −sin φ·sin λ, cos φ)
* Up = ( cos φ·cos λ, cos φ·sin λ, sin φ)
*
* @param lat - Geodetic latitude in degrees
* @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 sinPhi = Math.sin(phi), cosPhi = Math.cos(phi)
const sinLam = Math.sin(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]
return { east, north, up }
}
/**
* Convert a vector from ECEF frame to the local ENU frame at the observer.
* The ECEF input vector represents the displacement from observer to target.
*
* @param ecefDelta - ECEF displacement vector (target - observer), any units
* @param lat - Observer geodetic latitude in degrees
* @param lon - Observer longitude in degrees
* @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)
const dot = (a: Vec3, b: Vec3) => a[0]*b[0] + a[1]*b[1] + a[2]*b[2]
return [dot(ecefDelta, east), dot(ecefDelta, north), dot(ecefDelta, up)]
}
/**
* Convert a topocentric ENU vector to azimuth and altitude.
*
* Azimuth is measured from North clockwise: 0° = N, 90° = E, 180° = S, 270° = W.
* Altitude is degrees above the horizon (negative = below).
*
* @param enu - [east, north, up] components (any consistent units)
* @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
// 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 }
}
// ─── Topocentric parallax ─────────────────────────────────────────────────────
/**
* Compute the topocentric position of a body by subtracting the observer's
* geocentric position from the geocentric body position.
*
* This is the primary parallax correction. The Moon's horizontal parallax
* (~57 arcmin) makes this non-trivial — topocentric position differs from
* geocentric by up to the Moon's angular diameter.
*
* @param bodyGCRS - Geocentric body position in km (GCRS)
* @param observerGCRS - Observer position in km (GCRS)
* @returns Topocentric body position in km (GCRS)
*/
export function topocentricPosition(bodyGCRS: Vec3, observerGCRS: Vec3): Vec3 {
return [
bodyGCRS[0] - observerGCRS[0],
bodyGCRS[1] - observerGCRS[1],
bodyGCRS[2] - observerGCRS[2],
]
}
// ─── Full pipeline: GCRS → az/alt ────────────────────────────────────────────
/**
* Compute topocentric azimuth and altitude for a body, given its GCRS position
* and the observer's location.
*
* Pipeline:
* 1. body GCRS → body ITRS (via gcrsToItrs)
* 2. observer geodetic → observer ECEF (ITRS, meters → km)
* 3. delta_ITRS = body_ITRS - observer_ITRS
* 4. delta_ITRS → ENU via projection on local basis
* 5. ENU → az/alt
* 6. Apply Bennett refraction if not airless
*
* @param bodyGCRS - Geocentric body position in km (GCRS)
* @param observer - Observer location
* @param ts - Time scales for the epoch (needed for GCRS→ITRS rotation)
* @param airless - If true, return geometric altitude without refraction
* @returns Azimuth (degrees from N, clockwise) and altitude (degrees)
*/
export function computeAzAlt(
bodyGCRS: Vec3,
observer: Observer,
ts: TimeScales,
airless: boolean,
): AzAlt {
// 1. Convert body position from GCRS to ITRS (km)
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]
// 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],
]
// 4. Project onto local ENU basis at the observer's location
const enu = ecefToENU(delta, observer.lat, observer.lon)
// 5. Convert ENU to azimuth + altitude
const azAlt = enuToAzAlt(enu)
// 6. Refraction correction
if (!airless) {
azAlt.altitude = applyRefraction(
azAlt.altitude,
observer.pressure,
observer.temperature,
)
}
return azAlt
}
// ─── Atmospheric refraction ───────────────────────────────────────────────────
/**
* Bennett (1982) atmospheric refraction correction.
* Adds the refraction amount to the geometric (airless) altitude.
*
* Accurate to ~0.1 arcmin for altitudes > 5°; degrades below that.
* At 0° altitude, refraction ≈ 34 arcmin.
*
* Formula: R = cot(h + 7.31 / (h + 4.4)) / 60 [degrees]
* where h is the geometric altitude in degrees.
*
* Pressure/temperature correction:
* R_adj = R × (P / 1010) × (283 / (273 + T))
*
* @param altitudeDeg - Geometric (airless) altitude in degrees
* @param pressure - Atmospheric pressure in millibars (default 1013.25)
* @param temperature - Temperature in Celsius (default 15)
* @returns Refraction to add to the altitude, in degrees
*/
export function bennettRefraction(
altitudeDeg: number,
pressure = 1013.25,
temperature = 15,
): number {
// No refraction below the geometric horizon (Bennett formula diverges below ~−1°)
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
// Pressure and temperature correction
const corrected = R * (pressure / 1010) * (283 / (273 + temperature))
return Math.max(0, corrected)
}
/**
* Apply refraction correction to an airless altitude.
* Returns the apparent (observed) altitude.
*/
export function applyRefraction(
airlessAlt: number,
pressure = 1013.25,
temperature = 15,
): number {
return airlessAlt + bennettRefraction(airlessAlt, pressure, temperature)
}
/**
* Remove refraction from an apparent altitude to get the airless altitude.
* Iterative inversion of the Bennett formula — converges in 3 iterations.
*/
export function removeRefraction(
apparentAlt: number,
pressure = 1013.25,
temperature = 15,
): number {
// Start from the apparent altitude and iterate
let airless = apparentAlt
for (let i = 0; i < 4; i++) {
airless = apparentAlt - bennettRefraction(airless, pressure, temperature)
}
return airless
}

410
src/spk/index.ts Normal file
View file

@ -0,0 +1,410 @@
/**
* spk — DAF/SPK kernel reader and Chebyshev segment evaluator.
*
* JPL planetary ephemerides are distributed as SPK (Spacecraft and Planet Kernel)
* files using the DAF (Double Precision Array File) binary format. DE442S uses
* SPK data type 2 (Chebyshev position-only) for all planetary segments.
*
* DAF structure:
* - File record (1024 bytes): magic bytes 'DAF/SPK', encoding, summary params
* - Summary records: each 1024 bytes, linked list, describe each segment
* - Data records: Chebyshev coefficient arrays for each segment
*
* SPK Type 2 record layout (one record per time interval):
* [0] = MID (midpoint of interval, ET seconds past J2000)
* [1] = RADIUS (half-width of interval, seconds)
* [2..n+1] = Chebyshev coefficients for X
* [n+2..2n+1] = Chebyshev coefficients for Y
* [2n+2..3n+1] = Chebyshev coefficients for Z
* where n = polynomial degree (RSIZE = 3n + 2)
*
* References:
* NAIF SPK Required Reading (spk.req)
* NAIF DAF Required Reading (daf.req)
* SPICE source: SPKE02.f, DAFRA.f
*/
import type { SpkSegment, StateVector } from '../types.js'
import { chebyshevEvalWithDerivative } from '../math/index.js'
// ─── NAIF body IDs ────────────────────────────────────────────────────────────
/** NAIF integer body IDs used in DE442S segment chaining */
export const NAIF_IDS = {
SSB: 0, // Solar System Barycenter
MERCURY_BARYCENTER: 1,
VENUS_BARYCENTER: 2,
EMB: 3, // Earth-Moon Barycenter
MARS_BARYCENTER: 4,
JUPITER_BARYCENTER: 5,
SATURN_BARYCENTER: 6,
URANUS_BARYCENTER: 7,
NEPTUNE_BARYCENTER: 8,
PLUTO_BARYCENTER: 9,
SUN: 10,
MOON: 301,
EARTH: 399,
} as const
/** Frame code for ICRF/J2000 (the inertial reference frame used by DE442S) */
export const FRAME_J2000 = 1
/** DAF record size in bytes */
const DAF_RECORD_SIZE = 1024
const BYTES_PER_DOUBLE = 8
// ─── SPK Kernel ───────────────────────────────────────────────────────────────
/**
* A loaded SPK kernel with segment index.
*/
export class SpkKernel {
private readonly buffer: ArrayBuffer
private readonly segments: SpkSegment[]
private readonly index: Map<string, SpkSegment[]>
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()
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)
}
}
/** 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)
}
/** 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)
}
/**
* Compute the state vector (position + velocity) for a body relative to a center.
* 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)
}
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
}
private getChained(target: number, center: number, et: number): StateVector {
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)
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)
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)
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)
// Earth/SSB = EMB/SSB - Earth/EMB
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)
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}`)
}
getSegments(): ReadonlyArray<SpkSegment> {
return this.segments
}
}
// ─── DAF parsing ──────────────────────────────────────────────────────────────
function parseDafFileRecord(buffer: ArrayBuffer): {
nd: number; ni: number; fward: number; bward: number; free: number; le: boolean
} {
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)
if (nd < 1 || nd > 100) {
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)
return { nd, ni, fward, bward, free, le }
}
function parseSummaryRecords(
buffer: ArrayBuffer,
fward: number,
nd: number,
ni: number,
le: boolean,
): SpkSegment[] {
const dv = new DataView(buffer)
const segments: SpkSegment[] = []
const summaryBytes = nd * BYTES_PER_DOUBLE + ni * 4
let recordNum = fward
while (recordNum !== 0) {
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))
let offset = recOffset + 24 // skip 3 control doubles (24 bytes)
for (let i = 0; i < nSummaries; i++) {
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 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
segments.push({ target, center, frame, dataType, startET, endET, dataOffset, dataSize })
}
recordNum = Math.round(nextRecord)
}
return segments
}
// ─── Segment evaluation ───────────────────────────────────────────────────────
function evaluateSegment(
buffer: ArrayBuffer,
seg: SpkSegment,
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}`)
}
/**
* Evaluate a Type 2 SPK segment at the given ET.
*/
export function evaluateType2(
buffer: ArrayBuffer,
seg: SpkSegment,
et: number,
le = true,
): StateVector {
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)
// degree = (rsize - 2) / 3 (Type 2 stores 3 components)
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))
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 [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] }
}
/**
* Evaluate a Type 3 SPK segment at the given ET.
* Type 3 stores separate position and velocity Chebyshev fits.
*/
export function evaluateType3(
buffer: ArrayBuffer,
seg: SpkSegment,
et: number,
le = true,
): StateVector {
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 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))
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 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]
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)
for (let k = 0; k < n; k++) {
arr[k] = dv.getFloat64(offset + k * BYTES_PER_DOUBLE, le)
}
return arr
}
function subtractSV(a: StateVector, b: StateVector): StateVector {
return {
position: [
a.position[0] - b.position[0],
a.position[1] - b.position[1],
a.position[2] - b.position[2],
],
velocity: [
a.velocity[0] - b.velocity[0],
a.velocity[1] - b.velocity[1],
a.velocity[2] - b.velocity[2],
],
}
}
// ─── Leap-second kernel ───────────────────────────────────────────────────────
/**
* Parse a NAIF LSK (Leap Second Kernel) text file.
* Extracts DELTET/DELTA_AT pairs and converts to (JD_UTC, deltaAT) pairs.
*/
export function parseLsk(text: string): ReadonlyArray<readonly [number, number]> {
const results: [number, number][] = []
const match = text.match(/DELTET\/DELTA_AT\s*=\s*\(\s*([\s\S]*?)\)/m)
if (!match) return results
const block = match[1]
const months: Record<string, number> = {
JAN: 1, FEB: 2, MAR: 3, APR: 4, MAY: 5, JUN: 6,
JUL: 7, AUG: 8, SEP: 9, 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
while ((m = pairRe.exec(block)) !== null) {
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 jdNoon = day + Math.floor((153 * mo + 2) / 5) + 365 * y +
Math.floor(y / 4) - Math.floor(y / 100) + Math.floor(y / 400) - 32045
// Midnight = JD - 0.5
results.push([jdNoon - 0.5, deltaAT])
}
return results
}

264
src/time/index.ts Normal file
View file

@ -0,0 +1,264 @@
/**
* time — Time scale conversions and Julian Date arithmetic.
*
* JPL DE ephemerides use Barycentric Dynamical Time (TDB) as their time argument,
* stored as seconds past J2000.0 (2000 Jan 1, 12:00:00 TDB = JD 2451545.0 TDB).
*
* Conversion chain: UTC → TAI → TT → TDB
* TAI = UTC + ΔAT (ΔAT from leap-second table, integer seconds)
* TT = TAI + 32.184s (exact, by definition)
* TDB ≈ TT + 0.001658 * sin(g) + ... (sub-millisecond correction)
*
* For Earth-rotation-based quantities (hour angle, ERA), we also need UT1:
* UT1 = UTC + (UT1 - UTC) (from IERS Bulletin A, typically < ±0.9s)
* ΔT = TT - UT1 (historically large; ~69s in 2024)
*
* References:
* IERS Conventions (2010) — Chapter 5, time scales
* NAIF LSK kernel (naif0012.tls) — ΔAT table + TDB approximation constants
* Espenak & Meeus — ΔT polynomial expressions
*/
import type { TimeScales } from '../types.js'
// ─── Constants ────────────────────────────────────────────────────────────────
/** Julian Date of J2000.0 epoch (2000 Jan 1, 12:00 TT) */
export const J2000 = 2451545.0
/** TT - TAI offset in seconds (exact, by definition) */
export const TT_MINUS_TAI = 32.184
/** Seconds per day */
export const SECONDS_PER_DAY = 86400.0
/** Days per Julian century */
export const DAYS_PER_JULIAN_CENTURY = 36525.0
// ─── Leap-second table ────────────────────────────────────────────────────────
/**
* ΔAT table (TAI - UTC in seconds), chronological.
* Source: NAIF naif0012.tls. Update when new leap seconds are announced.
*
* Each entry: [JD(UTC) when step takes effect, new ΔAT value]
* The step applies for all UTC times >= the entry JD.
*/
export const LEAP_SECOND_TABLE: ReadonlyArray<readonly [number, number]> = [
[2441317.5, 10], // 1972 Jan 1
[2441499.5, 11], // 1972 Jul 1
[2441683.5, 12], // 1973 Jan 1
[2442048.5, 13], // 1974 Jan 1
[2442413.5, 14], // 1975 Jan 1
[2442778.5, 15], // 1976 Jan 1
[2443144.5, 16], // 1977 Jan 1
[2443509.5, 17], // 1978 Jan 1
[2443874.5, 18], // 1979 Jan 1
[2444239.5, 19], // 1980 Jan 1
[2444786.5, 20], // 1981 Jul 1
[2445151.5, 21], // 1982 Jul 1
[2445516.5, 22], // 1983 Jul 1
[2446247.5, 23], // 1985 Jul 1
[2447161.5, 24], // 1988 Jan 1
[2447892.5, 25], // 1990 Jan 1
[2448257.5, 26], // 1991 Jan 1
[2448804.5, 27], // 1992 Jul 1
[2449169.5, 28], // 1993 Jul 1
[2449534.5, 29], // 1994 Jul 1
[2450083.5, 30], // 1996 Jan 1
[2450630.5, 31], // 1997 Jul 1
[2451179.5, 32], // 1999 Jan 1
[2453736.5, 33], // 2006 Jan 1
[2454832.5, 34], // 2009 Jan 1
[2456109.5, 35], // 2012 Jul 1
[2457204.5, 36], // 2015 Jul 1
[2457754.5, 37], // 2017 Jan 1
] 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
for (const [jd, dat] of LEAP_SECOND_TABLE) {
if (jdUTC >= jd) deltaAT = dat
else break
}
return deltaAT
}
// ─── Julian Date ─────────────────────────────────────────────────────────────
/**
* Convert a JavaScript Date (UTC) to Julian Date in UTC.
* Uses the standard formula; valid for dates after the Gregorian reform.
*/
export function dateToJD(date: Date): number {
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)
}
/**
* Julian centuries from J2000.0 (in TT).
* Used as the standard argument for precession and nutation polynomials.
*/
export function jdTTtoT(jdTT: number): number {
return (jdTT - J2000) / DAYS_PER_JULIAN_CENTURY
}
// ─── Time scale conversions ───────────────────────────────────────────────────
/**
* Compute all relevant time scales for a given UTC Date.
*
* @param utc - Input time in UTC
* @param ut1utcOverride - UT1 - UTC in seconds (from IERS Bulletin A). Falls
* back to 0 when not provided, which introduces ~0.9s error in ERA at most.
* @param deltaTOverride - TT - UT1 in seconds. When provided, overrides the
* polynomial approximation. Provide this for maximum accuracy.
*/
export function computeTimeScales(
utc: Date,
ut1utcOverride?: number,
deltaTOverride?: number,
): TimeScales {
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
// TT → TDB (periodic correction, sub-millisecond)
const tdbCorrection = tdbMinusTT(jdTT) / SECONDS_PER_DAY
const jdTDB = jdTT + tdbCorrection
// UT1
let jdUT1: number
let deltaT: number
if (ut1utcOverride !== undefined) {
jdUT1 = jdUTC + ut1utcOverride / SECONDS_PER_DAY
deltaT = (jdTT - jdUT1) * SECONDS_PER_DAY
} else if (deltaTOverride !== undefined) {
deltaT = deltaTOverride
jdUT1 = jdTT - deltaT / SECONDS_PER_DAY
} else {
deltaT = deltaTPolynomial(jdTT)
jdUT1 = jdTT - deltaT / SECONDS_PER_DAY
}
return { utc, jdUTC, jdTT, jdTDB, jdUT1, deltaT, deltaAT }
}
/**
* Convert Julian Date in TT to ET seconds past J2000.0.
* ET (Ephemeris Time) is used as the time argument in SPK kernels.
* TDB and TT differ by at most ~1.7ms; applying the correction here
* gives the proper SPICE-compatible ET value.
*/
export function jdTTtoET(jdTT: number): number {
const tdbCorr = tdbMinusTT(jdTT)
return (jdTT - J2000) * SECONDS_PER_DAY + tdbCorr
}
/**
* Approximate TDB - TT in seconds.
* Uses the SPICE-consistent formulation from the LSK kernel constants.
*
* TDB - TT = 0.001658 * sin(g) + 0.000014 * sin(2g)
* where g = 357.53° + 0.9856003° * (JD_TT - 2451545.0) [mean anomaly of the Sun]
*
* Maximum error: ~30 microseconds (acceptable for crescent work).
*/
export function tdbMinusTT(jdTT: number): number {
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)
}
/**
* Delta-T polynomial: TT - UT1 in seconds.
* Uses Espenak & Meeus expressions, piecewise by year range.
*
* Reference: NASA Five Millennium Canon of Solar Eclipses, Espenak & Meeus (2009)
*/
export function deltaTPolynomial(jdTT: number): number {
// Convert JD to decimal year
const y = 2000 + (jdTT - J2000) / 365.25
if (y < -500) {
const u = (y - 1820) / 100
return -20 + 32 * u * u
} else if (y < 500) {
const u = y / 100
return (
10583.6 - 1014.41 * u + 33.78311 * u * u - 5.952053 * u * u * u -
0.1798452 * u ** 4 + 0.022174192 * u ** 5 + 0.0090316521 * u ** 6
)
} else if (y < 1600) {
const u = (y - 1000) / 100
return (
1574.2 - 556.01 * u + 71.23472 * u * u + 0.319781 * u ** 3 -
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
} 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
)
} else if (y < 1860) {
const t = y - 1800
return (
13.72 - 0.332447 * t + 0.0068612 * t * t + 0.0041116 * t ** 3 -
0.00037436 * t ** 4 + 0.0000121272 * t ** 5 -
0.0000001699 * t ** 6 + 0.000000000875 * t ** 7
)
} else if (y < 1900) {
const t = y - 1860
return (
7.62 + 0.5737 * t - 0.251754 * t * t + 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
)
} else if (y < 1941) {
const t = y - 1920
return 21.20 + 0.84493 * t - 0.076100 * 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
} else if (y < 1986) {
const t = y - 1975
return 45.45 + 1.067 * t - t * t / 260 - t ** 3 / 718
} else if (y < 2005) {
const t = y - 2000
return (
63.86 + 0.3345 * t - 0.060374 * t * t + 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
} else if (y < 2150) {
return -20 + 32 * ((y - 1820) / 100) ** 2 - 0.5628 * (2150 - y)
} else {
const u = (y - 1820) / 100
return -20 + 32 * u * u
}
}

377
src/types.ts Normal file
View file

@ -0,0 +1,377 @@
// ─── Primitive geometry ──────────────────────────────────────────────────────
/** 3-element position or velocity vector in km or km/s */
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
}
/** Azimuth + altitude in degrees */
export interface AzAlt {
/** Degrees from North, measured clockwise (0 = N, 90 = E, 180 = S, 270 = W) */
azimuth: number
/** Degrees above the horizon (negative = below) */
altitude: number
}
// ─── Time ────────────────────────────────────────────────────────────────────
/** All relevant time scale values for a single moment */
export interface TimeScales {
utc: Date
/** Julian Date in UTC */
jdUTC: number
/** Julian Date in Terrestrial Time (TT = TAI + 32.184s) */
jdTT: number
/** Julian Date in Barycentric Dynamical Time (used by JPL ephemerides) */
jdTDB: number
/** Julian Date in UT1 (Earth rotation time) */
jdUT1: number
/** TT - UT1 in seconds (delta-T) */
deltaT: number
/** TAI - UTC in seconds (leap seconds count) */
deltaAT: number
}
// ─── Observer ────────────────────────────────────────────────────────────────
/** Observer location and environmental parameters */
export interface Observer {
/** Geodetic latitude in degrees (north positive) */
lat: number
/** Longitude in degrees (east positive) */
lon: number
/** Height above WGS84 ellipsoid in meters */
elevation: number
/** Optional label for the location */
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
/**
* Override UT1 - UTC in seconds (from IERS Bulletin A).
* Takes precedence over deltaT when both are provided.
*/
ut1utc?: number
/** Atmospheric pressure in millibars (default 1013.25) */
pressure?: number
/** Ambient temperature in Celsius (default 15) */
temperature?: number
}
// ─── Crescent geometry ───────────────────────────────────────────────────────
/**
* The five geometric quantities used by all major crescent visibility criteria.
* All values computed at best time (T_b) unless noted.
*/
export interface CrescentGeometry {
/** Arc of light: topocentric Sun-Moon angular separation (elongation), degrees */
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
/**
* Relative azimuth: Sun azimuth minus Moon azimuth, normalized to [-180, 180], degrees.
* Positive = Moon north of Sun.
*/
DAZ: number
/**
* Topocentric crescent width in arc minutes.
* Used directly in Odeh's polynomial V expression.
*/
W: number
/** Moonset minus sunset in minutes. Negative = Moon sets before Sun (no sighting possible). */
lag: number
}
// ─── Yallop q-test ───────────────────────────────────────────────────────────
/** Yallop q-test visibility category (NAO Technical Note 69) */
export type YallopCategory = 'A' | 'B' | 'C' | 'D' | 'E' | 'F'
/**
* Published q thresholds (Yallop 1997, NAO TN 69):
* A: q > +0.216 — Easily visible to the naked eye
* B: q > -0.014 — Visible under perfect conditions
* C: q > -0.160 — May need optical aid to find; visible to naked eye
* D: q > -0.232 — Optical aid needed; will not be visible to naked eye
* E: q > -0.293 — Not visible even with telescope
* F: q <= -0.293 — Below Danjon limit (Moon too close to Sun)
*/
export const YALLOP_THRESHOLDS = {
A: 0.216,
B: -0.014,
C: -0.160,
D: -0.232,
E: -0.293,
} as const
export const YALLOP_DESCRIPTIONS: Record<YallopCategory, string> = {
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
/** Visibility category A through F */
category: YallopCategory
/** Human-readable interpretation */
description: string
/** True for categories A and B */
isVisibleNakedEye: boolean
/** True for categories C and D */
requiresOpticalAid: boolean
/** True for category F */
isBelowDanjonLimit: boolean
/** Topocentric crescent width W' used in the q formula, arc minutes */
Wprime: number
}
// ─── Odeh criterion ──────────────────────────────────────────────────────────
/** Odeh visibility zone (Experimental Astronomy 2006) */
export type OdehZone = 'A' | 'B' | 'C' | 'D'
/**
* Published V thresholds (Odeh 2006):
* A: V >= 5.65 — Visible with naked eye
* B: V >= 2.00 — Visible with optical aid; may be seen with naked eye
* C: V >= -0.96 — Visible with optical aid only
* D: V < -0.96 — Not visible even with optical aid
*/
export const ODEH_THRESHOLDS = {
A: 5.65,
B: 2.00,
C: -0.96,
} as const
export const ODEH_DESCRIPTIONS: Record<OdehZone, string> = {
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
/** Visibility zone A through D */
zone: OdehZone
/** Human-readable interpretation */
description: string
/** True for zone A */
isVisibleNakedEye: boolean
/** True for zones A and B */
isVisibleWithOpticalAid: boolean
}
// ─── Moon phase ──────────────────────────────────────────────────────────────
export type MoonPhaseName =
| '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
/** Illuminated fraction 0-100 (percent) */
illumination: number
/** Hours since last new moon */
age: number
/** Ecliptic longitude of the Moon minus the Sun, degrees [0, 360) */
elongationDeg: number
/** True when Moon is moving away from the Sun (illumination increasing) */
isWaxing: boolean
/** UTC date of the next new moon */
nextNewMoon: Date
/** UTC date of the next full moon */
nextFullMoon: Date
/** UTC date of the previous new moon */
prevNewMoon: Date
}
// ─── Event times ─────────────────────────────────────────────────────────────
export interface SunMoonEvents {
/** UTC time of sunset for the given date at the observer's location */
sunsetUTC: Date | null
/** UTC time of moonset for the given date at the observer's location */
moonsetUTC: Date | null
/** UTC time of sunrise */
sunriseUTC: Date | null
/** UTC time of moonrise */
moonriseUTC: Date | null
/** UTC time when civil twilight ends (Sun at -6°) */
civilTwilightEndUTC: Date | null
/** UTC time when nautical twilight ends (Sun at -12°) */
nauticalTwilightEndUTC: Date | null
/** UTC time when astronomical twilight ends (Sun at -18°) */
astronomicalTwilightEndUTC: Date | null
}
// ─── Full moon sighting report ────────────────────────────────────────────────
export interface MoonSightingReport {
/** Date for which the sighting report was computed */
date: Date
/** Observer location used */
observer: Observer
// Event times
sunsetUTC: Date | null
moonsetUTC: Date | null
/** Moonset minus sunset, in minutes. Null if either event is null. */
lagMinutes: number | null
/** Best observation time (Odeh/Yallop: T_s + 4/9 * Lag) */
bestTimeUTC: Date | null
/** Conservative observation window [bestTime - 20min, bestTime + 20min] */
bestTimeWindowUTC: [Date, Date] | null
// At best time
/** Topocentric Moon position at best time */
moonPosition: AzAlt | null
/** Topocentric Sun position at best time */
sunPosition: AzAlt | null
/** Moon illumination percent at best time */
illumination: number | null
/** Hours since conjunction (new moon) */
moonAge: number | null
// Crescent geometry at best time
geometry: CrescentGeometry | null
// Visibility criteria results
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
// Metadata
/** Source ephemeris used for this calculation */
ephemerisSource: 'DE442S' | 'approximate'
/** Whether the Moon is even above the horizon at best time */
moonAboveHorizon: boolean | null
/** Whether sighting is geometrically possible (lag > 0, Moon above horizon at best time) */
sightingPossible: boolean
}
// ─── Kernel configuration ─────────────────────────────────────────────────────
/**
* How to source a binary kernel file.
* 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-calc
export interface KernelConfig {
/** Planetary SPK kernel — defaults to de442s.bsp via auto-download */
planetary?: KernelSource
/** Leap-second kernel — defaults to naif0012.tls via auto-download */
leapSeconds?: KernelSource
/**
* Directory for the download cache.
* Defaults to ~/.cache/moon-calc on POSIX, %LOCALAPPDATA%\moon-calc on Windows.
*/
cacheDir?: string
/**
* SHA-256 checksum of de442s.bsp for download verification.
* Bundled default matches the NAIF distribution as of 2024.
*/
checksumOverride?: string
}
// ─── Top-level options ────────────────────────────────────────────────────────
export interface SightingOptions {
/** Kernel acquisition configuration. Defaults to auto-download. */
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'
}
// ─── WGS84 constants ─────────────────────────────────────────────────────────
/** WGS84 reference ellipsoid parameters */
export const WGS84 = {
/** Semi-major axis in meters */
a: 6378137.0,
/** Inverse flattening */
invF: 298.257223563,
/** Flattening */
f: 1 / 298.257223563,
/** Semi-minor axis in meters */
get b() { return this.a * (1 - this.f) },
/** First eccentricity squared */
get e2() { return 2 * this.f - this.f * this.f },
} 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
/** NAIF body ID of the center body */
center: number
/** Reference frame code */
frame: number
/** SPK data type (2 = Chebyshev position only, 3 = Chebyshev position + velocity) */
dataType: 2 | 3
/** Segment start time in ET seconds past J2000 */
startET: number
/** Segment end time in ET seconds past J2000 */
endET: number
/** Byte offset of the data array in the file */
dataOffset: number
/** Number of double-precision values in the data array */
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
/** Half-width of the record interval in seconds */
radius: number
/** Chebyshev coefficients for X, Y, Z [3][degree+1] */
coeffs: Float64Array[]
/** Degree of the polynomial */
degree: number
}

270
src/visibility/index.ts Normal file
View file

@ -0,0 +1,270 @@
/**
* visibility — Crescent visibility criteria: Yallop and Odeh.
*
* Both criteria transform the five geometric quantities (ARCL, ARCV, DAZ, W, Lag)
* into a single score that maps to a visibility category.
*
* Yallop (NAO Technical Note 69, 1997):
* q = (ARCV - (11.8371 - 6.3226*W' + 0.7319*W'^2 - 0.1018*W'^3)) / 10
* where W' is the topocentric crescent width in arc minutes.
* Categories: A (q > 0.216) through F (q <= -0.293)
*
* Odeh (Experimental Astronomy 2006):
* V = ARCV - arcv_min(W)
* where arcv_min(W) = 11.8371 - 6.3226*W + 0.7319*W^2 - 0.1018*W^3 (same polynomial as Yallop)
* Zones: A (V >= 5.65) through D (V < -0.96)
*
* The geometric quantities must be computed at best time T_b using AIRLESS
* (refraction-free) topocentric altitudes for both models.
*
* References:
* Yallop (1997), A Method for Predicting the First Sighting of the New Crescent Moon,
* NAO Technical Note No. 69, Royal Greenwich Observatory
* Odeh (2006), New Criterion for Lunar Crescent Visibility,
* Experimental Astronomy 18(1), 39-64
*/
import type {
CrescentGeometry,
YallopResult,
YallopCategory,
OdehResult,
OdehZone,
} 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'
// ─── Shared polynomial ────────────────────────────────────────────────────────
/**
* The polynomial ARCV minimum as a function of crescent width W (arc minutes).
* This expression appears in both Yallop (as the denominator basis) and Odeh
* (as arcv_min directly).
*
* arcv_min(W) = 11.8371 - 6.3226*W + 0.7319*W^2 - 0.1018*W^3
*
* Represents the minimum arc of vision required for detection as a function
* of crescent width, derived empirically from historical observations.
*
* @param W - Topocentric crescent width in arc minutes
* @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
}
// ─── Yallop q-test ────────────────────────────────────────────────────────────
/**
* Compute the Yallop q parameter.
*
* q = (ARCV - arcv_min(W')) / 10
*
* A positive q means the actual ARCV exceeds the minimum required, scaled
* so that thresholds A–F correspond to intuitive distance from the boundary.
*
* @param ARCV - Topocentric airless arc of vision in degrees
* @param Wprime - Topocentric crescent width W' in arc minutes (Yallop's topocentric form)
* @returns q parameter (continuous)
*/
export function computeYallopQ(ARCV: number, Wprime: number): number {
return (ARCV - arcvMinimum(Wprime)) / 10
}
/**
* Map a q value to the Yallop category (A–F).
*
* Thresholds (Yallop 1997 Table 1):
* A: q > +0.216
* B: q > -0.014
* C: q > -0.160
* D: q > -0.232
* E: q > -0.293
* 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'
}
/**
* Compute the full Yallop result from crescent geometry.
*
* @param geometry - CrescentGeometry (W is assumed to be Wprime for Yallop purposes)
* @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)
return {
q,
category,
description: YALLOP_DESCRIPTIONS[category],
isVisibleNakedEye: category === 'A' || category === 'B',
requiresOpticalAid: category === 'C' || category === 'D',
isBelowDanjonLimit: category === 'F',
Wprime,
}
}
// ─── Odeh criterion ───────────────────────────────────────────────────────────
/**
* Compute the Odeh V parameter.
*
* V = ARCV - arcv_min(W)
*
* Where W is the topocentric crescent width in arc minutes (Odeh's formulation
* uses W directly, not the Yallop W' correction).
*
* A positive V indicates the observed ARCV exceeds the minimum threshold for
* visibility at that crescent width.
*
* @param ARCV - Topocentric airless arc of vision in degrees
* @param W - Topocentric crescent width in arc minutes
*/
export function computeOdehV(ARCV: number, W: number): number {
return ARCV - arcvMinimum(W)
}
/**
* Map a V value to the Odeh zone (A–D).
*
* Thresholds (Odeh 2006 Table 1):
* A: V >= 5.65 — Visible with naked eye
* B: V >= 2.00 — Visible with optical aid; may be seen naked eye
* C: V >= -0.96 — Visible with optical aid only
* 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'
}
/**
* Compute the full Odeh result from crescent geometry.
* 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)
return {
V,
zone,
description: ODEH_DESCRIPTIONS[zone],
isVisibleNakedEye: zone === 'A',
isVisibleWithOpticalAid: zone === 'A' || zone === 'B',
}
}
// ─── Geometry computation ─────────────────────────────────────────────────────
/**
* Compute all five crescent geometry quantities (ARCL, ARCV, DAZ, W, Lag)
* at a given best time.
*
* All angular quantities use AIRLESS topocentric positions as required by
* both Yallop and Odeh criteria.
*
* @param moonAirlessAzAlt - Moon topocentric airless az/alt at best time
* @param sunAirlessAzAlt - Sun topocentric airless az/alt at best time
* @param moonGCRS - Topocentric Moon position vector (km) — for ARCL and W
* @param sunGCRS - Topocentric Sun position vector (km) — for ARCL
* @param sunsetUTC - UTC time of sunset
* @param moonsetUTC - UTC time of moonset
*/
export function computeCrescentGeometry(
moonAirlessAzAlt: { azimuth: number; altitude: number },
sunAirlessAzAlt: { azimuth: number; altitude: number },
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
// 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
// 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)
// W: topocentric crescent width in arc minutes
const { W } = computeCrescentWidth(moonGCRS, ARCL)
// lag: moonset minus sunset in minutes (negative = Moon sets before Sun)
const lag = (moonsetUTC.getTime() - sunsetUTC.getTime()) / 60000
return { ARCL, ARCV, DAZ, W, lag }
}
// ─── Guidance text ────────────────────────────────────────────────────────────
/**
* Generate human-readable sighting guidance based on the crescent report.
*
* @param yallop - Yallop result
* @param odeh - Odeh result
* @param moonAz - Moon azimuth at best time (degrees from North)
* @param moonAlt - Moon altitude at best time (degrees)
* @param bestTimeUTC - Best observation time
* @param lagMinutes - Lag in minutes
* @returns Guidance string for observers
*/
export function buildGuidanceText(
yallop: YallopResult,
odeh: OdehResult,
moonAz: number,
moonAlt: number,
bestTimeUTC: Date,
lagMinutes: number,
): string {
const direction = azimuthToCardinal(moonAz)
const timeStr = bestTimeUTC.toISOString().replace('T', ' ').replace(/\.\d+Z$/, ' UTC')
const lagStr = `${Math.round(lagMinutes)} min after sunset`
let visibility: string
if (yallop.isVisibleNakedEye && odeh.isVisibleNakedEye) {
visibility = 'should be visible to the naked eye'
} else if (odeh.isVisibleWithOpticalAid) {
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)'
} else {
visibility = 'is not expected to be visible even with optical aid'
}
return (
`Best time to look: ${timeStr} (${lagStr}). ` +
`Look ${direction} at ${Math.round(moonAlt)}° above the horizon. ` +
`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]
}

90
test-cjs.cjs Normal file
View file

@ -0,0 +1,90 @@
'use strict'
/**
* moon-calc CJS test suite
* Runs with: node test-cjs.cjs
* Verifies the CommonJS build is functional.
*/
const assert = require('node:assert/strict')
const {
YALLOP_THRESHOLDS,
YALLOP_DESCRIPTIONS,
ODEH_THRESHOLDS,
ODEH_DESCRIPTIONS,
WGS84,
getMoonPhase,
initKernels,
downloadKernels,
verifyKernels,
getMoonSightingReport,
getSunMoonEvents,
} = require('./dist/index.cjs')
let passed = 0
let failed = 0
function test(name, fn) {
try {
fn()
console.log(` [${name}]... PASS`)
passed++
} catch (err) {
console.error(` [${name}]... FAIL: ${err.message}`)
failed++
}
}
console.log('CJS compatibility:')
test('require() works', () => {
assert.ok(YALLOP_THRESHOLDS !== undefined)
})
test('YALLOP_THRESHOLDS.A is 0.216', () => {
assert.equal(YALLOP_THRESHOLDS.A, 0.216)
})
test('ODEH_THRESHOLDS.A is 5.65', () => {
assert.equal(ODEH_THRESHOLDS.A, 5.65)
})
test('WGS84.a is 6378137.0', () => {
assert.equal(WGS84.a, 6378137.0)
})
test('All API functions are exported', () => {
assert.equal(typeof getMoonPhase, 'function')
assert.equal(typeof initKernels, 'function')
assert.equal(typeof downloadKernels, 'function')
assert.equal(typeof verifyKernels, 'function')
assert.equal(typeof getMoonSightingReport, 'function')
assert.equal(typeof getSunMoonEvents, 'function')
})
console.log('\nCJS getMoonPhase:')
test('getMoonPhase returns valid phase', () => {
const valid = new Set([
'new-moon', 'waxing-crescent', 'first-quarter', 'waxing-gibbous',
'full-moon', 'waning-gibbous', 'last-quarter', 'waning-crescent',
])
const p = getMoonPhase(new Date('2025-03-14T12:00:00Z'))
assert.ok(valid.has(p.phase), `got: ${p.phase}`)
})
test('getMoonPhase illumination in [0, 100]', () => {
const p = getMoonPhase(new Date('2025-03-01T12:00:00Z'))
assert.ok(p.illumination >= 0 && p.illumination <= 100)
})
test('getMoonPhase near full moon has high illumination', () => {
const p = getMoonPhase(new Date('2025-03-14T12:00:00Z'))
assert.ok(p.illumination > 85, `illumination=${p.illumination.toFixed(1)}%`)
})
test('getMoonPhase Dates are Date objects', () => {
const p = getMoonPhase(new Date('2025-03-01T12:00:00Z'))
assert.ok(p.nextNewMoon instanceof Date)
assert.ok(p.prevNewMoon instanceof Date)
assert.ok(p.nextFullMoon instanceof Date)
})
console.log(`\n${passed + failed} tests: ${passed} passed, ${failed} failed`)
if (failed > 0) {
process.exit(1)
}

233
test.mjs Normal file
View file

@ -0,0 +1,233 @@
/**
* moon-calc ESM test suite
* Runs with: node test.mjs
* All tests use plain assert — no test framework.
*/
import assert from 'node:assert/strict'
import {
// Constants
YALLOP_THRESHOLDS,
YALLOP_DESCRIPTIONS,
ODEH_THRESHOLDS,
ODEH_DESCRIPTIONS,
WGS84,
// API
getMoonPhase,
initKernels,
downloadKernels,
verifyKernels,
getMoonSightingReport,
getSunMoonEvents,
} from './dist/index.mjs'
let passed = 0
let failed = 0
function test(name, fn) {
try {
fn()
console.log(` [${name}]... PASS`)
passed++
} catch (err) {
console.error(` [${name}]... FAIL: ${err.message}`)
failed++
}
}
// ─── Constants ────────────────────────────────────────────────────────────────
console.log('Constants:')
test('YALLOP_THRESHOLDS.A is 0.216', () => {
assert.equal(YALLOP_THRESHOLDS.A, 0.216)
})
test('YALLOP_THRESHOLDS.E is -0.293', () => {
assert.equal(YALLOP_THRESHOLDS.E, -0.293)
})
test('All Yallop thresholds are defined', () => {
for (const key of ['A', 'B', 'C', 'D', 'E']) {
assert.ok(typeof YALLOP_THRESHOLDS[key] === 'number', `${key} should be a number`)
}
})
test('Yallop thresholds descend A > B > C > D > E', () => {
assert.ok(YALLOP_THRESHOLDS.A > YALLOP_THRESHOLDS.B)
assert.ok(YALLOP_THRESHOLDS.B > YALLOP_THRESHOLDS.C)
assert.ok(YALLOP_THRESHOLDS.C > YALLOP_THRESHOLDS.D)
assert.ok(YALLOP_THRESHOLDS.D > YALLOP_THRESHOLDS.E)
})
test('ODEH_THRESHOLDS.A is 5.65', () => {
assert.equal(ODEH_THRESHOLDS.A, 5.65)
})
test('ODEH_THRESHOLDS.C is -0.96', () => {
assert.equal(ODEH_THRESHOLDS.C, -0.96)
})
test('Odeh thresholds descend A > B > C', () => {
assert.ok(ODEH_THRESHOLDS.A > ODEH_THRESHOLDS.B)
assert.ok(ODEH_THRESHOLDS.B > ODEH_THRESHOLDS.C)
})
test('WGS84.a is 6378137.0', () => {
assert.equal(WGS84.a, 6378137.0)
})
test('WGS84.invF is 298.257223563', () => {
assert.equal(WGS84.invF, 298.257223563)
})
test('WGS84.e2 is positive and < 1', () => {
assert.ok(WGS84.e2 > 0 && WGS84.e2 < 1, `e2=${WGS84.e2}`)
})
test('WGS84.b < WGS84.a (oblate spheroid)', () => {
assert.ok(WGS84.b < WGS84.a)
})
test('Yallop descriptions are non-empty strings', () => {
for (const cat of ['A', 'B', 'C', 'D', 'E', 'F']) {
assert.ok(typeof YALLOP_DESCRIPTIONS[cat] === 'string' && YALLOP_DESCRIPTIONS[cat].length > 0)
}
})
test('Odeh descriptions are non-empty strings', () => {
for (const zone of ['A', 'B', 'C', 'D']) {
assert.ok(typeof ODEH_DESCRIPTIONS[zone] === 'string' && ODEH_DESCRIPTIONS[zone].length > 0)
}
})
// ─── API function exports ──────────────────────────────────────────────────────
console.log('\nAPI exports:')
test('getMoonPhase is a function', () => {
assert.equal(typeof getMoonPhase, 'function')
})
test('initKernels is a function', () => {
assert.equal(typeof initKernels, 'function')
})
test('downloadKernels is a function', () => {
assert.equal(typeof downloadKernels, 'function')
})
test('verifyKernels is a function', () => {
assert.equal(typeof verifyKernels, 'function')
})
test('getMoonSightingReport is a function', () => {
assert.equal(typeof getMoonSightingReport, 'function')
})
test('getSunMoonEvents is a function', () => {
assert.equal(typeof getSunMoonEvents, 'function')
})
// ─── getMoonPhase (synchronous, no kernel) ─────────────────────────────────────
console.log('\ngetMoonPhase — structure:')
const VALID_PHASES = new Set([
'new-moon', 'waxing-crescent', 'first-quarter', 'waxing-gibbous',
'full-moon', 'waning-gibbous', 'last-quarter', 'waning-crescent',
])
// Test with a known reference date: 2025-03-01 UTC
// At this date the Moon was a waxing crescent (~2 days after new moon Feb 28)
const DATE_MARCH_1_2025 = new Date('2025-03-01T12:00:00Z')
const phase_march1 = getMoonPhase(DATE_MARCH_1_2025)
test('getMoonPhase returns an object', () => {
assert.ok(phase_march1 !== null && typeof phase_march1 === 'object')
})
test('getMoonPhase.phase is a valid phase name', () => {
assert.ok(VALID_PHASES.has(phase_march1.phase), `got: ${phase_march1.phase}`)
})
test('getMoonPhase.illumination is in [0, 100]', () => {
assert.ok(phase_march1.illumination >= 0 && phase_march1.illumination <= 100,
`illumination=${phase_march1.illumination}`)
})
test('getMoonPhase.age is >= 0', () => {
assert.ok(phase_march1.age >= 0, `age=${phase_march1.age}`)
})
test('getMoonPhase.elongationDeg is in [0, 180]', () => {
assert.ok(phase_march1.elongationDeg >= 0 && phase_march1.elongationDeg <= 180,
`elongationDeg=${phase_march1.elongationDeg}`)
})
test('getMoonPhase.isWaxing is a boolean', () => {
assert.equal(typeof phase_march1.isWaxing, 'boolean')
})
test('getMoonPhase.nextNewMoon is a Date', () => {
assert.ok(phase_march1.nextNewMoon instanceof Date)
})
test('getMoonPhase.prevNewMoon is a Date', () => {
assert.ok(phase_march1.prevNewMoon instanceof Date)
})
test('getMoonPhase.nextFullMoon is a Date', () => {
assert.ok(phase_march1.nextFullMoon instanceof Date)
})
test('getMoonPhase.prevNewMoon is before reference date', () => {
assert.ok(phase_march1.prevNewMoon < DATE_MARCH_1_2025,
`prevNewMoon=${phase_march1.prevNewMoon.toISOString()}`)
})
test('getMoonPhase.nextNewMoon is after prevNewMoon', () => {
assert.ok(phase_march1.nextNewMoon > phase_march1.prevNewMoon)
})
console.log('\ngetMoonPhase — phase boundaries:')
// 2025-03-14 was close to full moon (illumination should be high)
const DATE_FULL_MOON = new Date('2025-03-14T12:00:00Z')
const phase_full = getMoonPhase(DATE_FULL_MOON)
test('Near full moon: illumination > 85%', () => {
assert.ok(phase_full.illumination > 85,
`illumination at full moon=${phase_full.illumination.toFixed(1)}%`)
})
test('Near full moon: phase is full-moon or waxing/waning gibbous', () => {
const valid = new Set(['full-moon', 'waxing-gibbous', 'waning-gibbous'])
assert.ok(valid.has(phase_full.phase), `got: ${phase_full.phase}`)
})
test('Near full moon: elongation > 120°', () => {
assert.ok(phase_full.elongationDeg > 120, `elongation=${phase_full.elongationDeg}`)
})
// 2025-03-29 is close to new moon (illumination should be low)
const DATE_NEW_MOON = new Date('2025-03-29T12:00:00Z')
const phase_new = getMoonPhase(DATE_NEW_MOON)
test('Near new moon: illumination < 10%', () => {
assert.ok(phase_new.illumination < 10,
`illumination at new moon=${phase_new.illumination.toFixed(1)}%`)
})
test('Near new moon: elongation < 30°', () => {
assert.ok(phase_new.elongationDeg < 30, `elongation=${phase_new.elongationDeg}`)
})
console.log('\ngetMoonPhase — consistency:')
// Two dates: one clearly waxing, one clearly waning
const DATE_WAXING = new Date('2025-03-05T12:00:00Z') // ~7 days after new moon
const DATE_WANING = new Date('2025-03-20T12:00:00Z') // ~6 days after full moon
const phase_waxing = getMoonPhase(DATE_WAXING)
const phase_waning = getMoonPhase(DATE_WANING)
test('5 days after new moon: isWaxing = true', () => {
assert.equal(phase_waxing.isWaxing, true)
})
test('6 days after full moon: isWaxing = false', () => {
assert.equal(phase_waning.isWaxing, false)
})
test('getMoonPhase with default date (now) returns valid result', () => {
const nowPhase = getMoonPhase()
assert.ok(VALID_PHASES.has(nowPhase.phase))
assert.ok(nowPhase.illumination >= 0 && nowPhase.illumination <= 100)
})
// Synodic month duration check: nextNewMoon - prevNewMoon ≈ 29.53 days
test('Synodic month duration is ~29.5 days (±0.5)', () => {
const synodicMs = phase_march1.nextNewMoon.getTime() - phase_march1.prevNewMoon.getTime()
const synodicDays = synodicMs / 86400000
assert.ok(
synodicDays > 29.0 && synodicDays < 30.1,
`synodic month=${synodicDays.toFixed(2)} days`,
)
})
// ─── Summary ─────────────────────────────────────────────────────────────────
console.log(`\n${passed + failed} tests: ${passed} passed, ${failed} failed`)
if (failed > 0) {
process.exit(1)
}

18
tsconfig.json Normal file
View file

@ -0,0 +1,18 @@
{
"compilerOptions": {
"target": "ES2020",
"module": "ESNext",
"moduleResolution": "bundler",
"strict": true,
"forceConsistentCasingInFileNames": true,
"esModuleInterop": true,
"declaration": true,
"declarationMap": true,
"sourceMap": true,
"outDir": "dist",
"rootDir": "src",
"types": ["node"]
},
"include": ["src"],
"exclude": ["node_modules", "dist"]
}

36
tsup.config.ts Normal file
View file

@ -0,0 +1,36 @@
import { defineConfig } from 'tsup'
export default defineConfig([
// Library build (CJS + ESM)
{
entry: { index: 'src/index.ts' },
format: ['cjs', 'esm'],
dts: true,
clean: true,
outDir: 'dist',
splitting: false,
sourcemap: true,
target: 'es2020',
platform: 'node',
outExtension({ format }) {
return { js: format === 'cjs' ? '.cjs' : '.mjs' }
},
},
// CLI build (CJS only, with shebang)
{
entry: { 'cli/index': 'src/cli/index.ts' },
format: ['cjs'],
dts: false,
outDir: 'dist',
splitting: false,
sourcemap: false,
target: 'es2020',
platform: 'node',
banner: {
js: '#!/usr/bin/env node',
},
outExtension() {
return { js: '.cjs' }
},
},
])