Matlab simulation for planet motion

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]

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)

Answer 1

[James Tursa]

The initial condition for position and velocity need to be outside the loop, prior to loop entry.

[Niklas Kurz]

Thank you for pointing it out.

Answer 2

[KSSV]

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]

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