clf;
clear all;
%%%%%%%%%%% Equations of interest
%%% m1_new = m1 - dm2
%%% m1*v1 + dm2*v2 = m1_new*v1new

%%%%%%%%% Set Global Constants
global G
G = 6.67*(10^-11); % m^3/kg*s^2
%%%%%%%%% Set masses
ms1 = 1.3*(10^22); % kg, mass of pluto
ms2 = 3*(10^21); % kg, mass of charon
mass1(1) = ms1;
mass2(1) = ms2;
%%%%%%%%% Initial positions
rs1 = [20;0];
rs2 = [-10;0];
r_cm = (rs1*ms1 + rs2*ms2)/(ms1+ms2);

%%% Shift everything by center of mass
rs1 = rs1 - r_cm;
rs2 = rs2 - r_cm;

%%%%%%%%% Initial velocities
mu_12 = (ms1 * ms2)/(ms1 + ms2);
R12 = [rs1(1)-rs2(1);rs1(2)-rs2(2)];
omega = sqrt(((G*ms1*ms2)/(mu_12*(norm(R12)^3))));
vs1 = flipud(rs1.*omega);
vs1(1) = -vs1(1);
vs2 = flipud(rs2.*omega);
vs2(1) = -vs2(1);

%%%%%%%%% Set y vector
y = [rs1;vs1;rs2;vs2];
%%% Find Period
T = 2*pi*sqrt(((norm(rs1))^3)/(G*ms2));

%%%%%%%%%%%%%%%%%%% Runga Kutta Method %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dt = T/100;
t = 0;
tmax = T*50;

dms1 = ms2/(tmax/dt); % mass transfer rate
dms2 = -dms1;

i = 1;

figure(1)
xlabel('X axis')
ylabel('Y axis')
title('Binary Star System, Center of Mass Frame')

while t < tmax
    
    R = norm([(y(1)-y(5));(y(2)-y(6))]);
    E = .5*((ms1*ms2)/(ms1+ms2))*((norm([y(7);y(8)])).^2) - ((G*ms1*ms2)/R)
    
    tt(i) = t;
    yt(i,:) = y';
    
    theta1 = atan(y(2)/y(1));
    theta2 = atan(y(6)/y(5));

%     if (theta1 == theta2)||((theta1 == 0)&&(theta2 ==(2*pi)))||((theta2 == 0)&&(theta1 ==(2*pi)))
        ms1_new = ms1 + dms1;
        ms2_new = ms2 + dms2;
%         y(3) = ((ms1*y(3))+(dms1*y(7)))/ms1_new;
%         y(4) = ((ms1*y(4))+(dms1*y(8)))/ms1_new;
        y(7) = ((ms2*y(7))-(dms2*y(3)))/ms2_new;
        y(8) = ((ms2*y(8))-(dms2*y(4)))/ms2_new;
        ms1 = ms1_new;
        ms2 = ms2_new;
%         theta1
%         theta2
%         pause
%     end

    mass1(i+1) = ms1;
    mass2(i+1) = ms2;

    F1 = kepler_ode(t,ms1,ms2,y);
    F2 = kepler_ode(t+((dt/2)),ms1,ms2,y+((dt/2).*F1));
    F3 = kepler_ode(t+(dt/2),ms1,ms2,y+((dt/2).*F2));
    F4 = kepler_ode((t+dt),ms1,ms2,y+(dt.*F3));

    y = y + (F1 + F2.*2 + F3.*2 + F4).*(dt/6);

    i = i + 1;
    t = t + dt;

end
x1_cm = yt(:,1);
y1_cm = yt(:,2);
x2_cm = yt(:,5);
y2_cm = yt(:,6);

for i = 1:2:length(x1_cm)
    hold
    plot(x1_cm(i),y1_cm(i),'o','MarkerSize',ceil(50*((mass1(i))/(10^22))),'MarkerFaceColor','b','MarkerEdgeColor','b')
    plot(x2_cm(i),y2_cm(i),'o','MarkerSize',ceil(50*((mass2(i))/(10^22))),'MarkerFaceColor','r','MarkerEdgeColor','r')
    plot(x1_cm,y1_cm)
    plot(x2_cm,y2_cm,'r')
    saveas(gcf,strcat('test_movie_',num2str(999+((i+1)/2)),'.jpg'));
    fprintf('Frame %4d of %4d\n',j,20);
    pause(0.01)
end
% %%% Change reference frame to normal (Combine r' = r - r_cm for stars 1&2)
% x1 = ((ms1*(x1_cm-x2_cm))+((ms1+ms2)*x2_cm))/ms2;
% y1 = ((ms1*(y1_cm-y2_cm))+((ms1+ms2)*y2_cm))/ms2;
% x2 = (((ms1*(x1_cm-x2_cm))+((ms1+(2*ms2))*x2_cm))/ms2)-x1_cm;
% y2 = (((ms1*(y1_cm-y2_cm))+((ms1+(2*ms2))*y2_cm))/ms2)-y1_cm;
% figure(2)
% hold on
% xlabel('X axis')
% ylabel('Y axis')
% title('Binary Star System, Normal Frame')
% plot(x1,y1)
% plot(x2,y2,'r')
