package suncalc

import (
	"errors"
	"math"
	"time"
)

const (
	dayMs = 1000 * 60 * 60 * 24
	j1970 = 2440588
	j2000 = 2451545
	rad   = math.Pi / 180
	e     = rad * 23.4397
)

var ErrNoCivilTwilight = errors.New("civil twilight is not available for the selected date and location")

type CivilTwilight struct {
	Dawn time.Time
	Dusk time.Time
}

func GetCivilTwilight(date time.Time, lat, lng float64) (CivilTwilight, error) {
	if math.IsNaN(lat) || math.IsNaN(lng) || lat < -90 || lat > 90 || lng < -180 || lng > 180 {
		return CivilTwilight{}, errors.New("invalid coordinates")
	}

	d := toDays(date)
	lw := -lng * rad
	phi := lat * rad
	n := julianCycle(d, lw)
	ds := approxTransit(0, lw, n)
	M := solarMeanAnomaly(ds)
	L := eclipticLongitude(M)
	dec := declination(L, 0)
	noonJ := solarTransitJ(ds, M, L)

	duskJ, ok := getSetJ(-6*rad, lw, phi, dec, n, M, L)
	if !ok {
		return CivilTwilight{}, ErrNoCivilTwilight
	}
	dawnJ := 2*noonJ - duskJ

	return CivilTwilight{
		Dawn: fromJulian(dawnJ),
		Dusk: fromJulian(duskJ),
	}, nil
}

func toJulian(date time.Time) float64 {
	return float64(date.UTC().Unix())/86400 - 0.5 + j1970
}

func fromJulian(j float64) time.Time {
	return time.Unix(int64((j+0.5-float64(j1970))*86400), 0).UTC()
}

func toDays(date time.Time) float64 {
	return toJulian(date) - j2000
}

func rightAscension(l, b float64) float64 {
	return math.Atan2(math.Sin(l)*math.Cos(e)-math.Tan(b)*math.Sin(e), math.Cos(l))
}

func declination(l, b float64) float64 {
	return math.Asin(math.Sin(b)*math.Cos(e) + math.Cos(b)*math.Sin(e)*math.Sin(l))
}

func solarMeanAnomaly(d float64) float64 {
	return rad * (357.5291 + 0.98560028*d)
}

func eclipticLongitude(M float64) float64 {
	C := rad * (1.9148*math.Sin(M) + 0.02*math.Sin(2*M) + 0.0003*math.Sin(3*M))
	P := rad * 102.9372
	return M + C + P + math.Pi
}

func julianCycle(d, lw float64) float64 {
	return math.Round(d - 0.0009 - lw/(2*math.Pi))
}

func approxTransit(Ht, lw, n float64) float64 {
	return 0.0009 + (Ht+lw)/(2*math.Pi) + n
}

func solarTransitJ(ds, M, L float64) float64 {
	return j2000 + ds + 0.0053*math.Sin(M) - 0.0069*math.Sin(2*L)
}

func hourAngle(h, phi, d float64) float64 {
	cosH := (math.Sin(h) - math.Sin(phi)*math.Sin(d)) / (math.Cos(phi) * math.Cos(d))
	if cosH < -1 || cosH > 1 {
		return math.NaN()
	}
	return math.Acos(cosH)
}

func getSetJ(h, lw, phi, dec, n, M, L float64) (float64, bool) {
	w := hourAngle(h, phi, dec)
	if math.IsNaN(w) {
		return 0, false
	}
	a := approxTransit(w, lw, n)
	return solarTransitJ(a, M, L), true
}