summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--src/solar.py72
1 files changed, 22 insertions, 50 deletions
diff --git a/src/solar.py b/src/solar.py
index 34233be..aceacfc 100644
--- a/src/solar.py
+++ b/src/solar.py
@@ -463,36 +463,35 @@ def is_winter(latitude, t = None):
# TODO document
-def future_past_equinox(delta, t = None):
- epsilon = 0.000001
+def solar_prediction(delta, requested, fun, epsilon = 0.000001, span = 0.01, t = None):
t = julian_centuries() if t is None else t
t1 = t2 = t
- d1 = d0 = solar_declination(t)
+ v1 = v0 = fun(t)
while True:
- if abs(t2 - t) > 0.01:
+ if abs(t2 - t) > span:
return None
t2 += delta
- d2 = solar_declination(t2)
- if (d1 <= 0 <= d2) or ((0 >= d1 >= d2) and (0 <= d0)):
+ v2 = fun(t2)
+ if (v1 <= requested <= v2) or ((requested >= v1 >= v2) and (requested <= v0)):
break
- if (d1 >= 0 >= d2) or ((0 <= d1 <= d2) and (0 >= d0)):
+ if (v1 >= requested >= v2) or ((requested <= v1 <= v2) and (requested >= v0)):
break
t1 = t2
- d2 = d1
+ v2 = v1
for _itr in range(1000):
tm = (t1 + t2) / 2
- d1 = solar_declination(t1)
- d2 = solar_declination(t2)
- dm = solar_declination(tm)
- if abs(d1 - d2) < epsilon:
- return tm if abs(dm) < epsilon else None
- if d1 < d2:
- if 0 < dm:
+ v1 = fun(t1)
+ v2 = fun(t2)
+ vm = fun(tm)
+ if abs(v1 - v2) < epsilon:
+ return tm if abs(vm) < epsilon else None
+ if v1 < v2:
+ if requested < vm:
t2 = tm
else:
t1 = tm
- elif d1 > d2:
- if 0 > dm:
+ elif v1 > v2:
+ if requested > vm:
t2 = tm
else:
t1 = tm
@@ -500,6 +499,11 @@ def future_past_equinox(delta, t = None):
# TODO document
+def future_past_equinox(delta, t = None):
+ return solar_prediction(delta, 0, solar_declination, t = t)
+
+
+# TODO document
def future_equinox(t = None):
return future_past_equinox(0.01 / 2000, t)
@@ -520,39 +524,7 @@ def future_past_elevation(delta, latitude, longitude, elevation, t = None):
@param t:float? The time in Julian Centuries, `None` for the current time
@return :float The calculated time point, `None` if none were found within a year
'''
- epsilon = 0.000001
- t = julian_centuries() if t is None else t
- t1 = t2 = t
- e1 = e0 = solar_elevation(latitude, longitude, t)
- while True:
- if abs(t2 - t) > 0.01:
- return None
- t2 += delta
- e2 = solar_elevation(latitude, longitude, t2)
- if (e1 <= elevation <= e2) or ((elevation >= e1 >= e2) and (elevation <= e0)):
- break
- if (e1 >= elevation >= e2) or ((elevation <= e1 <= e2) and (elevation >= e0)):
- break
- t1 = t2
- e2 = e1
- for _itr in range(1000):
- tm = (t1 + t2) / 2
- e1 = solar_elevation(latitude, longitude, t1)
- e2 = solar_elevation(latitude, longitude, t2)
- em = solar_elevation(latitude, longitude, tm)
- if abs(e1 - e2) < epsilon:
- return tm if abs(em - elevation) < epsilon else None
- if e1 < e2:
- if elevation < em:
- t2 = tm
- else:
- t1 = tm
- elif e1 > e2:
- if elevation > em:
- t2 = tm
- else:
- t1 = tm
- return None
+ return solar_prediction(delta, elevation, lambda t : solar_elevation(latitude, longitude, t), t = t)
def future_elevation(latitude, longitude, elevation, t = None):