Matlab simulation for planet motion

4 views (last 30 days)
Niklas Kurz
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
  1 Comment
KSSV
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)

Sign in to comment.

Accepted Answer

James Tursa
James Tursa on 16 Mar 2022
The initial condition for position and velocity need to be outside the loop, prior to loop entry.

More Answers (1)

KSSV
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)
  1 Comment
Niklas Kurz
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

Sign in to comment.

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by