%N_random_planets_about_star
% Many-body gravity sim. Andy French Oct 2019.
% All masses in solar masses, distances in AU and times in years.

function N_random_planets_about_star
global d
dt = 0.01; fsize = 16; rmax = 10; N = 300; tmax = 20; rmin = 2;
d.stop = 0; d.quit = 0; M1 = 10; M2 = 3; a = 5; d.t = 0;

%N masses orbiting a central star of mass M1
M = 0.001*ones(1,N);
theta = 2*pi*rand(1,N);
r = rmax*rand(1,N) + rmin; r = repmat( r, [3,1] );
T = ( 1/sqrt(M1+M2) )*r.^(3/2);
x = r.*[ cos(theta); sin(theta); zeros(1,N) ];
v = [ -sin(theta) ; cos(theta); zeros(1,N) ]*2*pi.*r./T;

%Stars
[x1,x2,v1,v2,T12] = kepler( M1,M2,a );
M = [M1,M2,M];
x = [x1,x2,x]; v = [v1,v2,v];
N = N + 2;

%Run simulation
figure('renderer','zbuffer','keypressfcn',@keypressfunc,...
    'units','normalized','position',[0.2,0.2,0.5,0.5]);
vidObj = VideoWriter(['gravity ',strrep(datestr(now),':',''),'.avi']); open(vidObj);
grid on; hold on; box on; set(gca,'fontsize',fsize); colours = rand(N,3);
p = scatter3( x(1,:).',x(2,:).',x(3,:).',...
    [200;10*ones(N-2,1);100],colours,'o','filled');
view(2);
axis equal; xlim([-rmax-rmin,rmax+rmin]); ylim([-rmax-rmin,rmax+rmin]);
tit = title( 't = 0.00 years','fontsize',fsize );
while ( d.t<tmax ) && d.quit==0
    if d.stop==0
        [d.t,x,v,anew,a] = verlet(x,v,d.t,dt,M);
        set( p, 'xdata', x(1,:).', 'ydata', x(2,:).', 'zdata', x(3,:).' );
        scatter3( x(1,:).',x(2,:).',x(3,:).',ones(N,1),colours,'o','filled');
        set( tit, 'string', ['t = ',num2str(d.t,3),' years'] );
        writeVideo(vidObj,getframe)
        pause(dt);
        drawnow;
    else
        pause(1);
    end
end
close(vidObj);
close(gcf);

%%

%Verlet solver
function [tnew, xnew,vnew,anew,a] = verlet(x,v,t,dt,M)
a = newtonv(x,M);
xnew = x + v*dt + 0.5*a*(dt^2);
tnew = t + dt;
anew = newtonv(xnew,M);
vnew = v + 0.5*( a + anew )*dt;

%%

%Calculate acceleration via Newtonian gravity
function a = newton( x, M )
rmin = 0.1;      %Minimum inter-mass distance /AU for gravity effect
N = length(M);   %Number of masses (must be at least two)
a = zeros(3,N);  %Initialize acceleration
for i=1:N
    for j=1:N
        if i~=j
            %Determine acceleration of mass i due to mass j
            dx = x(:,i) - x(:,j);
            r = sqrt( dx(1)^2 + dx(2)^2 + dx(3)^2 );
            if r > rmin
                aj = -( 4*pi^2 )*M(j)*dx./(r.^3);
            else
                aj = [0;0;0];
            end
            a(:,i) = a(:,i) + aj;
        end
    end
end

%%

%Calculate acceleration via Newtonian gravity - vectorized version
function a = newtonv( x, M )
%Minimum inter-mass distance /AU for gravity effect
rmin = 0.1;

%Prepare mass matrix
N = length(M); MM = zeros(1,N,N);
M = repmat(M, [N,1] ); MM(1,:,:) = M(:,:); MM = repmat( MM, [3,1,1] );

%Find inter-mass distances /AU (NxN matrix)
r = distance(x,x); rr = zeros(1,N,N); rr(1,:,:) = r(:,:); 
rr = repmat(rr,[3,1,1]);

%Find displacement matrix (3xNxN)
xi = repmat(x,[1,1,N] ); xj = zeros(3,1,N); xj(:,1,:) = x(:,:);
xj = repmat(xj,[1,N,1]);

%Calculate acceleration
aa = -( 4*pi^2 )*MM.*( xi - xj)./(rr.^3); aa( rr < rmin ) = 0;
a = sum( aa, 3 );

%%

%Fast vectorized distance calculator for 3xN matrices a,b
% (a-b)^2 = a^2 +b^2 - 2*a*b
function d = distance(a,b)
aa=sum(a.*a,1); bb=sum(b.*b,1); ab=a'*b;
d = sqrt( abs(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab) );

%%

%Kepler
% Use Kepler's laws to determine circuar orbit initial conditions (x,v)
% and period T for a pair of objects separated by a
function [x1,x2,v1,v2,T] = kepler( M1,M2,a )
T = sqrt( (1/(M1+M2))*a^3 );
x1 = zeros(3,1); v1 = zeros(3,1);
x2 = zeros(3,1); v2 = zeros(3,1);
x1(1) = -M2*a/(M1+M2); x1(2) = 0; x1(3) = 0;
x2(1) = M1*a/(M1+M2); x2(2) = 0; x2(3) = 0;
v1(1) = 0;  v1(2) = 2*pi*x1(1)/T; v1(3) = 0;
v2(1) = 0;  v2(2) = 2*pi*x2(1)/T; v2(3) = 0;

%%

% Key press handler
function keypressfunc( fig,evnt )
global d
if strcmp(get(fig,'currentkey'),'s')==1
    d.stop = 1;
elseif strcmp(get(fig,'currentkey'),'c')==1
    d.stop = 0;
elseif strcmp(get(fig,'currentkey'),'q')==1
    d.stop = 1; d.quit = 1;
elseif strcmp(get(fig,'currentkey'),'p')==1
    fname = ['gravity t=',num2str(d.t(end))];
    fname = strrep(fname,'.','p');
    print( gcf, [fname,'.png'],'-dpng','-r300');
end

%End of code