% ### EXtwobodyB.m ###       11.18.16 C. Bergevin

% Example code to simulate two body motion via numeric integration.
% Modified from original source code by K. Berwick (re Giordano's book; see
% his original comments at bottom of code)
% NOTE: This version of the code simply plots everything at the end

% ---
% As per Giordano sec.4.1, assumptions made are:
% o there are two masses (m_S and m_E)
% o m_S >> m_E such that its motion can be neglected
% o by virtue of expressing things w/ units of AU, the consts. G and m_S
% (which always appears together as a product) can appear simply as 4*pi^2
% when using this choice of units. Proof:
%   - m_S= 1.989 × 10^30 kg
%   - G= 6.67408 × 10-11 m^3 kg^-1 s^-2
%   --> G*m_S = 1.3275 x 10^20 m^3 s^-2
%   - 1 AU = 1.5 x 10^11 m
%   - 1 yr ~ 3.15 x 10^7 s
%   --> 1 m^3 s^-2 = (2.94 x 10^-19) AU^3 yr^-2
%   ** --> G*m_S = 1.3275 x 10^20 m^3 s^-2 ~ 39 AU^3 yr^-2 ~ 4*pi^2 AU^3 yr^-2
% o code uses a 2nd order Runge-Kutta routine (as opposed to the simpler
% Euler method)

clear
% =============================================

Np= 10000;    % # of time steps to compute {5000}
dt = 0.001;     % time step [yrs] {0.002}
x0=1;            % initial x-position of planet [AU] {1}
y0=0;            %    "    y-position     "   [AU] {0}
vx0= 0;          %    "    x velocity     "   [AU/yr] {0}
vy0= 1.3*pi;          %    "    y velocity     "   [AU/yr] {2*pi}

% =============================================

% initialize variables
t=0;            
x(1)= x0;   y(1)= y0;   
vx(1)= vx0; vy(1)= vy0;

% +++++++++++
% loop over the timesteps
for nn = 1:Np-1;

        r= sqrt(x(nn)^2+ y(nn)^2);    % determine radius
        y_dash= y(nn) +0.5* vy(nn)* dt;
        vy_dash= vy(nn)- 0.5* (4* pi^2* y(nn)* dt)/(r^3);
        % update positions and new y velocity
        y_new=y(nn)+ vy_dash* dt;
        vy_new=vy(nn)-(4* pi^2* y_dash* dt)/(r^3);
        % Compute Runge Kutta values for the x equations
        x_dash= x(nn)+ 0.5* vx(nn)* dt;
        vx_dash= vx(nn) - 0.5* (4*pi^2* x(nn)* dt)/(r^3);
        % update positions using newly calculated velocity
        x_new= x(nn)+ vx_dash* dt; 
        vx_new= vx(nn)-(4*pi^2* x_dash* dt)/(r^3); 
      
        % Update x and y velocities with new velocities
        vx(nn+1)= vx_new;
        vy(nn+1)= vy_new;
        % Update x and y with new positions
        x(nn+1)= x_new;
        y(nn+1)= y_new;
  
 
end;

% +++++++++++
% Plot the Sun at the origin
figure(1); clf;
plot(0,0,'oy','MarkerSize',30, 'MarkerFaceColor','yellow');
xlim(1.2*[min(x) max(x)]);
ylim(1.2*[min(y) max(y)]);
xlabel('x_E (AU)'); ylabel('y_E (AU)'); hold on; grid on;
title('Two-body solutions (m_S and m_E; assumes m_S>>m_E and thus S movement is negligible')
plot(x,y,'k.','MarkerSize',8);
h1=plot(x(1),y(1),'gs','LineWidth',2,'MarkerSize',8);
h2=plot(x(end),y(end),'ro','LineWidth',2,'MarkerSize',8);
legend([h1 h2],'start point','end point');


% =============================================
% [original Berwick comments below]
% Planetary orbit using second order Runge-Kutta method.
% by Kevin Berwick,
% based on 'Computational Physics' book by N Giordano and H Nakanishi
% Section 4.2 %