prayertimes-go/pkg/diyanetcalc/provider.go

package diyanetcalc

import (
	"context"
	"errors"
	"fmt"
	"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 (Provider) SearchLocations(_ context.Context, _ string) ([]prayer.Location, error) {
	return nil, fmt.Errorf("failed to search locations: %w", errNotSupported)
}

func (Provider) Get(_ context.Context, _ string) (prayer.TimesResult, error) {
	return prayer.TimesResult{}, fmt.Errorf("failed to get prayer times by location id: %w", errNotSupported)
}

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)
	})

	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),
		}
	})

	return prayer.TimesResult{
		Location: prayer.Location{
			Latitude:  coords.Latitude,
			Longitude: coords.Longitude,
			Timezone:  fmt.Sprintf("UTC%+d", int(math.Round(coords.Longitude/degPerHour))),
		},
		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 }