%solar_system_verlet
% Solar system simulator using numerical solution of Newton's Law of
% Universal Gravitation by the Verlet method.
%
% LAST UPDATED by Andy French January 2017.

function solar_system_verlet

%Timestep /years
dt = 0.01;

%Number of years for simulation
tmax = 500;

%Planets
planets = {'Mercury','Venus','Earth','Mars',...
    'Jupiter','Saturn','Uranus','Neptune'};

%Orbital semi-major axis /AU (we will ignore any eccentricity here, and initial
%condictions shall be to result in circular orbits).
a = [0.387,...
    0.723,...
    1.000,...
    1.523,...
    5.202,...
    9.576,...
    19.293,...
    30.246
    ];

%Fontsize
fsize = 16;

%

%Compute length of time computations
N = ceil(tmax/dt) + 1;

%Initialize outputs
x = zeros( N,length(planets) );
y = x;
vx = x;
vy = x;
t = 0 : dt : (N-1)*dt;

%Circular orbit initial conditions
x(1,:) = a;
vy(1,:) = 2*pi./sqrt(a);

%Compute orbits via Verlet method
for n=1:N-1
    r3 = ( x(n,:).^2 + y(n,:).^2 ).^(3/2);
    x(n+1,:) = x(n,:) + vx(n,:)*dt - 0.5*(4*(pi^2)*(dt^2))*x(n,:)./r3;
    y(n+1,:) = y(n,:) + vy(n,:)*dt - 0.5*(4*(pi^2)*(dt^2))*y(n,:)./r3;
    rr3 = ( x(n+1,:).^2 + y(n+1,:).^2 ).^(3/2);
    vx(n+1,:) = vx(n,:) - 0.5*(4*pi^2)*dt*( x(n,:)./r3 + x(n+1,:)./rr3 ) ;
    vy(n+1,:) = vy(n,:) - 0.5*(4*pi^2)*dt*( y(n,:)./r3 + y(n+1,:)./rr3 )  ;   
end

%Plot orbits
figure('name','verlet orbit solar system','color',[1 1 1]);
plot(x,y);
legend(planets);
xlabel('x /AU','fontsize',fsize);
ylabel('y /AU','fontsize',fsize);
title('Solar system simulator using Verlet method','fontsize',fsize);
set(gca,'fontsize',fsize)
grid on;
axis equal
print( gcf, 'solar_system.png','-dpng','-r300' );

%End of code