% Simulates the circular orbit of a planet around a massive star
function circular_orbit_verlet

%Orbital parameters
a = 9.576; %Semi-major axis of orbit in AU
M = 2.0; %Mass of star in solar masses

%Verlet trajectory parameters
num_periods = 7; %Number of orbital periods for trajectory calculation
dt = 0.02; %Timestep in years
fsize = 18; %Fontsize

%

%Determine period in years via Kepler III
T = sqrt( (1/M)*a^3 );

%Set initial x,y,vx,vy,t parameters
x = a; y = 0; vx = 0; vy = 2*pi*sqrt( M/a ); t = 0;

%Determine initial accelerations from Newton II
[ax,ay] = gravity( x,y, 0, 0, M );

%Now determine orbit using Verlet method
n=1;
while t(end) < num_periods * T
    t(n+1) = t(n) + dt;
    x(n+1) = x(n) + vx(n)*dt + 0.5*ax*dt^2;
    y(n+1) = y(n) + vy(n)*dt + 0.5*ay*dt^2;
    ax_old = ax; ay_old = ay;
    [ax,ay] = gravity( x(n+1),y(n+1), 0, 0, M );
    vx(n+1) = vx(n) + 0.5*dt*( ax_old + ax );
    vy(n+1) = vy(n) + 0.5*dt*( ay_old + ay );
    n = n + 1;
end

%Plot trajectory
figure('color',[1 1 1],'name','star and planet','units','normalized',...
    'position',[0.2,0.2, 0.5, 0.5],'renderer','painters');
plot(x,y,'b'); hold on
plot( 0,0,'r*','markersize',80);
plot( 0,0,'ro','markersize',50,'markerfacecolor','r' );
p = plot(x(1),y(1),'bo','markersize',20,'markerfacecolor','b' );
axis equal; axis manual;
grid on; set(gca,'fontsize',fsize);
xlabel('x /AU'); ylabel('y /AU');
title({'Star and planet',...
    ['M=',num2str(M,3),' solar masses, T=',num2str(T,3),' years, a=',...
    num2str(a,3),'AU.']})
print( gcf,'circular orbit.png','-dpng','-r300');

%Animation
stop = 0; N = length(x); n=1;
while stop == 0
    %Update plot
    try
        set( p, 'xdata', x(n), 'ydata', y(n) ); drawnow; n = n+1;
    catch
        return
    end
    %Reset once N iterations have elapsed
    if n>N; n=1; end
end

%%

%Newton's law of gravitation. x,y accelerations of a mass at location x,y
%resulting from the force of gravity of mass M1 at location x1,y1.
%Units of distance are AU and time are years.
function [ax,ay] = gravity( x,y, x1, y1, M1 )
ax = -( 4*pi^2 )*M1*(x-x1)/ ( ( (x-x1)^2 + (y-y1)^2 )^(3/2) ) ;
ay = -( 4*pi^2 )*M1*(y-y1)/ ( ( (x-x1)^2 + (y-y1)^2 )^(3/2) ) ;

%End of code


