%Code from Leathers et al., 2004, NRL report
% function Updating Direction Cosines
function [alp,bet,gam] = chgdir(alpha, beta, gamma, gammas, phis)
% Determines the new direction cosines (alp,bet,gam) from old
% direction cosines (alpha, beta, gamma), polar scattering cosine
% (gammas), and scattering azimuthal angle (phis).
%
% Robert A Leathers & Trijntje Downes, 17 March 2000.
stheta = sqrt(1-gamma*gamma); % sin(theta) of previous photon direction
sthetas = sqrt(1-gammas*gammas); % sin(theta) of scattering angle
alphas = sthetas*cos(phis); % x-direction cosine of scattering
betas = sthetas*sin(phis); % y--direction cosine of scattering
% new photon directions...
if (stheta > 1e-12) % if initial direction not straight up or straight down
B = [alpha*gamma/stheta, -beta/stheta, alpha;...
beta*gamma/stheta, alpha/stheta, beta;...
-stheta, 0, gamma];
[alphas;betas;gammas];
newdir = B*[alphas;betas;gammas];
alp = newdir(1);
bet = newdir(2);
gam = newdir(3);
else % initial direction straight up or down (sin(theta)=0)
% then new direction cosines = +/- scattering direction cosines
s = sign(gamma);
alp = s*alphas;
gam = s*gammas;
bet = s*betas;
end