Some mistakes in "precessionlib.py"

[reza-ghazi]

Hello,

There are some mistakes in the precessionlib.py $\psi_A$ and $\omega_A$ series terms. I attached the picture shot from the table on (ESAA¹, page 216), highlighting the correct numbers that should be replaced with the current values in the file.

---

On the other hand, the Precession matrix is defined as:

$$ \mathbf{P}=\mathbf{R}_3\left(\chi_A\right) \mathbf{R}_1\left(-\omega_A\right) \mathbf{R}_3\left(-\psi_A\right) \mathbf{R}_1\left(\varepsilon_0\right) \rightarrow \text{(ESAA, page 220)} $$

But you have mentioned in the docstring that "R3(-chi_a) R1(omega_a) R3(psi_a) R1(-epsilon_0)"

---

Let me explain more here:

The equations describe the two precession matrices:

$\mathbf{P_1} = \mathbf{R}_3(\chi_A) \mathbf{R}_1(-\omega_A) \mathbf{R}_3(-\psi_A) \mathbf{R}_1(\varepsilon_0)$

$\mathbf{P_2} = \mathbf{R}_3(-\chi_A) \mathbf{R}_1(\omega_A) \mathbf{R}_3(\psi_A) \mathbf{R}_1(-\varepsilon_0)$

They are not the same, but they are closely related. Their relationship depends on the reference frames and the direction of transformation being considered.

Key Differences

1. Transformation Direction:

   • The first matrix is likely designed for transforming coordinates from an earlier epoch (e.g., J2000.0) to the epoch of date.
   • The second matrix is likely designed for transforming coordinates from the epoch of date back to an earlier epoch (e.g., J2000.0).

2. Rotational Signs:

   • The angles $\chi_A, \omega_A, \psi_A, \varepsilon_0$ represent precession parameters.
   • Switching the signs of these angles corresponds to inverting the direction of the rotation. For example:
     • $\mathbf{R}_3(\chi_A)$ rotates counterclockwise around the $z$-axis by $\chi_A$.
     • $\mathbf{R}_3(-\chi_A)$ rotates clockwise around the $z$-axis by the same angle.

3. Inverse Transformation:

   • The second matrix is the inverse of the first matrix:

$$ \mathbf{P}_2 = \mathbf{P}_1^{-1} $$

• This is because reversing the rotations' order and flipping the angles' signs results in the inverse transformation.

---

Mathematical Relationship

The relationship between the two matrices is:

$$ \mathbf{P}_2 = \mathbf{P}_1^{-1} $$

Where:

• $\mathbf{P}_1 = \mathbf{R}_3(\chi_A) \mathbf{R}_1(-\omega_A) \mathbf{R}_3(-\psi_A) \mathbf{R}_1(\varepsilon_0)$
• $\mathbf{P}_2 = \mathbf{R}_3(-\chi_A) \mathbf{R}_1(\omega_A) \mathbf{R}_3(\psi_A) \mathbf{R}_1(-\varepsilon_0)$

The inverse of a rotation matrix is equivalent to its transpose when the matrix is orthonormal (as all rotation matrices are):

$$ \mathbf{P}_1^{-1} = \mathbf{P}_1^T $$

Thus, the two matrices are inversely related, meaning applying one matrix followed by the other will return the original coordinates:

$$ \mathbf{P}_2 \mathbf{P}_1 = \mathbf{I} $$

---

Physical Interpretation

1. First Matrix ($\mathbf{P}_1$):
   • Transforms a position vector $\mathbf{r}_0$ from the mean equator and equinox of J2000.0 to the mean equator and equinox of date:

$$ \mathbf{r} = \mathbf{P}_1 \mathbf{r}_0 $$

2. Second Matrix ($\mathbf{P}_2$):
   • Transforms a position vector $\mathbf{r}$ from the mean equator and equinox of date back to the mean equator and equinox of J2000.0:

$$ \mathbf{r}_0 = \mathbf{P}_2 \mathbf{r} $$

---

Conclusion

• The two matrices represent inverse transformations.
• Both are valid but correspond to opposite transformation directions:
  • $\mathbf{P}_1$: Forward transformation (J2000.0 to epoch of date).
  • $\mathbf{P}_2$: Reverse transformation (epoch of date to J2000.0).
