Erroneous documentation of rotations in precessionlib.compute_precession?

[wedesoft]

The documentation in the precession method is:

# Compute elements of precession rotation matrix equivalent to
# R3(chi_a) R1(-omega_a) R3(-psi_a) R1(epsilon_0).

I implemented the matrix multiplication using separate rotation matrices and tested it with real values and it seems the documentation should rather say:

# Compute elements of precession rotation matrix equivalent to
# R3(-chi_a) R1(omega_a) R3(psi_a) R1(-epsilon_0).

Obviously this means also flipping the signs in the sinus and cosinus terms above and flipping the signs for the sinus terms in the rotation matrix below as well for consistency. Let me know if I should make a corresponding pull request.

[brandon-rhodes]

Doing a quick search, it looks like that comment was copied from the NOVAS library, that I used for the core astrometry calculations inside Skyfield:

https://github.com/brandon-rhodes/python-novas/blob/ac85c2b0577026c6a64a2af1564bbb699d1ecd06/Cdist/novas.c#L6455

Before making any change to the comment in Skyfield itself, I think I'd like to understand whether (a) it was already an incorrect comment in the NOVAS code, or (b) it was correct in the original code but I swapped signs as I copied the code into Skyfield and so it's now incorrect.

Also, it would be helpful to see your code that uses separate rotation matrices. Is it in a condition that you could share it as a small stand-alone script? Thanks!

[wedesoft]

Hi, Thanks for your quick response. The code is just a Clojure Gist at the moment. I tested it with tdb = 2460461.5 and I can reproduce the result of the Python code in Skyfield. Great work by the way :)

[wedesoft]

I get the following result with R3(-chi_a) R1(omega_a) R3(psi_a) R1(-epsilon_0) which matches skyfield:

[[0.9999822866051156, -0.005458997301346284, -0.0023718820519731975]
 [0.0054589973730140105, 0.999985099542067, -6.4438907311543E-6]
 [0.0023718818870265717, -6.504321302780447E-6, 0.9999971870630475]]

If I use R3(chi_a) R1(-omega_a) R3(-psi_a) R1(epsilon_0), I get:

[[0.9999822866051156, 0.005458997301346284, -0.0023718820519731975]
 [-0.0054589973730140105, 0.999985099542067, 6.4438907311543E-6]
 [0.0023718818870265717, 6.504321302780447E-6, 0.9999971870630475]]

Kind regards Jan

[brandon-rhodes]

We should check whether the rotations are defined in the same, or opposite directions, in the two libraries. Maybe you could take a simple vector like [1 2 3] and rotate it by a positive R1 of +45°, and compare the output in both Clojure and Skyfield?

[wedesoft]

Ok, using my definition of R1 in Clojure, I get:

(mulv (r1 (/ PI 4)) (vec3 1 2 3))
; [1.0, -0.7071067811865472, 3.5355339059327378]

[wedesoft]

Using skyfield.functions, I get the same:

import math
import numpy as np
from skyfield.functions import rot_x
np.dot(rot_x(math.pi/4), [1, 2, 3])
# array([ 1.        , -0.70710678,  3.53553391])

[wedesoft]

Using skyfield.functions.mxv instead of np.dot also returns the same.

[brandon-rhodes]

Thanks for those final checks; I've made the change, so hopefully the comment doesn't cause any further confusion for folks!

[wedesoft]

Ok, that's great!