% % Massive Objects Simulation
% % This code uses my collision function, a potential function
% % modified for Newton's Law of Gravity, and a ODE solver to simulate
% % many massive objects interacting.
% 
% % Henry Peterson 2016

clearvars;
close all;

% buff is the buffer for collisions used in collisions.m
% c is the constant used for calculating potential in potential.m
% L is the length of a side used in ode for dx, dy
global c L
c = 0.5;
L = 10;

%number of objects
nParts = 20;
%small time step
dt = 1e-2;
%mass of the planets
m = 1e-5;
%mass of the stars
M = 1e2;

%initialize coords and vels
x = zeros(1, nParts);
y = zeros(1, nParts);
vx = zeros(1, nParts);
vy = zeros(1, nParts);

% initial conditions - positions of stars
x(1) = L/2;
y(1) = L/4;
x(2) = L/2;
y(2) = 3*L/4;

%set up the masses 
ms = m*ones(1, nParts);
ms(1:2) = M;

%initial conditions of the planets - orbiting the starts
for i=3:nParts
    r = rand()*L/6 + L/15;
    t = rand()*2*pi;
    if i <= nParts/2
        V = potential(ms(i),ms(1),r);
        v0 = sqrt(V/ms(i));
        x(i) = r*cos(t) + x(1);
        y(i) = r*sin(t) + y(1);
        vx(i) = -v0*sin(t);
        vy(i) = v0*cos(t);
    else
        V = potential(ms(i),ms(2),r);
        v0 = sqrt(V/ms(i));
        x(i) = r*cos(t) + x(2);
        y(i) = r*sin(t) + y(2);
        vx(i) = -v0*sin(t);
        vy(i) = v0*cos(t);
    end
end

%how fast their systems are moving to the left/right
v_rel = 3;
vx(1) = v_rel;
vx(2) = -v_rel;
vx(3:nParts/2) = vx(3:nParts/2) + v_rel;
vx((nParts/2+1):end) = vx((nParts/2+1):end) - v_rel;

%keep track of frames when making a gif
frame = 1;

while 1
    
    %catch collisions
    [vx,vy] = collisions(x,y,vx,vy,ms);
    
    %solve ode
    A = [x; y; vx; vy];
    F1 = ode(A, ms);
    F2 = ode(A + F1*dt/2, ms);
    F3 = ode(A + F2*dt/2, ms);
    F4 = ode(A + F3*dt, ms);    
    A = A + dt/6*(F1 + 2*F2 + 2*F3 + F4);
    
    x = A(1,:);
    y = A(2,:);
    vx = A(3,:);
    vy = A(4,:);
    
    
    %plot the results
    color1 = [1,1,0];
    color2 = [0.5,1,0.75];
    color3 = [1,0,1];
    
    plot(x(3:nParts/2),y(3:nParts/2),'.',...
        'MarkerSize',8,'MarkerFaceColor',color2,'MarkerEdgeColor',color2);
    hold on
    plot(x((nParts/2+1):end),y((nParts/2+1):end),'.',...
        'MarkerSize',8,'MarkerFaceColor',color3,'MarkerEdgeColor',color3);
    p = plot(x(1:2),y(1:2),'.',...
        'MarkerSize',20,'MarkerFaceColor',color1,'MarkerEdgeColor',color1);
    hold off
    
    axis([-L/2, 3/2*L, -L/2, 3/2*L])
    set(gca,'Color','black','YTick',[],'XTick',[])
    set(gcf,'Color','white','InvertHardCopy','off')
    title('Massive Objects - Gravitational Forces')
    legend({'Planets of Star 1','Planets of Star 2','Stars'},'Location','southeast','Color','white')
    
    pause(dt);
    
    %if you wish to make a gif, use this code for saving images:
    print(gcf,['mass_',num2str(frame,'%04d'),'.png'],'-dpng','-r300')
    frame = frame + 1;
    %status for user
    if mod(frame,10) == 0
        fprintf('frame: %d\n',frame);
    end
    if frame >= 300
        break;
    end
    
    %if user closed the plot, break
    if ishandle(p) == 0
        break;
    end
    
end

%thanks for checking out my code!