I would like to use pyephem or skyfield1 to find the degree of the local horizon (starting from South = 0°) at which is located the highest point of the ecliptic.
Note that I am not trying to find the culmination of a planet. It could well be that, when looking South, there is no planet on the ecliptic by that time/latitude, but there will still be a highest point of the ecliptic. Can anyone help?
1Skyfield is preferred over PyEphem for new projects for several reasons
Edit:
Getting the exact value winds up being simple, but it was only later in the day that the answer occurred to me. All we need to do is ask for the altitude and azimuth of the Ecliptic’s south pole! The ecliptic itself will be highest above our horizon in exactly that direction, at an altitude that is exactly 90° above its pole.
Thus the answer is given by the coordinates of a single vector:
from skyfield import api, framelib
from skyfield.positionlib import Apparent
ts = api.load.timescale()
t = ts.utc(2021, 2, 16, 22, 52)
bluffton = api.Topos('40.74 N', '84.11 W')
f = framelib.ecliptic_frame
r = api.Distance([0, 0, -1])
v = api.Velocity([0, 0, 0])
p = Apparent.from_time_and_frame_vectors(t, f, r, v)
p.center = bluffton
alt, az, distance = p.altaz()
print('Altitude of highest point on ecliptic:', alt.degrees + 90.0)
print('Azimuth of highest point on ecliptic:', az.degrees)
The answer:
Altitude of highest point on ecliptic: 67.54977465183501
Azimuth of highest point on ecliptic: 162.61593086793783
— is far more exact, and faster to compute, than my first try. Enjoy!
Original answer:
While it’s likely that someone will soon jump in with a closed-form solution that uses spherical trigonometry — in which case it’s likely that Skyfield or PyEphem will only be used to determine the Earth orientation — here’s a quick way to get an answer within about a degree:
The result agrees closely to what I see if I open Stellarium, turn on the ecliptic, and choose a field star right near the point where the ecliptic reaches the highest point in the sky.
import numpy as np
from skyfield import api, framelib
from skyfield.positionlib import Apparent
𝜏 = api.tau
ts = api.load.timescale()
eph = api.load('de421.bsp')
bluffton = api.Topos('40.74 N', '84.11 W')
t = ts.utc(2021, 2, 16, 22, 52)
angle = np.arange(360) / 360.0 * 𝜏
zero = angle * 0.0
f = framelib.ecliptic_frame
d = api.Distance([np.sin(angle), np.cos(angle), zero])
v = api.Velocity([zero, zero, zero])
p = Apparent.from_time_and_frame_vectors(t, f, d, v)
p.center = bluffton
alt, az, distance = p.altaz()
i = np.argmax(alt.degrees) # Which of the 360 points has highest altitude?
print('Altitude of highest point on ecliptic:', alt.degrees[i])
print('Azimuth of highest point on ecliptic:', az.degrees[i])
The result:
Altitude of highest point on ecliptic: 67.5477569215633
Azimuth of highest point on ecliptic: 163.42529398930515
This is probably a computationally expensive enough approach that it won’t interest you once you or someone else does the spherical trigonometry to find an equation for the azimuth; but at the very least this might provide numbers to check possible formulae against.
