% Quick single particle gravity simulation
function elliptical_orbit_verlet

%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.9; %k = 1 means circular orbit

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

%

vy = k*2*pi*sqrt( M/a );  %Initial y velocity in AU per year
vx = 0; %Initial x velocity in AU per year
x = a; %Initial x position (from star) in AU
y = 0; %Initial y position (from star) in AU
t = 0; %Time in years

%Compute trajectory using Verlet method
N = ceil( num_periods*T/dt ); %Number of simulation steps
n = 1;
while n<N
    t(n+1) = t(n) + dt;   %Increment time
    r = sqrt( x(n)^2 + y(n)^2 );  %Distance from planet to star
    ax = -(4*pi^2)*x(n)/r^3 ;  %Acceleration in x direction
    ay = -(4*pi^2)*y(n)/r^3 ;  %Acceleration in y direction
    dx = vx(n)*dt +0.5*ax*(dt^2); %XUVAT between time steps in x
    dy = vy(n)*dt + 0.5*ay*(dt^2); %XUVAT between time steps in x
    x(n+1) = x(n) + dx;  %Update vector of x positions
    y(n+1) = y(n) + dy;  %Update vector of x positions
    r = sqrt( x(n+1)^2 + y(n+1)^2 );  %Distance from planet to star
    axx = -(4*pi^2)*x(n+1)/r^3 ; %Acceleration in x direction
    ayy = -(4*pi^2)*y(n+1)/r^3 ; %Acceleration in y direction
    vx(n+1) = vx(n) + 0.5*( axx+ax )*dt;  %Update vector of x velocities
    vy(n+1) = vy(n) + 0.5*( ayy+ay )*dt; %Update vector of y velocities
    n = n + 1;
end

%Plot trajectory
figure('name','gravity','numbertitle','off','color',[1 1 1],...
    'units','normalized','position',[0 0 1 1],'renderer','painters',...
    'units','normalized',...
    'position',[0.2,0.2,0.5,0.5] );
plot(x,y,'b'); 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);
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');
axis equal; stop = 0; n = 1;
print( gcf,'elliptical orbit.png','-dpng','-r300');
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


