function u = rotate(u0,ws)
%rotate.m Arbitrary direction vector deflection (random azimuth)
% 
% FileName      : rotate.m
% Creation      : 09/29/2016
% Modifications : September 29, 2016 Initial creation
% Author(s)     : Alex Bielajew
%
% Inputs        : u0 (vector, double): 
%                    = (u_0,v_0,w_0) initial (x,y,z)-axis direction cosines
%                    u = sin(theta_0)*cos(phi_0)
%                    v = sin(theta_0)*sin(phi_0)
%                    w = cos(theta_0)
%               : ws (scalar, double)
%                   = cos(Theta) cosine of the polar scattering angle, Theta
% Outputs
%               : u (vector, double)
%                    = (u,v,w) final (x,y,z)-axis direction cosines
%
% Purpose       : NERS 544 Monte Carlo Methods
%               : Deflects the particle direction into a random azimuthal angle 
%                    about it's  direction u0, by cos(Theta) = ws
%               : From (7:13) of the Monte Carlo book        
%               : Demonstration, production, part of the code library      

if (ws == 1)
   u = u0; % There is no scattering, return the original vector
   return
end

% General scattering case

[us,vs] = azimuthal(); % Unit normalized 2-D transverse direction cosines
                       % i.e. cos(Phi), sin(Phi)
      
% Convert 2-D transverse to 3-D transverse direction cosines
      
sinThetas = sqrt((1 - ws)*(1 + ws)); % This handles the near 1 situation for
                                     % ws with more precision than
                                     % sqrt(1 - ws^2))
us = sinThetas*us;     % i.e. us = sin(Theta)*cos(Phi)
vs = sinThetas*vs;     % i.e. vs = sin(Theta)*sin(Phi)
 
u2v2 = u0(1)^2 + u0(2)^2;    % i.e. u2v2 = sin(theta_0)^2
    
if (u2v2 == 0) % Particle's original direction was along the z-axis,
               % a common source condition
   u(1) = us;
   u(2) = vs;
   u(3) = ws;       % Bug found 10/24/2016! 
   u(3) = u0(3)*ws; % B10/24/2016 change. Now handles backscatter correctly
   return
end

sintheta  = sqrt(u2v2);       % i.e. sin(theta_0)
sinthetai = 1/sintheta;       % inverse used below
cosphi    = u0(1)*sinthetai;  % i.e. cos(phi)/sin(Theta)
sinphi    = u0(2)*sinthetai;  % i.e. sin(phi)

% Equation (7.13) from the book
u(1) = u0(1)*ws + us*u0(3)*cosphi - vs*sinphi;
u(2) = u0(2)*ws + us*u0(3)*sinphi + vs*cosphi;
u(3) = u0(3)*ws - us*sintheta;
return