function [r] = new_direction_3D(p,gamma,phi)
% angles in radians 

% initial direction vector
p_n = p./norm(p);
gamma = gamma/180*pi;
phi = phi/180*pi;

% trigonometric method
phi_rot = atan2(p_n(2),p_n(1));
theta_rot = acos(p_n(3));

% phi_rot = atan(p_n(2)/p_n(1));
% if p_n(1)<0 % taking into account the atan domain
%     phi_rot = pi+phi_rot;
% end

if theta_rot == 0 | theta_rot == pi % case for p||z => no need to rotate around z
    Rz = eye(3); % identity matrix
else
    Rz = [cos(phi_rot) -sin(phi_rot) 0; % rotation matrix around the z-axis
        sin(phi_rot) cos(phi_rot) 0;
        0 0 1];
end

Ry = [cos(theta_rot) 0 sin(theta_rot); % rotation matrix around the y-axis 
    0 1 0;
    -sin(theta_rot) 0 cos(theta_rot)];

r_scat = [sin(gamma)*cos(phi); sin(gamma)*sin(phi); cos(gamma)]; % in the rotated system
% r = r_scat'*Ry^(-1)*Rz^(-1); % the resulting direction in the original system
% r = r_scat'*(Rz*Ry)^(-1); %faster
r =((Rz*Ry)*r_scat)'; % trick without inversion! much faster 

%% equivalent coordinate form
%{
d = sqrt(p_n(2)^2+p_n(3)^2);

if d == 0 % no need to rotate
    Rx = eye(3); % identity matrix
else
    Rx = [1 0 0;
        0 p_n(3)/d p_n(2)/d;
        0 -p_n(2)/d p_n(3)/d];
end

Ry = [d 0 p_n(1); 
    0 1 0; 
    -p_n(1) 0 d];

% test p = Rx*Ry*z
Rx*Ry*z'
%}

end