• Their relationship is such that:

$$ \mathbf{P}_2 = \mathbf{P}_1^{-1} $$

Thank you,

Footnotes

1. Urban, S. E., & Seidelmann, P. K. (2013). Explanatory Supplement to the Astronomical
   Almanac (3rd ed.). University Science Books. ISBN: 978-1-891389-85-6. ↩

[brandon-rhodes]

@reza-ghazi — The formulae in precessionlib.py do not compute the precession of the ecliptic, as seems to be the case in the table you have shared, but the precession of the equator. The numbers Skyfield uses are the official numbers from IAU 2000 from this paper:

https://www.aanda.org/articles/aa/full/2003/48/aa4068/aa4068.html

Look in particular at formulae 37 and 39 from that paper, and I think you will see that they agree exactly with Skyfield's code.

It is interesting that you mention the comment — it used to say:

# R3(chi_a) R1(-omega_a) R3(-psi_a) R1(epsilon_0).

But then a contributor complained that the rotation matrices in the comment were all backwards. Here is the issue they opened:

#972

And so the comment got changed, because when I did a quick test, it did indeed look like rotation matrices worked in the direction that they claimed — was I wrong in letting them convince me, and I should not have made the change? If you can point out where that contributor made a mistake, I will be happy to reverse the change and return the comment to its original form! That way, Skyfield would again agree with the NOVAS library that the logic was adapted from, which has exactly the same comment.

[reza-ghazi]

@brandon-rhodes - I've read the #972; I think the following explanation makes it clear:

The definitions of the rotation matrices depend on the direction of rotation (clockwise vs. counterclockwise) and the convention (e.g., physics vs. mathematics vs. computer graphics). The definitions from Wikipedia and my earlier explanation are correct but represent different conventions.

---

1. Rotation Matrix Differences

Wikipedia's Definitions (resource of #972)

1. Rotation About the $x$-Axis:

$$ R_x(\theta) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos\theta & -\sin\theta \\ 0 & \sin\theta & \cos\theta \end{bmatrix} $$

2. Rotation About the $y$-Axis:

$$ R_y(\theta) = \begin{bmatrix} \cos\theta & 0 & \sin\theta \\ 0 & 1 & 0 \\ -\sin\theta & 0 & \cos\theta \end{bmatrix} $$

3. Rotation About the $z$-Axis:

$$ R_z(\theta) = \begin{bmatrix} \cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \end{bmatrix} $$

My Definitions

1. Rotation About the $x$-Axis:

$$ R_x(\phi) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos\phi & \sin\phi \\ 0 & -\sin\phi & \cos\phi \end{bmatrix} $$

2. Rotation About the $y$-Axis:

$$ R_y(\phi) = \begin{bmatrix} \cos\phi & 0 & -\sin\phi \\ 0 & 1 & 0 \\ \sin\phi & 0 & \cos\phi \end{bmatrix} $$

3. Rotation About the $z$-Axis:

$$ R_z(\phi) = \begin{bmatrix} \cos\phi & \sin\phi & 0 \\ -\sin\phi & \cos\phi & 0 \\ 0 & 0 & 1 \end{bmatrix} $$

---

2. Key Difference: Direction of Rotation

• Wikipedia uses a clockwise rotation convention when viewed along the positive axis ( e.g., positive $z$-axis for $R_z(\theta)$ ).

• My earlier explanation follows the counterclockwise rotation convention when viewed along the positive axis.

Why This Happens

• A clockwise rotation in a right-handed coordinate system is equivalent to a counterclockwise rotation with a negated angle.
• This means the two conventions are mathematically related:

$$ R_{\text{Wikipedia}}(\theta) = R_{\text{My Definition}}(-\phi) $$

---

3. Which One to Use?

The choice depends on our application:

• Physics and Mathematics: Counterclockwise is more common because it aligns with the right-hand rule.
• Computer Graphics: We often use Clockwise rotation because of screen-space conventions.
• Astronomy (Precession): Counterclockwise matches the celestial coordinate system's conventions.

---

4. Converting Between Conventions

If we want to convert from one to the other:

• Negate the angle:

$$ R_{\text{My Definition}}(\phi) = R_{\text{Wikipedia}}(-\phi) $$

