Earth orbiting around the sun using Euler

I am not fluent in matlab, nor in physics, but I am trying to make a simple simulation of the earth orbiting around the sun, though it doesn't seem to work. I want to use Euler to calculate the coordinates in 2-D. I thought the points would end up in a more or less circular plot, but I am just getting a straight line... I'm using nasa to get the data, and I want the sun to be in origo. After doing some reading online, I still can't seem to figure what I am doing wrong..

Here is what I've got so far:

%bodies = {'sun','earth','moon'};
mass = [1.989e8 5.972e2 7.349].*1e22;
 
% The positions (x,y) 
xx = [0 -7.461557189835504E-01];
yy = [0 6.451970837257494E-01];
 
% The velocities (vx,vy) 
ux= [0 -1.154054959303684E-02];
uy= [0 -1.307896092314608E-02];
 
G = 6.67e-11;
%M = mass(1)*mass(2);
M = mass(1);
x0 = xx(2); y0 = yy(2); vx0 = ux(2); vy0=uy(2); 
 
t0 = 0;
 
x(1)=x0;
y(1)=y0;
vx(1)=vx0;
vy(1)=vy0;
t(1) = t0;
 
fvx = @(x,y) -(G*M*x/((x^2+y^2)^(3/2)));
fvy = @(x,y) -(G*M*y/((x^2+y^2)^(3/2)));

steps = 1000;
dt = 1/steps;
 
for k=1:steps-1
    
    vx(k+1) = vx(k) + fvx(x(k),y(k))*dt;
    vy (k+1) = vy(k) + fvy(x(k),y(k))*dt; 
    
    x(k+1) = x(k)+vx(k+1)*dt;
    y(k+1) = y(k)+vy(k+1)*dt;
    
    t(k+1) = t(k) + dt;
    
    %plot(x(k+1),y(1+k),'-b');
    %drawnow;
    
end
 
figure
hold on
plot(x,y,'r')
drawnow;

Answer 1

[Walter Roberson]

The time unit for G is s^-2

Your maximum time is (steps-1)/steps which is going to be just less than 1. 1 what? Well since you do not apply any scaling to G we have to conclude "1 second"

Therefore you have carefully plotted orbits in thousands of a second for one second, and you are expecting visible curvature.

How far is your Earth from your Sun? Well the unit of G involves m^3 and your first object is 7e-1 away from the second, so we conclude that your simulated Earth is about 70 cm away from the center of mass of the Sun.

Unfortunately, if you concentrate the mass of the sun into that small of a radius then you would have a black hole and the radius of the black hole would be between 5 and 6 kilometres. Thus if you want to plot the orbit of the Earth 75 cm from the center of mass of the Sun, you either have to model conditions inside a black hole or else you have to reduce the effective mass of the Sun as much of it would be "above" the Earth and so pulling up on the Earth reducing the effective mass.

Answer 2

[Signe Helene Tveitaskog Vesterøy]

Aha, okay, I wasn't aware that my units didn't match! I've added the AU-data (which I wasn't aware of that I was missing) and kept the G as it was. If my calculations are correct I've fixed all the values to be in the SI-units, as I would prefer my orbit not to be in a black hole. Now my method seems to be orbiting around the sun, hurray!

I've also changed my timesteps, so I think this is more correct? :)

I am still wondering, why will it not plot inside the loop? It is'nt very essential, I'm fine with plotting the vectors x and y outside the loop, but why is'nt it working with plotting the x- and y-coordinates inside directly?

%bodies = {'sun','earth','moon'};
mass = [1.989e8 5.972e2 7.349].*1e22; %kg

AU = 149597870700; %meters
%The positions (x,y) 
xx = [0 -7.461557189835504E-01 ].*AU; %m
yy = [0 6.451970837257494E-01 ].*AU;  %m
 
%The velocities (vx,vy) %m/s
ux= [0 -1.154054959303684E-02].*AU/86400; %converting from m/day to m/s
uy= [0 -1.307896092314608E-02].*AU/86400; 
 
G = 6.67e-11; %m^3/(kg*s^2)
M = mass(1);
x0 = xx(2); y0 = yy(2); vx0 = ux(2); vy0=uy(2); 
 
t0 = 0;
 
x(1)=x0;
y(1)=y0;
vx(1)=vx0;
vy(1)=vy0;
t(1) = t0;

fvx = @(x,y) -(G*M*x/((x^2+y^2)^(3/2)));
fvy = @(x,y) -(G*M*y/((x^2+y^2)^(3/2)));
 
T = 20*365.25*24*60*60; % the time domain: 0 --> tf %10 years in seconds
%T = 365.25*24*60*60;
N = 500;
dt = T/N;

%plot(x(1),y(1),'r'); %first point
%axis([-2 2 -2 2].*1e11);
%hold on 

for k=1:N

    %vx(k+1) = vx(k)-((G*M)/((sqrt(x(k)^2+y(k)^2))^3))*x(k)*dt;
    vx(k+1) = vx(k) + fvx(x(k),y(k))*dt;
    
    %vy(k+1) = vy(k)-((G*M)/((sqrt(x(k)^2+y(k)^2))^3))*y(k)*dt;
    vy (k+1) = vy(k) + fvy(x(k),y(k))*dt; 
    
    x(k+1) = x(k)+vx(k+1)*dt;
    y(k+1) = y(k)+vy(k+1)*dt;
    
    t(k+1) = t(k) + dt;

    %plot(x(k+1),y(k+1),'r');
    %hold on
    %drawnow;

end
 
figure
plot(0,0,'o','LineWidth',25,'markersize',5,'color','k'); %sun
hold on
plot(x,y,'r') %earth orbiting around the sun using euler

[Walter Roberson]

    %plot(x(k+1),y(k+1),'r');

That would only plot one point at a time, and plot() only draws lines between plot points if there are at least two consecutive finite endpoints. If you had specified a marker such as 'r*-' then you would have seen the markers.

See also animatedline()