clc
clear all
close all
maxNumCompThreads(1)

% The argument to the function specifies the number N of particles in the
% system
[x1, x2, v1, v2, dt, q] = setup_system(100);
for t = 0:100
    %[r1, r2] = slow_loop(x1, x2);
    %[r1, r2] = fast1(x1, x2);
    [r1, r2] = fast2(x1, x2);
    
    % ?1 and ?2 are the x and y vector components respectively
    %
    % r: distance between target particle vs other particle
    % m: mass of target particle
    % M: mass of other particle
    % F: force on target particle
    % a: acceleration of target particle
    % Gravity in 2-D looks like F = m*M/r (where we've implicitly set G=1)
    % Acceleration of our particle is F=ma -> a=F/m
    % So substituting gravitational force into acceleration eqn.
    % a = M/r
    %
    % Here we make use of the fact that 1/r = r/r^2 to get the vector
    % components of each particle's position
    r_squared = r1.^2 + r2.^2 + 0.0001;
    r1 = r1./r_squared;
    r2 = r2./r_squared;
    
    % Compute new velocity based on old velocity plus change in velocity
    % after the force of all other particles acts on a particle for a short
    % time period delta-t or "dt" (each particles mass is in "q")
    v1 = v1 - r1*q*dt;
    v2 = v2 - r2*q*dt;

    % Our new position is our old position plus how far we've gone after
    % "dt" time travelling at a particular velocity
    x1 = x1 + v1*dt;
    x2 = x2 + v2*dt;
    
    if 1
        scatter(x1,x2,'.')
        axis([-1 1 -1 1])
        axis square
        drawnow
        if t==0; pause(10); end
    end
end

% Loop over each source and then over each target particle
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

% Answers below, no cheating :)
% =
% ==
% ===
% ====
% =====
% ======
% =======
% ========
% =========
% ========
% =======
% ======
% =====
% ====
% ===
% ==
% =

% We shape the source and target vectors such that they
% are implicitly expanded to calc all combinations
function [r1, r2] = fast1(x1, x2)
    r1 = x1' - x1;
    r2 = x2' - x2;
end

% We can explicitly call for a Binary Singleton Expansion
% This is useful for more complex binary operations
function [r1, r2] = fast2(x1, x2)
    r1 = bsxfun(@minus,x1,x1');
    r2 = bsxfun(@minus,x2,x2');
end

% Code for setup that isn't relevant to the intended point
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
