use quaternions

ctapipe/coordinates/ground_frames.py

@@ -108,9 +108,21 @@ def _get_shower_trans_matrix(azimuth, zenith, inverse=False):
     -------
     trans: 3x3 ndarray transformation matrix
     """
-    rot = Rotation.from_euler(
-        "zy", np.stack((azimuth.to_value(u.rad), -zenith.to_value(u.rad)), -1)
-    )
+    # We construct a quaternion for this rotation by multiplying the quarternion of 'z' rotation [0,0,sin(theta/2),cos(theta/2)] with the quarternion of 'y' rotation [0,sin(theta/2),0,cos(theta/2)].
+    # See Hamilton product at https://en.wikipedia.org/wiki/Quaternion#Multiplication_of_basis_elements
+    # By doing this we can see that 12 of 16 terms vanish. The other 4 construct the quaternion vector used here.
+
+    zen2 = -zenith / 2  # minus sign due to corsika coords
+    az2 = azimuth / 2
+    s_zen = np.sin(zen2)
+    c_zen = np.cos(zen2)
+    s_az = np.sin(az2)
+    c_az = np.cos(az2)
+
+    # this is the hamilton product of ('y'*'z')
+    quat = np.array([s_zen * s_az, s_zen * c_az, c_zen * s_az, c_zen * c_az]).T
+
+    rot = Rotation(quat)
     if inverse:
         rot = rot.inv()
     return rot.as_matrix()
@@ -150,8 +162,9 @@ def ground_to_tilted(ground_coord, tilted_frame):
     """
     xyz_grd = _get_xyz(ground_coord)
 
-    zenith = tilted_frame.pointing_direction.zen
-    azimuth = tilted_frame.pointing_direction.az
+    # convert to rad first and substract. Faster than .zen
+    zenith = np.pi / 2 - tilted_frame.pointing_direction.alt.to_value(u.rad)
+    azimuth = tilted_frame.pointing_direction.az.to_value(u.rad)
 
     rotation_matrix = _get_shower_trans_matrix(azimuth, zenith)
 
@@ -180,8 +193,8 @@ def tilted_to_ground(tilted_coord, ground_frame):
     """
     xyz_tilt = _get_xyz(tilted_coord)
 
-    zenith = tilted_coord.pointing_direction.zen
-    azimuth = tilted_coord.pointing_direction.az
+    zenith = np.pi / 2 - tilted_coord.pointing_direction.alt.to_value(u.rad)
+    azimuth = tilted_coord.pointing_direction.az.to_value(u.rad)
 
     rotation_matrix = _get_shower_trans_matrix(azimuth, zenith, inverse=True)
 
@@ -219,8 +232,8 @@ def project_to_ground(tilt_system):
     z_initial = ground_system.z.value
 
     trans = _get_shower_trans_matrix(
-        tilt_system.pointing_direction.az,
-        tilt_system.pointing_direction.alt,
+        tilt_system.pointing_direction.az.to_value(u.rad),
+        np.pi / 2 - tilt_system.pointing_direction.alt.to_value(u.rad),
     )
 
     x_projected = x_initial - trans[2][0] * z_initial / trans[2][2]