Error in ecliptic_latlon method of position objects

[JoshPaterson]

To recreate the error first initialize some objects:

from skyfield.api import Loader
from numpy import array

path = ''
load = Loader(path)
ts = load.timescale()
ephem = load('de430.bsp')
earth = ephem['earth']

Here are the results when I call the ecliptic_latlon method using a few different Time objects:

>>> earth.at(ts.utc(2017, 1)).ecliptic_latlon()
(<Angle -00deg 00' 40.3">, <Angle 100deg 17' 25.1">, <Distance 0.986598 au>)

>>> earth.at(ts.utc(2017, 1, array([1, 2]))).ecliptic_latlon()
ValueError: not enough values to unpack (expected 3, got 2)

>>> earth.at(ts.utc(2017, 1, array([1, 2, 3]))).ecliptic_latlon()
(<Angle 3 values from -38deg 49' 06.5" to +35deg 06' 55.5">,
 <Angle 3 values from 44deg 54' 02.9" to 227deg 39' 35.8">,
 <Distance [3.35845341e-01 1.67537362e+00 3.32415619e-04] au>)

>>> earth.at(ts.utc(2017, 1, array([1, 2, 3, 4]))).ecliptic_latlon()
ValueError: too many values to unpack (expected 3)

As far as I can tell it works when Time.tt is either a float or an array of size 3, but not otherwise. The traceback is the same for each of the above two exceptions:

  File "C:\Users\Josh\Miniconda3\lib\site-packages\skyfield\positionlib.py", line 248, in ecliptic_latlon
    d, lat, lon = to_polar(vector.au)

  File "C:\Users\Josh\Miniconda3\lib\site-packages\skyfield\functions.py", line 59, in to_polar
    x, y, z = xyz

[JoshPaterson]

It looks like this is the same error that's causing #207.

[brandon-rhodes]

It looks like it might involve @GammaSagittarii's code that tries to make the ecliptic coordinate system of different dates available. I'll go take a look and at least commit a test showing that Skyfield is broken until this gets fixed!

[brandon-rhodes]

Here's that commit, so that I can find the link again easily if I need to:

a0d2779

[brandon-rhodes]

I've just released Skyfield 1.9 with this fix. Thanks for the report!

[GammaSagittarii]

I am really sorry for this mess, I should have written more tests along with my code.

[brandon-rhodes]

@GammaSagittarii ­— I hope you are happy that you added a whole new feature, and even in tweaking your code I kept the essential bit you added, by knowing which of the Earth rotation attributes is essential to ecliptic coordinates! Getting coordinate operations written so that they work for both individual times and arrays of times is VERY fiddly, and it was only by drawing on years of experience doing it in Skyfield that I knew how to pivot the code.

Next time, we'll write a few more tests and raise the bar. But for now your feature is rescued and should give everyone good results!

[GammaSagittarii]

Of course I am happy, but I still felt embarrassed, @brandon-rhodes thank you for your encouraging words!

I think I still have this implementation of precession calculation for long time intervals - from this paper https://www.aanda.org/articles/aa/pdf/2011/10/aa17274-11.pdf .

I haven't made a pull request for it, because it is not working with arrays of times obviously :D .
It is more work than the ecliptic stuff I've done. I converted the code from SOFA implementation in c. (nice benefit are fixes for wrong numbers in the paper! )
I am not expert on licenses and stuff but since SOFA is not "pure" public domain, I think I would have to add some sort of notice like it is based on that function etc... (If I remember properly, basically it is public domain with attribution and not using the same names for functions)
Also since the calculations are more complex, I am guessing it might inflict performance regression on Skyfield if I just change precession function to this one. So this means I would have to introduce some sort of switch to Skyfield, for the precession model.
I don't know if this is valuable, but if you think I should proceed with this I will.

[brandon-rhodes]

I'm not sure whether we can port code from SOFA. I have always pulled code from USNO NOVAS and from the JPL SPICE toolkit instead, and I think it will be safest for Skyfield if we continue that tradition. Can you find any other implementations of that paper's code that aren't part of SOFA?

If not, then I'd recommend you release your Python version as a separate package. Skyfielders who wanted to use it would simply pip install your other package and tell Skyfield to use that as its precession algorithm. If there did turn out to be problems with permission to publish, you would simply take down the new package and destroy all copies of it, and Skyfield itself wouldn't be encumbered.

[GammaSagittarii]

OK, so it turns out there is already published code on PyPi (https://pypi.org/project/vondrak/) by @nebulousdog so I think publishing my version would be redundant. Unfortunately his code also contains functions from SOFA, and also needs to be modified to work for arrays of times if it was to be option for Skyfield.

I have not found any other implementations, aside from the original one in the paper that is in Fortran, but that one does not have any license stated explicitly, and it also uses a few sofa functions, although I think they can be avoided by using numpy, or Skyfield functions.

So with everything as it is, I think I will leave it be.

[brandon-rhodes]

Okay, that sounds like a reasonable approach for now — I'll let you know if I run across other implementations that are independent!