# Simulating earth rotating the sun with eulers method.

7 views (last 30 days)
Roble Noor on 9 Dec 2021
Commented: William Rose on 10 Dec 2021
Hi! I’m trying to make a simulation of earths rotation around the sun, and for some weird reason my angle does not change. I’m using a formula I found from the internet, but it seems to work for that person. I think its something wrong with my constants, because I just get a linear movement when I plot it.
%% Initial values
%Constants and initialising vectors
G=6.674*10^(-11); %gravitational constant
M=1.989e30; %Mass of sun
m=5.927e24; %Mass of earth
h=1; %timestep
t=[0:h:12];
%%t=[0:h:3.154e7];
distance=[];
speed=[];
angle=[];
angular_velocity=[];
x=[];
y=[];
%Initial conditions for speed
init_distance=1.496e11;
init_speed=0;
%Intial conditions for angle
init_angular_velocity=1.99e-7;
%% Settings first values
distance(1)=init_distance;
speed(1)=init_speed;
angle(1)=0;
angular_velocity(1)=init_angular_velocity;
dvdt=0;
dsdt=0;
x(1)=distance(1)*sin(angle(1));
y(1)=distance(1)*cos(angle(1));
%% Eulers Method
for i=2:length(t) %eulers method
dvdt=P(G,M,distance(i-1),angular_velocity(i-1));
speed(i)=speed(i-1)+dvdt*h;
distance(i)=distance(i-1)+speed(i)*h;
dsdt=S(angular_velocity(i-1),speed(i-1),distance(i-1));
angular_velocity(i)=angular_velocity(i-1)+dsdt*h;
angle(i)=angle(i-1)+angular_velocity(i)*h;
x(i)=distance(i)*sin(angle(i));
y(i)=distance(i)*cos(angle(i));
end
plot(0,0,'o','LineWidth',25,'markersize',5,'color','k'); %sun
hold on
%%plot(radius,t,'r') %earth orbiting around the sun using euler
plot(y,x);
end
end

William Rose on 9 Dec 2021
init_speed=0;
%...
init_angular_velocity=1.99e-7;
init_speed=0 is wrong. The Earth will be pulled straight toward the Sun if it has no initial tangential velocity. initial_angular_velocity has the correct value, in radians per second. Therefore your initial conditions are inconsistent.
The units for
end
do not make sense. On the left hand side, the units of velocity shoudl be m/s. The first term on the right hand side has units of m/(s^2), and the second term on the RHS has units of m^2/(s^2).
The units for
end
do not make sense. On the left hand side, the units for angle shoudl be non-dimensional. The right hand side has units of 1/(s^2).
##### 2 CommentsShowHide 1 older comment
William Rose on 9 Dec 2021
@Roble Noor, what units do you intend for h? This is another example where comments in your code are a good idea. Since G and initial_angular_velocity are specified with time units of seconds, h should also be in seconds. This means your simulation runs for 12 seconds. Did you mean 12 months? For 12 seconds, Earth's path is not discernably different from a straight line.
You define
but when you call the function, you call it by
dvdt=P(G,M,distance(i-1),angular_velocity(i-1));
which makes me wonder if P() is supposed to calculate velocity or acceleration. Right now it calculates neither, since the right hand side terms have inconsistent units.
x(i)=distance(i)*sin(angle(i));
y(i)=distance(i)*cos(angle(i));
which shows that distance equals distance from the origin (where the sun is). In the formulas above, you reverse the usual convention: x=r*cos(a), y=r*sin(a) is standard, but you have them reversed. That would not be a big problem if your convention were applied consistently throughout.
Your code also has the equations
speed(i)=speed(i-1)+dvdt*h;
distance(i)=distance(i-1)+speed(i)*h;
These are OK for one dimensional motion, but not for motion in two dimensions, which you are simulating. In uniform circular motion, speed is constant, but distance from center does not change. Your formula does not allow this to happen.
Overall, it seems like your equations are a mixed and inconsisent combinaiton of equations for straight line motion and circular motion. I recommend that you do it all in Cartesian coordinates. In other words, write the equaitons for the x- and y- components of acceleration, and the x- and y- components of velocity, and use those to update the x- and y- position, by Euler's method. If you want to do it in polar coordinates, i.e. using angle and angular velocity and angular acceleration, go for it, but do it right, which is tricky.

Jan on 9 Dec 2021
When I run your code, I do not see a linear movement:
%Constants and initialising vectors
G=6.674*10^(-11); %gravitational constant
M=1.989e30; %Mass of sun
m=5.927e24; %Mass of earth
h=1; %timestep
t=0:h:12;
%%t=[0:h:3.154e7];
distance=[];
speed=[];
angle=[];
angular_velocity=[];
x=[];
y=[];
%Initial conditions for speed
init_distance=1.496e11;
init_speed=0;
%Intial conditions for angle
init_angular_velocity=1.99e-7;
%% Settings first values
distance(1)=init_distance;
speed(1)=init_speed;
angle(1)=0;
angular_velocity(1)=init_angular_velocity;
dvdt=0;
dsdt=0;
x(1)=distance(1)*sin(angle(1));
y(1)=distance(1)*cos(angle(1));
%% Eulers Method
for i=2:length(t) %eulers method
dvdt=P(G,M,distance(i-1),angular_velocity(i-1));
speed(i)=speed(i-1)+dvdt*h;
distance(i)=distance(i-1)+speed(i)*h;
dsdt=S(angular_velocity(i-1),speed(i-1),distance(i-1));
angular_velocity(i)=angular_velocity(i-1)+dsdt*h;
angle(i)=angle(i-1)+angular_velocity(i)*h;
x(i)=distance(i)*sin(angle(i));
y(i)=distance(i)*cos(angle(i));
end
plot(0,0,'o','LineWidth',25,'markersize',5,'color','y'); %sun
hold on
plot(y,x); end
end
William Rose on 10 Dec 2021
@Roble Noor, you're welcome! Good luck with your work.

R2021a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!