Fix modulo op for negative values

[trufanov-nok]

use numpy's fmod instead of python's % to be equal to C fmod() in case of negative values. resolves https://github.com/skyfielders/python-skyfield/issues/873

[brandon-rhodes]

@trufanov-nok — Thanks for pointing this out! I will also want to write a test that prevents a contributor from breaking this again in the future. Could you share the specific coordinates that led you to discover the difference in how modulo was handled, so that I could use them in the test? Thanks!

[trufanov-nok]

The test case would be

from skyfield import nutationlib
nutationlib.equation_of_the_equinoxes_complimentary_terms(jd_tt=2457506.7901703566)

It's a date somewhere in 2016 if I racall right
It'll give you after patch a fa array: [ 4.63751774e+00 8.26363095e+00 2.16215241e+00 1.07606250e+01 -3.32755093e+00 4.30222197e+02 1.69882417e+02 1.04308925e+02 6.07306002e+01 9.24542295e+00 4.35558967e+00 6.70191591e+00 5.93431245e+00 3.97985269e-03]

There the fa[4] is negative and it's the only value that will be counted in a = ke1.dot(fa) because ke1 = [0 0 0 0 1 0 0 0 0 0 0 0 0 0].

[brandon-rhodes]

It's a date somewhere in 2016 if I recall right It'll give you after patch a fa array…

I get only a scalar value:

2.1371671320680886e-09

Why would the output from the equation_of_the_equinoxes_complimentary_terms(jd_tt=2457506.7901703566) switch to an array, when given a scalar input like 2457506.7901703566? Maybe I misunderstood your comment?

[trufanov-nok]

When I referred to fa array and a variable I meant local variables which values I can check via print command. The function result is scalar.

[brandon-rhodes]

@trufanov-nok — Thank you for explaining! I understand now that you were referring to the internal value fa that's inside the function. But for a test case, I'm going to need a value that's visible outside the function, because the only thing a test case can do is call equation_of_the_equinoxes_complimentary_terms() and look at the result. Do you have a date I can try calling it with, that will give a different result once your patch is applied? With the two changes in your patch applied, I still get 2.137167132068086e-09 as the return value. (But maybe I've applied them wrong, so please double-check me!)

(Also: could you double-check whether the first line you change in your patch makes a difference to the numeric outcome? I cannot tell quite why, in the stack of equations for fa[0] and fa[1] and so forth, you choose fa[4] alone to adjust.)

[trufanov-nok]

The fa[4] is chosen because python's modulo operation behaves differently only with negative arguments, and assuming time t is always positive the fa[4] is the only place that have a negative coefficient.

Indeed, I was comparing the local variables instead of function result. And now I see that effect on final result is almost neglectable small: for t = 2457506.7901703566

before patch: 2.1371671320680886e-09
after patch:  2.137167132068086e-09

[brandon-rhodes]

before patch: 2.1371671320680886e-09
after patch:  2.137167132068086e-09

Yes, that difference looks to be down at the level of the machine precision of a 64-bit floating point number. Unless this affects an astronomical computation to within its error bounds — and please re-open this PR if you find a case where it does! — I am inclined to keep using the built-in % Python operator, whose operation is much simpler to reason about — since something like '% 7.0' has a simple-to-explain output range [0.0,7.0), instead of fmod's rather crazy range (-7.0,7.0) where pairs of values like -4.0 and +3.0 kind-of mean exactly the same thing except expressed differently.