% Approximation of Two Body Motion
% David Bloore, August 2011

clear all;
clc;

% Universal Constants

G=3;

% case selector (each case is a set of initial conditions)

ini=input('Case?');
switch ini
    case 1
        m1=1;
        x1=-5;
        y1=-1;
        vx1=3;
        vy1=0;
        m2=4;
        x2=0;
        y2=0;
        vx2=0;
        vy2=0.1;
        t=0;
        endtime=3;
        i=10000;
    case 2
        m1=1;
        x1=-3;
        y1=-0.4;
        vx1=3;
        vy1=0;
        m2=4;
        x2=0;
        y2=0;
        vx2=-0.5;
        vy2=0;
        t=0;
        endtime=3;
        i=1000000;
    case 3
        m1=1;
        x1=-2;
        y1=0;
        vx1=3;
        vy1=1;
        m2=4;
        x2=2;
        y2=-1;
        vx2=1;
        vy2=1;
        t=0;
        endtime=7;
        i=10000000;
    case 3.2
        m1=1;
        x1=-2;
        y1=0;
        vx1=3;
        vy1=1;
        m2=4;
        x2=2;
        y2=-1;
        vx2=1;
        vy2=1;
        t=0;
        endtime=30;
        i=10000000;
    case 4
        m1=1;
        x1=-2;
        y1=0;
        vx1=3;
        vy1=1;
        m2=4;
        x2=2;
        y2=-1;
        vx2=1;
        vy2=1;
        t=0;
        endtime=100;
        i=10000000;
    case 5
        m1=1;
        x1=-20;
        y1=0;
        vx1=3;
        vy1=0;
        m2=300;
        x2=0;
        y2=0;
        vx2=0;
        vy2=3;
        t=0;
        endtime=30;
        i=1000000;
    case 6
        m1=1;
        x1=-20;
        y1=0;
        vx1=0;
        vy1=-3.5;
        m2=300;
        x2=0;
        y2=0;
        vx2=0;
        vy2=3;
        t=0;
        endtime=100;
        i=1000000;
    case 7
        m1=1;
        x1=-15;
        y1=0;
        vx1=0;
        vy1=-3;
        m2=300;
        x2=0;
        y2=0;
        vx2=0;
        vy2=0;
        t=0;
        endtime=100;
        i=1000000;
    case 7.1
        m1=1;
        x1=-20;
        y1=0;
        vx1=0;
        vy1=-3;
        m2=300;
        x2=0;
        y2=0;
        vx2=0;
        vy2=0;
        t=0;
        endtime=100;
        i=1000000;
    case 8
        m1=1;
        x1=-15;
        y1=0;
        vx1=0;
        vy1=-7.5;
        m2=300;
        x2=0;
        y2=0;
        vx2=0;
        vy2=0;
        t=0;
        endtime=100;
        i=1000000;
    case 9
        m1=1;
        x1=-15;
        y1=0;
        vx1=0;
        vy1=-2;
        m2=300;
        x2=0;
        y2=0;
        vx2=0;
        vy2=0;
        t=0;
        endtime=100;
        i=1000000;
    case 9.1
        m1=1;
        x1=-15;
        y1=0;
        vx1=0;
        vy1=-2;
        m2=300;
        x2=0;
        y2=0;
        vx2=0;
        vy2=0;
        t=0;
        endtime=100;
        i=10000000;
    otherwise
        clc;
        display('Cases available are:  1, 2, 3, 3.2, 4, 5, 6, 7, 7.1, 8, 9, and 9.1');
        break
end

% time interval

dt=(endtime)/i;

% generation of blank vectors to be filled whose elements are values at a
% specific interval

X1=zeros(1,i);
Y1=zeros(1,i);
X2=zeros(1,i);
Y2=zeros(1,i);
VX1=zeros(1,i);
VY1=zeros(1,i);
VX2=zeros(1,i);
VY2=zeros(1,i);

% assignment of initial condition to the first elements of parametrics

