% Quick single particle gravity simulation
function elliptical_orbit

%Orbital parameters
a = 5.202; %Semi-major axis of orbit in AU
T = 11.86; %Orbital period in years
M = 1.0; %Mass of star in solar masses
k = 0.6; %k = 1 means circular orbit

dt = 0.005; %Timestep in years. See 'precession error' if dt = 0.05
num_periods = 7; %Number of orbital periods
fsize = 18; %Fontsize

%

%Determine max simulation time /Yr from Kepler III
tmax = num_periods*T;

%Run solver. Star is the first object, planet is the second
x0 = [0,a]; y0 = [0,0]; z0 = [0,0]; 
vy0 = [ 0,k*2*pi*sqrt( M/a ) ]; vx0 = [0,0]; vz0 = [0,0];
[t,x,y,z,vx,vy,vz] = gravity( dt,tmax,x0,y0,z0,vx0,vy0,vz0,[M,0] );
x = x(:,2); y = y(:,2); vx = vx(:,2); vy = vy(:,2);

%Plot trajectory
figure('name','gravity','numbertitle','off','color',[1 1 1],...
    'units','normalized','position',[0 0 1 1],'renderer','painters',...
    'units','normalized',...
    'position',[0.05 0.05 0.85 0.85] );
plot(x,y,'r'); hold on; plot( 0,0,'r*','markersize',80); grid on;
plot( 0,0,'ro','markersize',50,'markerfacecolor','r' );
p = plot(x(1),y(1),'bo','markersize',20,'markerfacecolor','b' );
set(gca,'fontsize',fsize); axis equal; 
title({'Star and planet',...
    ['M=',num2str(M,3),' solar masses, T=',num2str(T,3),' years, a=',...
    num2str(a,3),'AU, k=',num2str(k)]});
xlabel('x /AU'); ylabel('y /AU');
stop = 0; n = 1; N = length(x);
while stop==0  %Animate tajectory
    n = n+1; if n>N; n = 1; end
    try
        set( p, 'xdata',x(n),'ydata',y(n) ); drawnow;
    catch
        return
    end
end

%End of code


