summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rwxr-xr-xsrc/__main__.py1
-rw-r--r--src/solar.py159
2 files changed, 160 insertions, 0 deletions
diff --git a/src/__main__.py b/src/__main__.py
index 1170aa8..400c675 100755
--- a/src/__main__.py
+++ b/src/__main__.py
@@ -20,6 +20,7 @@ import time
import signal
import datetime
+from solar import *
from curve import *
from colour import *
from monitor import *
diff --git a/src/solar.py b/src/solar.py
new file mode 100644
index 0000000..a0620bc
--- /dev/null
+++ b/src/solar.py
@@ -0,0 +1,159 @@
+#!/usr/bin/env python3
+
+# Copyright © 2014 Mattias Andrée (maandree@member.fsf.org)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Affero General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+import math
+import time
+
+
+def julian_day_to_epoch(t):
+ return (jd - 2440587.5) * 86400.0
+
+def epoch_to_julian_day(t):
+ return t / 86400.0 + 2440587.5
+
+def julian_day_to_julian_centuries(t):
+ return (t - 2451545.0) / 36525.0
+
+def julian_centuries_to_julian_day(t):
+ return t * 36525.0 + 2451545.0
+
+def epoch_to_julian_centuries(t):
+ return julian_day_to_julian_centuries(epoch_to_julian_day(t))
+
+def julian_centuries_to_epoch(t):
+ return julian_day_to_epoch(julian_centuries_to_julian_day(t))
+
+def epoch():
+ return time.time()
+
+def julian_day():
+ return epoch_to_julian_day(epoch())
+
+def julian_centuries():
+ return epoch_to_julian_centuries(epoch())
+
+def radians(deg):
+ return deg * math.pi / 180
+
+def degrees(rad):
+ return rad * 180 / math.pi
+
+def sun_geometric_mean_longitude(t):
+ return radians((0.0003032 * t ** 2 + 36000.76983 * t + 280.46646) % 360)
+
+def sun_geometric_mean_anomaly(t):
+ return radians(-0.0001537 * t ** 2 + 35999.05029 * t + 357.52911)
+
+def earth_orbit_eccentricity(t):
+ return -0.0000001267 * t ** 2 - 0.000042037 * t + 0.016708634
+
+def sun_equation_of_centre(t):
+ a = sun_geometric_mean_anomaly(t)
+ rc = math.sin(1 * a) * (-0.000014 * t ** 2 - 0.004817 * t + 1.914602)
+ rc += math.sin(2 * a) * (-0.000101 * t + 0.019993)
+ rc += math.sin(3 * a) * 0.000289
+ return radians(rc)
+
+def sun_real_longitude(t):
+ rc = sun_geometric_mean_longitude(t)
+ return rc + sun_equation_of_centre(t)
+
+def sun_apparent_longitude(t):
+ rc = degrees(sun_real_longitude(t)) - 0.00569
+ rc -= 0.00478 * math.sin(radians(-1934.136 * t + 125.04))
+ return radians(rc)
+
+def mean_ecliptic_obliquity(t):
+ rc = 0.001813 * t ** 3 - 0.00059 * t ** 2 - 46.815 * t + 21.448
+ rc = 26 + rc / 60
+ rc = 23 + rc / 60
+ return radians(rc)
+
+def corrected_mean_ecliptic_obliquity(t):
+ rc = -1934.136 * t + 125.04
+ rc = 0.00256 * math.cos(radians(rc))
+ rc += degrees(mean_ecliptic_obliquity(t))
+ return radians(rc)
+
+def solar_declination(t):
+ rc = math.sin(corrected_mean_ecliptic_obliquity(t))
+ rc *= math.sin(sun_apparent_longitude(t))
+ return math.asin(rc)
+
+def equation_of_time(t):
+ l = sun_geometric_mean_longitude(t)
+ e = earth_orbit_eccentricity(t)
+ m = sun_geometric_mean_anomaly(t)
+ y = corrected_mean_ecliptic_obliquity(t)
+ y = math.tan(y / 2) ** 2
+ rc = y * math.sin(2 * l)
+ rc += (4 * y * math.cos(2 * l) - 2) * e * math.sin(m)
+ rc -= 0.5 * y ** 2 * math.sin(4 * l)
+ rc -= 1.25 * e ** 2 * math.sin(2 * m)
+ return 4 * degrees(rc)
+
+def hour_angle_from_elevation(latitude, declinaton, elevation):
+ if elevation == 0:
+ return 0
+ rc = math.cos(abs(elevation))
+ rc -= math.sin(radians(latitude)) * math.sin(declinaton)
+ rc /= math.cos(radians(latitude)) * math.cos(declinaton)
+ rc = math.acos(rc)
+ return -rc if (rc < 0) == (elevation < 0) else rc;
+
+def elevation_from_hour_angle(latitude, declinaton, hour_angle):
+ rc = math.cos(radians(latitude))
+ rc *= math.cos(hour_angle) * math.cos(declinaton)
+ rc += math.sin(radians(latitude)) * math.sin(declinaton)
+ return math.asin(rc)
+
+def time_of_solar_noon(t, longitude):
+ t, rc = julian_centuries_to_julian_day(t), longitude
+ for (k, m) in ((-360, 0), (1440, -0.5)):
+ rc = julian_day_to_julian_centuries(t + m + rc / k)
+ rc = 720 - 4 * longitude - equation_of_time(rc)
+ return rc
+
+def time_of_solar_elevation(t, noon, latitude, longitude, elevation):
+ rc = noon
+ rc, et = solar_declination(rc), equation_of_time(rc)
+ rc = hour_angle_from_elevation(latitude, rc, elevation)
+ rc = 720 - 4 * (longitude + degrees(rc)) - et
+
+ rc = julian_day_to_julian_centuries(julian_centuries_to_julian_day(t) + rc / 1440)
+ rc, et = solar_declination(rc), equation_of_time(rc)
+ rc = hour_angle_from_elevation(latitude, rc, elevation)
+ rc = 720 - 4 * (longitude + degrees(rc)) - et
+ return rc
+
+def solar_elevation_from_time(t, latitude, longitude):
+ rc = julian_centuries_to_julian_day(t)
+ rc = (rc - float(int(rc + 0.5)) - 0.5) * 1440
+ rc = 720 - rc - equation_of_time(t)
+ rc = radians(rc / 4 - longitude)
+ return elevation_from_hour_angle(latitude, solar_declination(t), rc)
+
+def solar_elevation(latitude, longitude, t = None):
+ rc = julian_centuries() if t is None else t
+ rc = solar_elevation_from_time(rc, latitude, longitude)
+ return degrees(rc)
+
+def sun(latitude, longitude, t = None, low = -6.0, high = 3.0):
+ t = julian_centuries() if t is None else t
+ e = solar_elevation(latitude, longitude, t)
+ e = (e - low) / (high - low)
+ return min(max(0, e), 1)
+