clf
clear

% set constants: masses, gravity

global ms1 ms2 mp1 G

ms1 = 10^7;
ms2 = 10^6;
mp1 = 10^1;
G = 1;
dt = 19/100*10^-3;
t=10;

% set initial positions for sun, planet, asteroid

R=20;
rs1 = [0,0];
rs2 = [20,0];

rcm = (ms1.*rs1 + ms2.*rs2) ./ (ms1+ms2);

rs1 = rs1 - rcm;
rs2 = rs2 - rcm;
rp1 = [0,80];

% find omega (w), the angular velocity of the objects

mu = ms1*ms2/(ms1+ms2);

w = sqrt(G*ms1*ms2/(mu*R^3));

% set initial velocities for sun, planet, asteroid

vs1 = -w.*[0,-rs1(1)];
vs2 = w.*[0,rs2(1)];
vp1 = [-sqrt(G*(ms1+ms2)/80),0];

% compile matrix y

y = [ rs1(1) ; rs1(2) ; vs1(1) ; vs1(2) ; rs2(1) ; rs2(2) ; vs2(1) ; vs2(2) ; rp1(1) ; rp1(2) ; vp1(1) ; vp1(2) ];

% iterate

set(gca,'nextplot','replacechildren');
aviobj = avifile('real99','quality',100,'compression','Cinepak')

for i = 1:8000
    dydt = binaryandplanet_kepler(t,y(:,i));
    
    F1 = binaryandplanet_kepler(t,y(:,i));
    F2 = binaryandplanet_kepler(t+dt/2,y(:,i)+dt/2.*F1);
    F3 = binaryandplanet_kepler(t+dt/2,y(:,i)+dt/2.*F2);
    F4 = binaryandplanet_kepler(t+dt,y(:,i)+dt.*F3);
    
    y(:,i+1) = y(:,i) + dt/6.*(F1+2.*F2+2.*F3+F4);
    
    rs1x(i) = y(1,i);
    rs1y(i) = y(2,i);
    rs2x(i) = y(5,i);
    rs2y(i) = y(6,i);
    rp1x(i) = y(9,i);
    rp1y(i) = y(10,i);

    plot(rs1x(max(i-1,1):i-1),rs1y(max(i-1,1):i-1),'og',rs2x(max(i-1,1):i-1),rs2y(max(i-1,1):i-1),'ob','MarkerSize',20)
    hold on
    plot (rp1x(max(i-1,1):i-1),rp1y(max(i-1,1):i-1),'or','MarkerSize',10)
    hold on
    plot (rs1x(max(i-2000,1):i),rs1y(max(i-2000,1):i),'.g',rs2x(max(i-2000,1):i),rs2y(max(i-2000,1):i),'.b',rp1x(max(i-2000,1):i),rp1y(max(i-2000,1):i),'.r','MarkerSize',5)
    hold off
    axis ([-100 100 -100 100])
    pause (0.01)
    
    if mod(i,10) == 1
    F(i) = getframe;
    aviobj = addframe(aviobj,F(i));
    end
    
end

aviobj = close(aviobj);