1 files changed, 134 insertions, 0 deletions
diff --git a/calculator.go b/calculator.go
new file mode 100644
index 0000000..84ed8c8
--- /dev/null
+++ b/calculator.go
@@ -0,0 +1,134 @@
+// calculator.go - Weather calculations
+package main
+
+import (
+	"math"
+	"time"
+)
+
+// CalculateFeelsLike calculates apparent temperature
+func CalculateFeelsLike(temp, humidity, windSpeed float64) float64 {
+	if temp <= 10 && windSpeed > 1.34 {
+		// Wind chill formula for cold temperatures
+		return 13.12 + 0.6215*temp - 11.37*math.Pow(windSpeed, 0.16) + 0.3965*temp*math.Pow(windSpeed, 0.16)
+	} else if temp >= 27 {
+		// Simplified heat index for warm temperatures
+		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
+}