Merge pull request #2537 from cta-observatory/shower_axis_point

Add function to get point on shower axis in altaz

docs/changes/2537.feature.rst

@@ -0,0 +1,3 @@
+Add function ``ctapipe.coordinates.get_point_on_shower_axis``
+that computes a point on the shower axis in alt/az as seen
+from a telescope.

src/ctapipe/coordinates/__init__.py

@@ -20,7 +20,7 @@
 from .impact_distance import impact_distance, shower_impact_distance
 from .nominal_frame import NominalFrame
 from .telescope_frame import TelescopeFrame
-from .utils import altaz_to_righthanded_cartesian
+from .utils import altaz_to_righthanded_cartesian, get_point_on_shower_axis
 
 __all__ = [
     "TelescopeFrame",
@@ -35,6 +35,7 @@
     "altaz_to_righthanded_cartesian",
     "impact_distance",
     "shower_impact_distance",
+    "get_point_on_shower_axis",
 ]
 
 

src/ctapipe/coordinates/ground_frames.py

@@ -46,11 +46,15 @@ def _earthlocation_to_altaz(location, reference_location):
 
 
 class GroundFrame(BaseCoordinateFrame):
-    """Ground coordinate frame.  The ground coordinate frame is a simple
-    cartesian frame describing the 3 dimensional position of objects
-    compared to the array ground level in relation to the nomial
-    centre of the array.  Typically this frame will be used for
-    describing the position on telescopes and equipment.
+    """Ground coordinate frame.
+
+    The ground coordinate frame is a simple cartesian frame describing the
+    3 dimensional position of objects compared to the array ground level
+    in relation to the nominal center of the array.
+
+    Typically this frame will be used for describing the position of telescopes,
+    equipment and shower impact coordinates.
+
     In this frame x points north, y points west and z is meters above array center.
 
     Frame attributes: None

src/ctapipe/coordinates/tests/test_utils.py

@@ -0,0 +1,55 @@
+import astropy.units as u
+import numpy as np
+import pytest
+from astropy.coordinates import AltAz
+
+
+def test_point_on_shower_axis_far(subarray_prod5_paranal):
+    """Test for get_point_on_shower_axis"""
+    from ctapipe.coordinates import get_point_on_shower_axis
+
+    core_x = 50 * u.m
+    core_y = 100 * u.m
+    alt = 68 * u.deg
+    az = 30 * u.deg
+    # for a very large distance, should be identical to the shower direction
+    distance = 1000 * u.km
+
+    point = get_point_on_shower_axis(
+        core_x=core_x,
+        core_y=core_y,
+        alt=alt,
+        az=az,
+        telescope_position=subarray_prod5_paranal.tel_coords,
+        distance=distance,
+    )
+
+    np.testing.assert_allclose(point.alt, alt, rtol=1e-3)
+    np.testing.assert_allclose(point.az, az, rtol=1e-2)
+
+
+def test_single_telescope(subarray_prod5_paranal):
+    from ctapipe.coordinates import (
+        MissingFrameAttributeWarning,
+        get_point_on_shower_axis,
+    )
+
+    core_x = 50 * u.m
+    core_y = 100 * u.m
+    alt = 68 * u.deg
+    az = 30 * u.deg
+    distance = 10 * u.km
+
+    point = get_point_on_shower_axis(
+        core_x=core_x,
+        core_y=core_y,
+        alt=alt,
+        az=az,
+        telescope_position=subarray_prod5_paranal.tel_coords[0],
+        distance=distance,
+    )
+
+    source = AltAz(alt=alt, az=az)
+    # 10 km is around the shower maximum, should be around 1 degree from the source
+    with pytest.warns(MissingFrameAttributeWarning):
+        assert u.isclose(source.separation(point), 1.0 * u.deg, atol=0.1 * u.deg)

src/ctapipe/coordinates/utils.py

@@ -1,12 +1,21 @@
 import astropy.units as u
+import numpy as np
+from astropy.coordinates import AltAz
+from erfa.ufunc import p2s as cartesian_to_spherical
 from erfa.ufunc import s2p as spherical_to_cartesian
 
+from .ground_frames import _get_xyz
+
 __all__ = [
     "altaz_to_righthanded_cartesian",
+    "get_point_on_shower_axis",
 ]
 
 
-def altaz_to_righthanded_cartesian(alt, az):
+_zero_m = u.Quantity(0, u.m)
+
+
+def altaz_to_righthanded_cartesian(alt, az, distance=1.0):
     """Turns an alt/az coordinate into a 3d direction in a right-handed coordinate
     system.  This is because azimuths are in a left-handed system.
 
@@ -19,8 +28,55 @@ def altaz_to_righthanded_cartesian(alt, az):
     az: u.Quantity
         azimuth
     """
+    # this is an optimization, shaves of ca. 50 % compared to handing over units
+    # to the erfa function
     if hasattr(az, "unit"):
         az = az.to_value(u.rad)
         alt = alt.to_value(u.rad)
 
-    return spherical_to_cartesian(-az, alt, 1.0)
+    unit = None
+    if hasattr(distance, "unit"):
+        unit = distance.unit
+        distance = distance.value
+
+    result = spherical_to_cartesian(-az, alt, distance)
+    if unit is not None:
+        return u.Quantity(result, unit=unit, copy=False)
+    return result
+
+
+@u.quantity_input(core_x=u.m, core_y=u.m, alt=u.rad, az=u.rad, distance=u.km)
+def get_point_on_shower_axis(core_x, core_y, alt, az, telescope_position, distance):
+    """
+    Get a point on the shower axis in AltAz frame as seen by a telescope at the given position.
+
+    Parameters
+    ----------
+    core_x : u.Quantity[length]
+        Impact x-coordinate
+    core_y : u.Quantity[length]
+        Impact y-coordinate
+    alt : u.Quantity[angle]
+        Altitude of primary
+    az : u.Quantity[angle]
+        Azimuth of primary
+    telescope_position : GroundFrame
+        Telescope position
+    distance : u.Quantity[length]
+        Distance along the shower axis from the ground of the point returned.
+
+    Returns
+    -------
+    coord : AltAz
+        The AltAz coordinate of a point on the shower axis at the given distance
+        from the impact point.
+    """
+    impact = u.Quantity([core_x, core_y, _zero_m], unit=u.m)
+    # move up the shower axis by slant_distance
+    point = impact + altaz_to_righthanded_cartesian(alt=alt, az=az, distance=distance)
+
+    # offset by telescope positions and convert to sperical
+    # to get local AltAz for each telescope
+    cartesian = point[np.newaxis, :] - _get_xyz(telescope_position).T
+    lon, lat, _ = cartesian_to_spherical(cartesian)
+    return AltAz(alt=lat, az=-lon, copy=False)