clc
clear all
close all
maxNumCompThreads(1)

saved = [];
N = 2.^[0:7]*100;
for N=N
[x1, x2, v1, v2, dt, q] = setup_system(N);
tic
for t = 0:100
    %[r1, r2] = slow_loop(x1, x2);
    %[r1, r2] = fast1(x1, x2);
    [r1, r2] = fast2(x1, x2);
    
    r_squared = r1.^2 + r2.^2 + 0.0001;
    r1 = r1./r_squared;
    r2 = r2./r_squared;
    v1 = v1 - r1*q*dt;
    v2 = v2 - r2*q*dt;
    x1 = x1 + v1*dt;
    x2 = x2 + v2*dt;    
end
saved = [saved, toc]
end
p = polyfit(log(N),log(saved),1)

function [r1, r2] = slow_loop(x1, x2)
    for i = 1:N
        for j = 1:N
         r1(i,j) = x1(i) - x1(j);
         r2(i,j) = x2(i) - x2(j);
        end
    end
end

function [r1, r2] = fast1(x1, x2)
    r1 = x1' - x1;
    r2 = x2' - x2;
end

function [r1, r2] = fast2(x1, x2)
    r1 = bsxfun(@minus,x1,x1');
    r2 = bsxfun(@minus,x2,x2');
end

function [x1, x2, v1, v2, dt, q] = setup_system(N)
    dt = 0.001; % timestep
    vmag = 2*(log(N)/log(10)-1); % Magnitude of init velocity of particles
    rng(0, 'simdTwister')
    theta = rand(N, 1)*2*pi; % position particles randomly along circumference
    tt = 1./[(N-1)/4:5*(N-1)/4]'; % Init radial position of particles
    rr = cumsum((0.5/sum(tt))*tt); % Use tt to eveny spread particles
    x1 = rr.*cos(theta);
    x2 = rr.*sin(theta);
    q = rand(N, 1) + 1;  % Particle mass, we might choose to ignore this and set q=1
    q = ones(N,1);
    vr = x1.^2 + x2.^2;
    v1 = vmag*-x2./vr;
    v2 = vmag*x1./vr;
end
