%
%
%
%   What if some planets from out solar system were to orbit the 
%   binary star system of Alpha Centauri?
%   What conditions would it require to be stable?
%   What would it look like in various cases?
%
%
%

close all
clear all
clf

f=0;        %   for saving .png
iter=0;     %   for saving .png

solmass=1.98855e30;  % solar mass unit kg

global G mCentauriA mCentauriB mEarth mMars mJupiter mCentauri
mCentauriA    = 1.1*solmass; % kg
mCentauriB    = 0.907*solmass; % kg
mCentauri     = mCentauriA+mCentauriB;

mEarth        = 5.972e24; % kg
mMars         = 6.417e23; % kg
mJupiter      = 1.899e27; % kg   

G             = 6.67384e-11; % m^3/(kg*s)
AU  = 1.4960e11; % m

xCA = 0;
xCB = 35.6*AU;
rCA = [xCA,0];
rCB = [xCB,0];
rEarth = [AU+1e13,0];
rMars = [5*AU+2e13,0];
rJupiter = [5.2034*AU+.5e14,0];

rCMAB = (rCA*mCentauriA+rCB*mCentauriB)/(mCentauri);
rCM = (mCentauri*rCMAB+rEarth*mEarth+rMars*mMars+rJupiter*mJupiter)/(mCentauri+mEarth+mMars+mJupiter);



rCA = (rCA-rCMAB);    % distance to center of mass
rCB = (rCB-rCMAB);    % distance to center of mass
% rEarth = (rEarth-rCMAB);
% rMars = (rMars-rCMAB);
% rJupiter = (rJupiter-rCMAB);
RAB = abs(rCA-rCB); % distance from Centauri A to Centauri B


% if planets orbit center of mass (binary stars acts as one point)


t    =  0.0;
T = sqrt((4*pi^2*norm(rJupiter)^3)/(G*(mCentauri+mJupiter)));
n    = 1000; % number of steps
dt   =  T/n;

mewAB = (mCentauriA*mCentauriB)/(mCentauriA+mCentauriB);
wAB = sqrt(G*mCentauriA*mCentauriB./(mewAB*norm(RAB).^3));
vCA = wAB*rCA;
vCA(:,[1,2]) = vCA(:,[2,1]);
vCB = wAB*rCB;
vCB(:,[1,2]) = vCB(:,[2,1]);

mCentauri = mCentauriA+mCentauriB; %    Binary star combined mass


RCE = abs(rEarth-rCM);
mewCE = (mCentauri*mEarth)/(mCentauri+mEarth);
wCE = sqrt(G*mCentauri*mEarth./(mewCE*norm(rEarth).^3));
vEarth = wCE*rEarth;
vEarth(:,[1,2]) = vEarth(:,[2,1]);



mewCM = (mCentauri*mMars)/(mCentauri+mMars);
wCM = sqrt(G*mCentauri*mMars./(mewCM*norm(rMars).^3));
vMars = wCM*rMars;
vMars(:,[1,2]) = vMars(:,[2,1]);

mewCJ = (mCentauri*mJupiter)/(mCentauri+mJupiter);
wCJ = sqrt(G*mCentauri*mJupiter./(mewCJ*norm(rJupiter).^3));
vJupiter = wCJ*rJupiter;
vJupiter(:,[1,2]) = vJupiter(:,[2,1]);



Y = [rCA(1);rCA(2);vCA(1);vCA(2);rCB(1);rCB(2);vCB(1);vCB(2);rEarth(1);rEarth(2);vEarth(1);vEarth(2);rMars(1);rMars(2);vMars(1);vMars(2);rJupiter(1);rJupiter(2);vJupiter(1);vJupiter(2)];


tt = zeros(n+1,1);
Yt = zeros(n+1,size(Y,1));
tt(1)   = t;
Yt(1,:) = Y';



for i=1:n*1.3
    f1 = kepler_odef(t     , Y                );
    f2 = kepler_odef(t+dt/2, Y+f1 .*  (dt/2)  );
    f3 = kepler_odef(t+dt/2, Y+f2 .*  (dt/2)  );
    f4 = kepler_odef(t+dt  , Y+f3 .*   dt     );
   
    Y = Y + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    
 

    t = t+dt;
    tt(i+1)   = t;  % this records the time at each step
    Yt(i+1,:) = Y'; % this records the vector Y at each step --> Yt(:,1) has all x coordinates

    whitebg('k')
    
    subplot(1,2,1)
    plot(0,0,'X')   % plot center of mass
    hold on
    plot(Y(1,:),Y(2,:),'y.','MarkerSize',40)   %   Alpha Centauri A
    
    plot(Y(5,:),Y(6,:),'c.','MarkerSize',30)    %   Alpha Centauri B
    
    plot(Y(9,:),Y(10,:),'b.','MarkerSize',6)   %   Earth
    
    plot(Y(13,:),Y(14,:),'r.','MarkerSize',5)   %   Mars
       
    plot(Y(17,:),Y(18,:),'m.','MarkerSize',20)   %   Jupiter
    
    
    axis equal
    axis([-.6e14 .6e14 -.5e14 .5e14]) 
    
    hold off
    
    subplot(1,2,2)
    plot(0,0,'X')   % plot center of mass
    hold on
    plot(Y(1,:),Y(2,:),'y-')   %   Alpha Centauri A
    hold on
    plot(Y(5,:),Y(6,:),'c-')    %   Alpha Centauri B
    hold on
    plot(Y(9,:),Y(10,:),'b-')   %   Earth
    hold on
    plot(Y(13,:),Y(14,:),'r-')   %   Mars
    hold on
    plot(Y(17,:),Y(18,:),'m-')   %   Jupiter
    hold on

%     title('Kepler XY Plot ','fontsize',20);
    axis equal
     axis([-.6e14 .6e14 -.5e14 .5e14]) 
    pause(0.00001)
   
   
   
   iter=iter+1;
    if mod(iter+4,5)==0
        f=f+1;
        saveas(gcf, ['final_',sprintf('%03d',f),'.png'])
        pause(0.001)
    end

end