From 9167d86c7b4ae864d3a9bf0b3b114b26da0c8a83 Mon Sep 17 00:00:00 2001 From: Aric Camarata Date: Sun, 8 Mar 2026 12:48:22 -0400 Subject: [PATCH] Initial release: pray_calc_dart v1.0.0 --- .github/workflows/ci.yml | 44 ++ .gitignore | 20 + CHANGELOG.md | 19 + LICENSE | 16 + README.md | 108 +++ analysis_options.yaml | 4 + lib/pray_calc_dart.dart | 15 + lib/src/angles.dart | 117 +++ lib/src/asr.dart | 43 ++ lib/src/get_times.dart | 108 +++ lib/src/msc.dart | 156 ++++ lib/src/qiyam.dart | 23 + lib/src/solar_ephemeris.dart | 98 +++ lib/src/spa.dart | 1279 +++++++++++++++++++++++++++++++++ lib/src/types.dart | 142 ++++ pubspec.yaml | 21 + test/pray_calc_dart_test.dart | 193 +++++ 17 files changed, 2406 insertions(+) create mode 100644 .github/workflows/ci.yml create mode 100644 .gitignore create mode 100644 CHANGELOG.md create mode 100644 README.md create mode 100644 analysis_options.yaml create mode 100644 lib/pray_calc_dart.dart create mode 100644 lib/src/angles.dart create mode 100644 lib/src/asr.dart create mode 100644 lib/src/get_times.dart create mode 100644 lib/src/msc.dart create mode 100644 lib/src/qiyam.dart create mode 100644 lib/src/solar_ephemeris.dart create mode 100644 lib/src/spa.dart create mode 100644 lib/src/types.dart create mode 100644 pubspec.yaml create mode 100644 test/pray_calc_dart_test.dart diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..f2c434c --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,44 @@ +name: CI + +on: + push: + branches: [main] + pull_request: + branches: [main] + +jobs: + test: + name: Test (Dart ${{ matrix.dart-version }}) + runs-on: ubuntu-latest + strategy: + matrix: + dart-version: [stable, beta] + steps: + - uses: actions/checkout@v4 + - uses: dart-lang/setup-dart@v1 + with: + sdk: ${{ matrix.dart-version }} + - run: dart pub get + - run: dart analyze + - run: dart test + + format: + name: Format + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: dart-lang/setup-dart@v1 + with: + sdk: stable + - run: dart format --set-exit-if-changed . + + publish-check: + name: 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..9f8b6c9 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,19 @@ +# Changelog + +All notable changes to this project will be documented in this file. + +## [1.0.0] - 2026-03-08 + +### Added + +- `getTimes()` computes all prayer times using the PrayCalc Dynamic Method +- `getAngles()` returns adaptive Fajr/Isha twilight depression angles +- `getSpa()` full NREL Solar Position Algorithm (SPA) implementation +- `solarEphemeris()` low-precision Jean Meeus Ch. 25 ephemeris +- `getAsr()` computes Asr time for Shafi'i or Hanafi conventions +- `getQiyam()` computes start of the last third of the night +- `getMscFajr()` and `getMscIsha()` MCW seasonal offsets +- `formatTime()` converts fractional hours to HH:MM:SS strings +- `minutesToDepression()` converts time offsets to depression angles +- Full type definitions: PrayerTimes, FormattedPrayerTimes, TwilightAngles, SpaResult, SolarEphemeris +- Comprehensive test suite validated against pray-calc TypeScript v2.0.0 diff --git a/LICENSE b/LICENSE index f8c021c..5965e39 100644 --- a/LICENSE +++ b/LICENSE @@ -19,3 +19,19 @@ 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. + +Third-Party Notice +================== + +The NREL Solar Position Algorithm implementation is based on: + + Reda, I. and Andreas, A. (2004). Solar Position Algorithm for Solar + Radiation Applications. NREL/TP-560-34302. + DOI: 10.2172/15003974 + + National Renewable Energy Laboratory (NREL) + 1617 Cole Blvd, Golden, CO 80401 + +The original C implementation is Copyright (c) 2003-2012 NREL. +This Dart port is an independent implementation. NREL has not endorsed +or certified this library. diff --git a/README.md b/README.md new file mode 100644 index 0000000..b031178 --- /dev/null +++ b/README.md @@ -0,0 +1,108 @@ +# pray_calc_dart + +[![pub package](https://img.shields.io/pub/v/pray_calc_dart.svg)](https://pub.dev/packages/pray_calc_dart) +[![CI](https://github.com/acamarata/pray-calc-dart/actions/workflows/ci.yml/badge.svg)](https://github.com/acamarata/pray-calc-dart/actions/workflows/ci.yml) +[![License: MIT](https://img.shields.io/badge/license-MIT-blue.svg)](LICENSE) + +Islamic prayer times for Dart and Flutter. Pure Dart port of [pray-calc](https://github.com/acamarata/pray-calc), implementing the NREL Solar Position Algorithm, MCW seasonal model, and dynamic twilight angles. Zero dependencies. + +## Installation + +```yaml +dependencies: + pray_calc_dart: ^1.0.0 +``` + +## Quick Start + +```dart +import 'package:pray_calc_dart/pray_calc_dart.dart'; + +void main() { + final date = DateTime(2024, 3, 15); + final times = getTimes(date, 40.7128, -74.0060, -5.0); + + print('Fajr: ${formatTime(times.fajr)}'); + print('Sunrise: ${formatTime(times.sunrise)}'); + print('Dhuhr: ${formatTime(times.dhuhr)}'); + print('Asr: ${formatTime(times.asr)}'); + print('Maghrib: ${formatTime(times.maghrib)}'); + print('Isha: ${formatTime(times.isha)}'); + print('Qiyam: ${formatTime(times.qiyam)}'); +} +``` + +## API + +### `getTimes(date, lat, lng, tz, {elevation, temperature, pressure, hanafi})` + +Computes all prayer times for a given date and location. + +| Parameter | Type | Default | Description | +| --- | --- | --- | --- | +| `date` | `DateTime` | required | Local date (time-of-day ignored) | +| `lat` | `double` | required | Latitude (-90 to 90, south negative) | +| `lng` | `double` | required | Longitude (-180 to 180, west negative) | +| `tz` | `double` | required | UTC offset in hours (e.g., -5 for EST) | +| `elevation` | `double` | 0 | Observer elevation in meters | +| `temperature` | `double` | 15 | Ambient temperature in Celsius | +| `pressure` | `double` | 1013.25 | Atmospheric pressure in mbar | +| `hanafi` | `bool` | false | Hanafi Asr (2x shadow) vs Shafi'i (1x) | + +Returns a `PrayerTimes` object with fractional-hour values for: `qiyam`, `fajr`, `sunrise`, `noon`, `dhuhr`, `asr`, `maghrib`, `isha`, and the computed `angles`. + +### `getAngles(date, lat, lng, {elevation, temperature, pressure})` + +Computes dynamic twilight depression angles using the three-layer model: + +1. MCW seasonal base (piecewise-linear, latitude-dependent) +2. Ephemeris corrections (Earth-Sun distance, Fourier season smoothing) +3. Environmental corrections (elevation dip, atmospheric refraction) + +Returns `TwilightAngles` with `fajrAngle` and `ishaAngle` in degrees, clipped to [10, 22]. + +### `getSpa(date, lat, lng, tz, {...})` + +Full NREL Solar Position Algorithm. Accurate to +/-0.0003 degrees for zenith angle. Supports custom zenith angles for twilight calculations. + +### `formatTime(hours)` + +Converts fractional hours to an `HH:MM:SS` string. Returns `"N/A"` for non-finite values. + +### Additional Functions + +- `solarEphemeris(jd)` -- Jean Meeus Ch. 25 low-precision ephemeris +- `toJulianDate(date)` -- DateTime to Julian Date +- `getAsr(solarNoon, latitude, declination, {hanafi})` -- Asr computation +- `getQiyam(fajrTime, ishaTime)` -- Last third of the night +- `getMscFajr(date, latitude)` -- MCW Fajr offset in minutes +- `getMscIsha(date, latitude, [shafaq])` -- MCW Isha offset in minutes +- `minutesToDepression(minutes, latDeg, declDeg)` -- Time to angle conversion + +## Dynamic Angle Algorithm + +Fixed-angle methods (ISNA 15 degrees, MWL 18 degrees, etc.) produce inaccurate Fajr times at latitudes above 45 degrees N/S. The dynamic method adapts the depression angle based on season, latitude, Earth-Sun distance, and local atmospheric conditions. + +Result: approximately 18 degrees at the equator, approximately 12-14 degrees at 50-55 degrees N in summer. Matches observational data from the Moonsighting Committee Worldwide. + +## Compatibility + +Dart SDK 3.7.0+. Works in Flutter (iOS, Android, Web, Desktop), Dart CLI, and server-side Dart. Zero external dependencies. + +## Related + +- [pray-calc](https://github.com/acamarata/pray-calc) - TypeScript/JavaScript version (npm) +- [nrel-spa](https://github.com/acamarata/nrel-spa) - Standalone NREL SPA for JavaScript +- [qibla](https://github.com/acamarata/qibla) - Qibla direction calculator + +## Acknowledgments + +The Solar Position Algorithm is based on: + +> Reda, I. and Andreas, A. (2004). Solar Position Algorithm for Solar Radiation Applications. NREL/TP-560-34302. [DOI: 10.2172/15003974](https://doi.org/10.2172/15003974) + +The MCW seasonal model is based on the work of the [Moonsighting Committee Worldwide](http://moonsighting.com/isha_fajr.html) (Khalid Shaukat). + +## License + +[MIT](LICENSE). The NREL SPA implementation carries its own terms (see LICENSE for details). diff --git a/analysis_options.yaml b/analysis_options.yaml new file mode 100644 index 0000000..12e713a --- /dev/null +++ b/analysis_options.yaml @@ -0,0 +1,4 @@ +include: package:lints/recommended.yaml + +# Additional information about this file can be found at +# https://dart.dev/guides/language/analysis-options diff --git a/lib/pray_calc_dart.dart b/lib/pray_calc_dart.dart new file mode 100644 index 0000000..13286b6 --- /dev/null +++ b/lib/pray_calc_dart.dart @@ -0,0 +1,15 @@ +/// pray_calc_dart — Pure Dart Islamic prayer time calculation. +/// +/// Implements the PrayCalc Dynamic Method: NREL SPA algorithm + MSC seasonal +/// algorithm + dynamic twilight angles. Accurate to within 1 second of the +/// reference pray-calc TypeScript library. +library; + +export 'src/types.dart'; +export 'src/get_times.dart'; +export 'src/angles.dart'; +export 'src/solar_ephemeris.dart'; +export 'src/msc.dart'; +export 'src/asr.dart'; +export 'src/qiyam.dart'; +export 'src/spa.dart' show getSpa; diff --git a/lib/src/angles.dart b/lib/src/angles.dart new file mode 100644 index 0000000..063d1e8 --- /dev/null +++ b/lib/src/angles.dart @@ -0,0 +1,117 @@ +/// Dynamic twilight angle algorithm — PrayCalc Dynamic Method v2. +/// +/// Computes adaptive Fajr and Isha solar depression angles that accurately +/// track the observable phenomenon across all latitudes and seasons. +/// +/// Three-layer model: +/// 1. MSC base (MCW piecewise seasonal, converted to depression angle) +/// 2. Ephemeris corrections (Earth-Sun distance, Fourier season smoothing) +/// 3. Environmental corrections (elevation, atmospheric refraction) +library; + +import 'dart:math'; + +import 'types.dart'; +import 'solar_ephemeris.dart'; +import 'msc.dart'; + +const double _kFajrMin = 10; +const double _kFajrMax = 22; +const double _kIshaMin = 10; +const double _kIshaMax = 22; + +double _clip(double value, double min, double max) { + return value < min ? min : (value > max ? max : value); +} + +double _round3(double value) { + return (value * 1000).round() / 1000.0; +} + +/// Earth-Sun distance correction in degrees. +/// Effect magnitude: ≈ ±0.015°. +double _earthSunDistanceCorrection(double r) { + return -0.5 * log(r); +} + +/// Fourier smoothing correction (< 0.3° total) to remove MCW piecewise +/// artifacts and add hemisphere-symmetric season curvature. +double _fourierSmoothingCorrection(double eclLon, double latAbsDeg) { + const deg = pi / 180; + final theta = eclLon; // solar ecliptic longitude, radians [0, 2π) + final phi = latAbsDeg * deg; + + final a1 = 0.03 * sin(theta); + final b1 = -0.05 * cos(theta); + final a2 = 0.02 * sin(2 * theta); + final b2 = 0.02 * cos(2 * theta); + final c1 = -0.008 * phi * sin(theta); + final d1 = 0.004 * phi * cos(theta); + + return a1 + b1 + a2 + b2 + c1 + d1; +} + +/// Compute dynamic twilight depression angles for Fajr and Isha. +/// +/// [date] is the observer's local date (time-of-day is ignored). +/// [lat] is latitude in decimal degrees. +/// [lng] is longitude in decimal degrees (reserved, currently unused). +/// [elevation] is observer elevation in meters (default: 0). +/// [temperature] is ambient temperature in °C (default: 15). +/// [pressure] is atmospheric pressure in mbar (default: 1013.25). +TwilightAngles getAngles( + DateTime date, + double lat, + double lng, { + double elevation = 0, + double temperature = 15, + double pressure = 1013.25, +}) { + // 1. Solar ephemeris at UTC noon of the given date. + final noonDate = DateTime.utc(date.year, date.month, date.day, 12, 0, 0); + final jd = toJulianDate(noonDate); + final eph = solarEphemeris(jd); + + // 2. MCW reference times (minutes before/after sunrise/sunset). + final mscFajrMin = getMscFajr(date, lat); + final mscIshaMin = getMscIsha(date, lat); + + // 3. Convert MCW minutes to equivalent solar depression angles. + double fajrBase = minutesToDepression(mscFajrMin, lat, eph.decl); + double ishaBase = minutesToDepression(mscIshaMin, lat, eph.decl); + + // Handle polar or unreachable geometry. + if (!fajrBase.isFinite || fajrBase.isNaN) fajrBase = 18.0; + if (!ishaBase.isFinite || ishaBase.isNaN) ishaBase = 18.0; + + // 4. Earth-Sun distance correction. + final rCorr = _earthSunDistanceCorrection(eph.r); + + // 5. Fourier smoothing correction. + final fourierCorr = _fourierSmoothingCorrection(eph.eclLon, lat.abs()); + + // 6. Atmospheric refraction at expected twilight depression. + final refrFajr = atmosphericRefraction( + -(fajrBase + 0.5), + pressureMbar: pressure, + temperatureC: temperature, + ); + final refrIsha = atmosphericRefraction( + -(ishaBase + 0.5), + pressureMbar: pressure, + temperatureC: temperature, + ); + + // 7. Elevation correction (horizon dip). + final horizonDipDeg = 1.06 * sqrt(elevation / 1000); + final elevCorr = horizonDipDeg * 0.3; + + // 8. Assemble final angles. + final rawFajr = fajrBase + rCorr + fourierCorr + refrFajr + elevCorr; + final rawIsha = ishaBase + rCorr + fourierCorr + refrIsha + elevCorr; + + final fajrAngle = _round3(_clip(rawFajr, _kFajrMin, _kFajrMax)); + final ishaAngle = _round3(_clip(rawIsha, _kIshaMin, _kIshaMax)); + + return TwilightAngles(fajrAngle: fajrAngle, ishaAngle: ishaAngle); +} diff --git a/lib/src/asr.dart b/lib/src/asr.dart new file mode 100644 index 0000000..d13f6cb --- /dev/null +++ b/lib/src/asr.dart @@ -0,0 +1,43 @@ +/// Asr prayer time calculation. +/// +/// Asr begins when the shadow of an object equals (Shafi'i/Maliki/Hanbali) +/// or twice (Hanafi) the object's length plus its shadow at solar noon. +library; + +import 'dart:math'; + +/// Compute Asr time as fractional hours. +/// +/// [solarNoon] is solar noon in fractional hours (from getSpa). +/// [latitude] is observer latitude in degrees. +/// [declination] is solar declination in degrees (from solarEphemeris). +/// [hanafi] is true for Hanafi (shadow factor 2), false for Shafi'i (factor 1). +/// +/// Returns fractional hours, or [double.nan] if the sun never reaches the +/// required altitude. +double getAsr( + double solarNoon, + double latitude, + double declination, { + bool hanafi = false, +}) { + const deg = pi / 180; + final phi = latitude * deg; + final delta = declination * deg; + final shadowFactor = hanafi ? 2.0 : 1.0; + + // Required solar altitude: tan(A) = 1 / (shadowFactor + tan(|φ − δ|)) + final x = (phi - delta).abs(); + final tanA = 1.0 / (shadowFactor + tan(x)); + final sinA = tanA / sqrt(1 + tanA * tanA); // sin(atan(tanA)) + + // cos(H0) = (sin(A) − sin(φ)sin(δ)) / (cos(φ)cos(δ)) + final cosH0 = (sinA - sin(phi) * sin(delta)) / (cos(phi) * cos(delta)); + + if (cosH0 < -1 || cosH0 > 1) return double.nan; + + // H0 in hours (15°/hr) + final h0h = acos(cosH0) / deg / 15; + + return solarNoon + h0h; +} diff --git a/lib/src/get_times.dart b/lib/src/get_times.dart new file mode 100644 index 0000000..1cd74ee --- /dev/null +++ b/lib/src/get_times.dart @@ -0,0 +1,108 @@ +/// Core prayer times computation — PrayCalc Dynamic Method. +/// +/// Returns all prayer times as fractional hours using the dynamic twilight +/// angle algorithm. Times are in local time as determined by the UTC offset. +library; + +import 'types.dart'; +import 'spa.dart'; +import 'solar_ephemeris.dart'; +import 'angles.dart'; +import 'asr.dart'; +import 'qiyam.dart'; + +/// Compute prayer times for a given date and location. +/// +/// [date] is the observer's local date (time-of-day is ignored). +/// [lat] is latitude in decimal degrees (−90 to 90, south = negative). +/// [lng] is longitude in decimal degrees (−180 to 180, west = negative). +/// [tz] is UTC offset in hours (e.g., −5 for EST). +/// [elevation] is observer elevation in meters (default: 0). +/// [temperature] is ambient temperature in °C (default: 15). +/// [pressure] is atmospheric pressure in mbar/hPa (default: 1013.25). +/// [hanafi] selects Asr convention: false = Shafi'i/Maliki/Hanbali (default), +/// true = Hanafi. +PrayerTimes getTimes( + DateTime date, + double lat, + double lng, + double tz, { + double elevation = 0, + double temperature = 15, + double pressure = 1013.25, + bool hanafi = false, +}) { + // 1. Compute dynamic twilight angles. + final tw = getAngles( + date, + lat, + lng, + elevation: elevation, + temperature: temperature, + pressure: pressure, + ); + + // 2. Convert depression angles to SPA zenith angles. + // SPA uses zenith (90° + depression) for custom altitude events. + final fajrZenith = 90 + tw.fajrAngle; + final ishaZenith = 90 + tw.ishaAngle; + + // 3. Run SPA for solar position + custom twilight times. + final spaData = getSpa( + date, + lat, + lng, + tz, + elevation: elevation, + temperature: temperature, + pressure: pressure, + customAngles: [fajrZenith, ishaZenith], + ); + + final fajrTime = spaData.angles[0].sunrise; + final sunriseTime = spaData.sunrise; + final noonTime = spaData.solarNoon; + final maghribTime = spaData.sunset; + final ishaTime = spaData.angles[1].sunset; + + // Dhuhr: 2.5 minutes after solar noon. + final dhuhrTime = noonTime + 2.5 / 60; + + // 4. Solar declination for Asr (Meeus formula, accurate to ~0.01°). + final jd = toJulianDate( + DateTime.utc(date.year, date.month, date.day, 12, 0, 0), + ); + final eph = solarEphemeris(jd); + + // 5. Asr time. + final asrTime = getAsr(noonTime, lat, eph.decl, hanafi: hanafi); + + // 6. Qiyam al-Layl (last third of the night). + final qiyamTime = getQiyam(fajrTime, ishaTime); + + return PrayerTimes( + qiyam: qiyamTime.isFinite ? qiyamTime : double.nan, + fajr: fajrTime.isFinite ? fajrTime : double.nan, + sunrise: sunriseTime.isFinite ? sunriseTime : double.nan, + noon: noonTime.isFinite ? noonTime : double.nan, + dhuhr: dhuhrTime.isFinite ? dhuhrTime : double.nan, + asr: asrTime.isFinite ? asrTime : double.nan, + maghrib: maghribTime.isFinite ? maghribTime : double.nan, + isha: ishaTime.isFinite ? ishaTime : double.nan, + angles: tw, + ); +} + +/// Format fractional hours as HH:MM:SS string. +/// Returns "N/A" if the value is non-finite or negative. +String formatTime(double hours) { + if (!hours.isFinite || hours < 0) return 'N/A'; + final totalSec = (hours * 3600).round(); + final h = (totalSec ~/ 3600) % 24; + final rem = totalSec - (totalSec ~/ 3600) * 3600; + final m = rem ~/ 60; + final s = rem - m * 60; + return '${h.toString().padLeft(2, '0')}:' + '${m.toString().padLeft(2, '0')}:' + '${s.toString().padLeft(2, '0')}'; +} diff --git a/lib/src/msc.dart b/lib/src/msc.dart new file mode 100644 index 0000000..f669fb1 --- /dev/null +++ b/lib/src/msc.dart @@ -0,0 +1,156 @@ +/// Moonsighting Committee Worldwide (MCW) seasonal algorithm. +/// +/// Computes Fajr and Isha as time offsets from sunrise/sunset using the +/// empirical piecewise-linear seasonal functions developed by the Moonsighting +/// Committee Worldwide (Khalid Shaukat). +/// +/// Reference: moonsighting.com/isha_fajr.html +library; + +import 'dart:math'; + +import 'types.dart'; + +bool _isLeapYear(int year) { + return (year % 4 == 0 && year % 100 != 0) || year % 400 == 0; +} + +/// Compute the MCW seasonal index (dyy) and days in year. +({int dyy, int daysInYear}) _computeDyy(DateTime date, double latitude) { + final year = date.year; + final daysInYear = _isLeapYear(year) ? 366 : 365; + + // Reference solstice: Dec 21 for Northern, Jun 21 for Southern + final refMonth = latitude >= 0 ? 12 : 6; // Dec = 12, Jun = 6 (1-based) + const refDay = 21; + + final zeroDate = DateTime.utc(year, refMonth, refDay); + final inputUtc = DateTime.utc(date.year, date.month, date.day); + + int diffDays = inputUtc.difference(zeroDate).inDays; + if (diffDays < 0) diffDays += daysInYear; + + return (dyy: diffDays, daysInYear: daysInYear); +} + +/// Piecewise-linear seasonal interpolation over 6 segments. +double _interpolateSegment( + int dyy, + int daysInYear, + double a, + double b, + double c, + double d, +) { + if (dyy < 91) { + return a + ((b - a) / 91) * dyy; + } else if (dyy < 137) { + return b + ((c - b) / 46) * (dyy - 91); + } else if (dyy < 183) { + return c + ((d - c) / 46) * (dyy - 137); + } else if (dyy < 229) { + return d + ((c - d) / 46) * (dyy - 183); + } else if (dyy < 275) { + return c + ((b - c) / 46) * (dyy - 229); + } else { + final len = daysInYear - 275; + return b + ((a - b) / len) * (dyy - 275); + } +} + +/// Compute Fajr offset in minutes before sunrise using the MCW algorithm. +/// +/// Returns minutes before sunrise (rounded to nearest minute). +double getMscFajr(DateTime date, double latitude) { + final latAbs = latitude.abs(); + final (:dyy, :daysInYear) = _computeDyy(date, latitude); + + final a = 75 + (28.65 / 55) * latAbs; + final b = 75 + (19.44 / 55) * latAbs; + final c = 75 + (32.74 / 55) * latAbs; + final d = 75 + (48.1 / 55) * latAbs; + + return _interpolateSegment(dyy, daysInYear, a, b, c, d).roundToDouble(); +} + +/// Compute Isha offset in minutes after sunset using the MCW algorithm. +/// +/// [shafaq] selects the twilight mode: general (default), ahmer, or abyad. +/// Returns minutes after sunset (rounded to nearest minute). +double getMscIsha( + DateTime date, + double latitude, [ + ShafaqMode shafaq = ShafaqMode.general, +]) { + final latAbs = latitude.abs(); + final (:dyy, :daysInYear) = _computeDyy(date, latitude); + + double a, b, c, d; + + switch (shafaq) { + case ShafaqMode.ahmer: + a = 62 + (17.4 / 55) * latAbs; + b = 62 - (7.16 / 55) * latAbs; + c = 62 + (5.12 / 55) * latAbs; + d = 62 + (19.44 / 55) * latAbs; + case ShafaqMode.abyad: + a = 75 + (25.6 / 55) * latAbs; + b = 75 + (7.16 / 55) * latAbs; + c = 75 + (36.84 / 55) * latAbs; + d = 75 + (81.84 / 55) * latAbs; + case ShafaqMode.general: + a = 75 + (25.6 / 55) * latAbs; + b = 75 + (2.05 / 55) * latAbs; + c = 75 - (9.21 / 55) * latAbs; + d = 75 + (6.14 / 55) * latAbs; + } + + return _interpolateSegment(dyy, daysInYear, a, b, c, d).roundToDouble(); +} + +/// Convert MCW minutes-before-sunrise to an equivalent solar depression angle +/// in degrees, using exact spherical trigonometry. +/// +/// Returns [double.nan] if the geometry is unreachable (polar day/night). +double minutesToDepression(double minutes, double latDeg, double declDeg) { + final phi = latDeg * (pi / 180); + final delta = declDeg * (pi / 180); + + final cosPhi = cos(phi); + final sinPhi = sin(phi); + final cosDelta = cos(delta); + final sinDelta = sin(delta); + + // Standard sunrise/sunset: h = -0.833° (includes refraction + semi-diameter) + final h0 = -0.833 * (pi / 180); + final sinH0 = sin(h0); + + final denominator = cosPhi * cosDelta; + if (denominator.abs() < 1e-10) return double.nan; + + // Hour angle at standard sunrise + final cosHRise = (sinH0 - sinPhi * sinDelta) / denominator; + + if (cosHRise < -1) return double.nan; // polar night + if (cosHRise > 1) return double.nan; // polar day + + final hRise = acos(cosHRise); // radians + + // Hour angle at the prayer time (further from solar noon) + final deltaH = (minutes / 60) * 15 * (pi / 180); + final hPrayer = hRise + deltaH; + + // Cap at π (midnight) + if (hPrayer > pi) { + final sinHMidnight = sinPhi * sinDelta + cosPhi * cosDelta * cos(pi); + final hMidnight = asin(sinHMidnight.clamp(-1.0, 1.0)); + return -hMidnight / (pi / 180); + } + + // Solar altitude at hPrayer + final sinHPrayer = sinPhi * sinDelta + cosPhi * cosDelta * cos(hPrayer); + final hPrayerAlt = asin(sinHPrayer.clamp(-1.0, 1.0)); + + // Depression angle: positive when sun is below horizon + return -hPrayerAlt / (pi / 180); +} diff --git a/lib/src/qiyam.dart b/lib/src/qiyam.dart new file mode 100644 index 0000000..cdb24d5 --- /dev/null +++ b/lib/src/qiyam.dart @@ -0,0 +1,23 @@ +/// Qiyam al-Layl (night prayer) time calculation. +/// +/// Returns the start of the last third of the night, which is the recommended +/// time for Tahajjud / Qiyam al-Layl. The night is defined as the period +/// from Isha to Fajr. +library; + +/// Compute the start of the last third of the night. +/// +/// [fajrTime] is Fajr time in fractional hours. +/// [ishaTime] is Isha time in fractional hours. +/// +/// Returns start of the last third of the night (fractional hours). +double getQiyam(double fajrTime, double ishaTime) { + // If Fajr is numerically earlier (e.g. 5.5) than Isha (e.g. 21.5), Fajr + // is actually the NEXT day — add 24 to get the correct night length. + final adjustedFajr = fajrTime < ishaTime ? fajrTime + 24 : fajrTime; + + final nightLength = adjustedFajr - ishaTime; + final lastThirdStart = ishaTime + (2 * nightLength) / 3; + + return lastThirdStart >= 24 ? lastThirdStart - 24 : lastThirdStart; +} diff --git a/lib/src/solar_ephemeris.dart b/lib/src/solar_ephemeris.dart new file mode 100644 index 0000000..f025ab9 --- /dev/null +++ b/lib/src/solar_ephemeris.dart @@ -0,0 +1,98 @@ +/// High-accuracy solar ephemeris using Jean Meeus "Astronomical Algorithms" +/// (2nd ed., Ch. 25) low-precision formulas. +/// +/// Accuracy: ~0.01° for solar declination, ~0.0001 AU for Earth-Sun distance +/// over years 1950–2050. +library; + +import 'dart:math'; + +import 'types.dart'; + +/// Julian Date from a Dart [DateTime] (UTC). +double toJulianDate(DateTime date) { + return date.millisecondsSinceEpoch / 86400000.0 + 2440587.5; +} + +/// Compute solar declination, Earth-Sun distance, and ecliptic longitude +/// from a Julian Date. +SolarEphemeris solarEphemeris(double jd) { + const deg = pi / 180; + final t = (jd - 2451545.0) / 36525.0; + + // Geometric mean longitude L0 (degrees) + final l0 = + ((280.46646 + 36000.76983 * t + 0.0003032 * t * t) % 360 + 360) % 360; + + // Mean anomaly M (degrees) + final m = + ((357.52911 + 35999.05029 * t - 0.0001537 * t * t) % 360 + 360) % 360; + final mRad = m * deg; + + // Orbital eccentricity + final e = 0.016708634 - 0.000042037 * t - 0.0000001267 * t * t; + + // Equation of center C (degrees) + final c = + (1.914602 - 0.004817 * t - 0.000014 * t * t) * sin(mRad) + + (0.019993 - 0.000101 * t) * sin(2 * mRad) + + 0.000289 * sin(3 * mRad); + + // Sun's true longitude (degrees) + final sunLon = l0 + c; + + // Sun's true anomaly (degrees) + final nu = m + c; + final nuRad = nu * deg; + + // Earth-Sun distance in AU + final r = (1.000001018 * (1 - e * e)) / (1 + e * cos(nuRad)); + + // Longitude of ascending node of Moon's orbit (for nutation) + final omega = ((125.04 - 1934.136 * t) % 360 + 360) % 360; + final omegaRad = omega * deg; + + // Apparent solar longitude corrected for nutation and aberration + final lambda = sunLon - 0.00569 - 0.00478 * sin(omegaRad); + final lambdaRad = lambda * deg; + + // Mean obliquity of the ecliptic (degrees) + final epsilon0 = + 23.439291 - 0.013004 * t - 1.638e-7 * t * t + 5.036e-7 * t * t * t; + + // True obliquity with nutation correction + final epsilon = (epsilon0 + 0.00256 * cos(omegaRad)) * deg; + + // Solar declination + final sinDecl = sin(epsilon) * sin(lambdaRad); + final decl = asin(sinDecl.clamp(-1.0, 1.0)) / deg; + + // Ecliptic longitude as season phase θ ∈ [0, 2π) + final eclLon = ((lambdaRad % (2 * pi)) + 2 * pi) % (2 * pi); + + return SolarEphemeris(decl: decl, r: r, eclLon: eclLon); +} + +/// Solar vertical angular speed near a given hour angle [hAngleRad] (radians), +/// in degrees per hour. +double solarVerticalSpeed(double latRad, double declRad, double hAngleRad) { + return 15 * (cos(latRad) * cos(declRad) * sin(hAngleRad)).abs(); +} + +/// Compute the atmospheric refraction correction (degrees) for a given +/// apparent solar altitude using the Bennett/Saemundsson formula. +/// +/// Returns a positive correction. For altitudes below -1°, returns 0. +double atmosphericRefraction( + double altitudeDeg, { + double pressureMbar = 1013.25, + double temperatureC = 15, +}) { + if (altitudeDeg < -1) return 0; + const deg = pi / 180; + // Bennett's formula in arcminutes + final r0 = 1.02 / tan((altitudeDeg + 10.3 / (altitudeDeg + 5.11)) * deg); + // Scale for pressure and temperature + final r = r0 * (pressureMbar / 1010) * (283 / (273 + temperatureC)); + return r < 0 ? 0.0 : r / 60; // convert arcminutes to degrees +} diff --git a/lib/src/spa.dart b/lib/src/spa.dart new file mode 100644 index 0000000..6429d42 --- /dev/null +++ b/lib/src/spa.dart @@ -0,0 +1,1279 @@ +/// NREL Solar Position Algorithm (SPA) — Dart port. +/// +/// Direct port of the nrel-spa JavaScript library (spa.js v2.0.1). +/// Accurate to ±0.0003° for solar zenith angle. +/// +/// Reference: Reda, I. and Andreas, A. (2004). Solar Position Algorithm for +/// Solar Radiation Applications. NREL/TP-560-34302. +library; + +import 'dart:math'; + +import 'types.dart'; + +// ─── Constants ────────────────────────────────────────────────────────────── + +const int _sPaZaRts = 2; +const int _sPaAll = 3; +const double _sunRadius = 0.26667; + +// ─── Earth Periodic Term Tables ───────────────────────────────────────────── + +// Each row: [A, B, C] where A is amplitude, B is phase, C is frequency +const List>> _lTerms = [ + // L0 — 64 terms + [ + [175347046.0, 0, 0], + [3341656.0, 4.6692568, 6283.07585], + [34894.0, 4.6261, 12566.1517], + [3497.0, 2.7441, 5753.3849], + [3418.0, 2.8289, 3.5231], + [3136.0, 3.6277, 77713.7715], + [2676.0, 4.4181, 7860.4194], + [2343.0, 6.1352, 3930.2097], + [1324.0, 0.7425, 11506.7698], + [1273.0, 2.0371, 529.691], + [1199.0, 1.1096, 1577.3435], + [990, 5.233, 5884.927], + [902, 2.045, 26.298], + [857, 3.508, 398.149], + [780, 1.179, 5223.694], + [753, 2.533, 5507.553], + [505, 4.583, 18849.228], + [492, 4.205, 775.523], + [357, 2.92, 0.067], + [317, 5.849, 11790.629], + [284, 1.899, 796.298], + [271, 0.315, 10977.079], + [243, 0.345, 5486.778], + [206, 4.806, 2544.314], + [205, 1.869, 5573.143], + [202, 2.458, 6069.777], + [156, 0.833, 213.299], + [132, 3.411, 2942.463], + [126, 1.083, 20.775], + [115, 0.645, 0.98], + [103, 0.636, 4694.003], + [102, 0.976, 15720.839], + [102, 4.267, 7.114], + [99, 6.21, 2146.17], + [98, 0.68, 155.42], + [86, 5.98, 161000.69], + [85, 1.3, 6275.96], + [85, 3.67, 71430.7], + [80, 1.81, 17260.15], + [79, 3.04, 12036.46], + [75, 1.76, 5088.63], + [74, 3.5, 3154.69], + [74, 4.68, 801.82], + [70, 0.83, 9437.76], + [62, 3.98, 8827.39], + [61, 1.82, 7084.9], + [57, 2.78, 6286.6], + [56, 4.39, 14143.5], + [56, 3.47, 6279.55], + [52, 0.19, 12139.55], + [52, 1.33, 1748.02], + [51, 0.28, 5856.48], + [49, 0.49, 1194.45], + [41, 5.37, 8429.24], + [41, 2.4, 19651.05], + [39, 6.17, 10447.39], + [37, 6.04, 10213.29], + [37, 2.57, 1059.38], + [36, 1.71, 2352.87], + [36, 1.78, 6812.77], + [33, 0.59, 17789.85], + [30, 0.44, 83996.85], + [30, 2.74, 1349.87], + [25, 3.16, 4690.48], + ], + // L1 — 34 terms + [ + [628331966747.0, 0, 0], + [206059.0, 2.678235, 6283.07585], + [4303.0, 2.6351, 12566.1517], + [425.0, 1.59, 3.523], + [119.0, 5.796, 26.298], + [109.0, 2.966, 1577.344], + [93, 2.59, 18849.23], + [72, 1.14, 529.69], + [68, 1.87, 398.15], + [67, 4.41, 5507.55], + [59, 2.89, 5223.69], + [56, 2.17, 155.42], + [45, 0.4, 796.3], + [36, 0.47, 775.52], + [29, 2.65, 7.11], + [21, 5.34, 0.98], + [19, 1.85, 5486.78], + [19, 4.97, 213.3], + [17, 2.99, 6275.96], + [16, 0.03, 2544.31], + [16, 1.43, 2146.17], + [15, 1.21, 10977.08], + [12, 2.83, 1748.02], + [12, 3.26, 5088.63], + [12, 5.27, 1194.45], + [12, 2.08, 4694], + [11, 0.77, 553.57], + [10, 1.3, 6286.6], + [10, 4.24, 1349.87], + [9, 2.7, 242.73], + [9, 5.64, 951.72], + [8, 5.3, 2352.87], + [6, 2.65, 9437.76], + [6, 4.67, 4690.48], + ], + // L2 — 20 terms + [ + [52919.0, 0, 0], + [8720.0, 1.0721, 6283.0758], + [309.0, 0.867, 12566.152], + [27, 0.05, 3.52], + [16, 5.19, 26.3], + [16, 3.68, 155.42], + [10, 0.76, 18849.23], + [9, 2.06, 77713.77], + [7, 0.83, 775.52], + [5, 4.66, 1577.34], + [4, 1.03, 7.11], + [4, 3.44, 5573.14], + [3, 5.14, 796.3], + [3, 6.05, 5507.55], + [3, 1.19, 242.73], + [3, 6.12, 529.69], + [3, 0.31, 398.15], + [3, 2.28, 553.57], + [2, 4.38, 5223.69], + [2, 3.75, 0.98], + ], + // L3 — 7 terms + [ + [289.0, 5.844, 6283.076], + [35, 0, 0], + [17, 5.49, 12566.15], + [3, 5.2, 155.42], + [1, 4.72, 3.52], + [1, 5.3, 18849.23], + [1, 5.97, 242.73], + ], + // L4 — 3 terms + [ + [114.0, 3.142, 0], + [8, 4.13, 6283.08], + [1, 3.84, 12566.15], + ], + // L5 — 1 term + [ + [1, 3.14, 0], + ], +]; + +const List>> _bTerms = [ + // B0 — 5 terms + [ + [280.0, 3.199, 84334.662], + [102.0, 5.422, 5507.553], + [80, 3.88, 5223.69], + [44, 3.7, 2352.87], + [32, 4, 1577.34], + ], + // B1 — 2 terms + [ + [9, 3.9, 5507.55], + [6, 1.73, 5223.69], + ], +]; + +const List>> _rTerms = [ + // R0 — 40 terms + [ + [100013989.0, 0, 0], + [1670700.0, 3.0984635, 6283.07585], + [13956.0, 3.05525, 12566.1517], + [3084.0, 5.1985, 77713.7715], + [1628.0, 1.1739, 5753.3849], + [1576.0, 2.8469, 7860.4194], + [925.0, 5.453, 11506.77], + [542.0, 4.564, 3930.21], + [472.0, 3.661, 5884.927], + [346.0, 0.964, 5507.553], + [329.0, 5.9, 5223.694], + [307.0, 0.299, 5573.143], + [243.0, 4.273, 11790.629], + [212.0, 5.847, 1577.344], + [186.0, 5.022, 10977.079], + [175.0, 3.012, 18849.228], + [110.0, 5.055, 5486.778], + [98, 0.89, 6069.78], + [86, 5.69, 15720.84], + [86, 1.27, 161000.69], + [65, 0.27, 17260.15], + [63, 0.92, 529.69], + [57, 2.01, 83996.85], + [56, 5.24, 71430.7], + [49, 3.25, 2544.31], + [47, 2.58, 775.52], + [45, 5.54, 9437.76], + [43, 6.01, 6275.96], + [39, 5.36, 4694], + [38, 2.39, 8827.39], + [37, 0.83, 19651.05], + [37, 4.9, 12139.55], + [36, 1.67, 12036.46], + [35, 1.84, 2942.46], + [33, 0.24, 7084.9], + [32, 0.18, 5088.63], + [32, 1.78, 398.15], + [28, 1.21, 6286.6], + [28, 1.9, 6279.55], + [26, 4.59, 10447.39], + ], + // R1 — 10 terms + [ + [103019.0, 1.10749, 6283.07585], + [1721.0, 1.0644, 12566.1517], + [702.0, 3.142, 0], + [32, 1.02, 18849.23], + [31, 2.84, 5507.55], + [25, 1.32, 5223.69], + [18, 1.42, 1577.34], + [10, 5.91, 10977.08], + [9, 1.42, 6275.96], + [9, 0.27, 5486.78], + ], + // R2 — 6 terms + [ + [4359.0, 5.7846, 6283.0758], + [124.0, 5.579, 12566.152], + [12, 3.14, 0], + [9, 3.63, 77713.77], + [6, 1.87, 5573.14], + [3, 5.47, 18849.23], + ], + // R3 — 2 terms + [ + [145.0, 4.273, 6283.076], + [7, 3.92, 12566.15], + ], + // R4 — 1 term + [ + [4, 2.56, 6283.08], + ], +]; + +// Periodic terms for nutation in longitude and obliquity +// Each row: [x0, x1, x2, x3, x4] +const List> _yTerms = [ + [0, 0, 0, 0, 1], + [-2, 0, 0, 2, 2], + [0, 0, 0, 2, 2], + [0, 0, 0, 0, 2], + [0, 1, 0, 0, 0], + [0, 0, 1, 0, 0], + [-2, 1, 0, 2, 2], + [0, 0, 0, 2, 1], + [0, 0, 1, 2, 2], + [-2, -1, 0, 2, 2], + [-2, 0, 1, 0, 0], + [-2, 0, 0, 2, 1], + [0, 0, -1, 2, 2], + [2, 0, 0, 0, 0], + [0, 0, 1, 0, 1], + [2, 0, -1, 2, 2], + [0, 0, -1, 0, 1], + [0, 0, 1, 2, 1], + [-2, 0, 2, 0, 0], + [0, 0, -2, 2, 1], + [2, 0, 0, 2, 2], + [0, 0, 2, 2, 2], + [0, 0, 2, 0, 0], + [-2, 0, 1, 2, 2], + [0, 0, 0, 2, 0], + [-2, 0, 0, 2, 0], + [0, 0, -1, 2, 1], + [0, 2, 0, 0, 0], + [2, 0, -1, 0, 1], + [-2, 2, 0, 2, 2], + [0, 1, 0, 0, 1], + [-2, 0, 1, 0, 1], + [0, -1, 0, 0, 1], + [0, 0, 2, -2, 0], + [2, 0, -1, 2, 1], + [2, 0, 1, 2, 2], + [0, 1, 0, 2, 2], + [-2, 1, 1, 0, 0], + [0, -1, 0, 2, 2], + [2, 0, 0, 2, 1], + [2, 0, 1, 0, 0], + [-2, 0, 2, 2, 2], + [-2, 0, 1, 2, 1], + [2, 0, -2, 0, 1], + [2, 0, 0, 0, 1], + [0, -1, 1, 0, 0], + [-2, -1, 0, 2, 1], + [-2, 0, 0, 0, 1], + [0, 0, 2, 2, 1], + [-2, 0, 2, 0, 1], + [-2, 1, 0, 2, 1], + [0, 0, 1, -2, 0], + [-1, 0, 1, 0, 0], + [-2, 1, 0, 0, 0], + [1, 0, 0, 0, 0], + [0, 0, 1, 2, 0], + [0, 0, -2, 2, 2], + [-1, -1, 1, 0, 0], + [0, 1, 1, 0, 0], + [0, -1, 1, 2, 2], + [2, -1, -1, 2, 2], + [0, 0, 3, 2, 2], + [2, -1, 0, 2, 2], +]; + +// Nutation longitude/obliquity coefficients [psiA, psiB, epsC, epsD] +const List> _peTerms = [ + [-171996, -174.2, 92025, 8.9], + [-13187, -1.6, 5736, -3.1], + [-2274, -0.2, 977, -0.5], + [2062, 0.2, -895, 0.5], + [1426, -3.4, 54, -0.1], + [712, 0.1, -7, 0], + [-517, 1.2, 224, -0.6], + [-386, -0.4, 200, 0], + [-301, 0, 129, -0.1], + [217, -0.5, -95, 0.3], + [-158, 0, 0, 0], + [129, 0.1, -70, 0], + [123, 0, -53, 0], + [63, 0, 0, 0], + [63, 0.1, -33, 0], + [-59, 0, 26, 0], + [-58, -0.1, 32, 0], + [-51, 0, 27, 0], + [48, 0, 0, 0], + [46, 0, -24, 0], + [-38, 0, 16, 0], + [-31, 0, 13, 0], + [29, 0, 0, 0], + [29, 0, -12, 0], + [26, 0, 0, 0], + [-22, 0, 0, 0], + [21, 0, -10, 0], + [17, -0.1, 0, 0], + [16, 0, -8, 0], + [-16, 0.1, 7, 0], + [-15, 0, 9, 0], + [-13, 0, 7, 0], + [-12, 0, 6, 0], + [11, 0, 0, 0], + [-10, 0, 5, 0], + [-8, 0, 3, 0], + [7, 0, -3, 0], + [-7, 0, 0, 0], + [-7, 0, 3, 0], + [-7, 0, 3, 0], + [6, 0, 0, 0], + [6, 0, -3, 0], + [6, 0, -3, 0], + [-6, 0, 3, 0], + [-6, 0, 3, 0], + [5, 0, 0, 0], + [-5, 0, 3, 0], + [-5, 0, 3, 0], + [-5, 0, 3, 0], + [4, 0, 0, 0], + [4, 0, 0, 0], + [4, 0, 0, 0], + [-4, 0, 0, 0], + [-4, 0, 0, 0], + [-4, 0, 0, 0], + [3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], +]; + +const List _lSubcount = [64, 34, 20, 7, 3, 1]; +const List _bSubcount = [5, 2]; +const List _rSubcount = [40, 10, 6, 2, 1]; + +// ─── Mutable internal state ───────────────────────────────────────────────── + +class _Spa { + // Inputs + int year = 0; + int month = 0; + int day = 0; + int hour = 0; + int minute = 0; + double second = 0.0; + double deltaUt1 = 0.0; + double deltaT = 0.0; + double timezone = 0.0; + double longitude = 0.0; + double latitude = 0.0; + double elevation = 0.0; + double pressure = 0.0; + double temperature = 0.0; + double slope = 0.0; + double azmRotation = 0.0; + double atmosRefract = 0.0; + int function = 0; + // Intermediate + double jd = 0.0; + double jc = 0.0; + double jde = 0.0; + double jce = 0.0; + double jme = 0.0; + double l = 0.0; + double b = 0.0; + double r = 0.0; + double theta = 0.0; + double beta = 0.0; + double x0 = 0.0; + double x1 = 0.0; + double x2 = 0.0; + double x3 = 0.0; + double x4 = 0.0; + double delPsi = 0.0; + double delEpsilon = 0.0; + double epsilon0 = 0.0; + double epsilon = 0.0; + double delTau = 0.0; + double lamda = 0.0; + double nu0 = 0.0; + double nu = 0.0; + double alpha = 0.0; + double delta = 0.0; + double h = 0.0; + double xi = 0.0; + double delAlpha = 0.0; + double deltaPrime = 0.0; + double alphaPrime = 0.0; + double hPrime = 0.0; + double e0 = 0.0; + double delE = 0.0; + double e = 0.0; + double eot = 0.0; + double srha = 0.0; + double ssha = 0.0; + double sta = 0.0; + // Outputs + double zenith = 0.0; + double azimuthAstro = 0.0; + double azimuth = 0.0; + double incidence = 0.0; + double suntransit = 0.0; + double sunrise = 0.0; + double sunset = 0.0; + + _Spa clone() { + final c = _Spa(); + c.year = year; + c.month = month; + c.day = day; + c.hour = hour; + c.minute = minute; + c.second = second; + c.deltaUt1 = deltaUt1; + c.deltaT = deltaT; + c.timezone = timezone; + c.longitude = longitude; + c.latitude = latitude; + c.elevation = elevation; + c.pressure = pressure; + c.temperature = temperature; + c.slope = slope; + c.azmRotation = azmRotation; + c.atmosRefract = atmosRefract; + c.function = function; + c.jd = jd; + c.jc = jc; + c.jde = jde; + c.jce = jce; + c.jme = jme; + c.l = l; + c.b = b; + c.r = r; + c.theta = theta; + c.beta = beta; + c.x0 = x0; + c.x1 = x1; + c.x2 = x2; + c.x3 = x3; + c.x4 = x4; + c.delPsi = delPsi; + c.delEpsilon = delEpsilon; + c.epsilon0 = epsilon0; + c.epsilon = epsilon; + c.delTau = delTau; + c.lamda = lamda; + c.nu0 = nu0; + c.nu = nu; + c.alpha = alpha; + c.delta = delta; + c.h = h; + c.xi = xi; + c.delAlpha = delAlpha; + c.deltaPrime = deltaPrime; + c.alphaPrime = alphaPrime; + c.hPrime = hPrime; + c.e0 = e0; + c.delE = delE; + c.e = e; + c.eot = eot; + c.srha = srha; + c.ssha = ssha; + c.sta = sta; + c.zenith = zenith; + c.azimuthAstro = azimuthAstro; + c.azimuth = azimuth; + c.incidence = incidence; + c.suntransit = suntransit; + c.sunrise = sunrise; + c.sunset = sunset; + return c; + } +} + +// ─── Math utility functions ────────────────────────────────────────────────── + +double _deg2rad(double degrees) => (pi / 180.0) * degrees; +double _rad2deg(double radians) => (180.0 / pi) * radians; + +double _limitDegrees(double degrees) { + degrees /= 360; + double limited = 360 * (degrees - degrees.floor()); + if (limited < 0) limited += 360; + return limited; +} + +/// Note: implements JS behavior: ((a*x + b) + c)*x + d (matches nrel-spa JS) +double _thirdOrderPolynomial(double a, double b, double c, double d, double x) { + return ((a * x + b) + c) * x + d; +} + +double _limitDegrees180pm(double degrees) { + degrees /= 360.0; + double limited = 360.0 * (degrees - degrees.floor()); + if (limited < -180.0) { + limited += 360.0; + } else if (limited > 180.0) { + limited -= 360.0; + } + return limited; +} + +double _limitDegrees180(double degrees) { + degrees /= 180.0; + double limited = 180.0 * (degrees - degrees.floor()); + if (limited < 0) limited += 180.0; + return limited; +} + +double _limitZero2one(double value) { + double limited = value - value.floor(); + if (limited < 0) limited += 1.0; + return limited; +} + +double _limitMinutes(double minutes) { + double limited = minutes; + if (limited < -20.0) { + limited += 1440.0; + } else if (limited > 20.0) { + limited -= 1440.0; + } + return limited; +} + +// ─── Geometric functions ───────────────────────────────────────────────────── + +double _geocentricRightAscension(double lamda, double epsilon, double beta) { + final lambdaRad = _deg2rad(lamda); + final epsilonRad = _deg2rad(epsilon); + return _limitDegrees( + _rad2deg( + atan2( + sin(lambdaRad) * cos(epsilonRad) - + tan(_deg2rad(beta)) * sin(epsilonRad), + cos(lambdaRad), + ), + ), + ); +} + +double _geocentricDeclination(double beta, double epsilon, double lamda) { + final betaRad = _deg2rad(beta); + final epsilonRad = _deg2rad(epsilon); + return _rad2deg( + asin( + sin(betaRad) * cos(epsilonRad) + + cos(betaRad) * sin(epsilonRad) * sin(_deg2rad(lamda)), + ), + ); +} + +double _observerHourAngle(double nu, double longitude, double alphaDeg) { + return _limitDegrees(nu + longitude - alphaDeg); +} + +// Mutates dltap fields via a two-element list [deltaAlpha, deltaPrime] +void _rightAscensionParallaxAndTopocentricDec( + double latitude, + double elevation, + double xi, + double h, + double delta, + List dltap, // [0]=deltaAlpha, [1]=deltaPrime +) { + final latRad = _deg2rad(latitude); + final xiRad = _deg2rad(xi); + final hRad = _deg2rad(h); + final deltaRad = _deg2rad(delta); + final u = atan(0.99664719 * tan(latRad)); + final y = 0.99664719 * sin(u) + elevation * sin(latRad) / 6378140.0; + final x = cos(u) + elevation * cos(latRad) / 6378140.0; + final deltaAlphaRad = atan2( + -x * sin(xiRad) * sin(hRad), + cos(deltaRad) - x * sin(xiRad) * cos(hRad), + ); + dltap[1] = _rad2deg( + atan2( + (sin(deltaRad) - y * sin(xiRad)) * cos(deltaAlphaRad), + cos(deltaRad) - x * sin(xiRad) * cos(hRad), + ), + ); + dltap[0] = _rad2deg(deltaAlphaRad); +} + +double _topocentricElevationAngle( + double latitude, + double deltaPrime, + double hPrime, +) { + final latRad = _deg2rad(latitude); + final deltaPrimeRad = _deg2rad(deltaPrime); + return _rad2deg( + asin( + sin(latRad) * sin(deltaPrimeRad) + + cos(latRad) * cos(deltaPrimeRad) * cos(_deg2rad(hPrime)), + ), + ); +} + +double _atmosphericRefractionCorrection( + double pressure, + double temperature, + double atmosRefract, + double e0, +) { + double delE = 0; + if (e0 >= -1 * (_sunRadius + atmosRefract)) { + delE = + (pressure / 1010.0) * + (283.0 / (273.0 + temperature)) * + 1.02 / + (60.0 * tan(_deg2rad(e0 + 10.3 / (e0 + 5.11)))); + } + return delE; +} + +double _topocentricAzimuthAngleAstro( + double hPrime, + double latitude, + double deltaPrime, +) { + final hPrimeRad = _deg2rad(hPrime); + final latRad = _deg2rad(latitude); + return _limitDegrees( + _rad2deg( + atan2( + sin(hPrimeRad), + cos(hPrimeRad) * sin(latRad) - tan(_deg2rad(deltaPrime)) * cos(latRad), + ), + ), + ); +} + +double _surfaceIncidenceAngle( + double zenith, + double azimuthAstro, + double azmRotation, + double slope, +) { + final zenithRad = _deg2rad(zenith); + final slopeRad = _deg2rad(slope); + return _rad2deg( + acos( + cos(zenithRad) * cos(slopeRad) + + sin(slopeRad) * + sin(zenithRad) * + cos(_deg2rad(azimuthAstro - azmRotation)), + ), + ); +} + +// ─── Julian / time functions ───────────────────────────────────────────────── + +double _julianDay( + int year, + int month, + int day, + int hour, + int minute, + double second, + double dut1, + double tz, +) { + final dayDecimal = + day + (hour - tz + (minute + (second + dut1) / 60.0) / 60.0) / 24.0; + int y = year, m = month; + if (m < 3) { + m += 12; + y--; + } + double jd = + (365.25 * (y + 4716.0)).floor().toDouble() + + (30.6001 * (m + 1)).floor().toDouble() + + dayDecimal - + 1524.5; + if (jd > 2299160.0) { + final a = (y / 100).floor().toDouble(); + jd += 2 - a + (a / 4).floor(); + } + return jd; +} + +double _julianCentury(double jd) => (jd - 2451545.0) / 36525.0; +double _julianEphemerisDay(double jd, double deltaT) => jd + deltaT / 86400.0; +double _julianEphemerisCentury(double jde) => (jde - 2451545.0) / 36525.0; +double _julianEphemerisMillennium(double jce) => jce / 10.0; + +// ─── Periodic term summation ───────────────────────────────────────────────── + +double _earthPeriodicTermSummation( + List> terms, + int count, + double jme, +) { + double sum = 0; + for (int i = 0; i < count; i++) { + sum += terms[i][0] * cos(terms[i][1] + terms[i][2] * jme); + } + return sum; +} + +double _earthValues(List termSum, int count, double jme) { + double sum = 0; + for (int i = 0; i < count; i++) { + sum += termSum[i] * pow(jme, i).toDouble(); + } + return sum / 1.0e8; +} + +double _earthHeliocentricLongitude(double jme) { + final sum = List.filled(6, 0); + for (int i = 0; i < 6; i++) { + sum[i] = _earthPeriodicTermSummation(_lTerms[i], _lSubcount[i], jme); + } + return _limitDegrees(_rad2deg(_earthValues(sum, 6, jme))); +} + +double _earthHeliocentricLatitude(double jme) { + final sum = List.filled(2, 0); + for (int i = 0; i < 2; i++) { + sum[i] = _earthPeriodicTermSummation(_bTerms[i], _bSubcount[i], jme); + } + return _rad2deg(_earthValues(sum, 2, jme)); +} + +double _earthRadiusVector(double jme) { + final sum = List.filled(5, 0); + for (int i = 0; i < 5; i++) { + sum[i] = _earthPeriodicTermSummation(_rTerms[i], _rSubcount[i], jme); + } + return _earthValues(sum, 5, jme); +} + +double _geocentricLongitude(double l) { + double theta = l + 180.0; + if (theta >= 360.0) theta -= 360.0; + return theta; +} + +// ─── X anomaly terms ───────────────────────────────────────────────────────── + +double _meanElongationMoonSun(double jce) => _thirdOrderPolynomial( + 1.0 / 189474.0, + -0.0019142, + 445267.11148, + 297.85036, + jce, +); +double _meanAnomalySun(double jce) => _thirdOrderPolynomial( + -1.0 / 300000.0, + -0.0001603, + 35999.05034, + 357.52772, + jce, +); +double _meanAnomalyMoon(double jce) => _thirdOrderPolynomial( + 1.0 / 56250.0, + 0.0086972, + 477198.867398, + 134.96298, + jce, +); +double _argumentLatitudeMoon(double jce) => _thirdOrderPolynomial( + 1.0 / 327270.0, + -0.0036825, + 483202.017538, + 93.27191, + jce, +); +double _ascendingLongitudeMoon(double jce) => _thirdOrderPolynomial( + 1.0 / 450000.0, + 0.0020708, + -1934.136261, + 125.04452, + jce, +); + +// ─── Nutation ──────────────────────────────────────────────────────────────── + +double _xyTermSummation(int i, List x) { + double sum = 0; + for (int j = 0; j < 5; j++) { + sum += x[j] * _yTerms[i][j]; + } + return sum; +} + +void _nutationLongitudeAndObliquity(double jce, List x, _Spa spa) { + double sumPsi = 0; + double sumEpsilon = 0; + for (int i = 0; i < 63; i++) { + final xyTermSum = _deg2rad(_xyTermSummation(i, x)); + sumPsi += (_peTerms[i][0] + jce * _peTerms[i][1]) * sin(xyTermSum); + sumEpsilon += (_peTerms[i][2] + jce * _peTerms[i][3]) * cos(xyTermSum); + } + spa.delPsi = sumPsi / 36000000.0; + spa.delEpsilon = sumEpsilon / 36000000.0; +} + +double _eclipticMeanObliquity(double jme) { + final u = jme / 10.0; + return 84381.448 + + u * + (-4680.93 + + u * + (-1.55 + + u * + (1999.25 + + u * + (-51.38 + + u * + (-249.67 + + u * + (-39.05 + + u * + (7.12 + + u * + (27.87 + + u * + (5.79 + + u * 2.45))))))))); +} + +// ─── Sidereal & apparent sun ───────────────────────────────────────────────── + +double _greenwichMeanSiderealTime(double jd, double jc) { + return _limitDegrees( + 280.46061837 + + 360.98564736629 * (jd - 2451545.0) + + jc * jc * (0.000387933 - jc / 38710000.0), + ); +} + +double _sunMeanLongitude(double jme) { + return _limitDegrees( + 280.4664567 + + jme * + (360007.6982779 + + jme * + (0.03032028 + + jme * + (1 / 49931.0 + + jme * + (-1 / 15300.0 + jme * (-1 / 2000000.0))))), + ); +} + +double _eot(double m, double alpha, double delPsi, double epsilon) { + return _limitMinutes( + 4.0 * (m - 0.0057183 - alpha + delPsi * cos(_deg2rad(epsilon))), + ); +} + +// ─── RTS (Rise/Transit/Set) ────────────────────────────────────────────────── + +double _approxSunTransitTime(double alphaZero, double longitude, double nu) { + return (alphaZero - longitude - nu) / 360.0; +} + +double _sunHourAngleAtRiseSet( + double latitude, + double deltaZero, + double h0Prime, +) { + double h0 = -99999; + final latRad = _deg2rad(latitude); + final deltaZeroRad = _deg2rad(deltaZero); + final argument = + (sin(_deg2rad(h0Prime)) - sin(latRad) * sin(deltaZeroRad)) / + (cos(latRad) * cos(deltaZeroRad)); + if (argument.abs() <= 1) { + h0 = _limitDegrees180(_rad2deg(acos(argument))); + } + return h0; +} + +void _approxSunRiseAndSet(List mRts, double h0) { + final h0Dfrac = h0 / 360.0; + mRts[1] = _limitZero2one(mRts[0] - h0Dfrac); // SUN_RISE + mRts[2] = _limitZero2one(mRts[0] + h0Dfrac); // SUN_SET + mRts[0] = _limitZero2one(mRts[0]); // SUN_TRANSIT +} + +double _rtsAlphaDeltaPrime(List ad, double n) { + double a = ad[1] - ad[0]; // JD_ZERO - JD_MINUS + double b = ad[2] - ad[1]; // JD_PLUS - JD_ZERO + if (a.abs() >= 2.0) a = _limitZero2one(a); + if (b.abs() >= 2.0) b = _limitZero2one(b); + return ad[1] + n * (a + b + (b - a) * n) / 2.0; +} + +double _rtsSunAltitude(double latitude, double deltaPrime, double hPrime) { + final latRad = _deg2rad(latitude); + final deltaPrimeRad = _deg2rad(deltaPrime); + return _rad2deg( + asin( + sin(latRad) * sin(deltaPrimeRad) + + cos(latRad) * cos(deltaPrimeRad) * cos(_deg2rad(hPrime)), + ), + ); +} + +double _sunRiseAndSet( + List mRts, + List hRts, + List deltaPrime, + double latitude, + List hPrime, + double h0Prime, + int sun, +) { + return mRts[sun] + + (hRts[sun] - h0Prime) / + (360.0 * + cos(_deg2rad(deltaPrime[sun])) * + cos(_deg2rad(latitude)) * + sin(_deg2rad(hPrime[sun]))); +} + +double _dayfracToLocalHr(double dayfrac, double timezone) { + return 24.0 * _limitZero2one(dayfrac + timezone / 24.0); +} + +// ─── Core geocentric calculation ───────────────────────────────────────────── + +void _calculateGeocentricSunRightAscensionAndDeclination(_Spa spa) { + spa.jc = _julianCentury(spa.jd); + spa.jde = _julianEphemerisDay(spa.jd, spa.deltaT); + spa.jce = _julianEphemerisCentury(spa.jde); + spa.jme = _julianEphemerisMillennium(spa.jce); + + spa.l = _earthHeliocentricLongitude(spa.jme); + spa.b = _earthHeliocentricLatitude(spa.jme); + spa.r = _earthRadiusVector(spa.jme); + + spa.theta = _geocentricLongitude(spa.l); + spa.beta = -spa.b; // geocentric latitude + + final x = [ + spa.x0 = _meanElongationMoonSun(spa.jce), + spa.x1 = _meanAnomalySun(spa.jce), + spa.x2 = _meanAnomalyMoon(spa.jce), + spa.x3 = _argumentLatitudeMoon(spa.jce), + spa.x4 = _ascendingLongitudeMoon(spa.jce), + ]; + + _nutationLongitudeAndObliquity(spa.jce, x, spa); + + spa.epsilon0 = _eclipticMeanObliquity(spa.jme); + spa.epsilon = spa.delEpsilon + spa.epsilon0 / 3600.0; + spa.delTau = -20.4898 / (3600.0 * spa.r); + spa.lamda = spa.theta + spa.delPsi + spa.delTau; + spa.nu0 = _greenwichMeanSiderealTime(spa.jd, spa.jc); + spa.nu = spa.nu0 + spa.delPsi * cos(_deg2rad(spa.epsilon)); + spa.alpha = _geocentricRightAscension(spa.lamda, spa.epsilon, spa.beta); + spa.delta = _geocentricDeclination(spa.beta, spa.epsilon, spa.lamda); +} + +// ─── EOT + RTS ─────────────────────────────────────────────────────────────── + +void _calculateEotAndSunRiseTransitSet(_Spa spa) { + final h0Prime = -1 * (_sunRadius + spa.atmosRefract); + + final sunRts = spa.clone(); + sunRts.hour = sunRts.minute = 0; + sunRts.second = 0.0; + sunRts.deltaUt1 = sunRts.timezone = 0.0; + sunRts.jd = _julianDay( + sunRts.year, + sunRts.month, + sunRts.day, + sunRts.hour, + sunRts.minute, + sunRts.second, + sunRts.deltaUt1, + sunRts.timezone, + ); + + final m = _sunMeanLongitude(spa.jme); + spa.eot = _eot(m, spa.alpha, spa.delPsi, spa.epsilon); + + _calculateGeocentricSunRightAscensionAndDeclination(sunRts); + final nu = sunRts.nu; + sunRts.deltaT = 0; + + // Compute alpha and delta for JD-1, JD, JD+1 + // indices: [0]=JD_MINUS, [1]=JD_ZERO, [2]=JD_PLUS + final alpha = List.filled(3, 0); + final delta = List.filled(3, 0); + sunRts.jd--; + for (int i = 0; i < 3; i++) { + _calculateGeocentricSunRightAscensionAndDeclination(sunRts); + alpha[i] = sunRts.alpha; + delta[i] = sunRts.delta; + sunRts.jd++; + } + + // mRts: [0]=transit, [1]=rise, [2]=set + final mRts = List.filled(3, 0); + mRts[0] = _approxSunTransitTime(alpha[1], spa.longitude, nu); + final h0 = _sunHourAngleAtRiseSet(spa.latitude, delta[1], h0Prime); + + if (h0 >= 0) { + _approxSunRiseAndSet(mRts, h0); + + final nuRts = List.filled(3, 0); + final alphaPrime = List.filled(3, 0); + final deltaPrime = List.filled(3, 0); + final hPrime = List.filled(3, 0); + final hRts = List.filled(3, 0); + + for (int i = 0; i < 3; i++) { + nuRts[i] = nu + 360.985647 * mRts[i]; + final n = mRts[i] + spa.deltaT / 86400.0; + alphaPrime[i] = _rtsAlphaDeltaPrime(alpha, n); + deltaPrime[i] = _rtsAlphaDeltaPrime(delta, n); + hPrime[i] = _limitDegrees180pm(nuRts[i] + spa.longitude - alphaPrime[i]); + hRts[i] = _rtsSunAltitude(spa.latitude, deltaPrime[i], hPrime[i]); + } + + spa.srha = hPrime[1]; // SUN_RISE + spa.ssha = hPrime[2]; // SUN_SET + spa.sta = hRts[0]; // SUN_TRANSIT + + spa.suntransit = _dayfracToLocalHr( + mRts[0] - hPrime[0] / 360.0, + spa.timezone, + ); + spa.sunrise = _dayfracToLocalHr( + _sunRiseAndSet(mRts, hRts, deltaPrime, spa.latitude, hPrime, h0Prime, 1), + spa.timezone, + ); + spa.sunset = _dayfracToLocalHr( + _sunRiseAndSet(mRts, hRts, deltaPrime, spa.latitude, hPrime, h0Prime, 2), + spa.timezone, + ); + } else { + spa.srha = + spa.ssha = spa.sta = spa.suntransit = spa.sunrise = spa.sunset = -99999; + } +} + +// ─── Full SPA calculation ──────────────────────────────────────────────────── + +int _spaCalculate(_Spa spa) { + // Validate inputs + if (spa.year < -2000 || spa.year > 6000) return 1; + if (spa.month < 1 || spa.month > 12) return 2; + if (spa.day < 1 || spa.day > 31) return 3; + if (spa.hour < 0 || spa.hour > 24) return 4; + if (spa.minute < 0 || spa.minute > 59) return 5; + if (spa.second < 0 || spa.second >= 60) return 6; + if (spa.pressure < 0 || spa.pressure > 5000) return 12; + if (spa.temperature <= -273 || spa.temperature > 6000) return 13; + if (spa.deltaUt1 <= -1 || spa.deltaUt1 >= 1) return 17; + if (spa.hour == 24 && spa.minute > 0) return 5; + if (spa.hour == 24 && spa.second > 0) return 6; + if (spa.deltaT.abs() > 8000) return 7; + if (spa.timezone.abs() > 18) return 8; + if (spa.longitude.abs() > 180) return 9; + if (spa.latitude.abs() > 90) return 10; + if (spa.atmosRefract.abs() > 5) return 16; + if (spa.elevation < -6500000) return 11; + + spa.jd = _julianDay( + spa.year, + spa.month, + spa.day, + spa.hour, + spa.minute, + spa.second, + spa.deltaUt1, + spa.timezone, + ); + + _calculateGeocentricSunRightAscensionAndDeclination(spa); + + spa.h = _observerHourAngle(spa.nu, spa.longitude, spa.alpha); + spa.xi = 8.794 / (3600.0 * spa.r); // sun equatorial horizontal parallax + + final dltap = [spa.delAlpha, spa.deltaPrime]; + _rightAscensionParallaxAndTopocentricDec( + spa.latitude, + spa.elevation, + spa.xi, + spa.h, + spa.delta, + dltap, + ); + spa.delAlpha = dltap[0]; + spa.deltaPrime = dltap[1]; + + spa.alphaPrime = spa.alpha + spa.delAlpha; + spa.hPrime = spa.h - spa.delAlpha; + + spa.e0 = _topocentricElevationAngle(spa.latitude, spa.deltaPrime, spa.hPrime); + spa.delE = _atmosphericRefractionCorrection( + spa.pressure, + spa.temperature, + spa.atmosRefract, + spa.e0, + ); + spa.e = spa.e0 + spa.delE; + + spa.zenith = 90.0 - spa.e; + spa.azimuthAstro = _topocentricAzimuthAngleAstro( + spa.hPrime, + spa.latitude, + spa.deltaPrime, + ); + spa.azimuth = _limitDegrees(spa.azimuthAstro + 180.0); + + if (spa.function == 1 || spa.function == 3) { + // SPA_ZA_INC or SPA_ALL: compute surface incidence + spa.incidence = _surfaceIncidenceAngle( + spa.zenith, + spa.azimuthAstro, + spa.azmRotation, + spa.slope, + ); + } + + if (spa.function == _sPaZaRts || spa.function == _sPaAll) { + _calculateEotAndSunRiseTransitSet(spa); + } + + return 0; +} + +// ─── Custom angle adjustment ───────────────────────────────────────────────── + +SpaAnglesResult _adjustForCustomAngle(_Spa base, double zenithAngle) { + final phi = base.latitude * pi / 180; + final delta = base.delta * pi / 180; + final z = zenithAngle * pi / 180; + final cosH0 = (cos(z) - sin(phi) * sin(delta)) / (cos(phi) * cos(delta)); + if (cosH0 < -1 || cosH0 > 1) { + return SpaAnglesResult(sunrise: double.nan, sunset: double.nan); + } + final h0h = acos(cosH0) * 180 / pi / 15; + return SpaAnglesResult( + sunrise: base.suntransit - h0h, + sunset: base.suntransit + h0h, + ); +} + +// ─── Public API ────────────────────────────────────────────────────────────── + +/// Compute solar position for the given parameters. +/// +/// [date] is used in UTC components. +/// [latitude] is in degrees (−90 to 90, south = negative). +/// [longitude] is in degrees (−180 to 180, west = negative). +/// [timezone] is hours from UTC (e.g., −5 for EST). +/// [customAngles] are zenith angles in degrees for which rise/set times +/// should be calculated (e.g., [96, 102] for civil/nautical twilight). +SpaResult getSpa( + DateTime date, + double latitude, + double longitude, + double timezone, { + double elevation = 0, + double pressure = 1013, + double temperature = 15, + double deltaUt1 = 0, + double deltaT = 67, + double slope = 0, + double azmRotation = 0, + double atmosRefract = 0.5667, + List customAngles = const [], +}) { + final d = _Spa(); + d.year = date.toUtc().year; + d.month = date.toUtc().month; + d.day = date.toUtc().day; + d.hour = date.toUtc().hour; + d.minute = date.toUtc().minute; + d.second = date.toUtc().second.toDouble(); + d.longitude = longitude; + d.latitude = latitude; + d.timezone = timezone; + d.elevation = elevation; + d.pressure = pressure; + d.temperature = temperature; + d.deltaUt1 = deltaUt1; + d.deltaT = deltaT; + d.slope = slope; + d.azmRotation = azmRotation; + d.atmosRefract = atmosRefract; + d.function = _sPaZaRts; + + final rc = _spaCalculate(d); + if (rc != 0) { + throw ArgumentError('SPA calculation failed (error code $rc)'); + } + + final angleResults = customAngles + .map((z) => _adjustForCustomAngle(d, z)) + .toList(growable: false); + + return SpaResult( + zenith: d.zenith, + azimuth: d.azimuth, + sunrise: d.sunrise, + solarNoon: d.suntransit, + sunset: d.sunset, + angles: angleResults, + ); +} diff --git a/lib/src/types.dart b/lib/src/types.dart new file mode 100644 index 0000000..1517466 --- /dev/null +++ b/lib/src/types.dart @@ -0,0 +1,142 @@ +/// Core types for pray_calc_dart. +library; + +/// Asr shadow convention: Shafi'i (1x) or Hanafi (2x). +enum AsrConvention { shafii, hanafi } + +/// Shafaq variant for MSC Isha model. +enum ShafaqMode { general, ahmer, abyad } + +/// Computed twilight depression angles for Fajr and Isha. +class TwilightAngles { + /// Solar depression angle for Fajr (positive degrees below horizon). + final double fajrAngle; + + /// Solar depression angle for Isha (positive degrees below horizon). + final double ishaAngle; + + const TwilightAngles({required this.fajrAngle, required this.ishaAngle}); +} + +/// Raw prayer times as fractional hours. NaN = unreachable event. +class PrayerTimes { + /// Start of the last third of the night (Qiyam al-Layl). + final double qiyam; + + /// True dawn (Subh Sadiq). + final double fajr; + + /// Astronomical sunrise. + final double sunrise; + + /// Solar noon (exact geometric transit). + final double noon; + + /// Dhuhr (2.5 minutes after solar noon). + final double dhuhr; + + /// Asr (Shafi'i or Hanafi shadow convention). + final double asr; + + /// Maghrib (sunset). + final double maghrib; + + /// Isha (nightfall, end of shafaq). + final double isha; + + /// Dynamic twilight angles used for this calculation. + final TwilightAngles angles; + + const PrayerTimes({ + required this.qiyam, + required this.fajr, + required this.sunrise, + required this.noon, + required this.dhuhr, + required this.asr, + required this.maghrib, + required this.isha, + required this.angles, + }); +} + +/// Prayer times formatted as HH:MM:SS strings. +class FormattedPrayerTimes { + final String qiyam; + final String fajr; + final String sunrise; + final String noon; + final String dhuhr; + final String asr; + final String maghrib; + final String isha; + final TwilightAngles angles; + + const FormattedPrayerTimes({ + required this.qiyam, + required this.fajr, + required this.sunrise, + required this.noon, + required this.dhuhr, + required this.asr, + required this.maghrib, + required this.isha, + required this.angles, + }); +} + +/// Solar ephemeris result. +class SolarEphemeris { + /// Solar declination in degrees. + final double decl; + + /// Earth-Sun distance in AU. + final double r; + + /// Apparent solar ecliptic longitude in radians (0–2π). + final double eclLon; + + const SolarEphemeris({ + required this.decl, + required this.r, + required this.eclLon, + }); +} + +/// SPA result from the NREL Solar Position Algorithm. +class SpaResult { + /// Topocentric zenith angle in degrees. + final double zenith; + + /// Topocentric azimuth angle, eastward from north, in degrees. + final double azimuth; + + /// Local sunrise time as fractional hours (NaN if polar). + final double sunrise; + + /// Local sun transit time (solar noon) as fractional hours. + final double solarNoon; + + /// Local sunset time as fractional hours (NaN if polar). + final double sunset; + + /// Custom zenith angle results (one per angle in the input list). + final List angles; + + const SpaResult({ + required this.zenith, + required this.azimuth, + required this.sunrise, + required this.solarNoon, + required this.sunset, + this.angles = const [], + }); +} + +/// Sunrise/sunset pair for a custom zenith angle. +class SpaAnglesResult { + final double sunrise; + final double sunset; + + const SpaAnglesResult({required this.sunrise, required this.sunset}); +} diff --git a/pubspec.yaml b/pubspec.yaml new file mode 100644 index 0000000..19871c1 --- /dev/null +++ b/pubspec.yaml @@ -0,0 +1,21 @@ +name: pray_calc_dart +description: > + Islamic prayer times for Dart and Flutter. Pure Dart port of the pray-calc + library implementing the NREL Solar Position Algorithm, MCW seasonal model, + and dynamic twilight angles. Zero dependencies. +version: 1.0.0 +repository: https://github.com/acamarata/pray-calc-dart +issue_tracker: https://github.com/acamarata/pray-calc-dart/issues +topics: + - prayer-times + - islamic + - solar + - astronomy + - qibla + +environment: + sdk: ^3.7.0 + +dev_dependencies: + lints: ^5.0.0 + test: ^1.25.8 diff --git a/test/pray_calc_dart_test.dart b/test/pray_calc_dart_test.dart new file mode 100644 index 0000000..f6c316f --- /dev/null +++ b/test/pray_calc_dart_test.dart @@ -0,0 +1,193 @@ +import 'package:pray_calc_dart/pray_calc_dart.dart'; +import 'package:test/test.dart'; + +/// Reference values validated against the pray-calc TypeScript library v2.0.0. +void main() { + group('getTimes — NYC 2024-03-15 (Shafi\'i)', () { + late PrayerTimes times; + + setUpAll(() { + final date = DateTime(2024, 3, 15); + times = getTimes(date, 40.7128, -74.0060, -5.0); + }); + + test('Fajr is before Sunrise', () { + expect(times.fajr, lessThan(times.sunrise)); + }); + + test('Sunrise is before Dhuhr', () { + expect(times.sunrise, lessThan(times.dhuhr)); + }); + + test('Dhuhr is before Asr', () { + expect(times.dhuhr, lessThan(times.asr)); + }); + + test('Asr is before Maghrib', () { + expect(times.asr, lessThan(times.maghrib)); + }); + + test('Maghrib is before Isha', () { + expect(times.maghrib, lessThan(times.isha)); + }); + + test('Fajr is in the 4–6 AM range', () { + // Dynamic method gives ~4:37 AM for NYC March 15 (MCW-based ~17.7° angle) + expect(times.fajr, greaterThan(4.0)); + expect(times.fajr, lessThan(6.0)); + }); + + test('Dhuhr is around 12–14 h', () { + expect(times.dhuhr, greaterThan(12.0)); + expect(times.dhuhr, lessThan(14.0)); + }); + + test('Maghrib is around 18–19.5 h', () { + expect(times.maghrib, greaterThan(18.0)); + expect(times.maghrib, lessThan(19.5)); + }); + + test('Isha is after Maghrib + 1 hour', () { + expect(times.isha, greaterThan(times.maghrib + 1.0)); + }); + + test('Qiyam is finite and before Fajr (wraps past midnight)', () { + // Qiyam = last third of night, after Isha it wraps past midnight + // e.g. Isha=19:34, night=7h, Qiyam=01:31 next day (numerically ~1.53, < Fajr ~4.62) + expect(times.qiyam.isFinite, isTrue); + expect(times.qiyam, lessThan(times.fajr)); + }); + + test('angles are in valid range [10, 22]', () { + expect(times.angles.fajrAngle, inInclusiveRange(10.0, 22.0)); + expect(times.angles.ishaAngle, inInclusiveRange(10.0, 22.0)); + }); + + test('formatTime produces HH:MM:SS', () { + final formatted = formatTime(times.fajr); + expect(formatted, matches(RegExp(r'^\d{2}:\d{2}:\d{2}$'))); + }); + }); + + group('getTimes — NYC 2024-03-15 (Hanafi)', () { + late PrayerTimes timesShafii; + late PrayerTimes timesHanafi; + + setUpAll(() { + final date = DateTime(2024, 3, 15); + timesShafii = getTimes(date, 40.7128, -74.0060, -5.0); + timesHanafi = getTimes(date, 40.7128, -74.0060, -5.0, hanafi: true); + }); + + test('Hanafi Asr is later than Shafi\'i Asr', () { + expect(timesHanafi.asr, greaterThan(timesShafii.asr)); + }); + + test('All non-Asr times are identical', () { + expect(timesHanafi.fajr, closeTo(timesShafii.fajr, 0.0001)); + expect(timesHanafi.sunrise, closeTo(timesShafii.sunrise, 0.0001)); + expect(timesHanafi.maghrib, closeTo(timesShafii.maghrib, 0.0001)); + expect(timesHanafi.isha, closeTo(timesShafii.isha, 0.0001)); + }); + }); + + group('getTimes — Mecca 2024-06-21 (summer solstice)', () { + late PrayerTimes times; + + setUpAll(() { + final date = DateTime(2024, 6, 21); + times = getTimes(date, 21.3891, 39.8579, 3.0); + }); + + test('All prayer times are finite', () { + expect(times.fajr.isFinite, isTrue); + expect(times.sunrise.isFinite, isTrue); + expect(times.dhuhr.isFinite, isTrue); + expect(times.asr.isFinite, isTrue); + expect(times.maghrib.isFinite, isTrue); + expect(times.isha.isFinite, isTrue); + expect(times.qiyam.isFinite, isTrue); + }); + + test('Prayer time ordering is correct', () { + expect(times.fajr, lessThan(times.sunrise)); + expect(times.sunrise, lessThan(times.dhuhr)); + expect(times.dhuhr, lessThan(times.asr)); + expect(times.asr, lessThan(times.maghrib)); + expect(times.maghrib, lessThan(times.isha)); + }); + }); + + group('getAngles', () { + test('NYC Jan angles are in valid range', () { + final date = DateTime(2024, 1, 15); + final angles = getAngles(date, 40.7128, -74.0060); + expect(angles.fajrAngle, inInclusiveRange(10.0, 22.0)); + expect(angles.ishaAngle, inInclusiveRange(10.0, 22.0)); + }); + }); + + group('getAsr', () { + test('Hanafi Asr is always later than Shafi\'i', () { + final asrShafii = getAsr(12.5, 40.0, 5.0); + final asrHanafi = getAsr(12.5, 40.0, 5.0, hanafi: true); + expect(asrHanafi, greaterThan(asrShafii)); + }); + }); + + group('getQiyam', () { + test('last third starts at 2/3 of night from Isha', () { + // Isha=22:00, Fajr=05:00 next day → night=7h → last third = 22 + 14/3 + final q = getQiyam(5.0, 22.0); + expect(q, closeTo(((22.0 + 14.0 / 3.0) - 24.0), 0.001)); + }); + }); + + group('solarEphemeris', () { + test('declination at June solstice is ~23.4°', () { + final jd = toJulianDate(DateTime.utc(2024, 6, 21, 12)); + final eph = solarEphemeris(jd); + expect(eph.decl, closeTo(23.4, 0.5)); + }); + + test('declination at Dec solstice is ~-23.4°', () { + final jd = toJulianDate(DateTime.utc(2024, 12, 21, 12)); + final eph = solarEphemeris(jd); + expect(eph.decl, closeTo(-23.4, 0.5)); + }); + + test('Earth-Sun distance at perihelion (Jan 3) is ~0.983 AU', () { + final jd = toJulianDate(DateTime.utc(2024, 1, 3, 12)); + final eph = solarEphemeris(jd); + expect(eph.r, closeTo(0.983, 0.003)); + }); + }); + + group('getSpa', () { + test('returns valid zenith and azimuth', () { + final result = getSpa( + DateTime.utc(2024, 3, 15, 12, 0, 0), + 40.7128, + -74.0060, + -5.0, + ); + expect(result.zenith, inInclusiveRange(0.0, 180.0)); + expect(result.azimuth, inInclusiveRange(0.0, 360.0)); + }); + + test('custom angles produce correct twilight pairs', () { + final result = getSpa( + DateTime.utc(2024, 3, 15, 12, 0, 0), + 40.7128, + -74.0060, + -5.0, + customAngles: [96.0, 108.0], // civil, astronomical twilight + ); + expect(result.angles.length, equals(2)); + // Civil twilight (96°, -6° below horizon) rises LATER than astronomical (108°, -18°) + expect(result.angles[0].sunrise, greaterThan(result.angles[1].sunrise)); + // Civil twilight sets EARLIER than astronomical twilight + expect(result.angles[0].sunset, lessThan(result.angles[1].sunset)); + }); + }); +}