Simulating earth rotating the sun with eulers method.

7 views (last 30 days)
Roble Noor
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_angle=deg2rad(23.5);
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);
function vel = P(Gravity,Mass,radius,ang_velocity)
vel=(radius*(ang_velocity^2) - (Gravity*Mass)/radius);
end
function angle= S(ang_velocity,velocity,radius)
angle=-((2*velocity*ang_velocity)/radius);
end

Accepted Answer

William Rose
William Rose on 9 Dec 2021
In your code,
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.
Add comments to your code shwowing units, for every constant you define, and for formulas. This will help you find errors in your formulas, such as the following.
The units for
function vel = P(Gravity,Mass,radius,ang_velocity)
vel=(radius*(ang_velocity^2) - (Gravity*Mass)/radius);
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
function angle= S(ang_velocity,velocity,radius)
angle=-((2*velocity*ang_velocity)/radius);
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 Comments
William Rose
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
function vel = P(Gravity,Mass,radius,ang_velocity)
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.
Your code has equations
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.

Sign in to comment.

More Answers (1)

Jan
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_angle=deg2rad(23.5);
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);
function vel = P(Gravity,Mass,radius,ang_velocity)
vel=(radius*(ang_velocity^2) - (Gravity*Mass)/radius);
end
function angle= S(ang_velocity,velocity,radius)
angle=-((2*velocity*ang_velocity)/radius);
end
  4 Comments

Sign in to comment.

Tags

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!