X1(1)=x1;
Y1(1)=y1;
X2(1)=x2;
Y2(1)=y2;
VX1(1)=vx1;
VY1(1)=vy1;
VX2(1)=vx2;
VY2(1)=vy2;

% computation loop for x1,y1,x2,y2 as functions of t

c=1;

while c<i
    
    % calculation of acceleration for interval
    
    % distance between particles, distance from one to two when x1<x2 is +
    
    dx=X2(c)-X1(c);
    dy=Y2(c)-Y1(c);
    r=(dx^2+dy^2)^(1/2);
    
    % magnitude of acceleration
    
    a1=G*m2/r^2;
    a2=G*m1/r^2;
    
    % x,y components of acceleration for each mass
    
    ax1=a1*(dx/r);
    ay1=a1*(dy/r);
    ax2=a2*(-dx/r);
    ay2=a2*(-dy/r);
    
    % calculation and assignment of new positions
    
    X1(c+1)=X1(c)+VX1(c)*dt+(1/2)*ax1*dt^2;
    Y1(c+1)=Y1(c)+VY1(c)*dt+(1/2)*ay1*dt^2;
    X2(c+1)=X2(c)+VX2(c)*dt+(1/2)*ax2*dt^2;
    Y2(c+1)=Y2(c)+VY2(c)*dt+(1/2)*ay2*dt^2;
    
    % calculation and assignment of new velocities
    
    VX1(c+1)=VX1(c)+ax1*dt;
    VY1(c+1)=VY1(c)+ay1*dt;
    VX2(c+1)=VX2(c)+ax2*dt;
    VY2(c+1)=VY2(c)+ay2*dt;
    
    % increase counter by one before beginning next loop
    
    c=c+1;
end

% position plot

subplot(1,2,1);
posit=plot(X1,Y1,X2,Y2,x1,y1,'o',x2,y2,'o');
% legend(['m1=' sprintf('%d',m1)],['m2=' sprintf('%d',m2)],'m1 x0,y0','m2 x0,y0','Location','best');
title(['Position(t) for ' sprintf('%d',t) ' \leq t \leq ' sprintf('%d',endtime) ' , with dt=' sprintf('%d',dt)]);
xlabel('x');
ylabel('y');
daspect([1 1 1]);
set(gca,'PlotBoxAspectRatio',[1 1 1]);
grid on;

% velocity plot

subplot(1,2,2);
velo=plot(VX1,VY1,VX2,VY2,vx1,vy1,'o',vx2,vy2,'o');
% legend(['m1=' sprintf('%d',m1)],['m2=' sprintf('%d',m2)],'m1 vx0,vy0','m2 vx0,vy0','Location','best');
title(['Position(t) for ' sprintf('%d',t) ' \leq t \leq ' sprintf('%d',endtime) ' , with dt=' sprintf('%d',dt)]);
xlabel('vx');
ylabel('vy');
daspect([1 1 1]);
set(gca,'PlotBoxAspectRatio',[1 1 1]);
grid on;

set(gcf,'Position',[0 356 1440 900]);

for rd=1:i/300:i
    
    % rd must be an integer, because it is the index bounding the plot inputs
    
    rd=floor(rd);
    
    % reset elements of x,y,vx,vy assigned to plot, and redraw
    
    set(posit(1),'XData',X1(1:rd));
    set(posit(1),'YData',Y1(1:rd));
    set(posit(2),'XData',X2(1:rd));
    set(posit(2),'YData',Y2(1:rd));
    set(posit(3),'XData',X1(rd));
    set(posit(3),'YData',Y1(rd));
    set(posit(4),'XData',X2(rd));
    set(posit(4),'YData',Y2(rd));
    
    set(velo(1),'XData',VX1(1:rd));
    set(velo(1),'YData',VY1(1:rd));
    set(velo(2),'XData',VX2(1:rd));
    set(velo(2),'YData',VY2(1:rd));
    set(velo(3),'XData',VX1(rd));
    set(velo(3),'YData',VY1(rd));
    set(velo(4),'XData',VX2(rd));
    set(velo(4),'YData',VY2(rd));
    
    drawnow
end