# Matlab simulation for planet motion

4 views (last 30 days)
Niklas Kurz on 16 Mar 2022
Commented: Niklas Kurz on 17 Mar 2022
There were some attemps simulating planetary motion already, but I think mine is straightforward by solving and updating position via with Euler Cromers method:
t = 0;
while t < 10
pos1 = [1 2 3];
pos2 = [4 5 6];
m1 = 1;
m2 = 2;
G = 1;
r1 = pos1-pos2;
r2 = pos2-pos1;
F1 = G*m1*m2/norm(r1).^2.*r1/norm(r1);
F2 = G*m1*m2/norm(r2).^2.*r2/norm(r2);
dt = 0.1;
p1 = [0 100 0];
p2 = [0 100 0];
p1 = p1+F1.*dt;
p2 = p2+F2.*dt;
pos1 = pos1+p1/m1;
pos2 = pos2+p2/m2;
t = t+dt;
hold all;
plot3(pos1(1),pos1(2),pos1(3),'rx')
plot3(pos2(1),pos2(2),pos2(3),'bx')
end
However I don't really receive a plot of multiple data points, just 2 crosses remaining stationary. Also I get a 2-D plot even though I reverted to plot3
KSSV on 16 Mar 2022
You can change it to 3D using view.
plot3(pos1(1),pos1(2),pos1(3),'rx')
plot3(pos2(1),pos2(2),pos2(3),'bx')
view(3)

James Tursa on 16 Mar 2022
The initial condition for position and velocity need to be outside the loop, prior to loop entry.
Niklas Kurz on 17 Mar 2022
Thank you for pointing it out.

KSSV on 16 Mar 2022
t = 0;
m1 = 1;
m2 = 2;
G = 1;
pos01 = [1 2 3];
pos02 = [4 5 6];
pos1 = zeros([],3) ;
pos2 = zeros([],3) ;
iter = 0 ;
while t < 10
iter = iter+1 ;
r1 = pos01-pos02;
r2 = pos02-pos01;
F1 = G*m1*m2/norm(r1).^2.*r1/norm(r1);
F2 = G*m1*m2/norm(r2).^2.*r2/norm(r2);
dt = 0.1;
p1 = [0 100 0];
p2 = [0 100 0];
p1 = p1+F1.*dt;
p2 = p2+F2.*dt;
pos1(iter,:) = pos01+p1/m1;
pos2(iter,:) = pos02+p2/m2;
pos01 = pos1(iter,:) ;
pos02 = pos2(iter,:) ;
t = t+dt;
end
figure
hold on
plot3(pos1(:,1),pos1(:,2),pos1(:,3),'rx')
plot3(pos2(:,1),pos2(:,2),pos2(:,3),'bx')
view(3) Niklas Kurz on 17 Mar 2022
This looks much better to me regarding number of points. Nevertheless there is still something weird with the coding going around. Even if it should equal the code that I assimilated from Glowscript:
G = 1
star = sphere(pos = vector(0,0,0), radius = 0.2, color = color.yellow, mass = 1000, momentum = vector(0,0 ,0), make_trail=True)
plan = sphere(pos = vector(1,0,0), radius = 0.5, color = color.blue , mass = 1 , momentum = vector(0,30,0), make_trail=True)
while (True):
rate(500)
r_star = star.pos - plan.pos
r_plan = plan.pos - star.pos
star.force = -G*star.mass*plan.mass/(mag(r_star)**2)*(r_star)/mag(r_star)
plan.force = -G*star.mass*plan.mass/(mag(r_plan)**2)*(r_plan)/mag(r_plan)
star.momentum = star.momentum + star.force * dt
plan.momentum = plan.momentum + plan.force * dt
star.pos = star.pos + star.momentum/star.mass * dt
plan.pos = plan.pos + plan.momentum/plan.mass * dt
t = t+dt