commit fbb5bb5dc3c4c2411088f1aabefa61d9231416a6 Author: Aric Camarata Date: Sun Mar 8 13:07:24 2026 -0400 Initial release: moon_sighting v1.0.0 diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..4c0728e --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,33 @@ +name: CI + +on: + push: + branches: [main] + pull_request: + branches: [main] + +jobs: + test: + runs-on: ubuntu-latest + strategy: + matrix: + sdk: [stable, beta] + steps: + - uses: actions/checkout@v4 + - uses: dart-lang/setup-dart@v1 + with: + sdk: ${{ matrix.sdk }} + - run: dart pub get + - run: dart analyze --fatal-infos + - run: dart test + - run: dart format --set-exit-if-changed . + + publish-check: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: dart-lang/setup-dart@v1 + with: + sdk: stable + - run: dart pub get + - run: dart pub publish --dry-run diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..d736e86 --- /dev/null +++ b/.gitignore @@ -0,0 +1,20 @@ +.dart_tool/ +.packages +build/ +pubspec.lock +doc/api/ +*.iml +.idea/ +.DS_Store +.claude/ +.env +.env.* +.vscode/* +.codex/ +.cursor/ +.aider/ +.aider.chat.history.md +.continue/ +.windsurf/ +.gemini/ +.codeium/ diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..f4018e3 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,28 @@ +# Changelog + +All notable changes to this project will be documented in this file. + +The format is based on [Keep a Changelog](https://keepachangelog.com/). + +## [1.0.0] - 2026-03-08 + +### Added + +- `getMoonPhase()` returns phase name, emoji symbol, illumination percentage, age in hours, elongation, waxing/waning state, and next/previous new moon and full moon dates +- `getMoonPosition()` returns topocentric azimuth, altitude (with Bennett refraction), geocentric distance, and parallactic angle for any observer location +- `getMoonIllumination()` returns illumination fraction, phase cycle position, bright limb position angle, and waxing/waning indicator +- `getMoonVisibilityEstimate()` computes Odeh crescent visibility criterion (V parameter, zones A through D) from approximate Meeus positions at a given observation time +- `getMoon()` combined snapshot bundling all four functions in a single call +- `arcvMinimum()` the Odeh 2006 polynomial for minimum arc of vision as a function of crescent width +- Meeus Ch. 25 Sun position (geocentric, < 0.01 deg accuracy) +- Meeus Ch. 47 Moon position with 30 longitude terms and 20 latitude terms (< 0.3 deg accuracy) +- Meeus Ch. 49 new moon and full moon time estimation (within ~2 hours) +- Full delta-T polynomial (Espenak & Meeus, piecewise by year) +- Leap-second table through 2017 Jan 1 (37 seconds) +- WGS84 geodetic to ECEF conversion +- Bennett (1982) atmospheric refraction with pressure/temperature correction +- Simplified GCRS to ITRS frame rotation via Earth Rotation Angle (sufficient for lite-mode ~1 deg accuracy) +- Yallop q-test with categories A through F (NAO Technical Note 69) +- Odeh V-parameter with zones A through D (Experimental Astronomy 2006) +- Input validation for latitude, longitude ranges +- Complete type system: `MoonPhaseResult`, `MoonPosition`, `MoonIlluminationResult`, `MoonVisibilityEstimate`, `MoonSnapshot`, `YallopResult`, `OdehResult`, `CrescentGeometry` diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..f8c021c --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 index 0000000..284a4c7 --- /dev/null +++ b/README.md @@ -0,0 +1,149 @@ +# moon_sighting + +[![pub package](https://img.shields.io/pub/v/moon_sighting.svg)](https://pub.dev/packages/moon_sighting) +[![License: MIT](https://img.shields.io/badge/license-MIT-blue.svg)](LICENSE) + +Lunar crescent visibility for Dart and Flutter. Computes moon phase, topocentric position, illumination, and Yallop/Odeh crescent visibility criteria using Meeus algorithms. Zero dependencies. + +## Installation + +```yaml +dependencies: + moon_sighting: ^1.0.0 +``` + +## Quick Start + +```dart +import 'package:moon_sighting/moon_sighting.dart'; + +// Current moon phase +final phase = getMoonPhase(); +print('${phase.phaseName} ${phase.phaseSymbol}'); +print('Illumination: ${phase.illumination.toStringAsFixed(1)}%'); +print('Age: ${phase.age.toStringAsFixed(1)} hours'); + +// Moon position for an observer +final pos = getMoonPosition(DateTime.now(), 51.5074, -0.1278, elevation: 10); +print('Azimuth: ${pos.azimuth.toStringAsFixed(1)}'); +print('Altitude: ${pos.altitude.toStringAsFixed(1)}'); +print('Distance: ${pos.distance.toStringAsFixed(0)} km'); + +// Crescent visibility estimate (pass a post-sunset time) +final vis = getMoonVisibilityEstimate( + DateTime.utc(2025, 3, 31, 18, 30), + 21.4225, 39.8262, // Mecca +); +print('Zone: ${vis.zone.label}'); // A through D +print('Visible naked eye: ${vis.isVisibleNakedEye}'); +``` + +## API + +### getMoonPhase([DateTime? date]) + +Returns `MoonPhaseResult` with: + +| Field | Type | Description | +| --- | --- | --- | +| `phase` | `MoonPhaseName` | Enum: newMoon, waxingCrescent, firstQuarter, etc. | +| `phaseName` | `String` | Human-readable name | +| `phaseSymbol` | `String` | Moon emoji | +| `illumination` | `double` | Percent illuminated (0-100) | +| `age` | `double` | Hours since last new moon | +| `elongationDeg` | `double` | Moon-Sun elongation in degrees | +| `isWaxing` | `bool` | True when illumination is increasing | +| `nextNewMoon` | `DateTime` | Next new moon (UTC) | +| `nextFullMoon` | `DateTime` | Next full moon (UTC) | +| `prevNewMoon` | `DateTime` | Previous new moon (UTC) | + +### getMoonPosition(DateTime? date, double lat, double lon, {double elevation = 0}) + +Returns `MoonPosition` with: + +| Field | Type | Description | +| --- | --- | --- | +| `azimuth` | `double` | Degrees from North, clockwise | +| `altitude` | `double` | Degrees above horizon (refraction applied) | +| `distance` | `double` | Earth-Moon distance in km | +| `parallacticAngle` | `double` | Parallactic angle in radians | + +### getMoonIllumination([DateTime? date]) + +Returns `MoonIlluminationResult` with: + +| Field | Type | Description | +| --- | --- | --- | +| `fraction` | `double` | Illuminated fraction (0-1) | +| `phase` | `double` | Phase cycle (0=new, 0.25=Q1, 0.5=full, 0.75=Q3) | +| `angle` | `double` | Bright limb position angle (radians) | +| `isWaxing` | `bool` | True when waxing | + +### getMoonVisibilityEstimate(DateTime? date, double lat, double lon, {double elevation = 0}) + +Returns `MoonVisibilityEstimate` with: + +| Field | Type | Description | +| --- | --- | --- | +| `v` | `double` | Odeh V parameter | +| `zone` | `OdehZone` | Visibility zone (a through d) | +| `description` | `String` | Zone description | +| `isVisibleNakedEye` | `bool` | True for zone A | +| `isVisibleWithOpticalAid` | `bool` | True for zones A and B | +| `arcl` | `double` | Sun-Moon elongation (degrees) | +| `arcv` | `double` | Arc of vision (degrees) | +| `w` | `double` | Crescent width (arc minutes) | +| `moonAboveHorizon` | `bool` | Moon above horizon at given time | + +### getMoon(DateTime? date, double lat, double lon, {double elevation = 0}) + +Returns `MoonSnapshot` combining all four results: `phase`, `position`, `illumination`, `visibility`. + +## Visibility Criteria + +The Odeh criterion (2006) classifies crescent visibility into four zones: + +| Zone | V threshold | Meaning | +| --- | --- | --- | +| 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 | + +The Yallop q-test (1997) uses six categories A through F with different threshold semantics. Both criteria use the same underlying polynomial for the minimum arc of vision as a function of crescent width. + +## Accuracy + +All positions use Meeus (1998) approximations: + +- **Sun:** < 0.01 deg ecliptic longitude error +- **Moon:** < 0.3 deg longitude, < 0.2 deg latitude +- **Distance:** ~300 km error (~0.08%) +- **New/full moon times:** within ~2 hours + +The GCRS-to-topocentric conversion uses a simplified Earth rotation (ERA only, no nutation/precession), adding ~1 deg systematic error. This is acceptable for phase displays, illumination calculations, and approximate visibility estimates. + +For high-precision crescent sighting reports with DE442S ephemeris accuracy, see the TypeScript [moon-sighting](https://github.com/acamarata/moon-sighting) package. + +## Lite Mode vs Full Mode + +This Dart package implements lite mode only. It covers moon phase, position, illumination, and quick visibility estimates using analytical Meeus algorithms with no external data files. + +The full-mode features (DE442S JPL ephemeris, sub-arcsecond accuracy, rise/set event finding, best-time optimization, full sighting reports) are available in the TypeScript package. + +## Compatibility + +- Dart SDK >= 3.7.0 +- Works in Flutter, Dart CLI, and server-side Dart +- Zero runtime dependencies + +## Related Packages + +- [moon-sighting](https://github.com/acamarata/moon-sighting) (TypeScript) - Full-accuracy lunar crescent visibility with DE442S ephemeris +- [nrel-spa](https://github.com/acamarata/nrel-spa) (TypeScript) - Pure JS NREL Solar Position Algorithm +- [pray-calc](https://github.com/acamarata/pray-calc) (JavaScript) - Islamic prayer time calculation +- [luxon-hijri](https://github.com/acamarata/luxon-hijri) (TypeScript) - Hijri/Gregorian calendar conversion + +## License + +MIT diff --git a/analysis_options.yaml b/analysis_options.yaml new file mode 100644 index 0000000..572dd23 --- /dev/null +++ b/analysis_options.yaml @@ -0,0 +1 @@ +include: package:lints/recommended.yaml diff --git a/lib/moon_sighting.dart b/lib/moon_sighting.dart new file mode 100644 index 0000000..a8f9914 --- /dev/null +++ b/lib/moon_sighting.dart @@ -0,0 +1,16 @@ +/// Lunar crescent visibility for Dart and Flutter. +/// +/// Moon phase, position, illumination, and Yallop/Odeh visibility criteria +/// using Meeus algorithms. Zero dependencies. +/// +/// Five public functions: +/// - [getMoonPhase] - phase name, illumination, age, next events +/// - [getMoonPosition] - topocentric azimuth, altitude, distance +/// - [getMoonIllumination] - illumination fraction, phase cycle, bright limb +/// - [getMoonVisibilityEstimate] - Odeh crescent visibility estimate +/// - [getMoon] - combined snapshot of all four +library; + +export 'src/types.dart'; +export 'src/api.dart'; +export 'src/visibility.dart' show arcvMinimum; diff --git a/lib/src/api.dart b/lib/src/api.dart new file mode 100644 index 0000000..e554490 --- /dev/null +++ b/lib/src/api.dart @@ -0,0 +1,304 @@ +/// User-facing API: 5 lite-mode functions. +/// +/// All functions work without a kernel (Meeus-based approximations). +library; + +import 'dart:math' as math; + +import 'math.dart'; +import 'time.dart'; +import 'types.dart'; +import 'bodies.dart'; +import 'observer.dart'; +import 'visibility.dart'; + +// ── Phase display lookup ───────────────────────────────────────────────────── + +const _phaseDisplay = { + MoonPhaseName.newMoon: (name: 'New Moon', symbol: '\u{1F311}'), + MoonPhaseName.waxingCrescent: (name: 'Waxing Crescent', symbol: '\u{1F312}'), + MoonPhaseName.firstQuarter: (name: 'First Quarter', symbol: '\u{1F313}'), + MoonPhaseName.waxingGibbous: (name: 'Waxing Gibbous', symbol: '\u{1F314}'), + MoonPhaseName.fullMoon: (name: 'Full Moon', symbol: '\u{1F315}'), + MoonPhaseName.waningGibbous: (name: 'Waning Gibbous', symbol: '\u{1F316}'), + MoonPhaseName.lastQuarter: (name: 'Last Quarter', symbol: '\u{1F317}'), + MoonPhaseName.waningCrescent: (name: 'Waning Crescent', symbol: '\u{1F318}'), +}; + +/// Map elongation and waxing direction to a named phase. +MoonPhaseName _elongationToPhase(double elongationDeg, bool isWaxing) { + if (elongationDeg < 5) return MoonPhaseName.newMoon; + if (elongationDeg > 175) return MoonPhaseName.fullMoon; + if (elongationDeg < 85) { + return isWaxing + ? MoonPhaseName.waxingCrescent + : MoonPhaseName.waningCrescent; + } + if (elongationDeg < 95) { + return isWaxing ? MoonPhaseName.firstQuarter : MoonPhaseName.lastQuarter; + } + return isWaxing ? MoonPhaseName.waxingGibbous : MoonPhaseName.waningGibbous; +} + +// ── Validation ─────────────────────────────────────────────────────────────── + +void _validateLatitude(double lat, String label) { + if (!lat.isFinite || lat < -90 || lat > 90) { + throw RangeError('$label: latitude must be in [-90, 90], got $lat'); + } +} + +void _validateLongitude(double lon, String label) { + if (!lon.isFinite || lon < -180 || lon > 180) { + throw RangeError('$label: longitude must be in [-180, 180], got $lon'); + } +} + +// ── getMoonPhase ───────────────────────────────────────────────────────────── + +/// Compute the Moon's current phase, illumination, and next phase times. +/// +/// Works without a kernel (uses Meeus approximation). +/// +/// ```dart +/// final phase = getMoonPhase(); +/// print(phase.phaseName); // 'Waxing Crescent' +/// print(phase.illumination); // 14.3 (percent) +/// print(phase.nextFullMoon); // DateTime +/// ``` +MoonPhaseResult getMoonPhase([DateTime? date]) { + final d = date ?? DateTime.now(); + final ts = computeTimeScales(d); + final (:moonGCRS, :sunGCRS) = getMoonSunApproximate(ts.jdTT); + + final illum = computeIllumination(moonGCRS, sunGCRS); + final illuminationPct = illum.illumination * 100; + + final prevNewMoonJD = nearestNewMoon(ts.jdTT - 15); + final age = (ts.jdTT - prevNewMoonJD) * 24; + + final phaseKey = _elongationToPhase(illum.elongationDeg, illum.isWaxing); + final display = _phaseDisplay[phaseKey]!; + + final nextNewMoonJD = nearestNewMoon(ts.jdTT + 15); + final nextFullMoonJD = nearestFullMoon(ts.jdTT); + + return MoonPhaseResult( + phase: phaseKey, + phaseName: display.name, + phaseSymbol: display.symbol, + illumination: illuminationPct, + age: age, + elongationDeg: illum.elongationDeg, + isWaxing: illum.isWaxing, + nextNewMoon: jdToDate(nextNewMoonJD), + nextFullMoon: jdToDate(nextFullMoonJD), + prevNewMoon: jdToDate(prevNewMoonJD), + ); +} + +// ── getMoonPosition ────────────────────────────────────────────────────────── + +/// Compute the Moon's topocentric position for an observer. +/// +/// Works without a kernel (Meeus Ch. 47 approximation). +/// Accuracy: azimuth/altitude ~0.3 deg, distance ~300 km. +/// +/// ```dart +/// final pos = getMoonPosition(DateTime.now(), 51.5, -0.1); +/// print('${pos.azimuth}, ${pos.altitude}'); +/// ``` +MoonPosition getMoonPosition( + DateTime? date, + double lat, + double lon, { + double elevation = 0, +}) { + final d = date ?? DateTime.now(); + _validateLatitude(lat, 'getMoonPosition'); + _validateLongitude(lon, 'getMoonPosition'); + + final ts = computeTimeScales(d); + final (:moonGCRS, sunGCRS: _) = getMoonSunApproximate(ts.jdTT); + + final azAlt = computeAzAlt(moonGCRS, lat, lon, elevation, ts); + final distance = vnorm(moonGCRS); + + // Equatorial coordinates for parallactic angle + final raMoon = math.atan2(moonGCRS.$2, moonGCRS.$1); + final decMoon = math.asin((moonGCRS.$3 / distance).clamp(-1.0, 1.0)); + + final era = computeERA(ts.jdUT1); + final ha = era + lon * deg2rad - raMoon; + + final parallacticAngle = math.atan2( + math.sin(ha), + math.cos(lat * deg2rad) * math.tan(decMoon) - + math.sin(lat * deg2rad) * math.cos(ha), + ); + + return MoonPosition( + azimuth: azAlt.azimuth, + altitude: azAlt.altitude, + distance: distance, + parallacticAngle: parallacticAngle, + ); +} + +// ── getMoonIllumination ────────────────────────────────────────────────────── + +/// Compute the Moon's illumination fraction, phase cycle, and bright limb angle. +/// +/// Works without a kernel (Meeus Ch. 47/48 approximation). +/// Accuracy: fraction ~0.5%, phase fraction ~0.003. +/// +/// ```dart +/// final illum = getMoonIllumination(); +/// print(illum.fraction); // e.g. 0.43 +/// ``` +MoonIlluminationResult getMoonIllumination([DateTime? date]) { + final d = date ?? DateTime.now(); + final ts = computeTimeScales(d); + final (:moonGCRS, :sunGCRS) = getMoonSunApproximate(ts.jdTT); + + final illum = computeIllumination(moonGCRS, sunGCRS); + + // Phase fraction: 0 = new, 0.25 = first quarter, 0.5 = full, 0.75 = last quarter + final phase = + illum.isWaxing + ? illum.elongationDeg / 360 + : 1 - illum.elongationDeg / 360; + + // Bright limb position angle + final moonDist = vnorm(moonGCRS); + final sunDist = vnorm(sunGCRS); + final raMoon = math.atan2(moonGCRS.$2, moonGCRS.$1); + final decMoon = math.asin((moonGCRS.$3 / moonDist).clamp(-1.0, 1.0)); + final raSun = math.atan2(sunGCRS.$2, sunGCRS.$1); + final decSun = math.asin((sunGCRS.$3 / sunDist).clamp(-1.0, 1.0)); + + final dRA = raSun - raMoon; + final angle = math.atan2( + math.cos(decSun) * math.sin(dRA), + math.sin(decSun) * math.cos(decMoon) - + math.cos(decSun) * math.sin(decMoon) * math.cos(dRA), + ); + + return MoonIlluminationResult( + fraction: illum.illumination, + phase: phase, + angle: angle, + isWaxing: illum.isWaxing, + ); +} + +// ── getMoonVisibilityEstimate ──────────────────────────────────────────────── + +/// Quick kernel-free crescent visibility estimate using the Odeh criterion. +/// +/// Computes approximate crescent geometry from Meeus positions at the given +/// observation time and applies the Odeh V-parameter formula. +/// +/// ```dart +/// final est = getMoonVisibilityEstimate(obsTime, 21.42, 39.83); +/// print(est.zone.label); // 'A' through 'D' +/// ``` +MoonVisibilityEstimate getMoonVisibilityEstimate( + DateTime? date, + double lat, + double lon, { + double elevation = 0, +}) { + final d = date ?? DateTime.now(); + _validateLatitude(lat, 'getMoonVisibilityEstimate'); + _validateLongitude(lon, 'getMoonVisibilityEstimate'); + + final ts = computeTimeScales(d); + final (:moonGCRS, :sunGCRS) = getMoonSunApproximate(ts.jdTT); + + // Airless positions (no refraction) + final moonAirless = computeAzAlt( + moonGCRS, + lat, + lon, + elevation, + ts, + airless: true, + ); + final sunAirless = computeAzAlt( + sunGCRS, + lat, + lon, + elevation, + ts, + airless: true, + ); + + // ARCL = elongation + final illumData = computeIllumination(moonGCRS, sunGCRS); + final arcl = illumData.elongationDeg; + + // ARCV = Moon airless alt minus Sun airless alt + final arcv = moonAirless.altitude - sunAirless.altitude; + + // Topocentric Moon vector for crescent width + final obsECEF = geodeticToECEF(lat, lon, elevation); + final obsITRS = (obsECEF.$1 / 1000, obsECEF.$2 / 1000, obsECEF.$3 / 1000); + final obsGCRS = itrsToGcrsSimple(obsITRS, ts); + final moonTopo = vsub(moonGCRS, obsGCRS); + + final (:w, wprime: _) = computeCrescentWidth(moonTopo, arcl); + + final v = arcv - arcvMinimum(w); + + final zone = + v >= odehThresholds[OdehZone.a]! + ? OdehZone.a + : v >= odehThresholds[OdehZone.b]! + ? OdehZone.b + : v >= odehThresholds[OdehZone.c]! + ? OdehZone.c + : OdehZone.d; + + return MoonVisibilityEstimate( + v: v, + zone: zone, + description: odehDescriptions[zone]!, + isVisibleNakedEye: zone == OdehZone.a, + isVisibleWithOpticalAid: zone == OdehZone.a || zone == OdehZone.b, + arcl: arcl, + arcv: arcv, + w: w, + moonAboveHorizon: moonAirless.altitude > 0, + ); +} + +// ── getMoon ────────────────────────────────────────────────────────────────── + +/// Combined kernel-free moon snapshot for a time and location. +/// +/// Calls [getMoonPhase], [getMoonPosition], [getMoonIllumination], and +/// [getMoonVisibilityEstimate] in a single request. +/// +/// ```dart +/// final moon = getMoon(DateTime.now(), 51.5, -0.1); +/// print(moon.phase.phaseName); +/// print(moon.visibility.zone.label); +/// ``` +MoonSnapshot getMoon( + DateTime? date, + double lat, + double lon, { + double elevation = 0, +}) { + final d = date ?? DateTime.now(); + _validateLatitude(lat, 'getMoon'); + _validateLongitude(lon, 'getMoon'); + + return MoonSnapshot( + phase: getMoonPhase(d), + position: getMoonPosition(d, lat, lon, elevation: elevation), + illumination: getMoonIllumination(d), + visibility: getMoonVisibilityEstimate(d, lat, lon, elevation: elevation), + ); +} diff --git a/lib/src/bodies.dart b/lib/src/bodies.dart new file mode 100644 index 0000000..2620bdc --- /dev/null +++ b/lib/src/bodies.dart @@ -0,0 +1,398 @@ +/// Moon and Sun positions using Meeus Ch. 25/47, illumination geometry. +/// +/// References: +/// Meeus, J. (1998). Astronomical Algorithms, 2nd ed. Willmann-Bell. +/// Odeh, M. (2006). New Criterion for Lunar Crescent Visibility. +library; + +import 'dart:math' as math; + +import 'math.dart'; +import 'time.dart'; + +/// AU in km. +const double _auKm = 149597870.7; + +/// Mean radius of the Moon in km (IAU 2015). +const double moonRadiusKm = 1737.4; + +/// Low-accuracy Sun and Moon positions using Meeus Ch. 25/47. +/// +/// Error budget: +/// Sun: < 0.01 deg in ecliptic longitude +/// Moon: < 0.3 deg in ecliptic longitude, < 0.2 deg in latitude +({Vec3 moonGCRS, Vec3 sunGCRS}) getMoonSunApproximate(double jdTT) { + final t = (jdTT - j2000) / daysPerJulianCentury; + + // Sun (Meeus Ch. 25) + final l0 = 280.46646 + 36000.76983 * t + 0.0003032 * t * t; + final mSun = 357.52911 + 35999.05029 * t - 0.0001537 * t * t; + final mSunRad = (mSun % 360) * deg2rad; + final eSun = 0.016708634 - 0.000042037 * t - 0.0000001267 * t * t; + + final c = + (1.914602 - 0.004817 * t - 0.000014 * t * t) * math.sin(mSunRad) + + (0.019993 - 0.000101 * t) * math.sin(2 * mSunRad) + + 0.000289 * math.sin(3 * mSunRad); + + final sunLonDeg = l0 + c; + final nuRad = mSunRad + c * deg2rad; + final rAU = (1.000001018 * (1 - eSun * eSun)) / (1 + eSun * math.cos(nuRad)); + final rKm = rAU * _auKm; + + final omega = (125.04 - 1934.136 * t) * deg2rad; + final sunLonApp = sunLonDeg - 0.00569 - 0.00478 * math.sin(omega); + final sunLonRad = sunLonApp * deg2rad; + + final eps = + (23.439291111 - + 0.013004167 * t - + 0.0000001638 * t * t + + 0.0000005036 * t * t * t) * + deg2rad; + + final sunGCRS = ( + rKm * math.cos(sunLonRad), + rKm * math.sin(sunLonRad) * math.cos(eps), + rKm * math.sin(sunLonRad) * math.sin(eps), + ); + + // Moon (Meeus Ch. 47) + final lp = + 218.3164477 + + 481267.88123421 * t - + 0.0015786 * t * t + + (t * t * t) / 538841 - + (t * t * t * t) / 65194000; + final d = + 297.8501921 + + 445267.1114034 * t - + 0.0018819 * t * t + + (t * t * t) / 545868 - + (t * t * t * t) / 113065000; + final m = + 357.5291092 + + 35999.0502909 * t - + 0.0001536 * t * t + + (t * t * t) / 24490000; + final mp = + 134.9633964 + + 477198.8675055 * t + + 0.0087414 * t * t + + (t * t * t) / 69699 - + (t * t * t * t) / 14712000; + final f = + 93.272095 + + 483202.0175233 * t - + 0.0036539 * t * t - + (t * t * t) / 3526000 + + (t * t * t * t) / 863310000; + + final a1 = (119.75 + 131.849 * t) * deg2rad; + final a2 = (53.09 + 479264.29 * t) * deg2rad; + final a3 = (313.45 + 481266.484 * t) * deg2rad; + + final dR = (d % 360) * deg2rad; + final mR = (m % 360) * deg2rad; + final mpR = (mp % 360) * deg2rad; + final fR = (f % 360) * deg2rad; + + final e = 1 - 0.002516 * t - 0.0000074 * t * t; + + // Longitude and distance: 30 terms from Meeus Table 47.A + // [d, m, mp, f, Sl (1e-6 deg), Sr (1e-3 km)] + const ld = <(int, int, int, int, int, int)>[ + (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), + ]; + + var sl = 0.0; + var sr = 0.0; + for (final (di, mi, mpi, fi, slc, src) in ld) { + final arg = di * dR + mi * mR + mpi * mpR + fi * fR; + final eCorr = + mi.abs() == 2 + ? e * e + : mi.abs() == 1 + ? e + : 1.0; + sl += slc * eCorr * math.sin(arg); + sr += src * eCorr * math.cos(arg); + } + + sl += + 3958 * math.sin(a1) + + 1962 * math.sin((lp - f) * deg2rad) + + 318 * math.sin(a2); + + // Latitude: 20 terms from Meeus Table 47.B + const fb = <(int, int, int, int, int)>[ + (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), + ]; + + var sb = 0.0; + for (final (di, mi, mpi, fi, sbc) in fb) { + final arg = di * dR + mi * mR + mpi * mpR + fi * fR; + final eCorr = + mi.abs() == 2 + ? e * e + : mi.abs() == 1 + ? e + : 1.0; + sb += sbc * eCorr * math.sin(arg); + } + + sb += + -2235 * math.sin(lp * deg2rad) + + 382 * math.sin(a3) + + 175 * math.sin(a1 - fR) + + 175 * math.sin(a1 + fR) + + 127 * math.sin((lp - mp) * deg2rad) - + 115 * math.sin((lp + mp) * deg2rad); + + final moonLonDeg = lp + sl * 1e-6; + final moonLatDeg = sb * 1e-6; + final moonDistKm = 385000.56 + sr * 0.001; + + final moonLonRad = moonLonDeg * deg2rad; + final moonLatRad = moonLatDeg * deg2rad; + + final moonGCRS = ( + moonDistKm * math.cos(moonLatRad) * math.cos(moonLonRad), + moonDistKm * + (math.cos(eps) * math.cos(moonLatRad) * math.sin(moonLonRad) - + math.sin(eps) * math.sin(moonLatRad)), + moonDistKm * + (math.sin(eps) * math.cos(moonLatRad) * math.sin(moonLonRad) + + math.cos(eps) * math.sin(moonLatRad)), + ); + + return (moonGCRS: moonGCRS, sunGCRS: sunGCRS); +} + +/// Estimate the nearest new moon JD using Meeus Ch. 49. +/// Accurate to within ~2 hours. +double nearestNewMoon(double jdTT) { + final y = 2000.0 + (jdTT - j2000) / 365.25; + final k = ((y - 2000.0) * 12.3685).round().toDouble(); + final t = k / 1236.85; + + var jde = + 2451550.09766 + + 29.530588861 * k + + 0.00015437 * t * t - + 0.00000015 * t * t * t + + 0.00000000073 * t * t * t * t; + + final mArg = + (2.5534 + 29.1053567 * k - 0.0000014 * t * t - 0.00000011 * t * t * t) * + deg2rad; + final mpArg = + (201.5643 + + 385.81693528 * k + + 0.0107582 * t * t + + 0.00001238 * t * t * t) * + deg2rad; + final fc = + (160.7108 + + 390.67050284 * k - + 0.0016118 * t * t - + 0.00000227 * t * t * t) * + deg2rad; + final om = + (124.7746 - 1.56375588 * k + 0.0020672 * t * t + 0.00000215 * t * t * t) * + deg2rad; + final eCorr = 1 - 0.002516 * t - 0.0000074 * t * t; + + jde += + -0.4072 * math.sin(mpArg) + + 0.17241 * eCorr * math.sin(mArg) + + 0.01608 * math.sin(2 * mpArg) + + 0.01039 * math.sin(2 * fc) + + 0.00739 * eCorr * math.sin(mpArg - mArg) - + 0.00514 * eCorr * math.sin(mpArg + mArg) + + 0.00208 * eCorr * eCorr * math.sin(2 * mArg) - + 0.00111 * math.sin(mpArg - 2 * fc) - + 0.00057 * math.sin(mpArg + 2 * fc) + + 0.00056 * eCorr * math.sin(2 * mpArg + mArg) - + 0.00042 * math.sin(3 * mpArg) + + 0.00042 * eCorr * math.sin(mArg + 2 * fc) + + 0.00038 * eCorr * math.sin(mArg - 2 * fc) - + 0.00024 * eCorr * math.sin(2 * mpArg - mArg) - + 0.00017 * math.sin(om) - + 0.00007 * math.sin(mpArg + 2 * mArg) + + 0.00004 * math.sin(2 * mpArg - 2 * fc) + + 0.00004 * math.sin(3 * mArg) + + 0.00003 * math.sin(mpArg + mArg - 2 * fc) + + 0.00003 * math.sin(2 * mpArg + 2 * fc) - + 0.00003 * math.sin(mpArg + mArg + 2 * fc) + + 0.00003 * math.sin(mpArg - mArg + 2 * fc) - + 0.00002 * math.sin(mpArg - mArg - 2 * fc) - + 0.00002 * math.sin(3 * mpArg + mArg) + + 0.00002 * math.sin(4 * mpArg); + + return jde; +} + +/// Estimate the nearest full moon JD using Meeus Ch. 49. +double nearestFullMoon(double jdTT) { + final y = 2000 + (jdTT - j2000) / 365.25; + final kBase = ((y - 2000.0) * 12.3685).round().toDouble(); + final k1 = kBase - 0.5; + final k2 = kBase + 0.5; + final jde1 = _fullMoonJDE(k1); + final jde2 = _fullMoonJDE(k2); + final d1 = (jde1 - jdTT).abs(); + final d2 = (jde2 - jdTT).abs(); + return d1 < d2 ? jde1 : jde2; +} + +double _fullMoonJDE(double k) { + final t = k / 1236.85; + var jde = + 2451550.09766 + + 29.530588861 * k + + 0.00015437 * t * t - + 0.00000015 * t * t * t + + 0.00000000073 * t * t * t * t; + + final mArg = (2.5534 + 29.1053567 * k - 0.0000014 * t * t) * deg2rad; + final mpArg = (201.5643 + 385.81693528 * k + 0.0107582 * t * t) * deg2rad; + final fc = (160.7108 + 390.67050284 * k - 0.0016118 * t * t) * deg2rad; + final om = (124.7746 - 1.56375588 * k + 0.0020672 * t * t) * deg2rad; + final eCorr = 1 - 0.002516 * t - 0.0000074 * t * t; + + jde += + -0.40614 * math.sin(mpArg) + + 0.17302 * eCorr * math.sin(mArg) + + 0.01614 * math.sin(2 * mpArg) + + 0.01043 * math.sin(2 * fc) + + 0.00734 * eCorr * math.sin(mpArg - mArg) - + 0.00515 * eCorr * math.sin(mpArg + mArg) + + 0.00209 * eCorr * eCorr * math.sin(2 * mArg) - + 0.00111 * math.sin(mpArg - 2 * fc) - + 0.00057 * math.sin(mpArg + 2 * fc) + + 0.00056 * eCorr * math.sin(2 * mpArg + mArg) - + 0.00042 * math.sin(3 * mpArg) + + 0.00042 * eCorr * math.sin(mArg + 2 * fc) + + 0.00038 * eCorr * math.sin(mArg - 2 * fc) - + 0.00024 * eCorr * math.sin(2 * mpArg - mArg) - + 0.00017 * math.sin(om) - + 0.00007 * math.sin(mpArg + 2 * mArg) + + 0.00004 * math.sin(2 * mpArg - 2 * fc) + + 0.00004 * math.sin(3 * mArg) + + 0.00003 * math.sin(mpArg + mArg - 2 * fc) + + 0.00003 * math.sin(2 * mpArg + 2 * fc) - + 0.00003 * math.sin(mpArg + mArg + 2 * fc) + + 0.00003 * math.sin(mpArg - mArg + 2 * fc) - + 0.00002 * math.sin(mpArg - mArg - 2 * fc) - + 0.00002 * math.sin(3 * mpArg + mArg) + + 0.00002 * math.sin(4 * mpArg); + + return jde; +} + +/// Compute Moon illumination quantities from geocentric positions. +/// +/// Returns illumination fraction [0-1], phase angle (deg), elongation (deg), +/// and whether the Moon is waxing. +({ + double illumination, + double phaseAngleDeg, + double elongationDeg, + bool isWaxing, +}) +computeIllumination(Vec3 moonGCRS, Vec3 sunGCRS) { + final rMoon = vnorm(moonGCRS); + final rSun = vnorm(sunGCRS); + + final cosElong = (vdot(moonGCRS, sunGCRS) / (rMoon * rSun)).clamp(-1.0, 1.0); + final elongationDeg = math.acos(cosElong) / deg2rad; + + final moonToSun = ( + sunGCRS.$1 - moonGCRS.$1, + sunGCRS.$2 - moonGCRS.$2, + sunGCRS.$3 - moonGCRS.$3, + ); + final moonToEarth = (-moonGCRS.$1, -moonGCRS.$2, -moonGCRS.$3); + final rMoonToSun = vnorm(moonToSun); + + final cosPhase = (vdot(moonToEarth, moonToSun) / (rMoon * rMoonToSun)).clamp( + -1.0, + 1.0, + ); + final phaseAngleDeg = math.acos(cosPhase) / deg2rad; + + final illumination = (1 + math.cos(phaseAngleDeg * deg2rad)) / 2; + + final crossZ = sunGCRS.$1 * moonGCRS.$2 - sunGCRS.$2 * moonGCRS.$1; + final isWaxing = crossZ > 0; + + return ( + illumination: illumination, + phaseAngleDeg: phaseAngleDeg, + elongationDeg: elongationDeg, + isWaxing: isWaxing, + ); +} + +/// Compute topocentric crescent width W in arc minutes. +/// +/// W = SD * (1 - cos(ARCL)) where SD is the topocentric semi-diameter. +({double w, double wprime}) computeCrescentWidth( + Vec3 moonTopoVec, + double arcl, +) { + final rMoon = vnorm(moonTopoVec); + final sdArcmin = (math.atan(moonRadiusKm / rMoon) / deg2rad) * 60; + final arclRad = arcl * deg2rad; + final w = sdArcmin * (1 - math.cos(arclRad)); + return (w: w, wprime: w); +} diff --git a/lib/src/math.dart b/lib/src/math.dart new file mode 100644 index 0000000..1ec1100 --- /dev/null +++ b/lib/src/math.dart @@ -0,0 +1,65 @@ +/// Vector math and angle utilities. +/// +/// All computation is pure (no I/O, no state). +library; + +import 'dart:math' as math; + +/// Degrees to radians conversion factor. +const double deg2rad = math.pi / 180; + +/// Radians to degrees conversion factor. +const double rad2deg = 180 / math.pi; + +/// A three-element vector (position or velocity). +typedef Vec3 = (double, double, double); + +/// Add two 3-vectors. +Vec3 vadd(Vec3 a, Vec3 b) => (a.$1 + b.$1, a.$2 + b.$2, a.$3 + b.$3); + +/// Subtract b from a. +Vec3 vsub(Vec3 a, Vec3 b) => (a.$1 - b.$1, a.$2 - b.$2, a.$3 - b.$3); + +/// Scale a 3-vector. +Vec3 vscale(Vec3 a, double s) => (a.$1 * s, a.$2 * s, a.$3 * s); + +/// Dot product. +double vdot(Vec3 a, Vec3 b) => a.$1 * b.$1 + a.$2 * b.$2 + a.$3 * b.$3; + +/// Euclidean norm. +double vnorm(Vec3 a) => math.sqrt(vdot(a, a)); + +/// Cross product. +Vec3 vcross(Vec3 a, Vec3 b) => ( + a.$2 * b.$3 - a.$3 * b.$2, + a.$3 * b.$1 - a.$1 * b.$3, + a.$1 * b.$2 - a.$2 * b.$1, +); + +/// Unit vector (normalized). Throws if zero vector. +Vec3 vunit(Vec3 a) { + final n = vnorm(a); + if (n == 0) throw RangeError('Cannot normalize a zero vector'); + return vscale(a, 1 / n); +} + +/// Angular separation between two direction vectors in radians. +double angularSep(Vec3 a, Vec3 b) { + final cosAngle = (vdot(vunit(a), vunit(b))).clamp(-1.0, 1.0); + return math.acos(cosAngle); +} + +/// Normalize an angle to [0, 2*pi). +double mod2pi(double angle) { + const twoPi = 2 * math.pi; + return ((angle % twoPi) + twoPi) % twoPi; +} + +/// Normalize an angle in degrees to [0, 360). +double mod360(double deg) => ((deg % 360) + 360) % 360; + +/// Normalize an angle in degrees to [-180, 180). +double normalizeDeg180(double deg) { + deg = mod360(deg); + return deg >= 180 ? deg - 360 : deg; +} diff --git a/lib/src/observer.dart b/lib/src/observer.dart new file mode 100644 index 0000000..368a3e7 --- /dev/null +++ b/lib/src/observer.dart @@ -0,0 +1,171 @@ +/// Geodetic transforms, topocentric ENU, atmospheric refraction, az/alt. +/// +/// Observer position chain: +/// geodetic (lat, lon, elev) -> ECEF (ITRS) -> GCRS -> topocentric ENU -> az/alt +/// +/// References: +/// WGS84: NIMA TR8350.2, 3rd edition +/// Bennett (1982), Astronomical Refraction +/// Saemundsson (1986), Sky & Telescope 72, 70 +library; + +import 'dart:math' as math; + +import 'math.dart'; +import 'time.dart'; +import 'types.dart'; + +/// Convert geodetic coordinates to ECEF position in meters. +Vec3 geodeticToECEF(double lat, double lon, double elev) { + final phi = lat * math.pi / 180; + final lam = lon * math.pi / 180; + final sinPhi = math.sin(phi); + final cosPhi = math.cos(phi); + final 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, + ); +} + +/// Compute ENU basis vectors at a geodetic location. +({Vec3 east, Vec3 north, Vec3 up}) computeENUBasis(double lat, double lon) { + final phi = lat * math.pi / 180; + final lam = lon * math.pi / 180; + final sinPhi = math.sin(phi); + final cosPhi = math.cos(phi); + final sinLam = math.sin(lam); + final cosLam = math.cos(lam); + + return ( + east: (-sinLam, cosLam, 0.0), + north: (-sinPhi * cosLam, -sinPhi * sinLam, cosPhi), + up: (cosPhi * cosLam, cosPhi * sinLam, sinPhi), + ); +} + +/// Convert ECEF displacement to local ENU. +Vec3 ecefToENU(Vec3 ecefDelta, double lat, double lon) { + final basis = computeENUBasis(lat, lon); + return ( + vdot(ecefDelta, basis.east), + vdot(ecefDelta, basis.north), + vdot(ecefDelta, basis.up), + ); +} + +/// Convert ENU vector to azimuth and altitude. +AzAlt enuToAzAlt(Vec3 enu) { + final e = enu.$1; + final n = enu.$2; + final u = enu.$3; + final horiz = math.sqrt(e * e + n * n); + final altitude = math.atan2(u, horiz) * 180 / math.pi; + var azimuth = math.atan2(e, n) * 180 / math.pi; + if (azimuth < 0) azimuth += 360; + return AzAlt(azimuth: azimuth, altitude: altitude); +} + +/// Compute Earth Rotation Angle from JD in UT1. +double computeERA(double jdUT1) { + final du = jdUT1 - 2451545.0; + final era = 2 * math.pi * (0.779057273264 + 1.0027378119113546 * du); + return ((era % (2 * math.pi)) + 2 * math.pi) % (2 * math.pi); +} + +/// Simplified GCRS to ITRS rotation (Earth rotation only, no nutation/precession). +/// +/// For lite mode: rotate GCRS by -GMST around the z-axis. +/// This is accurate to ~1 deg, sufficient for phase/visibility estimates. +Vec3 gcrsToItrsSimple(Vec3 gcrs, TimeScales ts) { + // Use ERA (Earth Rotation Angle) from UT1 as the rotation angle + final era = computeERA(ts.jdUT1); + final cosR = math.cos(era); + final sinR = math.sin(era); + // Rotation around z-axis by +ERA (GCRS -> ITRS) + return ( + cosR * gcrs.$1 + sinR * gcrs.$2, + -sinR * gcrs.$1 + cosR * gcrs.$2, + gcrs.$3, + ); +} + +/// Simplified ITRS to GCRS rotation (inverse of [gcrsToItrsSimple]). +Vec3 itrsToGcrsSimple(Vec3 itrs, TimeScales ts) { + final era = computeERA(ts.jdUT1); + final cosR = math.cos(era); + final sinR = math.sin(era); + return ( + cosR * itrs.$1 - sinR * itrs.$2, + sinR * itrs.$1 + cosR * itrs.$2, + itrs.$3, + ); +} + +/// Compute topocentric azimuth and altitude for a body. +/// +/// Pipeline: +/// 1. body GCRS -> body ITRS (simplified Earth rotation) +/// 2. observer geodetic -> observer ECEF (ITRS, m -> km) +/// 3. delta_ITRS = body_ITRS - observer_ITRS +/// 4. delta_ITRS -> ENU +/// 5. ENU -> az/alt +/// 6. Apply Bennett refraction if not airless +AzAlt computeAzAlt( + Vec3 bodyGCRS, + double lat, + double lon, + double elevation, + TimeScales ts, { + bool airless = false, + double pressure = 1013.25, + double temperature = 15, +}) { + final bodyITRS = gcrsToItrsSimple(bodyGCRS, ts); + + final obsECEF = geodeticToECEF(lat, lon, elevation); + final obsITRS = (obsECEF.$1 / 1000, obsECEF.$2 / 1000, obsECEF.$3 / 1000); + + final delta = vsub(bodyITRS, obsITRS); + final enu = ecefToENU(delta, lat, lon); + final azAlt = enuToAzAlt(enu); + + if (!airless) { + azAlt.altitude = applyRefraction(azAlt.altitude, pressure, temperature); + } + + return azAlt; +} + +/// Bennett (1982) atmospheric refraction correction in degrees. +/// +/// Formula: R = cot(h + 7.31/(h + 4.4)) / 60 +/// with pressure/temperature correction. +double bennettRefraction( + double altitudeDeg, { + double pressure = 1013.25, + double temperature = 15, +}) { + if (altitudeDeg < -1) return 0; + final h = altitudeDeg; + final argDeg = h + 7.31 / (h + 4.4); + final argRad = argDeg * math.pi / 180; + final r = 1 / (math.tan(argRad) * 60); + final corrected = r * (pressure / 1010) * (283 / (273 + temperature)); + return math.max(0, corrected); +} + +/// Apply refraction correction to an airless altitude. +double applyRefraction( + double airlessAlt, [ + double pressure = 1013.25, + double temperature = 15, +]) { + return airlessAlt + + bennettRefraction( + airlessAlt, + pressure: pressure, + temperature: temperature, + ); +} diff --git a/lib/src/time.dart b/lib/src/time.dart new file mode 100644 index 0000000..3679359 --- /dev/null +++ b/lib/src/time.dart @@ -0,0 +1,262 @@ +/// Julian Date, time scale conversions, and delta-T. +/// +/// Conversion chain: UTC -> TAI -> TT -> TDB +/// TAI = UTC + delta_AT (from leap-second table, integer seconds) +/// TT = TAI + 32.184s (exact, by definition) +/// TDB ~ TT + 0.001658*sin(g) + ... (sub-millisecond correction) +/// +/// References: +/// IERS Conventions (2010), Chapter 5 +/// NAIF LSK kernel (naif0012.tls) +/// Espenak & Meeus, delta-T polynomial expressions +library; + +import 'dart:math' as math; + +/// Julian Date of J2000.0 epoch (2000 Jan 1, 12:00 TT). +const double j2000 = 2451545.0; + +/// TT - TAI offset in seconds (exact, by definition). +const double ttMinusTai = 32.184; + +/// Seconds per day. +const double secondsPerDay = 86400.0; + +/// Days per Julian century. +const double daysPerJulianCentury = 36525.0; + +/// All relevant time scale values for a single moment. +class TimeScales { + /// The original UTC [DateTime]. + final DateTime utc; + + /// Julian Date in UTC. + final double jdUTC; + + /// Julian Date in Terrestrial Time. + final double jdTT; + + /// Julian Date in Barycentric Dynamical Time. + final double jdTDB; + + /// Julian Date in UT1. + final double jdUT1; + + /// TT - UT1 in seconds (delta-T). + final double deltaT; + + /// TAI - UTC in seconds (leap seconds count). + final double deltaAT; + + const TimeScales({ + required this.utc, + required this.jdUTC, + required this.jdTT, + required this.jdTDB, + required this.jdUT1, + required this.deltaT, + required this.deltaAT, + }); +} + +/// Leap-second table: [JD(UTC), delta_AT]. +/// Source: NAIF naif0012.tls. +const List<(double, double)> leapSecondTable = [ + (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 +]; + +/// Get the current leap second count (TAI - UTC) for a given JD in UTC. +double getDeltaAT(double jdUTC) { + double deltaAT = 10; + for (final (jd, dat) in leapSecondTable) { + if (jdUTC >= jd) { + deltaAT = dat; + } else { + break; + } + } + return deltaAT; +} + +/// Convert a [DateTime] (UTC) to Julian Date in UTC. +double dateToJD(DateTime date) { + return date.toUtc().millisecondsSinceEpoch / 86400000.0 + 2440587.5; +} + +/// Convert a Julian Date in UTC to a [DateTime]. +DateTime jdToDate(double jd) { + final ms = ((jd - 2440587.5) * 86400000.0).round(); + return DateTime.fromMillisecondsSinceEpoch(ms, isUtc: true); +} + +/// Julian centuries from J2000.0 (in TT). +double jdTTtoT(double jdTT) => (jdTT - j2000) / daysPerJulianCentury; + +/// Compute all relevant time scales for a given UTC [DateTime]. +TimeScales computeTimeScales( + DateTime utc, { + double? ut1utcOverride, + double? deltaTOverride, +}) { + final jdUTC = dateToJD(utc); + final deltaAT = getDeltaAT(jdUTC); + + final jdTAI = jdUTC + deltaAT / secondsPerDay; + final jdTT = jdTAI + ttMinusTai / secondsPerDay; + + final tdbCorrection = tdbMinusTT(jdTT) / secondsPerDay; + final jdTDB = jdTT + tdbCorrection; + + double jdUT1; + double deltaT; + + if (ut1utcOverride != null) { + jdUT1 = jdUTC + ut1utcOverride / secondsPerDay; + deltaT = (jdTT - jdUT1) * secondsPerDay; + } else if (deltaTOverride != null) { + deltaT = deltaTOverride; + jdUT1 = jdTT - deltaT / secondsPerDay; + } else { + deltaT = deltaTPolynomial(jdTT); + jdUT1 = jdTT - deltaT / secondsPerDay; + } + + return TimeScales( + utc: utc, + jdUTC: jdUTC, + jdTT: jdTT, + jdTDB: jdTDB, + jdUT1: jdUT1, + deltaT: deltaT, + deltaAT: deltaAT, + ); +} + +/// Approximate TDB - TT in seconds. +/// +/// TDB - TT = 0.001658 * sin(g) + 0.000014 * sin(2g) +/// where g = 357.53 + 0.9856003 * (JD_TT - 2451545.0) +double tdbMinusTT(double jdTT) { + final d = jdTT - j2000; + final gDeg = 357.53 + 0.9856003 * d; + final 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 (2009). +double deltaTPolynomial(double jdTT) { + final y = 2000 + (jdTT - j2000) / 365.25; + + if (y < -500) { + final u = (y - 1820) / 100; + return -20 + 32 * u * u; + } else if (y < 500) { + final u = y / 100; + return 10583.6 - + 1014.41 * u + + 33.78311 * u * u - + 5.952053 * u * u * u - + 0.1798452 * math.pow(u, 4) + + 0.022174192 * math.pow(u, 5) + + 0.0090316521 * math.pow(u, 6); + } else if (y < 1600) { + final u = (y - 1000) / 100; + return 1574.2 - + 556.01 * u + + 71.23472 * u * u + + 0.319781 * math.pow(u, 3) - + 0.8503463 * math.pow(u, 4) - + 0.005050998 * math.pow(u, 5) + + 0.0083572073 * math.pow(u, 6); + } else if (y < 1700) { + final t = y - 1600; + return 120 - 0.9808 * t - 0.01532 * t * t + math.pow(t, 3) / 7129; + } else if (y < 1800) { + final t = y - 1700; + return 8.83 + + 0.1603 * t - + 0.0059285 * t * t + + 0.00013336 * math.pow(t, 3) - + math.pow(t, 4) / 1174000; + } else if (y < 1860) { + final t = y - 1800; + return 13.72 - + 0.332447 * t + + 0.0068612 * t * t + + 0.0041116 * math.pow(t, 3) - + 0.00037436 * math.pow(t, 4) + + 0.0000121272 * math.pow(t, 5) - + 0.0000001699 * math.pow(t, 6) + + 0.000000000875 * math.pow(t, 7); + } else if (y < 1900) { + final t = y - 1860; + return 7.62 + + 0.5737 * t - + 0.251754 * t * t + + 0.01680668 * math.pow(t, 3) - + 0.0004473624 * math.pow(t, 4) + + math.pow(t, 5) / 233174; + } else if (y < 1920) { + final t = y - 1900; + return -2.79 + + 1.494119 * t - + 0.0598939 * t * t + + 0.0061966 * math.pow(t, 3) - + 0.000197 * math.pow(t, 4); + } else if (y < 1941) { + final t = y - 1920; + return 21.2 + 0.84493 * t - 0.0761 * t * t + 0.0020936 * math.pow(t, 3); + } else if (y < 1961) { + final t = y - 1950; + return 29.07 + 0.407 * t - (t * t) / 233 + math.pow(t, 3) / 2547; + } else if (y < 1986) { + final t = y - 1975; + return 45.45 + 1.067 * t - (t * t) / 260 - math.pow(t, 3) / 718; + } else if (y < 2005) { + final t = y - 2000; + return 63.86 + + 0.3345 * t - + 0.060374 * t * t + + 0.0017275 * math.pow(t, 3) + + 0.000651814 * math.pow(t, 4) + + 0.00002373599 * math.pow(t, 5); + } else if (y < 2050) { + final t = y - 2000; + return 62.92 + 0.32217 * t + 0.005589 * t * t; + } else if (y < 2150) { + return -20 + 32 * math.pow((y - 1820) / 100, 2) - 0.5628 * (2150 - y); + } else { + final u = (y - 1820) / 100; + return -20 + 32 * u * u; + } +} diff --git a/lib/src/types.dart b/lib/src/types.dart new file mode 100644 index 0000000..1dedf80 --- /dev/null +++ b/lib/src/types.dart @@ -0,0 +1,370 @@ +/// All result types for the moon_sighting package. +library; + +/// Azimuth + altitude in degrees. +class AzAlt { + /// Degrees from North, measured clockwise (0=N, 90=E, 180=S, 270=W). + double azimuth; + + /// Degrees above the horizon (negative = below). + double altitude; + + AzAlt({required this.azimuth, required this.altitude}); +} + +/// Topocentric moon position. +/// +/// Computed via Meeus Ch. 47 (no kernel required). +/// Accuracy: azimuth/altitude ~0.3 deg, distance ~300 km. +class MoonPosition { + /// Azimuth in degrees from North, clockwise. + final double azimuth; + + /// Apparent altitude in degrees (atmospheric refraction applied). + final double altitude; + + /// Distance from Earth center to Moon center, km. + final double distance; + + /// Parallactic angle in radians. + final double parallacticAngle; + + const MoonPosition({ + required this.azimuth, + required this.altitude, + required this.distance, + required this.parallacticAngle, + }); +} + +/// Moon illumination result. +/// +/// Computed via Meeus Ch. 47/48 (no kernel required). +class MoonIlluminationResult { + /// Illuminated fraction 0 (new) to 1 (full). + final double fraction; + + /// Phase cycle fraction in [0, 1): + /// 0 = new, 0.25 = first quarter, 0.5 = full, 0.75 = last quarter. + final double phase; + + /// Position angle of the bright limb midpoint in radians. + final double angle; + + /// True while elongation is increasing (new toward full). + final bool isWaxing; + + const MoonIlluminationResult({ + required this.fraction, + required this.phase, + required this.angle, + required this.isWaxing, + }); +} + +/// Named moon phase identifiers. +enum MoonPhaseName { + newMoon('new-moon'), + waxingCrescent('waxing-crescent'), + firstQuarter('first-quarter'), + waxingGibbous('waxing-gibbous'), + fullMoon('full-moon'), + waningGibbous('waning-gibbous'), + lastQuarter('last-quarter'), + waningCrescent('waning-crescent'); + + final String id; + const MoonPhaseName(this.id); +} + +/// Moon phase result from [getMoonPhase]. +class MoonPhaseResult { + /// Named phase. + final MoonPhaseName phase; + + /// Human-readable name, e.g. "Waxing Crescent". + final String phaseName; + + /// Moon phase emoji symbol. + final String phaseSymbol; + + /// Illuminated fraction 0-100 (percent). + final double illumination; + + /// Hours since last new moon. + final double age; + + /// Ecliptic longitude of Moon minus Sun, degrees [0, 360). + final double elongationDeg; + + /// True when illumination is increasing. + final bool isWaxing; + + /// UTC date of the next new moon. + final DateTime nextNewMoon; + + /// UTC date of the next full moon. + final DateTime nextFullMoon; + + /// UTC date of the previous new moon. + final DateTime prevNewMoon; + + const MoonPhaseResult({ + required this.phase, + required this.phaseName, + required this.phaseSymbol, + required this.illumination, + required this.age, + required this.elongationDeg, + required this.isWaxing, + required this.nextNewMoon, + required this.nextFullMoon, + required this.prevNewMoon, + }); +} + +/// Yallop q-test visibility category (NAO Technical Note 69). +enum YallopCategory { + /// Easily visible to the naked eye. + a('A'), + + /// Visible under perfect conditions. + b('B'), + + /// May need optical aid to find; naked eye possible. + c('C'), + + /// Optical aid needed; naked eye not possible. + d('D'), + + /// Not visible even with telescope under good conditions. + e('E'), + + /// Below Danjon limit. + f('F'); + + final String label; + const YallopCategory(this.label); +} + +/// Yallop q-test result. +class YallopResult { + /// The continuous q parameter (higher = more visible). + final double q; + + /// Visibility category A through F. + final YallopCategory category; + + /// Human-readable interpretation. + final String description; + + /// True for categories A and B. + final bool isVisibleNakedEye; + + /// True for categories C and D. + final bool requiresOpticalAid; + + /// True for category F. + final bool isBelowDanjonLimit; + + /// Topocentric crescent width W' used in the q formula, arc minutes. + final double wprime; + + const YallopResult({ + required this.q, + required this.category, + required this.description, + required this.isVisibleNakedEye, + required this.requiresOpticalAid, + required this.isBelowDanjonLimit, + required this.wprime, + }); +} + +/// Odeh visibility zone (Experimental Astronomy 2006). +enum OdehZone { + /// Visible with naked eye. + a('A'), + + /// Visible with optical aid; may be seen with naked eye. + b('B'), + + /// Visible with optical aid only. + c('C'), + + /// Not visible even with optical aid. + d('D'); + + final String label; + const OdehZone(this.label); +} + +/// Odeh criterion result. +class OdehResult { + /// Continuous V parameter: V = ARCV - f(W). Positive = exceeds threshold. + final double v; + + /// Visibility zone A through D. + final OdehZone zone; + + /// Human-readable interpretation. + final String description; + + /// True for zone A. + final bool isVisibleNakedEye; + + /// True for zones A and B. + final bool isVisibleWithOpticalAid; + + const OdehResult({ + required this.v, + required this.zone, + required this.description, + required this.isVisibleNakedEye, + required this.isVisibleWithOpticalAid, + }); +} + +/// Crescent geometry quantities. +class CrescentGeometry { + /// Arc of light: Sun-Moon angular separation (elongation), degrees. + final double arcl; + + /// Arc of vision: Moon airless alt minus Sun airless alt, degrees. + final double arcv; + + /// Relative azimuth: Sun az minus Moon az, degrees. + final double daz; + + /// Topocentric crescent width in arc minutes. + final double w; + + const CrescentGeometry({ + required this.arcl, + required this.arcv, + required this.daz, + required this.w, + }); +} + +/// Kernel-free Odeh-based crescent visibility estimate. +class MoonVisibilityEstimate { + /// Odeh V parameter. + final double v; + + /// Visibility zone A through D. + final OdehZone zone; + + /// Human-readable zone description. + final String description; + + /// True for zone A. + final bool isVisibleNakedEye; + + /// True for zones A and B. + final bool isVisibleWithOpticalAid; + + /// Arc of light (Sun-Moon elongation) in degrees. + final double arcl; + + /// Arc of vision (Moon airless alt minus Sun airless alt) in degrees. + final double arcv; + + /// Topocentric crescent width in arc minutes. + final double w; + + /// True when Moon is above the horizon at the given time. + final bool moonAboveHorizon; + + /// Always true: computed via Meeus approximation. + final bool isApproximate = true; + + MoonVisibilityEstimate({ + required this.v, + required this.zone, + required this.description, + required this.isVisibleNakedEye, + required this.isVisibleWithOpticalAid, + required this.arcl, + required this.arcv, + required this.w, + required this.moonAboveHorizon, + }); +} + +/// Combined kernel-free moon snapshot. +class MoonSnapshot { + /// Phase name, illumination, age, and next events. + final MoonPhaseResult phase; + + /// Topocentric az/alt, distance, parallactic angle. + final MoonPosition position; + + /// Illumination fraction, phase cycle, bright limb angle. + final MoonIlluminationResult illumination; + + /// Quick Odeh-based crescent visibility estimate. + final MoonVisibilityEstimate visibility; + + const MoonSnapshot({ + required this.phase, + required this.position, + required this.illumination, + required this.visibility, + }); +} + +/// Published Yallop q thresholds. +const Map yallopThresholds = { + YallopCategory.a: 0.216, + YallopCategory.b: -0.014, + YallopCategory.c: -0.16, + YallopCategory.d: -0.232, + YallopCategory.e: -0.293, +}; + +/// Yallop category descriptions. +const Map yallopDescriptions = { + YallopCategory.a: 'Easily visible to the naked eye', + YallopCategory.b: 'Visible under perfect conditions', + YallopCategory.c: 'May need optical aid to find; naked eye possible', + YallopCategory.d: 'Optical aid needed; naked eye not possible', + YallopCategory.e: 'Not visible even with telescope under good conditions', + YallopCategory.f: 'Below Danjon limit — crescent cannot form', +}; + +/// Published Odeh V thresholds. +const Map odehThresholds = { + OdehZone.a: 5.65, + OdehZone.b: 2.0, + OdehZone.c: -0.96, +}; + +/// Odeh zone descriptions. +const Map odehDescriptions = { + OdehZone.a: 'Visible with naked eye', + OdehZone.b: + 'Visible with optical aid; may be seen with naked eye under excellent conditions', + OdehZone.c: 'Visible with optical aid only', + OdehZone.d: 'Not visible even with optical aid', +}; + +/// WGS84 reference ellipsoid parameters. +class WGS84 { + WGS84._(); + + /// Semi-major axis in meters. + static const double a = 6378137.0; + + /// Inverse flattening. + static const double invF = 298.257223563; + + /// Flattening. + static const double f = 1 / invF; + + /// Semi-minor axis in meters. + static final double b = a * (1 - f); + + /// First eccentricity squared. + static final double e2 = 2 * f - f * f; +} diff --git a/lib/src/visibility.dart b/lib/src/visibility.dart new file mode 100644 index 0000000..98da6be --- /dev/null +++ b/lib/src/visibility.dart @@ -0,0 +1,106 @@ +/// Crescent visibility criteria: Yallop and Odeh. +/// +/// Both criteria transform the geometric quantities (ARCL, ARCV, DAZ, W) +/// into a single score that maps to a visibility category. +/// +/// References: +/// Yallop (1997), NAO Technical Note No. 69 +/// Odeh (2006), Experimental Astronomy 18(1), 39-64 +library; + +import 'math.dart'; +import 'bodies.dart'; +import 'types.dart'; + +/// The polynomial ARCV minimum as a function of crescent width W (arc minutes). +/// +/// 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, derived +/// empirically from historical observations. +double arcvMinimum(double w) { + return 11.8371 - 6.3226 * w + 0.7319 * w * w - 0.1018 * w * w * w; +} + +/// Compute the Yallop q parameter. +/// +/// q = (ARCV - arcv_min(W')) / 10 +double computeYallopQ(double arcv, double wprime) { + return (arcv - arcvMinimum(wprime)) / 10; +} + +/// Map a q value to Yallop category. +YallopCategory yallopCategoryFromQ(double q) { + if (q > yallopThresholds[YallopCategory.a]!) return YallopCategory.a; + if (q > yallopThresholds[YallopCategory.b]!) return YallopCategory.b; + if (q > yallopThresholds[YallopCategory.c]!) return YallopCategory.c; + if (q > yallopThresholds[YallopCategory.d]!) return YallopCategory.d; + if (q > yallopThresholds[YallopCategory.e]!) return YallopCategory.e; + return YallopCategory.f; +} + +/// Compute the full Yallop result from crescent geometry. +YallopResult computeYallop(CrescentGeometry geometry, double wprime) { + final q = computeYallopQ(geometry.arcv, wprime); + final category = yallopCategoryFromQ(q); + + return YallopResult( + q: q, + category: category, + description: yallopDescriptions[category]!, + isVisibleNakedEye: + category == YallopCategory.a || category == YallopCategory.b, + requiresOpticalAid: + category == YallopCategory.c || category == YallopCategory.d, + isBelowDanjonLimit: category == YallopCategory.f, + wprime: wprime, + ); +} + +/// Compute the Odeh V parameter. +/// +/// V = ARCV - arcv_min(W) +double computeOdehV(double arcv, double w) { + return arcv - arcvMinimum(w); +} + +/// Map a V value to the Odeh zone. +OdehZone odehZoneFromV(double v) { + if (v >= odehThresholds[OdehZone.a]!) return OdehZone.a; + if (v >= odehThresholds[OdehZone.b]!) return OdehZone.b; + if (v >= odehThresholds[OdehZone.c]!) return OdehZone.c; + return OdehZone.d; +} + +/// Compute the full Odeh result from crescent geometry. +OdehResult computeOdeh(CrescentGeometry geometry) { + final v = computeOdehV(geometry.arcv, geometry.w); + final zone = odehZoneFromV(v); + + return OdehResult( + v: v, + zone: zone, + description: odehDescriptions[zone]!, + isVisibleNakedEye: zone == OdehZone.a, + isVisibleWithOpticalAid: zone == OdehZone.a || zone == OdehZone.b, + ); +} + +/// Compute crescent geometry quantities from airless topocentric positions. +CrescentGeometry computeCrescentGeometry({ + required AzAlt moonAirless, + required AzAlt sunAirless, + required Vec3 moonTopo, + required Vec3 sunTopo, +}) { + final arcv = moonAirless.altitude - sunAirless.altitude; + + var daz = sunAirless.azimuth - moonAirless.azimuth; + if (daz > 180) daz -= 360; + if (daz < -180) daz += 360; + + final arcl = angularSep(moonTopo, sunTopo) * (180 / 3.141592653589793); + final (:w, wprime: _) = computeCrescentWidth(moonTopo, arcl); + + return CrescentGeometry(arcl: arcl, arcv: arcv, daz: daz, w: w); +} diff --git a/pubspec.yaml b/pubspec.yaml new file mode 100644 index 0000000..b56e350 --- /dev/null +++ b/pubspec.yaml @@ -0,0 +1,21 @@ +name: moon_sighting +description: > + Lunar crescent visibility for Dart and Flutter. Moon phase, position, + illumination, and Yallop/Odeh visibility criteria using Meeus algorithms. + Zero dependencies. +version: 1.0.0 +repository: https://github.com/acamarata/moon-sighting-dart +issue_tracker: https://github.com/acamarata/moon-sighting-dart/issues +topics: + - moon + - lunar + - crescent + - astronomy + - islamic + +environment: + sdk: ^3.7.0 + +dev_dependencies: + lints: ^5.0.0 + test: ^1.25.8 diff --git a/test/moon_sighting_test.dart b/test/moon_sighting_test.dart new file mode 100644 index 0000000..542c50f --- /dev/null +++ b/test/moon_sighting_test.dart @@ -0,0 +1,441 @@ +import 'package:moon_sighting/moon_sighting.dart'; +import 'package:test/test.dart'; + +void main() { + // ── Constants ──────────────────────────────────────────────────────────── + + group('Constants', () { + test('Yallop threshold A is 0.216', () { + expect(yallopThresholds[YallopCategory.a], equals(0.216)); + }); + + test('Yallop threshold E is -0.293', () { + expect(yallopThresholds[YallopCategory.e], equals(-0.293)); + }); + + test('All Yallop thresholds are defined', () { + for (final cat in [ + YallopCategory.a, + YallopCategory.b, + YallopCategory.c, + YallopCategory.d, + YallopCategory.e, + ]) { + expect(yallopThresholds[cat], isNotNull); + } + }); + + test('Yallop thresholds descend A > B > C > D > E', () { + expect( + yallopThresholds[YallopCategory.a]!, + greaterThan(yallopThresholds[YallopCategory.b]!), + ); + expect( + yallopThresholds[YallopCategory.b]!, + greaterThan(yallopThresholds[YallopCategory.c]!), + ); + expect( + yallopThresholds[YallopCategory.c]!, + greaterThan(yallopThresholds[YallopCategory.d]!), + ); + expect( + yallopThresholds[YallopCategory.d]!, + greaterThan(yallopThresholds[YallopCategory.e]!), + ); + }); + + test('Odeh threshold A is 5.65', () { + expect(odehThresholds[OdehZone.a], equals(5.65)); + }); + + test('Odeh threshold C is -0.96', () { + expect(odehThresholds[OdehZone.c], equals(-0.96)); + }); + + test('Odeh thresholds descend A > B > C', () { + expect( + odehThresholds[OdehZone.a]!, + greaterThan(odehThresholds[OdehZone.b]!), + ); + expect( + odehThresholds[OdehZone.b]!, + greaterThan(odehThresholds[OdehZone.c]!), + ); + }); + + test('WGS84.a is 6378137.0', () { + expect(WGS84.a, equals(6378137.0)); + }); + + test('WGS84.e2 is positive and < 1', () { + expect(WGS84.e2, greaterThan(0)); + expect(WGS84.e2, lessThan(1)); + }); + + test('WGS84.b < WGS84.a (oblate spheroid)', () { + expect(WGS84.b, lessThan(WGS84.a)); + }); + + test('Yallop descriptions are non-empty strings', () { + for (final cat in YallopCategory.values) { + expect(yallopDescriptions[cat], isNotEmpty); + } + }); + + test('Odeh descriptions are non-empty strings', () { + for (final zone in OdehZone.values) { + expect(odehDescriptions[zone], isNotEmpty); + } + }); + }); + + // ── getMoonPhase ─────────────────────────────────────────────────────── + + final dateMarch1 = DateTime.utc(2025, 3, 1, 12); + final phaseMarch1 = getMoonPhase(dateMarch1); + + group('getMoonPhase structure', () { + test('illumination is in [0, 100]', () { + expect(phaseMarch1.illumination, greaterThanOrEqualTo(0)); + expect(phaseMarch1.illumination, lessThanOrEqualTo(100)); + }); + + test('age is >= 0', () { + expect(phaseMarch1.age, greaterThanOrEqualTo(0)); + }); + + test('elongationDeg is in [0, 180]', () { + expect(phaseMarch1.elongationDeg, greaterThanOrEqualTo(0)); + expect(phaseMarch1.elongationDeg, lessThanOrEqualTo(180)); + }); + + test('nextNewMoon is a DateTime', () { + expect(phaseMarch1.nextNewMoon, isA()); + }); + + test('prevNewMoon is before reference date', () { + expect(phaseMarch1.prevNewMoon.isBefore(dateMarch1), isTrue); + }); + + test('nextNewMoon is after prevNewMoon', () { + expect(phaseMarch1.nextNewMoon.isAfter(phaseMarch1.prevNewMoon), isTrue); + }); + }); + + group('getMoonPhase phase boundaries', () { + test('near full moon: illumination > 85%', () { + final phaseFull = getMoonPhase(DateTime.utc(2025, 3, 14, 12)); + expect(phaseFull.illumination, greaterThan(85)); + }); + + test('near full moon: elongation > 120 deg', () { + final phaseFull = getMoonPhase(DateTime.utc(2025, 3, 14, 12)); + expect(phaseFull.elongationDeg, greaterThan(120)); + }); + + test('near new moon: illumination < 10%', () { + final phaseNew = getMoonPhase(DateTime.utc(2025, 3, 29, 12)); + expect(phaseNew.illumination, lessThan(10)); + }); + + test('near new moon: elongation < 30 deg', () { + final phaseNew = getMoonPhase(DateTime.utc(2025, 3, 29, 12)); + expect(phaseNew.elongationDeg, lessThan(30)); + }); + }); + + group('getMoonPhase consistency', () { + test('5 days after new moon: isWaxing = true', () { + expect(getMoonPhase(DateTime.utc(2025, 3, 5, 12)).isWaxing, isTrue); + }); + + test('6 days after full moon: isWaxing = false', () { + expect(getMoonPhase(DateTime.utc(2025, 3, 20, 12)).isWaxing, isFalse); + }); + + test('default date (now) returns valid result', () { + final nowPhase = getMoonPhase(); + expect(nowPhase.illumination, greaterThanOrEqualTo(0)); + expect(nowPhase.illumination, lessThanOrEqualTo(100)); + }); + + test('synodic month duration is ~29.5 days', () { + final synodicMs = + phaseMarch1.nextNewMoon + .difference(phaseMarch1.prevNewMoon) + .inMilliseconds; + final synodicDays = synodicMs / 86400000; + expect(synodicDays, greaterThan(29.0)); + expect(synodicDays, lessThan(30.1)); + }); + }); + + group('getMoonPhase phaseName and phaseSymbol', () { + test('waxing crescent: phaseName is Waxing Crescent', () { + final p = getMoonPhase(DateTime.utc(2025, 3, 5, 12)); + expect(p.phaseName, equals('Waxing Crescent')); + }); + + test('phaseName is a valid string', () { + final validNames = { + 'New Moon', + 'Waxing Crescent', + 'First Quarter', + 'Waxing Gibbous', + 'Full Moon', + 'Waning Gibbous', + 'Last Quarter', + 'Waning Crescent', + }; + final p = getMoonPhase(dateMarch1); + expect(validNames.contains(p.phaseName), isTrue); + }); + }); + + // ── getMoonPosition ──────────────────────────────────────────────────── + + group('getMoonPosition', () { + final moonPos = getMoonPosition( + DateTime.utc(2025, 3, 14, 20), + 51.5074, + -0.1278, + elevation: 10, + ); + + test('azimuth in [0, 360)', () { + expect(moonPos.azimuth, greaterThanOrEqualTo(0)); + expect(moonPos.azimuth, lessThan(360)); + }); + + test('altitude in [-90, 90]', () { + expect(moonPos.altitude, greaterThanOrEqualTo(-90)); + expect(moonPos.altitude, lessThanOrEqualTo(90)); + }); + + test('distance in lunar orbit range [356000, 407000] km', () { + expect(moonPos.distance, greaterThanOrEqualTo(356000)); + expect(moonPos.distance, lessThanOrEqualTo(407000)); + }); + + test('parallacticAngle is finite', () { + expect(moonPos.parallacticAngle.isFinite, isTrue); + }); + + test('default date returns valid result', () { + final pos = getMoonPosition(null, 21.4225, 39.8262); + expect(pos.azimuth, greaterThanOrEqualTo(0)); + expect(pos.azimuth, lessThan(360)); + expect(pos.distance, greaterThan(350000)); + expect(pos.distance, lessThan(410000)); + }); + }); + + // ── getMoonIllumination ──────────────────────────────────────────────── + + group('getMoonIllumination', () { + final illumFull = getMoonIllumination(DateTime.utc(2025, 3, 14, 12)); + final illumNew = getMoonIllumination(DateTime.utc(2025, 3, 29, 12)); + final illumWaxing = getMoonIllumination(DateTime.utc(2025, 3, 5, 12)); + + test('near full moon: fraction > 0.85', () { + expect(illumFull.fraction, greaterThan(0.85)); + }); + + test('near full moon: phase close to 0.5', () { + expect(illumFull.phase, greaterThan(0.4)); + expect(illumFull.phase, lessThan(0.6)); + }); + + test('near new moon: fraction < 0.05', () { + expect(illumNew.fraction, lessThan(0.05)); + }); + + test('near new moon: phase close to 0 or 1', () { + expect(illumNew.phase < 0.08 || illumNew.phase > 0.92, isTrue); + }); + + test('waxing: isWaxing = true', () { + expect(illumWaxing.isWaxing, isTrue); + }); + + test('fraction in [0, 1]', () { + expect(illumFull.fraction, greaterThanOrEqualTo(0)); + expect(illumFull.fraction, lessThanOrEqualTo(1)); + expect(illumNew.fraction, greaterThanOrEqualTo(0)); + expect(illumNew.fraction, lessThanOrEqualTo(1)); + }); + + test('phase in [0, 1)', () { + expect(illumFull.phase, greaterThanOrEqualTo(0)); + expect(illumFull.phase, lessThan(1)); + }); + + test('angle is finite', () { + expect(illumFull.angle.isFinite, isTrue); + }); + + test('default date returns valid result', () { + final illum = getMoonIllumination(); + expect(illum.fraction, greaterThanOrEqualTo(0)); + expect(illum.fraction, lessThanOrEqualTo(1)); + }); + }); + + // ── getMoonVisibilityEstimate ────────────────────────────────────────── + + group('getMoonVisibilityEstimate', () { + final vis = getMoonVisibilityEstimate( + DateTime.utc(2025, 3, 2, 18, 30), + 51.5074, + -0.1278, + elevation: 10, + ); + + test('zone is A, B, C, or D', () { + expect(OdehZone.values.contains(vis.zone), isTrue); + }); + + test('V is finite', () { + expect(vis.v.isFinite, isTrue); + }); + + test('ARCL is in [0, 180]', () { + expect(vis.arcl, greaterThanOrEqualTo(0)); + expect(vis.arcl, lessThanOrEqualTo(180)); + }); + + test('W >= 0', () { + expect(vis.w, greaterThanOrEqualTo(0)); + }); + + test('isApproximate is true', () { + expect(vis.isApproximate, isTrue); + }); + + test('isVisibleNakedEye matches zone A', () { + expect(vis.isVisibleNakedEye, equals(vis.zone == OdehZone.a)); + }); + + test('isVisibleWithOpticalAid matches zone A or B', () { + expect( + vis.isVisibleWithOpticalAid, + equals(vis.zone == OdehZone.a || vis.zone == OdehZone.b), + ); + }); + + test('near new moon: zone is D or C', () { + final nearNew = getMoonVisibilityEstimate( + DateTime.utc(2025, 3, 29, 18), + 21.4225, + 39.8262, + ); + expect(nearNew.zone == OdehZone.c || nearNew.zone == OdehZone.d, isTrue); + }); + }); + + // ── getMoon ──────────────────────────────────────────────────────────── + + group('getMoon', () { + final moon = getMoon( + DateTime.utc(2025, 3, 5, 20), + 51.5074, + -0.1278, + elevation: 10, + ); + + test('returns object with phase, position, illumination, visibility', () { + expect(moon.phase, isNotNull); + expect(moon.position, isNotNull); + expect(moon.illumination, isNotNull); + expect(moon.visibility, isNotNull); + }); + + test('phase is consistent with getMoonPhase standalone', () { + final standalone = getMoonPhase(DateTime.utc(2025, 3, 5, 20)); + expect(moon.phase.phase, equals(standalone.phase)); + expect(moon.phase.phaseName, equals(standalone.phaseName)); + }); + + test('illumination.isWaxing matches phase.isWaxing', () { + expect(moon.illumination.isWaxing, equals(moon.phase.isWaxing)); + }); + + test('visibility.isApproximate is true', () { + expect(moon.visibility.isApproximate, isTrue); + }); + + test('position has valid azimuth and altitude', () { + expect(moon.position.azimuth, greaterThanOrEqualTo(0)); + expect(moon.position.azimuth, lessThan(360)); + expect(moon.position.altitude, greaterThanOrEqualTo(-90)); + expect(moon.position.altitude, lessThanOrEqualTo(90)); + }); + }); + + // ── Input validation ────────────────────────────────────────────────── + + group('Input validation', () { + test('getMoonPosition rejects latitude out of range', () { + expect(() => getMoonPosition(DateTime.now(), 91, 0), throwsRangeError); + }); + + test('getMoonPosition rejects longitude out of range', () { + expect(() => getMoonPosition(DateTime.now(), 0, 181), throwsRangeError); + }); + + test('getMoonPosition rejects NaN latitude', () { + expect( + () => getMoonPosition(DateTime.now(), double.nan, 0), + throwsRangeError, + ); + }); + + test('getMoonVisibilityEstimate rejects invalid coordinates', () { + expect( + () => getMoonVisibilityEstimate(DateTime.now(), -91, 0), + throwsRangeError, + ); + }); + + test('getMoon rejects invalid coordinates', () { + expect(() => getMoon(DateTime.now(), 0, 200), throwsRangeError); + }); + }); + + // ── arcvMinimum polynomial ───────────────────────────────────────────── + + group('arcvMinimum polynomial', () { + test('arcvMinimum(0) = 11.8371', () { + expect(arcvMinimum(0), closeTo(11.8371, 0.0001)); + }); + + test('arcvMinimum decreases with moderate W', () { + expect(arcvMinimum(1), lessThan(arcvMinimum(0))); + }); + }); + + // ── nearestNewMoon accuracy ──────────────────────────────────────────── + + group('nearestNewMoon accuracy', () { + test('2024-04-08 solar eclipse near new moon', () { + // Known new moon: 2024-04-08 ~18:21 UTC + final knownNewMoon = DateTime.utc(2024, 4, 8, 18, 21); + final phase = getMoonPhase(knownNewMoon); + // Near new moon: illumination should be very low + expect(phase.illumination, lessThan(1.0)); + expect(phase.elongationDeg, lessThan(5)); + }); + }); + + // ── distanceKm range ────────────────────────────────────────────────── + + group('distanceKm', () { + test('distance stays within lunar orbit bounds over a month', () { + for (var day = 1; day <= 30; day++) { + final pos = getMoonPosition(DateTime.utc(2025, 3, day, 12), 0, 0); + expect(pos.distance, greaterThan(355000)); + expect(pos.distance, lessThan(410000)); + } + }); + }); +}