skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.42k stars 212 forks source link

Sun geocentric true coordinates #695

Closed bsics closed 2 years ago

bsics commented 2 years ago

Hi, I implemented VSOP87 and ELP2000-82B to calculate sun and moon geocentric positions, and wanted to validate my implementations against JPL ephemerides. Out of convenience, I used skyfield. Everything was perfect for the Moon, but with the Sun positions I had consistently about 20 arcsec difference in the true coordinates between VSOP87 and skyfield. It was bad, because the estimated VSOP87 accuracy is ~2 arcsec

Here is the difference I get for 700+ positions over the span of 7 years. dpos_sun

Interestingly enough, the difference is gone when I remove the aberation (ligth time) correction in my code: dpos_sun_noaberration

For that I went straight to Horizons. I found that the Horizons web application results are coherent with my VSOP87 code and not with skyfield.

For exemple, skyfield gives:

>>> planets = load('de440.bsp')
    ts = load.timescale()
    y, mo, d = (2012, 9, 27.084351851851853)
    t = ts.utc(y, mo, d)
    psun = planets['Sun']
    pearth = planets['earth']
    pos_sun = pearth.at(t).observe(psun)
    lat_sun, lon_sun, distance_sun = pos_sun.frame_latlon(ecliptic_frame)
>>>    lon_sun.degrees
184.38185171462558
>>>    lat_sun.degrees
0.0001164730838003623

My VSOP87 code gives:

>>>lon_vsop
184.37626251244774

And Horizons gives:

Date_________JDTT     R.A.___(ICRF)___DEC  R.A._(a-appar)_DEC.     ObsEcLon    ObsEcLat
***************************************************************************************
$$SOE
2456197.585129456     183.85449  -1.66898  184.01641  -1.73901  184.3761625   0.0001113
$$EOE
***************************************************************************************

I use skyfield v1.35 installed from anaconda, the conda_forge channel. The same conclusions with ecliptic or equatorial coordiantes.

brandon-rhodes commented 2 years ago

If you want Skyfield to apply aberration, you need to call .apparent() on the astrometric position returned by .observe(). Could you try that and see if it changes the result?

(Oh, also: I've edited your question so that it uses GitHub's Markdown formatting on your code; double-check to make sure you like how I applied their formatting. Thanks!)

bsics commented 2 years ago

That is great. This solves the systematic discrepancies I observed. For the date I used above:

>>>    lon_sun.degrees
184.37617577276538

Thank you!

brandon-rhodes commented 2 years ago

I'm glad it agrees more closely now! Here's the documentation in case you at any point want to read more:

https://rhodesmill.org/skyfield/positions.html#barycentric-astrometric-apparent