1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
|
// calculator.go - Weather calculations
package main
import (
"math"
"time"
)
// CalculateFeelsLike calculates apparent temperature
// windSpeed is expected in km/h
func CalculateFeelsLike(temp, humidity, windSpeed float64) float64 {
if temp <= 10 && windSpeed > 4.8 {
// Wind chill formula (Environment Canada, wind in km/h)
vPow := math.Pow(windSpeed, 0.16)
return 13.12 + 0.6215*temp - 11.37*vPow + 0.3965*temp*vPow
} else if temp >= 27 {
// Simplified heat index for humid heat (Steadman's approximation)
e := humidity / 100 * 6.105 * math.Exp(17.27*temp/(237.7+temp))
return temp + 0.33*e - 0.70*windSpeed - 4.00
}
return temp
}
// SunCalculator calculates sunrise and sunset times
type SunCalculator struct {
latitude float64
longitude float64
utcOffset float64
}
// NewSunCalculator creates a new sun calculator
func NewSunCalculator(latitude, longitude, utcOffset float64) *SunCalculator {
return &SunCalculator{
latitude: latitude,
longitude: longitude,
utcOffset: utcOffset,
}
}
// Calculate computes sunrise and sunset for a given date
func (sc *SunCalculator) Calculate(date time.Time) (sunrise, sunset time.Time, err error) {
julianDay := sc.dateToJulianDay(date)
julianCentury := (julianDay - 2451545) / 36525.0
geomMeanLongSun := math.Mod(280.46646+julianCentury*(36000.76983+julianCentury*0.0003032), 360)
geomMeanAnomSun := 357.52911 + julianCentury*(35999.05029-0.0001537*julianCentury)
eccentEarthOrbit := 0.016708634 - julianCentury*(0.000042037+0.0000001267*julianCentury)
sunEqOfCtr := math.Sin(degToRad(geomMeanAnomSun))*(1.914602-julianCentury*(0.004817+0.000014*julianCentury)) +
math.Sin(degToRad(2*geomMeanAnomSun))*(0.019993-0.000101*julianCentury) +
math.Sin(degToRad(3*geomMeanAnomSun))*0.000289
sunTrueLong := geomMeanLongSun + sunEqOfCtr
sunAppLong := sunTrueLong - 0.00569 - 0.00478*math.Sin(degToRad(125.04-1934.136*julianCentury))
meanObliqEcliptic := 23 + (26+(21.448-julianCentury*(46.815+julianCentury*(0.00059-julianCentury*0.001813)))/60)/60
obliqCorr := meanObliqEcliptic + 0.00256*math.Cos(degToRad(125.04-1934.136*julianCentury))
sunDeclin := radToDeg(math.Asin(math.Sin(degToRad(obliqCorr)) * math.Sin(degToRad(sunAppLong))))
varY := math.Tan(degToRad(obliqCorr/2)) * math.Tan(degToRad(obliqCorr/2))
eqOfTime := 4 * radToDeg(varY*math.Sin(2*degToRad(geomMeanLongSun))-
2*eccentEarthOrbit*math.Sin(degToRad(geomMeanAnomSun))+
4*eccentEarthOrbit*varY*math.Sin(degToRad(geomMeanAnomSun))*math.Cos(2*degToRad(geomMeanLongSun))-
0.5*varY*varY*math.Sin(4*degToRad(geomMeanLongSun))-
1.25*eccentEarthOrbit*eccentEarthOrbit*math.Sin(2*degToRad(geomMeanAnomSun)))
haSunrise := radToDeg(math.Acos(math.Cos(degToRad(90.833))/(math.Cos(degToRad(sc.latitude))*math.Cos(degToRad(sunDeclin))) -
math.Tan(degToRad(sc.latitude))*math.Tan(degToRad(sunDeclin))))
solarNoon := (720 - 4*sc.longitude - eqOfTime + sc.utcOffset*60) / 1440
sunriseTime := solarNoon - haSunrise*4/1440
sunsetTime := solarNoon + haSunrise*4/1440
return sc.julianDayToDate(julianDay + sunriseTime), sc.julianDayToDate(julianDay + sunsetTime), nil
}
func (sc *SunCalculator) dateToJulianDay(date time.Time) float64 {
year := float64(date.Year())
month := float64(date.Month())
day := float64(date.Day())
if month <= 2 {
year -= 1
month += 12
}
a := math.Floor(year / 100)
b := 2 - a + math.Floor(a/4)
return math.Floor(365.25*(year+4716)) + math.Floor(30.6001*(month+1)) + day + b - 1524.5
}
func (sc *SunCalculator) julianDayToDate(julianDay float64) time.Time {
julianDay += 0.5
z := math.Floor(julianDay)
f := julianDay - z
var a float64
if z >= 2299161 {
alpha := math.Floor((z - 1867216.25) / 36524.25)
a = z + 1 + alpha - math.Floor(alpha/4)
} else {
a = z
}
b := a + 1524
c := math.Floor((b - 122.1) / 365.25)
d := math.Floor(365.25 * c)
e := math.Floor((b - d) / 30.6001)
day := b - d - math.Floor(30.6001*e) + f
month := e - 1
if month > 12 {
month -= 12
}
year := c - 4715
if month > 2 {
year -= 1
}
hour := int(math.Floor((day - math.Floor(day)) * 24))
minute := int(math.Floor(((day-math.Floor(day))*24 - float64(hour)) * 60))
second := int(math.Floor((((day-math.Floor(day))*24-float64(hour))*60 - float64(minute)) * 60))
return time.Date(int(year), time.Month(month), int(math.Floor(day)), hour, minute, second, 0, time.UTC).
Add(time.Hour * time.Duration(sc.utcOffset))
}
func degToRad(deg float64) float64 {
return deg * math.Pi / 180
}
func radToDeg(rad float64) float64 {
return rad * 180 / math.Pi
}
|
