Rutherford Scattering Simulation

rutherford_noprint.m

echo off;
h = 0.02; % RK4 time step
n = 1000; % number of random test particles
maxR = 10;
function A = calc(X)
  rval = norm([X(1); X(2)]).^(-3); % r^(-3)
  W = [
  0 0 1 0;
  0 0 0 1;
  rval 0 0 0;
  0 rval 0 0;
  ];
  A = W * X;
end
function y = randomY()
  maxY = 4;
  % generates random distance from x-axis
  % according to cylindrical distribution.
  y = maxY * sqrt(rand); % inverse function of antiderivate of distribution.
end
arr = [];
for i = 1:n
  init = [-3; randomY(); 4; 0];
  vec = init;
  while( norm([vec(1); vec(2)]) < maxR )
    k1 = h * calc(vec);
    k2 = h * calc(vec + 0.5 * k1);
    k3 = h * calc(vec + 0.5 * k2);
    k4 = h * calc(vec + k3);
    vec = vec + (1/6) * (k1 + 2 * k2 + 2 * k3 + k4);
  end
  arr = [arr, atan(vec(4)/vec(3))];
end
hist(arr, 40);