109 lines
2.5 KiB
Go
109 lines
2.5 KiB
Go
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
|
|
}
|