Pyephem 15 minutes discrepancy when manually computing equinoxes and solstices

Question

the code below tries to compute manually the first equinox of 2019. It returns

('d1=', 2019/3/20 21:43:48) ('d2=', 2019/3/20 21:43:49) 2019/3/20 21:58:31

that is, a discrepancy of 15 minutes with the real equinox. Is this normal? Did I forget something? The problem also occurs with the solstices, and also if I used the integrated newton method. Could it have something to do with the epoch of computation?

Thanks,

Dennis

import ephem
sun = ephem.Sun()

# computing Spring equinox:
d1 = ephem.Date('2019/03/15')
d2 = ephem.Date('2019/03/25')
a=ephem.degrees('180.0')

for i in range(20):
  #middle date
  d3=(d1+d2)/2
  sun.compute(d3)
  if sun.hlon>a:
      d2=d3
  else:
      d1=d3

print("d1=",ephem.Date(d1))
print("d2=",ephem.Date(d2))
d1 = ephem.next_equinox('2019')
print(d1)

You can click on the [edit](https://stackoverflow.com/posts/54208690/edit) link below your post to update your post :) — Gino Mempin, Jan 16 '19 at 00:42
The problem also occurs when using Newton's method within pyephem: — Dennis, Jan 16 '19 at 13:58

score 1 · Answer 1 · answered Jan 19 '19 at 05:09

It looks like the difference is because PyEphem's underlying astronomy library always measures heliocentric longitude relative to the coordinates of J2000, which by the date you are asking about is noticeably different from the coordinates-of-date which are used to define the equinox.

Try running this as your compute step:

sun.compute(d3, epoch=d3)

and then look for when sun.ra is zero degrees; the result should be the equinox. I'll see about getting the PyEphem Quick Reference updated to note that heliocentric coordinates don't seem to pay attention to the epoch= parameter.

score 0 · Answer 2 · answered Jan 20 '19 at 16:27

Many thanks, Brandon, this is very helpful and I am finally getting the correct value! In fact, it seems that the equinoxes are defined by the right ascension being equal to 0h, 6h, 12h, 18h, and not the heliocentric longitude being 0, 90, 180, 270. There is a slight difference between ra and hlon, when you run the code below. But this leads to another question. The Wikipedia page https://en.wikipedia.org/wiki/Equinox says that the equinoxes are defined by the longitude being 0 or 180. So who is correct?

import ephem
sun = ephem.Sun()

d1 = ephem.Date('2019/03/15')
d2 = ephem.Date('2019/03/25')
a=ephem.degrees('0.0') # or 90, or 180, or 270

def spring_equinox(date):
  sun.compute(date)
  return ephem.degrees(sun.ra - a).znorm

d = ephem.newton(spring_equinox, d1, d2)
print(ephem.Date(d))
print sun.ra
print sun.hlon

(Before someone else points it out: this might have fit better as an update to your main question, maybe by editing and having it in a section called "Follow up:".) The RA and longitude should be 0 at the same time. The difference is that PyEphem is willing to give you RA in the coordinate system of the date you asked about, but longitude only in J2000 coordinates. Try reading about "precession", maybe on the Wikipedia, to understand why the coordinate systems shift where 0 is every year. — Brandon Rhodes, Jan 20 '19 at 21:16
Thanks. I know about precession, it's just that I didn't realize that RA and longitude can use different epochs! Many thanks again. — Dennis, Jan 21 '19 at 10:49

Pyephem 15 minutes discrepancy when manually computing equinoxes and solstices

2 Answers2