import numpy as np
# copied from https://github.com/cta-observatory/cta-lstchain/blob/v0.10.7/lstchain/reco/disp.py#L16
# to avoid dependency on entire lstchain
[docs]
def disp(cog_x, cog_y, src_x, src_y, hillas_psi):
"""
Compute the disp parameters
Parameters
----------
cog_x: `numpy.ndarray` or float
cog_y: `numpy.ndarray` or float
src_x: `numpy.ndarray` or float
src_y: `numpy.ndarray` or float
hillas_psi: `numpy.ndarray` or float
Returns
-------
(disp_dx, disp_dy, disp_norm, disp_angle, disp_sign):
disp_dx: 'astropy.units.m`
disp_dy: 'astropy.units.m`
disp_norm: 'astropy.units.m`
disp_angle: 'astropy.units.rad`
disp_sign: `numpy.ndarray`
"""
disp_dx = src_x - cog_x
disp_dy = src_y - cog_y
disp_norm = disp_dx * np.cos(hillas_psi) + disp_dy * np.sin(hillas_psi)
disp_sign = np.sign(disp_norm)
disp_norm = np.abs(disp_norm)
# disp_sign : indicates in which direction, "positive" or "negative", we must move along the
# reconstructed image axis (with direction defined by the versor cos(hillas_psi), sin(hillas_psi))
# we must move from cog_x, cog_y to get closest to the true direction (src_x, src_y)
if hasattr(disp_dx, "__len__"):
disp_angle = np.arctan(disp_dy / disp_dx)
disp_angle[disp_dx == 0] = np.pi / 2.0 * np.sign(disp_dy[disp_dx == 0])
else:
if disp_dx == 0:
disp_angle = np.pi / 2.0 * np.sign(disp_dy)
else:
disp_angle = np.arctan(disp_dy / disp_dx)
return disp_dx, disp_dy, disp_norm, disp_angle, disp_sign