refactor(pkg/diyanetcalc): Encapsulate the algo + add adapter
parent
9cfa7f2dcf
commit
82f252f32f
@ -0,0 +1,244 @@
|
|||||||
|
package diyanetcalc
|
||||||
|
|
||||||
|
import (
|
||||||
|
"iter"
|
||||||
|
"math"
|
||||||
|
"time"
|
||||||
|
|
||||||
|
"prayertimes/pkg/hijricalendar"
|
||||||
|
|
||||||
|
"github.com/samber/lo"
|
||||||
|
)
|
||||||
|
|
||||||
|
const (
|
||||||
|
daysToGenerate = 30
|
||||||
|
degPerHour = 15.0
|
||||||
|
j2000 = 2451545.0 // Julian Date of J2000.0 epoch (Jan 1.5, 2000)
|
||||||
|
)
|
||||||
|
|
||||||
|
// Diyanet angular criteria (post-1983 reform)
|
||||||
|
//
|
||||||
|
// imsakAngle: Sun 18° below horizon = start of astronomical twilight (Fajr).
|
||||||
|
// Pre-1983 used -19° + temkin buffer; now -18° with zero temkin.
|
||||||
|
const imsakAngle = -18.0
|
||||||
|
|
||||||
|
// ishaAngle: Sun 17° below horizon = shafaq al-ahmar (red twilight) gone.
|
||||||
|
const ishaAngle = -17.0
|
||||||
|
|
||||||
|
// sunAngle: combined refraction (~0.567°) + solar semi-diameter (~0.267°)
|
||||||
|
// correction applied at sunrise/sunset. Equivalent to 50 arcmin.
|
||||||
|
const sunAngle = -0.833
|
||||||
|
|
||||||
|
// Temkin - precautionary time buffers (minutes, post-1983 standardized values)
|
||||||
|
//
|
||||||
|
// Ensures a single published time remains valid across the full geographical
|
||||||
|
// extent of a city. Pre-1983 values were 10-20 min; the reform moderated them.
|
||||||
|
//
|
||||||
|
// Imsak: 0 - no buffer; avoids starting Fajr too early / fast too late
|
||||||
|
// Sunrise: -7 - subtracted, ensuring the Sun has fully cleared the horizon
|
||||||
|
// Dhuhr: +5 - Sun has clearly begun its descent
|
||||||
|
// Asr: +4 - accounts for local elevation and horizon obstacles
|
||||||
|
// Maghrib: +7 - Sun has completely set before breaking fast
|
||||||
|
// Isha: 0 - no buffer needed at this twilight stage
|
||||||
|
var temkin = struct{ Imsak, Sunrise, Dhuhr, Asr, Maghrib, Isha float64 }{
|
||||||
|
Imsak: 0, Sunrise: -7, Dhuhr: 5, Asr: 4, Maghrib: 7, Isha: 0,
|
||||||
|
}
|
||||||
|
|
||||||
|
type CalculateParams struct {
|
||||||
|
Latitude float64
|
||||||
|
Longitude float64
|
||||||
|
StartingDay time.Time
|
||||||
|
}
|
||||||
|
|
||||||
|
type Times struct {
|
||||||
|
Date time.Time
|
||||||
|
DateHijri string
|
||||||
|
Fajr time.Time
|
||||||
|
Sunrise time.Time
|
||||||
|
Dhuhr time.Time
|
||||||
|
Asr time.Time
|
||||||
|
Sunset time.Time
|
||||||
|
Maghrib time.Time
|
||||||
|
Isha time.Time
|
||||||
|
}
|
||||||
|
|
||||||
|
func CalculatePrayerTimes(params CalculateParams) iter.Seq[Times] {
|
||||||
|
startingDay := params.StartingDay.UTC().Truncate(24 * time.Hour)
|
||||||
|
if startingDay.IsZero() {
|
||||||
|
startingDay = time.Now().UTC().Truncate(24 * time.Hour)
|
||||||
|
}
|
||||||
|
|
||||||
|
return func(yield func(Times) bool) {
|
||||||
|
for day := startingDay; ; day = day.AddDate(0, 0, 1) {
|
||||||
|
c := prayerTimes(params.Latitude, params.Longitude, day)
|
||||||
|
if !yield(Times{
|
||||||
|
Date: day,
|
||||||
|
DateHijri: hijricalendar.ToISODate(day),
|
||||||
|
Fajr: lo.FromPtr(c.Imsak),
|
||||||
|
Sunrise: lo.FromPtr(c.Sunrise),
|
||||||
|
Dhuhr: lo.FromPtr(c.Dhuhr),
|
||||||
|
Asr: lo.FromPtr(c.Asr),
|
||||||
|
Sunset: lo.FromPtr(c.Sunset),
|
||||||
|
Maghrib: lo.FromPtr(c.Maghrib),
|
||||||
|
Isha: lo.FromPtr(c.Isha),
|
||||||
|
}) {
|
||||||
|
return
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// julianDay converts a calendar date to a Julian Day Number.
|
||||||
|
//
|
||||||
|
// JDN is a continuous day count from Jan 1, 4713 BC, used in astronomy to
|
||||||
|
// avoid calendar-system ambiguities. January/February are treated as months
|
||||||
|
// 13/14 of the prior year. The Gregorian correction b = 2 - a + a/4 accounts
|
||||||
|
// for century-year leap-day rules introduced in 1582.
|
||||||
|
func julianDay(d time.Time) float64 {
|
||||||
|
y, m, day := d.Date()
|
||||||
|
if m <= 2 {
|
||||||
|
y--
|
||||||
|
m += 12
|
||||||
|
}
|
||||||
|
a := y / 100
|
||||||
|
b := 2 - a + a/4
|
||||||
|
return math.Floor(365.25*float64(y+4716)) +
|
||||||
|
math.Floor(30.6001*float64(m+1)) +
|
||||||
|
float64(day+b) - 1524.5
|
||||||
|
}
|
||||||
|
|
||||||
|
// sunParams computes solar declination and equation of time via Meeus Ch.25
|
||||||
|
// (~0.0003° accuracy - ~30x better than the USNO 6-term approximation).
|
||||||
|
//
|
||||||
|
// Returns:
|
||||||
|
// - delta: solar declination in degrees (Sun's angular distance north/south
|
||||||
|
// of the celestial equator; drives seasonal day length and noon altitude).
|
||||||
|
// - eot: equation of time in minutes (difference between apparent solar time
|
||||||
|
// and mean solar time; caused by orbital eccentricity + axial tilt; +-16 min).
|
||||||
|
//
|
||||||
|
// Source: Jean Meeus, "Astronomical Algorithms" 2nd ed., Chapter 25.
|
||||||
|
//
|
||||||
|
// Variables (degrees unless noted):
|
||||||
|
//
|
||||||
|
// T - Julian centuries since J2000.0
|
||||||
|
// L0 - geometric mean longitude of the Sun (Meeus eq. 25.2)
|
||||||
|
// M - mean anomaly (Meeus eq. 25.3)
|
||||||
|
// e - orbital eccentricity (Meeus eq. 25.4)
|
||||||
|
// C - equation of centre: true - mean anomaly (Meeus eq. 25.4)
|
||||||
|
// omega - Moon's ascending node longitude, used for nutation (Meeus eq. 25.11)
|
||||||
|
// lam - apparent longitude: true lon + nutation - aberration (Meeus eq. 25.9)
|
||||||
|
// eps - true obliquity of the ecliptic incl. nutation (Meeus eq. 22.2 / 25.8)
|
||||||
|
// EoT - Spencer/Meeus y-series (Meeus p.185), accurate to ~0.5 s
|
||||||
|
func sunParams(jd float64) (delta, eot float64) {
|
||||||
|
T := (jd - j2000) / 36525.0
|
||||||
|
|
||||||
|
L0 := math.Mod(280.46646+T*(36000.76983+T*0.0003032), 360)
|
||||||
|
M := math.Mod(357.52911+T*(35999.05029-T*0.0001537), 360)
|
||||||
|
Mr := rad(M)
|
||||||
|
e := 0.016708634 - T*(0.000042037+T*0.0000001267)
|
||||||
|
|
||||||
|
// Equation of centre: corrects uniform circular -> true elliptical motion.
|
||||||
|
C := (1.914602-T*(0.004817+T*0.000014))*math.Sin(Mr) +
|
||||||
|
(0.019993-T*0.000101)*math.Sin(2*Mr) +
|
||||||
|
0.000289*math.Sin(3*Mr)
|
||||||
|
|
||||||
|
omega := 125.04 - 1934.136*T // Moon's ascending node: drives nutation
|
||||||
|
|
||||||
|
// Apparent longitude: add nutation, subtract aberration (-0.00569°).
|
||||||
|
lam := rad(L0 + C - 0.00569 - 0.00478*math.Sin(rad(omega)))
|
||||||
|
|
||||||
|
// True obliquity: Laskar (1986) mean obliquity + nutation in obliquity.
|
||||||
|
eps0 := 84381.448 - T*(46.8150+T*(0.00059-T*0.001813)) // arcseconds
|
||||||
|
eps := rad(eps0/3600 + 0.00256*math.Cos(rad(omega)))
|
||||||
|
|
||||||
|
delta = deg(math.Asin(math.Sin(eps) * math.Sin(lam)))
|
||||||
|
|
||||||
|
// EoT via y-series (Spencer 1971 / Meeus p.185).
|
||||||
|
// y = tan²(e/2); multiply degrees by 4 to get minutes (1° = 4 min).
|
||||||
|
y, L0r := math.Pow(math.Tan(eps/2), 2), rad(L0)
|
||||||
|
eot = deg(y*math.Sin(2*L0r)-
|
||||||
|
2*e*math.Sin(Mr)+
|
||||||
|
4*e*y*math.Sin(Mr)*math.Cos(2*L0r)-
|
||||||
|
0.5*y*y*math.Sin(4*L0r)-
|
||||||
|
1.25*e*e*math.Sin(2*Mr)) * 4
|
||||||
|
|
||||||
|
return delta, eot
|
||||||
|
}
|
||||||
|
|
||||||
|
// hourAngle solves for the hour angle H (hours) at which the Sun reaches
|
||||||
|
// the given altitude, using the spherical law of cosines:
|
||||||
|
//
|
||||||
|
// cos H = (sin a - sin phi*sin delta) / (cos phi*cos delta)
|
||||||
|
//
|
||||||
|
// Returns (0, false) when |cos H| > 1 - the Sun never reaches that altitude
|
||||||
|
// (midnight sun / polar night). Diyanet resolves these via the Takdir method.
|
||||||
|
func hourAngle(altDeg, lat, delta float64) (float64, bool) {
|
||||||
|
cosH := (math.Sin(rad(altDeg)) - math.Sin(rad(lat))*math.Sin(rad(delta))) /
|
||||||
|
(math.Cos(rad(lat)) * math.Cos(rad(delta)))
|
||||||
|
if math.Abs(cosH) > 1 {
|
||||||
|
return 0, false
|
||||||
|
}
|
||||||
|
return deg(math.Acos(cosH)) / degPerHour, true
|
||||||
|
}
|
||||||
|
|
||||||
|
// asrAltitude returns the solar altitude at which Asr-i Avval begins.
|
||||||
|
//
|
||||||
|
// Diyanet (majority school): Asr starts when shadow length = object height +
|
||||||
|
// its shortest noon shadow (fey-i zeval). Shadow factor = 1 (Hanafi uses 2).
|
||||||
|
//
|
||||||
|
// cot a = 1 + tan|phi - delta| -> a = atan(1 / (1 + tan|phi - delta|))
|
||||||
|
func asrAltitude(lat, delta float64) float64 {
|
||||||
|
return deg(math.Atan(1 / (1 + math.Tan(rad(math.Abs(lat-delta))))))
|
||||||
|
}
|
||||||
|
|
||||||
|
type computedTimes struct {
|
||||||
|
Imsak, Sunrise, Dhuhr, Asr, Sunset, Maghrib, Isha *time.Time
|
||||||
|
}
|
||||||
|
|
||||||
|
// prayerTimes computes all Diyanet prayer times for the given date, returned
|
||||||
|
// as UTC-aware values. Solar noon is the central reference:
|
||||||
|
//
|
||||||
|
// T_noon(UTC) = 12 - lambda/15 - EoT/60
|
||||||
|
//
|
||||||
|
// Morning times (Imsak, Sunrise) = noon - H + temkin
|
||||||
|
// Afternoon/evening times = noon + H + temkin
|
||||||
|
func prayerTimes(lat, lon float64, d time.Time) computedTimes {
|
||||||
|
delta, eot := sunParams(julianDay(d))
|
||||||
|
noon := 12 - lon/degPerHour - eot/60
|
||||||
|
|
||||||
|
// offset converts a noon-relative hour angle to a UTC *time.Time,
|
||||||
|
// applying the given temkin (minutes). Returns nil for polar night/day.
|
||||||
|
offset := func(h float64, ok bool, sign int, tk float64) *time.Time {
|
||||||
|
if !ok {
|
||||||
|
return nil
|
||||||
|
}
|
||||||
|
t := utcTime(d, noon+float64(sign)*h+tk/60)
|
||||||
|
return &t
|
||||||
|
}
|
||||||
|
|
||||||
|
hSun, okSun := hourAngle(sunAngle, lat, delta)
|
||||||
|
hAsr, okAsr := hourAngle(asrAltitude(lat, delta), lat, delta)
|
||||||
|
hImsak, okImsak := hourAngle(imsakAngle, lat, delta)
|
||||||
|
hIsha, okIsha := hourAngle(ishaAngle, lat, delta)
|
||||||
|
|
||||||
|
tDhuhr := utcTime(d, noon+temkin.Dhuhr/60)
|
||||||
|
|
||||||
|
return computedTimes{
|
||||||
|
Imsak: offset(hImsak, okImsak, -1, temkin.Imsak),
|
||||||
|
Sunrise: offset(hSun, okSun, -1, temkin.Sunrise),
|
||||||
|
Dhuhr: &tDhuhr,
|
||||||
|
Asr: offset(hAsr, okAsr, +1, temkin.Asr),
|
||||||
|
Sunset: offset(hSun, okSun, +1, 0), // geometric sunset, no temkin
|
||||||
|
Maghrib: offset(hSun, okSun, +1, temkin.Maghrib),
|
||||||
|
Isha: offset(hIsha, okIsha, +1, temkin.Isha),
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// utcTime converts decimal hours (e.g. 10.5 = 10:30) to a UTC time.Time.
|
||||||
|
func utcTime(d time.Time, hours float64) time.Time {
|
||||||
|
base := time.Date(d.Year(), d.Month(), d.Day(), 0, 0, 0, 0, time.UTC)
|
||||||
|
return base.Add(time.Duration(hours * float64(time.Hour)))
|
||||||
|
}
|
||||||
|
|
||||||
|
func rad(d float64) float64 { return d * math.Pi / 180 }
|
||||||
|
func deg(r float64) float64 { return r * 180 / math.Pi }
|
||||||
Loading…
Reference in New Issue