%random_m_and_v_orbits
% 2D gravity simulation. Point masses of random sizes and initial velocities
% are arratcted to each other using Newtonian gravity. A Verlet solver is
% used to determine subsequent motion. 
%
% Unrealistic accelerations are 'overcome' by introducing a local repulsion
% force.
%
% LAST UPDATED by Andy French January 2017

function random_m_and_v_orbits

%Number of masses
Nm = 50;

%Mass of mass #1  - which is initially stationary at (0,0) / solar masses
M1 = 2;

%Mass of mass #2  - which is at (2*R0max,-2*R0max) / solar masses
M2 = 0.5;

%Minimum and maximum mass (in solar masses)
mm = 1e-3;
Mmin = 0.9*mm;
Mmax = 1.2*mm;

%Initial radius of mass cloud /AU
R0min = 2;
R0max = 10;

%Repulsion exponent i.e. law varies as r^-re
re = 6;

%Timestep /years
dt = 0.01;

%Simulation duration /years
tmax = 20;

%Fix axis limits
xlimits = [-2*R0max,2*R0max];
ylimits = [-2*R0max,2*R0max];

%Fontsize for graphs
fsize = 18;

%Markersize for animation
msize = 5;

%

%Generate initial x and y coordinates of masses
r0 = R0min + ( R0max-R0min )*rand(1,Nm);
theta = 2*pi*rand(1,Nm);
x0 = [0,2*R0max,r0.*cos(theta)];
y0 = [0,-2*R0max,r0.*sin(theta)];

%Generate masses
M = Mmin + rand(1,Nm)*( Mmax - Mmin );
M = [M1,M2,M];
Nm = Nm + 2;

%Compute length of time computations
N = ceil(tmax/dt) + 1;

%Initialize kinematic variables
x = zeros( N,Nm );
y = zeros( N,Nm );
vx = zeros( N,Nm );
vy = zeros( N,Nm );
t = 0 : dt : (N-1)*dt;

%Circular orbit initial conditions
x(1,:) = x0;
y(1,:) = y0;
vy(1,:) = [0,2*pi*sqrt(M1)./sqrt(2*sqrt(2)*R0max),2*pi*cos(theta)*sqrt(M1)./sqrt(r0)];
vx(1,:) = [0,0,-2*pi*sin(theta)*sqrt(M1)./sqrt(r0)];

%Compute orbits via Verlet method
disp(' ')
disp(' Computing gravity simulation using Verlet method ...')
for n=1:N-1
    [ax0,ay0] = newton( x(n,:), y(n,:), M, re );
    x(n+1,:) = x(n,:) + vx(n,:)*dt + 0.5*ax0*(dt^2);
    y(n+1,:) = y(n,:) + vy(n,:)*dt + 0.5*ay0*(dt^2);
    [ax1,ay1] = newton( x(n+1,:), y(n+1,:), M, re );
    vx(n+1,:) = vx(n,:) + 0.5*dt*(  ax0 + ax1 );
    vy(n+1,:) = vy(n,:) + 0.5*dt*(  ay0 + ay1 );
    disp( [' Step ',num2str(n),' of ',num2str(N-1),' complete.'] );
end

%Save gravity data
save gravity

%Plot orbits
figure('name','verlet gravity sim','color',[1 1 1]);
plot(x,y);
xlabel('x /AU','fontsize',fsize);
ylabel('y /AU','fontsize',fsize);
title(['Gravity simulator using Verlet method: ',num2str(Nm),' masses'],'fontsize',fsize);
set(gca,'fontsize',fsize)
grid on;
axis equal
xlim(xlimits);
ylim(ylimits);
print( gcf, 'gravity.png','-dpng','-r300' );

%Animate orbits
cla;
p = plot( x(1,3:end).',y(1,3:end).','.','markersize',msize );
hold on;
p1 = plot( x(1,1),y(1,1),'r.','markersize',msize );
p2 = plot( x(1,2),y(1,2),'g.','markersize',msize );
xlim(xlimits);
ylim(ylimits);
grid on;
xlabel('x /AU','fontsize',fsize);
ylabel('y /AU','fontsize',fsize);
title('Gravity simulator using Verlet method','fontsize',fsize);
set(gca,'fontsize',fsize)
axis manual
stop=0;
while stop == 0
    for n=2:N
        set( p, 'xdata', x(n,3:end).', 'ydata', y(n,3:end).' );
        set( p1, 'xdata', x(n,1), 'ydata', y(n,1) );
        set( p2, 'xdata', x(n,2), 'ydata', y(n,2) );
        drawnow
        pause(1/100);
    end
end

%%

%Determine accelerations based upon Newton's law of Universal Gravitation
function [ax,ay] = newton( x, y, M, re )

%Determine number of masses
Nm = length(M);

%Compute an NxN matrix of mass separations d. Use the result
% | a - b |^2 = |a|^2 + |b|^2 - 2*a.b to compute this
%in an efficient vectorized fashion.
xi = repmat( x.',[1,Nm] );
yi = repmat( y.',[1,Nm] );
xj = repmat( x,[Nm,1] );
yj = repmat( y,[Nm,1] );
d = xi.^2 + yi.^2 + xj.^2 + yj.^2 - 2*( xi.*xj + yi.*yj ) ;
d = sqrt(d);

%Determine accelerations using Newtonian gravity. Make sure i=j
%'interactions' are not included in the summation.
m = repmat( M.', [1,Nm] );
ax = -4*(pi^2)*m.*( xj - xi).*( d.^-3 - d.^-(re+1) );
ax(1:Nm+1:end) = 0;
ax = sum(ax,1);
ay = -4*(pi^2)*m.*( yj - yi).*( d.^-3 - d.^-(re+1) );
ay(1:Nm+1:end) = 0;
ay = sum(ay,1);

%End of code