
% kepler.m solves Kepler problem
% using Runge-Kutta 4nd order
%
% Equation:
%
% m d^2r/dt^2 - m*r*(dphi/dt)^2 + dV(r)/dr = 0; V(r) = -a/r
%   d^2phi/dt^2 + 2/r*(dr/dt)*(dphi/dt) = 0

% Initial condition:    r(0)   =r0, dr/dt = 0
%                       phi(0) = 0, dphi/dt(0) = phid0
%
% Example:
%
% [t,x,y] = kepler(1,1,0.01,10*628);
% 
% Note: (x,y) are Cartesian coordinates
%

function [t,x,y] = kepler(r0,phid0,dt,nSteps)

% Initialize vectors

t = zeros(1,nSteps);

x1 = zeros(1,nSteps); % saves r
x2 = zeros(1,nSteps); % saves dr/dt
x3 = zeros(1,nSteps); % saves phi
x4 = zeros(1,nSteps); % saves dphi/dt

tRun = 0;
xRun = [r0,0,0,phid0];

fdt = 1.0/6.0*dt;

for j=1:nSteps
	
        % compute time step using 4-th order Runge-Kutta

        k1 = rkFunc(tRun,xRun);
	    k2 = rkFunc(tRun+0.5*dt,xRun+0.5*dt*k1);
        k3 = rkFunc(tRun+0.5*dt,xRun+0.5*dt*k2);
        k4 = rkFunc(tRun+dt,xRun+dt*k3);
	    
        xRun = xRun + fdt*(k1+2*k2+2*k3+k4);
        tRun = tRun + dt;

        % save computed results

        t(j) = tRun;
	    x1(j) = xRun(1);
        x2(j) = xRun(2);
        x3(j) = xRun(3);
        x4(j) = xRun(4);
end

% x1 contains r(t)   and x2 contains dr/dt
% x3 contains phi(t) and x4 contains dphi/dt

% transformation to Cartesian coordinates

x = x1.*cos(x3);
y = x1.*sin(x3); 

% DONE
end

% System of ODEs 

function xPrime = rkFunc(t,x)

% Parameters:
%
% m: reduced mass
% a: constant in V(r) = -a/r

m =  % Your Code here
a =  % Your Code here
ep = % Your Code here

xPrime = zeros(1,length(x));

xPrime(1) =  % Your Code here
xPrime(2) =  % Your Code here
xPrime(3) =  % Your Code here
xPrime(4) =  % Your Code here

end
