function earthOrbit
clear;
clc;
hold on;
%Physical constants
m_sun = 1.989E30;           %kg
m_earth = 5.972E24;         %kg
r_earth = 149.59787e6;      %km  
u_earth = 30;               %km/s

%Initial conditions
X = [r_earth;0];
U = [0; u_earth]

%simulation parameters
steps = 50; %200 400 800
tfinal = 3*365*24*60*60; %s 
dt = tfinal/steps;

tarr = (1:steps)*dt;
uarr = zeros(2,length(tarr));
xarr = zeros(2,length(tarr));
uarr(:,1) = U;
xarr(:,1) = X;
for step = 2:steps
    uarr(:,step) = uarr(:,step-1) + dt * newton(xarr(:,step-1),m_earth,m_sun)/m_earth;
    %xarr(:,step) = xarr(:,step-1) + dt * uarr(:,step-1); %Explicit Euler
    xarr(:,step) = xarr(:,step-1) + dt * uarr(:,step-0); %Symplectic Euler
end
%scatter(xarr(1,:),xarr(2,:))
comet(xarr(1,:),xarr(2,:))
circle(0,0,r_earth);
fprintf('Done\n');

end
function circle(x,y,r)
%x and y are the coordinates of the center of the circle
%r is the radius of the circle
%0.01 is the angle step, bigger values will draw the circle faster but
%you might notice imperfections (not very smooth)
ang=0:0.01:2*pi; 
xp=r*cos(ang);
yp=r*sin(ang);
plot(x+xp,y+yp);
end
function F=newton(Xji,mi,mj)
%6.674×10?11 N?m2/kg2
    G = 6.673e-11;       %N · m2/kg2=(kg*m/s^2)*m^2/kg^2 (pick up 3 km/m scales)
    G_scale = G/((10^3)^3);
    eji = Xji/norm(Xji);
    F = -eji*(G_scale*mi*mj/(norm(Xji)^2));
end