function w = rutherford(a)

% Samples the Rutherford distribution for a given screening parameter a

% Inputs
%%%%%%%%
%  double
%     a:  Screening parameter

% Outputs
%%%%%%%%%
%  double
%     w: Direction cosine of scattering

if (a == 0) % straight ahead scattering
   w = 1;
   return
elseif (a == inf) % Isotropic scattering
   w = 1 - 2*rand;;
   return
else % General case
    w = 1 + a - a*(2 + a)./(2*rand + a);
end