---

5. Conclusion

• Both definitions are correct, but they represent rotations in opposite directions.
• In practice, choose the convention that aligns with our field or application. Counterclockwise (as I explained earlier) is standard for astronomy and celestial mechanics.

[reza-ghazi]

For your information, I checked the resource ¹ of The ESAA and found that your values for $\psi_A$ and $\omega_A$ are correct; strangely, they have made a mistake in producing the same table in ESAA.

---

I also checked the Expressions for IAU 2000 precession quantities again, and I found that the matrix multiplication is precisely the same as I mentioned.

Footnotes

1. Hilton, J. L., N. Capitaine, J. Chapront, J. M. Ferrandiz, A. Fienga, T. Fukushima, J. Getino, P. Mathews, J.-L. Simon, M. Soffel, J. Vondrak, P. Wallace, and J. Williams (2006). Report of the International Astronomical Union Division I Working Group on Precession and the Ecliptic. Celestial Mechanics and Dynamical Astronomy 94, 351-367. ↩

[brandon-rhodes]

For your information, I checked the resource [^1] of The ESAA and found that your values for $\psi_A$ and $\omega_A$ are correct; strangely, they have made a mistake in producing the same table in ESAA.

Thank you for double-checking another resource! I am happy to learn that Skyfield's code is correct in this instance.

In practice, choose the convention that aligns with our field or application. Counterclockwise (as I explained earlier) is standard for astronomy and celestial mechanics.

Thank you very much for such a full and thorough explanation! I had no idea that there were two competing rotation directions in use by different fields. To me, the traditional trigonometric interpretation of a rotation starting on the $x$-axis and swinging up towards the $y$-axis is so completely fundamental (image from Wikipedia):

— that I would never even have imagined defining $R_3(\theta)$ such that positive $\theta$ swings the vector in the opposite direction, down towards negative $y$. But here it is, right in my copy of The Explanatory Supplement to the Astronomical Almanac (which, alas, I don't think I owned yet at the time; though even had I owned it, it might never have occurred to me to flip through looking for its matrix rotations):

Those are exactly the matrices you suggested above, and the ones to which the NOVAS library would have been referring in the comment in precessionlib.py.

(One confusion: I am not sure exactly what you mean by "aligns with the right-hand rule", so could you explain that idea a bit more? When I curl the fingers of my right hand with my thumb sticking straight up, the fingers curl counter-clockwise from the $x$-axis towards the $y$-axis, just as on the standard $xy$-plane. Which, incidentally, is also the direction the Earth rotates on its $z$-axis, and the direction in which the Earth revolves around the Solar System's $z$-axis. So I am mystified that astronomy chose the opposite rotation direction for $R3(\theta)$, and would appreciate a bit more intuition about why that direction seemed natural.)

I will plan to do the following:

• I will probably add functions R1(), R2(), and R3() to functions.py that perform the opposite rotations from rot_x() and its friends. (Alas, I can't ever change the definitions of rot_x() and its friends, as that would break user code that might already call them.)
• Though I'm unlikely to go rewriting all of Skyfield right now to switch to the new functions, they will immediately be available for users, and I'll plan on gradually moving to them in the future so that Skyfield's internals better match the Explanatory Supplement.
• As soon as the new functions are defined, I will go ahead and revert the change I made in Erroneous documentation of rotations in precessionlib.compute_precession? #972 so that the comment returns to its original form that I cut-and-pasted from the NOVAS library. The next time a question about it comes up, I'll be able to refer to the new definitions for R1() etc in functions.py as Skyfield's answer for what those rotations mean.

Once those changes have landed, will that resolve this issue? Or can you think of any other tweaks that will be necessary before closing it?

[reza-ghazi]

Frankly, the best tweak is updating the current codes. However, users can also negate the angles and get the correct value. Recall that:

$$ R_{\text{Current Skyfield Definition}}(\theta) = R_{\text{Correct Definition}}(-\phi) $$

However, a deep understanding of the underlying concepts is needed. But at least it can be a tweak.

[brandon-rhodes]

A slight surprise as I sat down and submitted the code: the Supplement only reverses the rotation around the $x$ and $z$ axes. But for its $y$ rotation, it follows the standard right-hand rule! I wonder what the history is there. 🙂