Simulating earth rotating the sun with eulers method.

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

Answer 1

[William Rose]

@Roble Noor,

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).

[Roble Noor]

As the first speed is the speed of the Earth in the direction of the Sun, not the speed in the direction of the orbit, it should be right i think? I found the formulas from this link https://evgenii.com/blog/earth-orbit-simulation/

[William Rose]

@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.

Answer 2

[Jan]

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

[Roble Noor]

Yeah, i changed the the value t for some testing, but when i increase it to simulate an actual year, it becomes a linear model. When running the code underneath(which i have done only small changes to) it becomes linear.

%% 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=50; %timestep



t=[0:h:3.154e7];

distance=[];
speed=[];


angle=[];
angular_velocity=[];


x=[];
y=[];



%Initial conditions for speed
init_distance=1.496e11;
init_speed=29.78e3;

%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(x,y);




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

[William Rose]

@Roble Noor,

Here's the output with a 3 day long time step, simulated for 1 year.

[Roble Noor]

Perfect! Thank you so much! really appreciated

[William Rose]

@Roble Noor, you're welcome! Good luck with your work.