function [p_reflect, s_reflect] = getPlaneReflection(p_in, s_in, p_plane, n_plane)
% getPlaneReflection -
%
% Syntax: [p_reflect, s_reflect] = getPlaneReflection(p_in, s_in, p_plane, n_plane)
%
% Inputs:
%    - p_in, s_in, p_plane, n_plane -
%
% Outputs:
%    - p_reflect, s_reflect -

p_reflect = p_in + s_in*dot(p_plane - p_in, n_plane)/dot(s_in, n_plane);
s_reflect = s_in-2*n_plane*dot(s_in, n_plane);
s_reflect = s_reflect./norm(s_reflect);

end