From 9cfa7f2dcf1707d77349dbfd4628a3a208a00149 Mon Sep 17 00:00:00 2001 From: Abdussamet Kocak Date: Sun, 22 Feb 2026 00:46:28 +0300 Subject: [PATCH] fix: Use better algo (Meeus Ch.25) for estimating sun position when calculating prayer times --- go.mod | 1 + go.sum | 2 + pkg/diyanetcalc/provider.go | 363 ++++++++++++++---------------------- 3 files changed, 144 insertions(+), 222 deletions(-) diff --git a/go.mod b/go.mod index 4ca52ae..cc7384c 100644 --- a/go.mod +++ b/go.mod @@ -31,6 +31,7 @@ require ( github.com/quic-go/quic-go v0.57.1 // indirect github.com/refraction-networking/utls v1.8.2 // indirect github.com/remyoudompheng/bigfft v0.0.0-20230129092748-24d4a6f8daec // indirect + github.com/samber/lo v1.52.0 // indirect github.com/tinylib/msgp v1.6.3 // indirect github.com/valyala/bytebufferpool v1.0.0 // indirect github.com/valyala/fasthttp v1.69.0 // indirect diff --git a/go.sum b/go.sum index 88a937f..d71a2fd 100644 --- a/go.sum +++ b/go.sum @@ -73,6 +73,8 @@ github.com/rogpeppe/go-internal v1.10.0/go.mod h1:UQnix2H7Ngw/k4C5ijL5+65zddjncj github.com/rs/xid v1.4.0/go.mod h1:trrq9SKmegXys3aeAKXMUTdJsYXVwGY3RLcfgqegfbg= github.com/rs/zerolog v1.29.0 h1:Zes4hju04hjbvkVkOhdl2HpZa+0PmVwigmo8XoORE5w= github.com/rs/zerolog v1.29.0/go.mod h1:NILgTygv/Uej1ra5XxGf82ZFSLk58MFGAUS2o6usyD0= +github.com/samber/lo v1.52.0 h1:Rvi+3BFHES3A8meP33VPAxiBZX/Aws5RxrschYGjomw= +github.com/samber/lo v1.52.0/go.mod h1:4+MXEGsJzbKGaUEQFKBq2xtfuznW9oz/WrgyzMzRoM0= github.com/shamaton/msgpack/v3 v3.0.0 h1:xl40uxWkSpwBCSTvS5wyXvJRsC6AcVcYeox9PspKiZg= github.com/shamaton/msgpack/v3 v3.0.0/go.mod h1:DcQG8jrdrQCIxr3HlMYkiXdMhK+KfN2CitkyzsQV4uc= github.com/stretchr/objx v0.1.0/go.mod h1:HFkY916IF+rwdDfMAkV7OtwuqBVzrE8GR6GFx+wExME= diff --git a/pkg/diyanetcalc/provider.go b/pkg/diyanetcalc/provider.go index 70918b2..2cff32d 100644 --- a/pkg/diyanetcalc/provider.go +++ b/pkg/diyanetcalc/provider.go @@ -9,77 +9,49 @@ import ( "prayertimes/pkg/hijricalendar" "prayertimes/pkg/prayer" + + "github.com/samber/lo" ) const ( daysToGenerate = 30 degPerHour = 15.0 + j2000 = 2451545.0 // Julian Date of J2000.0 epoch (Jan 1.5, 2000) ) var errNotSupported = errors.New("not supported in calculation provider") -// Diyanet prayer times calculator. -// -// Implements the Turkish Presidency of Religious Affairs (Diyanet) methodology, -// standardized in the 1983 reform. Prayer times are indexed to the Sun's apparent -// altitude angle at the observer's location, solved via spherical trigonometry. - -// --------------------------------------------------------------------------- // Diyanet angular criteria (post-1983 reform) -// --------------------------------------------------------------------------- - -// Imsak (Fajr): Sun is 18deg below the horizon = start of astronomical twilight. -// Pre-1983 used -19deg plus a temkin buffer; now -18deg with zero temkin. +// +// 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 -// Isha (Yatsi): Sun is 17deg below the horizon = shafaq al-ahmar (red twilight) -// has fully disappeared from the western sky. +// ishaAngle: Sun 17° below horizon = shafaq al-ahmar (red twilight) gone. const ishaAngle = -17.0 -// Sunrise / Sunset (Tulu / Gurup): geometric horizon alone is insufficient. -// Two physical corrections are combined into a single -0.833deg value (50 arcmin): -// - Atmospheric refraction: ~0.567deg — air bends sunlight over the horizon -// - Solar semi-diameter: ~0.267deg — Sun is "up" when its upper limb clears +// 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) // -// Temkin ensures a single published time remains valid across the full -// geographical extent of a city (highest peak to lowest valley, east to west). -// Pre-1983 values were often 10–20 min; the reform moderated them. -// -// Imsak: 0 min — no buffer; avoids starting Fajr too early / fast too late -// Sunrise: -7 min — subtracted, ensuring the Sun has fully cleared the horizon -// Dhuhr: +5 min — added, ensuring the Sun has clearly begun its descent -// Asr: +4 min — accounts for local elevation and horizon obstacles -// Maghrib:+7 min — ensures the Sun has completely set before breaking fast -// Isha: 0 min — no buffer needed at this twilight stage +// 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. // -// --------------------------------------------------------------------------- -type temkinMinutes struct { - Imsak float64 - Sunrise float64 - Dhuhr float64 - Asr float64 - Maghrib float64 - Isha float64 -} - -var temkin = temkinMinutes{ - Imsak: 0, - Sunrise: -7, - Dhuhr: 5, - Asr: 4, - Maghrib: 7, - Isha: 0, +// 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 Provider struct{} -func New() Provider { - return Provider{} -} +func New() Provider { return Provider{} } func (Provider) SearchLocations(_ context.Context, _ string) ([]prayer.Location, error) { return nil, fmt.Errorf("failed to search locations: %w", errNotSupported) @@ -90,241 +62,188 @@ func (Provider) Get(_ context.Context, _ string) (prayer.TimesResult, error) { } func (Provider) GetByCoords(_ context.Context, coords prayer.Coordinates) (prayer.TimesResult, error) { - offset := estimateUTCOffsetHours(coords.Longitude) todayUTC := time.Now().UTC().Truncate(24 * time.Hour) - results := make([]prayer.Times, 0, daysToGenerate) - for i := 0; i < daysToGenerate; i++ { - day := todayUTC.AddDate(0, 0, i) - calculated := prayerTimes(coords.Latitude, coords.Longitude, day) + days := lo.Map(lo.Range(daysToGenerate), func(i, _ int) time.Time { + return todayUTC.AddDate(0, 0, i) + }) - results = append(results, prayer.Times{ + times := lo.Map(days, func(day time.Time, _ int) prayer.Times { + c := prayerTimes(coords.Latitude, coords.Longitude, day) + return prayer.Times{ Date: day, DateHijri: hijricalendar.ToISODate(day), - Fajr: derefOrZero(calculated.Imsak), - Sunrise: derefOrZero(calculated.Sunrise), - Dhuhr: derefOrZero(calculated.Dhuhr), - Asr: derefOrZero(calculated.Asr), - Sunset: derefOrZero(calculated.Sunset), - Maghrib: derefOrZero(calculated.Maghrib), - Isha: derefOrZero(calculated.Isha), - }) - } + 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 prayer.TimesResult{ Location: prayer.Location{ Latitude: coords.Latitude, Longitude: coords.Longitude, - Timezone: formatUTCOffsetTimezone(offset), + Timezone: fmt.Sprintf("UTC%+d", int(math.Round(coords.Longitude/degPerHour))), }, - Times: results, + Times: times, }, nil } -// Convert a calendar date to a Julian Day Number. +// 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. The -1524.5 offset shifts the epoch to -// noon UT, Jan 1, 4713 BC (the standard astronomical Julian Date epoch). -// The Gregorian calendar correction term b = 2 - a + a//4 accounts for -// the century-year leap-day rules introduced in 1582. +// 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 { - // January/February treated as months 13/14 of the prior year y-- m += 12 } a := y / 100 - b := 2 - a + a/4 // Gregorian correction - - return math.Floor(365.25*float64(y+4716)) + math.Floor(30.6001*float64(m+1)) + float64(day+b) - 1524.5 + 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 } -// Compute solar declination and equation of time for a given Julian Day. +// sunParams computes solar declination and equation of time via Meeus Ch.25 +// (~0.0003° accuracy — ~30× better than the USNO 6-term approximation). // // Returns: -// -// delta — solar declination in degrees: the Sun's angular distance -// north (+) or south (-) of the celestial equator. Controls -// seasonal day length and the Sun's maximum altitude. -// eot — equation of time in minutes: difference between apparent solar -// time (sundial) and mean solar time (clock). Caused by Earth's -// elliptical orbit and axial tilt; ranges roughly ±16 min. -// -// Algorithm uses low-precision USNO solar coordinates (~0.01deg accuracy): -// -// d — days since J2000.0 epoch (Jan 1.5, 2000 = JD 2451545.0) -// g — mean anomaly: Sun's angular position in its elliptical orbit -// q — mean longitude of the Sun -// L — ecliptic longitude: corrected for orbital eccentricity via the -// equation of centre (1.915deg*sin g + 0.020deg*sin 2g) -// e — obliquity of the ecliptic: Earth's axial tilt (~23.44deg, slowly -// decreasing at 0.00000036deg/day) -func sunParams(jd float64) (delta float64, eot float64) { - d := jd - 2451545.0 // days since J2000.0 - - g := toRadians(357.529 + 0.98560028*d) // mean anomaly - q := toRadians(280.459 + 0.98564736*d) // mean longitude - - // Ecliptic longitude: equation of centre adds up to ~1.9deg correction for - // the difference between uniform circular and actual elliptical motion. - L := q + toRadians(1.915*math.Sin(g)+0.020*math.Sin(2*g)) - e := toRadians(23.439 - 0.00000036*d) // obliquity of ecliptic - - // Declination: project ecliptic longitude onto the celestial equator. - delta = toDegrees(math.Asin(math.Sin(e) * math.Sin(L))) - - // Right ascension in hours (atan2 handles all four quadrants correctly). - ra := toDegrees(math.Atan2(math.Cos(e)*math.Sin(L), math.Cos(L))) / degPerHour - - // EoT = mean sun hour angle minus apparent sun RA, normalized to ±30 min. - // round() removes the large integer offset (q accumulates many full rotations) - // before converting to minutes; without it the raw difference is ~600 hours. - diff := toDegrees(q)/degPerHour - ra - eot = (diff - math.Round(diff)) * 60 +// - 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²(ε/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 } -// Solve for the hour angle H (hours) at which the Sun reaches a given altitude. -// -// Derived from the spherical law of cosines for the astronomical triangle: -// -// sin(a) = sin(phi)*sin(delta) + cos(phi)*cos(delta)*cos(H) +// hourAngle solves for the hour angle H (hours) at which the Sun reaches +// the given altitude, using the spherical law of cosines: // -// Rearranged: +// cos H = (sin a − sin φ·sin δ) / (cos φ·cos δ) // -// cos(H) = (sin(a) − sin(phi)*sin(delta)) / (cos(phi)*cos(delta)) -// -// H is converted from degrees to hours by dividing by 15 (360deg/24h = 15deg/h). -// Returns false when |cos H| > 1, i.e. the Sun never reaches that altitude -// (midnight sun or polar night) — Diyanet handles these with the Takdir method. -func hourAngle(altitudeDeg, lat, delta float64) (hours float64, ok bool) { - cosH := (math.Sin(toRadians(altitudeDeg)) - math.Sin(toRadians(lat))*math.Sin(toRadians(delta))) / - (math.Cos(toRadians(lat)) * math.Cos(toRadians(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 toDegrees(math.Acos(cosH)) / degPerHour, true + return deg(math.Acos(cosH)) / degPerHour, true } -// Compute the solar altitude at which Asr begins (Asr-i Avval / First Asr). +// asrAltitude returns the solar altitude at which Asr-i Avval begins. // -// Diyanet follows the majority-school definition: Asr starts when an object's -// shadow length equals the object's height plus its shortest noon shadow (fey-i zeval). -// The shadow factor is 1 (Asr-i Avval; Hanafi uses 2 for Asr-i Sani). +// Diyanet (majority school): Asr starts when shadow length = object height + +// its shortest noon shadow (fey-i zeval). Shadow factor = 1 (Hanafi uses 2). // -// The required altitude satisfies: cot(a) = 1 + tan(|phi − delta|) -// which is: a = atan(1 / (1 + tan(|phi − delta|))) -// where |phi − delta| is the Sun's angular distance from the zenith at solar noon. +// cot a = 1 + tan|φ − δ| → a = atan(1 / (1 + tan|φ − δ|)) func asrAltitude(lat, delta float64) float64 { - return toDegrees(math.Atan(1.0 / (1.0 + math.Tan(toRadians(math.Abs(lat-delta)))))) -} - -// Convert a decimal hour value (e.g. 10.5 = 10:30) to a UTC-aware datetime. -func decimalHoursToUTC(hours float64, d time.Time) time.Time { - dayUTC := time.Date(d.Year(), d.Month(), d.Day(), 0, 0, 0, 0, time.UTC) - return dayUTC.Add(time.Duration(hours * float64(time.Hour))) + return deg(math.Atan(1 / (1 + math.Tan(rad(math.Abs(lat-delta)))))) } type computedTimes struct { - Imsak *time.Time - Sunrise *time.Time - Dhuhr *time.Time - Asr *time.Time - Sunset *time.Time - Maghrib *time.Time - Isha *time.Time + Imsak, Sunrise, Dhuhr, Asr, Sunset, Maghrib, Isha *time.Time } -// Compute Diyanet prayer times, returning UTC-aware datetimes. +// prayerTimes computes all Diyanet prayer times for the given date, returned +// as UTC-aware values. Solar noon is the central reference: // -// Solar noon (Dhuhr) is the central reference. All other times are offsets: +// T_noon(UTC) = 12 − λ/15 − EoT/60 // -// Morning times (Imsak, Sunrise): noon − H + temkin -// Afternoon/evening times (Asr, Maghrib, Isha): noon + H + temkin -// -// Internally computes solar noon at UTC (tz=0), so results are in UTC. -// -// Solar noon formula: T_noon = 12 + TZ − lambda/15 − EoT/60 -// -// lambda/15 converts longitude to hours (15deg/h) -// EoT corrects the gap between mean solar time and apparent solar time +// Morning times (Imsak, Sunrise) = noon − H + temkin +// Afternoon/evening times = noon + H + temkin func prayerTimes(lat, lon float64, d time.Time) computedTimes { - jd := julianDay(d) - delta, eot := sunParams(jd) + delta, eot := sunParams(julianDay(d)) + noon := 12 - lon/degPerHour - eot/60 - // UTC solar noon: tz=0, so T_noon = 12 − lambda/15 − EoT/60 - noonUTC := 12 - lon/degPerHour - eot/60.0 - - compute := func(base, angle float64, sign int, temkinMinutes float64) *time.Time { - h, ok := hourAngle(angle, lat, delta) + // 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 := decimalHoursToUTC(base+float64(sign)*h+temkinMinutes/60.0, d) + t := utcTime(d, noon+float64(sign)*h+tk/60) return &t } - hAsr, hasAsr := hourAngle(asrAltitude(lat, delta), lat, delta) - var asr *time.Time - if hasAsr { - t := decimalHoursToUTC(noonUTC+hAsr+temkin.Asr/60.0, d) - asr = &t - } - - // Sunset for output is the geometric sunset without temkin offset. - geometricSunset := func(base, angle float64) *time.Time { - h, ok := hourAngle(angle, lat, delta) - if !ok { - return nil - } - t := decimalHoursToUTC(base+h, d) - 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 := decimalHoursToUTC(noonUTC+temkin.Dhuhr/60.0, d) + tDhuhr := utcTime(d, noon+temkin.Dhuhr/60) return computedTimes{ - Imsak: compute(noonUTC, imsakAngle, -1, temkin.Imsak), - Sunrise: compute(noonUTC, sunAngle, -1, temkin.Sunrise), + Imsak: offset(hImsak, okImsak, -1, temkin.Imsak), + Sunrise: offset(hSun, okSun, -1, temkin.Sunrise), Dhuhr: &tDhuhr, - Asr: asr, - Sunset: geometricSunset(noonUTC, sunAngle), - Maghrib: compute(noonUTC, sunAngle, +1, temkin.Maghrib), - Isha: compute(noonUTC, ishaAngle, +1, temkin.Isha), + 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), } } -func derefOrZero(dt *time.Time) time.Time { - if dt == nil { - return time.Time{} - } - return dt.UTC() -} - -func formatUTCOffsetTimezone(offset float64) string { - return fmt.Sprintf("UTC%+d", int(offset)) -} - -func estimateUTCOffsetHours(longitude float64) float64 { - offset := math.Round(longitude / degPerHour) - if offset < -12 { - return -12 - } - if offset > 14 { - return 14 - } - return offset +// 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 toRadians(deg float64) float64 { - return deg * math.Pi / 180.0 -} - -func toDegrees(rad float64) float64 { - return rad * 180.0 / math.Pi -} +func rad(d float64) float64 { return d * math.Pi / 180 } +func deg(r float64) float64 { return r * 180 / math.Pi }