Clarify documentation of RotateAndTiltWrapper's used of galaxy/sky PA and inclination and test explicitly that the docs are correct

Author: jobovy

galpy/potential/RotateAndTiltWrapperPotential.py

@@ -32,7 +32,13 @@ class RotateAndTiltWrapperPotential(WrapperPotential):
 
     * A rotation around the original z-axis (`galaxy_pa`), the `inclination`, and a rotation around the new z axis (`sky_pa`).
 
-    The second option allows one to specify the inclination and sky position angle (measured from North) in the usual manner in extragalactic observations.
+    The second option allows one to specify the inclination and sky position angle (measured from North) in the usual manner in extragalactic observations, that is, the relation between the observed coordinates :math:`(x',y',z')` and the intrinsic coordinates :math:`(x,y,z)` is given by
+
+    .. math::
+
+        \\begin{pmatrix} x \\\\ y \\\\ z \\end{pmatrix} = \\begin{pmatrix} \\phantom{-}\\cos(\\text{PA}_\\text{gal}) & \\sin(\\text{PA}_\\text{gal}) & 0 \\\\ -\\sin(\\text{PA}_\\text{gal}) & \\cos(\\text{PA}_\\text{gal}) & 0 \\\\ 0 & 0 & 1 \\end{pmatrix} \\begin{pmatrix} 1 & 0 & 0 \\\\ 0 & \\phantom{-}\\cos i & \\sin i \\\\ 0 & -\\sin i & \\cos i \\end{pmatrix} \\begin{pmatrix} \\sin(\\text{PA}_\\text{sky}) & -\\cos(\\text{PA}_\\text{sky}) & 0 \\\\ \\cos(\\text{PA}_\\text{sky}) & \\phantom{-}\\sin(\\text{PA}_\\text{sky}) & 0 \\\\ 0 & 0 & 1 \\end{pmatrix} \\begin{pmatrix} x' \\\\ y' \\\\ z' \\end{pmatrix}\\,,
+
+
     A final `offset` option allows one to apply a static offset in Cartesian coordinate space to be applied to the potential following the rotation and tilt.
     """
 

tests/test_potential.py

@@ -6236,6 +6236,61 @@ def test_RotateAndTiltWrapper():
     ), "Wrapped + Offset and Internally rotated NFW potentials do not match when evaluated at the same point"
 
 
+def test_RotateAndTiltWrapper_pa_inclination_rotation_matrix():
+    # Test that the formula for the rotation matrix given in the documentation agrees with the code
+    def rotation_matrix_docs(galaxy_pa, inclination, sky_pa):
+        return numpy.dot(
+            numpy.array(
+                [
+                    [numpy.cos(galaxy_pa), numpy.sin(galaxy_pa), 0.0],
+                    [-numpy.sin(galaxy_pa), numpy.cos(galaxy_pa), 0.0],
+                    [0.0, 0.0, 1.0],
+                ]
+            ),
+            numpy.dot(
+                numpy.array(
+                    [
+                        [1.0, 0.0, 0.0],
+                        [0.0, numpy.cos(inclination), numpy.sin(inclination)],
+                        [0.0, -numpy.sin(inclination), numpy.cos(inclination)],
+                    ]
+                ),
+                numpy.array(
+                    [
+                        [numpy.sin(sky_pa), -numpy.cos(sky_pa), 0.0],
+                        [numpy.cos(sky_pa), numpy.sin(sky_pa), 0.0],
+                        [0.0, 0.0, 1],
+                    ]
+                ),
+            ),
+        )
+
+    galaxy_pa, inclination, sky_pa = 0.3, -0.2, 0.1
+    rtwp = potential.RotateAndTiltWrapperPotential(
+        pot=potential.MWPotential2014,
+        galaxy_pa=galaxy_pa,
+        inclination=inclination,
+        sky_pa=sky_pa,
+    )
+    assert numpy.all(
+        numpy.fabs(rtwp._rot - rotation_matrix_docs(galaxy_pa, inclination, sky_pa))
+        < 1e-10
+    ), "Rotation matrix in RotateAndTiltWrapperPotential does not agree with the formula in the documentation"
+
+    galaxy_pa, inclination, sky_pa = -0.3, 1.2, 2.1
+    rtwp = potential.RotateAndTiltWrapperPotential(
+        pot=potential.MWPotential2014,
+        galaxy_pa=galaxy_pa,
+        inclination=inclination,
+        sky_pa=sky_pa,
+    )
+    assert numpy.all(
+        numpy.fabs(rtwp._rot - rotation_matrix_docs(galaxy_pa, inclination, sky_pa))
+        < 1e-10
+    ), "Rotation matrix in RotateAndTiltWrapperPotential does not agree with the formula in the documentation"
+    return None
+
+
 def test_integration_RotateAndTiltWrapper():
     ## test a quick orbit integration to hit the C code (also test pure python)
     # two potentials, one offset