diff --git a/pkg/diyanetcalc/provider.go b/pkg/diyanetcalc/provider.go index 2cff32d..b8bb60f 100644 --- a/pkg/diyanetcalc/provider.go +++ b/pkg/diyanetcalc/provider.go @@ -7,51 +7,16 @@ import ( "math" "time" - "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 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 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) @@ -62,26 +27,30 @@ func (Provider) Get(_ context.Context, _ string) (prayer.TimesResult, error) { } func (Provider) GetByCoords(_ context.Context, coords prayer.Coordinates) (prayer.TimesResult, error) { - todayUTC := time.Now().UTC().Truncate(24 * time.Hour) - - days := lo.Map(lo.Range(daysToGenerate), func(i, _ int) time.Time { - return todayUTC.AddDate(0, 0, i) + seq := CalculatePrayerTimes(CalculateParams{ + Latitude: coords.Latitude, + Longitude: coords.Longitude, + StartingDay: time.Now().UTC().Truncate(24 * time.Hour), }) - 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: 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), + times := make([]prayer.Times, 0, daysToGenerate) + for item := range seq { + times = append(times, prayer.Times{ + Date: item.Date, + DateHijri: item.DateHijri, + Fajr: item.Fajr, + Sunrise: item.Sunrise, + Dhuhr: item.Dhuhr, + Asr: item.Asr, + Sunset: item.Sunset, + Maghrib: item.Maghrib, + Isha: item.Isha, + }) + + if len(times) == daysToGenerate { + break } - }) + } return prayer.TimesResult{ Location: prayer.Location{ @@ -92,158 +61,3 @@ func (Provider) GetByCoords(_ context.Context, coords prayer.Coordinates) (praye Times: times, }, nil } - -// 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 — ~30× 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²(ε/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 φ·sin δ) / (cos φ·cos δ) -// -// 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|φ − δ| → a = atan(1 / (1 + tan|φ − δ|)) -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 − λ/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 } diff --git a/pkg/diyanetcalc/times.go b/pkg/diyanetcalc/times.go new file mode 100644 index 0000000..3c36f72 --- /dev/null +++ b/pkg/diyanetcalc/times.go @@ -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 }