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[^1], 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} = \mathbf{R}_3(\chi_A) \mathbf{R}_1(-\omega_A) \mathbf{R}_3(-\psi_A) \mathbf{R}_1(\varepsilon_0)$

$\mathbf{P} = \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} $$

[^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.

Thank you,

[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:

https://github.com/skyfielders/python-skyfield/issues/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 [^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.

---

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.

[^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.