diff --git a/docs/changes/2537.feature.rst b/docs/changes/2537.feature.rst new file mode 100644 index 00000000000..2804cc272a8 --- /dev/null +++ b/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. diff --git a/src/ctapipe/coordinates/__init__.py b/src/ctapipe/coordinates/__init__.py index 388dd3e43a5..f975b41e927 100644 --- a/src/ctapipe/coordinates/__init__.py +++ b/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", ] diff --git a/src/ctapipe/coordinates/ground_frames.py b/src/ctapipe/coordinates/ground_frames.py index e0b8b8c69e8..77a6e8a57b4 100644 --- a/src/ctapipe/coordinates/ground_frames.py +++ b/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 diff --git a/src/ctapipe/coordinates/tests/test_utils.py b/src/ctapipe/coordinates/tests/test_utils.py new file mode 100644 index 00000000000..1ebd419c679 --- /dev/null +++ b/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) diff --git a/src/ctapipe/coordinates/utils.py b/src/ctapipe/coordinates/utils.py index 3d46c64c73b..33c766aef9c 100644 --- a/src/ctapipe/coordinates/utils.py +++ b/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)