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
137
138
139
140
141
142
|
"""
Astronomical functions
"""
from math import pi, radians, tau
try:
# As of Python 3.9, included in the stdlib
from zoneinfo import ZoneInfo
except ImportError: # pre-3.9
from backports.zoneinfo import ZoneInfo
from skyfield.api import Loader
from skyfield import almanac
from skyfield_data import get_skyfield_data_path
from datetime import date, timedelta
# Parameter for the newton method to converge towards the closest solution
# to the function. By default it'll be an approximation of a 10th of a second.
hour = 1 / 24
minute = hour / 60
second = minute / 60
newton_precision = second / 10
def calculate_equinoxes(year, timezone='UTC'):
""" calculate equinox with time zone """
tz = ZoneInfo(timezone)
load = Loader(get_skyfield_data_path())
ts = load.timescale()
planets = load('de421.bsp')
t0 = ts.utc(year, 1, 1)
t1 = ts.utc(year, 12, 31)
datetimes, _ = almanac.find_discrete(t0, t1, almanac.seasons(planets))
vernal_equinox = datetimes[0].astimezone(tz).date()
autumn_equinox = datetimes[2].astimezone(tz).date()
return vernal_equinox, autumn_equinox
def get_current_longitude(current_date, earth, sun):
"""
Return the ecliptic longitude, in radians.
"""
astrometric = earth.at(current_date).observe(sun)
latitude, longitude, _ = astrometric.ecliptic_latlon(epoch='date')
return longitude.radians
def newton(f, x0, x1, precision=newton_precision,
**func_kwargs):
"""Return an x-value at which the given function reaches zero.
Stops and declares victory once the x-value is within ``precision``
of the solution, which defaults to a half-second of clock time.
"""
f0, f1 = f(x0, **func_kwargs), f(x1, **func_kwargs)
while f1 and abs(x1 - x0) > precision and f1 != f0:
new_x1 = x1 + (x1 - x0) / (f0 / f1 - 1)
x0, x1 = x1, new_x1
f0, f1 = f1, f(x1, **func_kwargs)
return x1
def newton_angle_function(t, ts, target_angle, body1, body2):
"""
Compute the longitude of body2 relative to body1
In our case, it's Earth & Sun, but it could be used as any other
combination of solar system planets/bodies.
"""
# We've got a float which is the `tt`
sky_tt = ts.tt_jd(t)
longitude = get_current_longitude(sky_tt, body1, body2)
result = target_angle - longitude
if result > pi:
result = result - pi
if result < -pi:
result = result + pi
return result
def solar_term(year, degrees, timezone='UTC'):
"""
Returns the date of the solar term for the given longitude
and the given year.
Solar terms are used for Chinese and Taiwanese holidays
(e.g. Qingming Festival in Taiwan).
More information:
- https://en.wikipedia.org/wiki/Solar_term
- https://en.wikipedia.org/wiki/Qingming
This function is adapted from the following topic:
https://answers.launchpad.net/pyephem/+question/110832
"""
# Target angle as radians
target_angle = radians(degrees)
load = Loader(get_skyfield_data_path())
planets = load('de421.bsp')
earth = planets['earth']
sun = planets['sun']
ts = load.timescale()
tz = ZoneInfo(timezone)
jan_first = ts.utc(date(year, 1, 1))
current_longitude = get_current_longitude(jan_first, earth, sun)
# Find approximately the right time of year.
difference = (target_angle - current_longitude) % tau
# Here we have an approximation of the number of julian days to go
date_delta = 365.25 * difference / tau
# convert to "tt" and reconvert it back to a Time object
t0 = ts.tt_jd(jan_first.tt + date_delta)
# Using datetimes to compute the next step date
t0_plus_one_minute = t0.utc_datetime() + timedelta(minutes=1)
# Back to Skyfield Time objects
t0_plus_one_minute = ts.utc(t0_plus_one_minute)
# Julian day for the starting date
t0 = t0.tt
# Adding one minute to have a second boundary
t0_plus_one_minute = t0_plus_one_minute.tt
# Newton method to converge towards the target angle
t = newton(
newton_angle_function, t0, t0_plus_one_minute,
precision=newton_precision,
# Passed as kwargs to the angle function
ts=ts,
target_angle=target_angle,
body1=earth,
body2=sun,
)
# Here we have a float to convert to julian days.
t = ts.tt_jd(t)
# To convert to datetime
t = t.utc_datetime()
# Convert in the timezone
result = t.astimezone(tz)
return result.date()
|
