brunobord
commented
5 years ago I think I'm going a bit crazy. I think I've managed to convert from the "ephem" function to a "skyfield" function, but unfortunately, my function fails at the year 2009, for example, while the "ephem" function passes all the tests.
Here it is (please forgive its eventual clunky look, it's the result of hours of try/fail attempts):
def solar_term_skyfield(year, degrees, timezone='UTC'):
# Target angle as Astropy unit
target_angle = (degrees * u.degree).to('rad').value
load = Loader(get_skyfield_data_path())
planets = load('de421.bsp')
earth = planets['earth']
sun = planets['sun']
ts = load.timescale()
tz = pytz.timezone(timezone)
jan_first = ts.utc(date(year, 1, 1))
current_longitude = current_longitude_skyfield(jan_first, earth, sun)
# Find approximately the right time of year.
difference = (target_angle - current_longitude) % twopi
# Here we have an approximation of the number of julian days to go
date_delta = 365.25 * difference / twopi
# convert to "tt" and reconvert it back to a Time object
t0 = ts.tt_jd(jan_first.tt + date_delta)
def f(t):
# We've got a float which is the `tt`
sky_tt = ts.tt_jd(t)
longitude = current_longitude_skyfield(sky_tt, earth, sun)
result = target_angle - longitude
if result > pi:
result = result - pi
elif result < -pi:
result = result + pi
return result
# 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(f, t0, t0_plus_one_minute, newton_precision_ephem)
# 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()The function current_longitude_skyfield is, in my opinion, the key, here. Using the same units (radians), what ephem gives is (slightly) different from what my function generates.
def current_longitude_skyfield(current_date, earth, sun):
"""
Return the ecliptic longitude, in radians.
"""
astrometric = earth.at(current_date).observe(sun)
latitude, longitude, _ = astrometric.ecliptic_latlon()
return longitude.radiansI thought this was the equivalent of the following code:
def current_longitude_ephem(year):
# Find out the sun's current longitude.
sun = ephem.Sun(ephem.Date(str(year)))
return sun.hlong - piThe angle is almost the same, but it's not. I've got hours of delta when calculating the solar term using my method, and when you're close to midnight, the consequence is here: I'm wrong.
Unfortunately, I'm unable to find what's wrong in this function. I genuinely thought it was the right way to calculate the ecliptic angle at a given date...
If anyone can help me, please, I'll be delighted.
GammaSagittarii
brandon-rhodes
The library
workalendaris using(py)ephemto compute solar terms[^1], used in Hong Kong, Taiwan and potentially other calendars. The code used to calculate the Solar term was based on https://answers.launchpad.net/pyephem/+question/110832For the record, here it is:
Unfortunately, I've tried to convert this function using Skyfield and as I felt it was almost ready, it started to fall apart.
Here are some example test functions that the new function should be able to pass:
(please note that the "ephem" function returns a 10/10 card)
References:
[^1]: Solar terms - https://en.wikipedia.org/wiki/Solar_term