MATLAB Prob with MATLAB code (Planetary motion)

  • Thread starter cloudy14
  • Start date
Hi,

I am writing a code for planetary motion using ODE with runge-kutta.

Currently it has only two body, however, I could not get the initial values to fit in properly, the bodies will not move in circular orbital motion. They are moving towards each other.
at the v=[10 0;-10 0], whenever I try to add a value for the vy, the second body will dissapear completely from the screen. Help please!!

Code:
function orbit
global X v n G m       % Make these params. global.

n=2;    %No. of planets
G = 4*pi^2; % Grav. const
m=[4 4];  %m=[400 1.23e-2]; 
X=[0 0;50 0]; %initial positions
v=[10 0;-10 0]; %initail velocities %17.9979
state=zeros(n,4); %Empty the state memory
for k=1:n, %Forming the state matrix from the previous values
    state(k,:)=[X(k,1) X(k,2) v(k,1) v(k,2)];
end

nStep = 2000; %number of steps
tau = 0.01; %time step (yr)
time = 0;

for iStep=1:nStep
    for i=1:n
        state(i,:) = rk4(state(i,:),time,tau,@dEqs,G,i);
        %update new values of X and v
        X(i,:)=[state(i,1) state(i,2)];
        v(i,:)=[state(i,3) state(i,4)];
    end
time = time + tau;
xpos1(iStep)=X(1,1);
ypos1(iStep)=X(1,2);
xpos2(iStep)=X(2,1);
ypos2(iStep)=X(2,2);
%xpos3(iStep)=X(3,1);
%ypos3(iStep)=X(3,2);
clf('reset');
hold on;
axis on
axis equal
axis([-150 150 -150 150]);
plot(X(1,1),X(1,2),'bo','MarkerFaceColor','r','MarkerSize',12);
plot(X(2,1),X(2,2),'bo','MarkerFaceColor','y','MarkerSize',12);
%plot(X(3,1),X(3,2),'bo','MarkerFaceColor','r','MarkerSize',12);
plot(xpos1(1:iStep),ypos1(1:iStep));
plot(xpos2(1:iStep),ypos2(1:iStep));
%plot(xpos3(1:iStep),ypos3(1:iStep));

pause(0.001)


%  clf reset
%    set(gcf,'numbertitle','off','name',' n-body')
%    n = size(state,1);
%    s = max(sqrt(diag(state*state')));
%    if n <= 3, s = 2*s; end
%    axis([-s s -s s])
%    axis square

end


function F = dEqs(x,t,G,num) % First-order differential
global n m X

a = zeros(1,2); 
tiny=1e-20;
    for j=1:n,
        if j~=num;
        a(1) =G* (a(1)+ m(j)* (X(j,1)-x(1,1))/norm((X(j,1)-x(1,1))+tiny)^3 );%x component
        a(2) =G* (a(2)+ m(j)* (X(j,2)-x(1,2))/norm((X(j,2)-x(1,2))+tiny)^3 );%x component        
        end
    end

F = [x(1,3) x(1,4) a(1) a(2)]; 

% Runge-Kutta integrator (4th order)
function xout = rk4(x,t,tau,derivsRK,param,count)
half_tau = 0.5*tau;
F1 = feval(derivsRK,x,t,param,count);
t_half = t + half_tau;
xtemp = x + half_tau*F1;
F2 = feval(derivsRK,xtemp,t_half,param,count);
xtemp = x + half_tau*F2;
F3 = feval(derivsRK,xtemp,t_half,param,count);
t_full = t + tau;
xtemp = x + tau*F3;
F4 = feval(derivsRK,xtemp,t_full,param,count);
xout = x + tau/6.*(F1 + F4 + 2.*(F2+F3));
return;
 

Want to reply to this thread?

"Prob with MATLAB code (Planetary motion)" You must log in or register to reply here.

Physics Forums Values

We Value Quality
• Topics based on mainstream science
• Proper English grammar and spelling
We Value Civility
• Positive and compassionate attitudes
• Patience while debating
We Value Productivity
• Disciplined to remain on-topic
• Recognition of own weaknesses
• Solo and co-op problem solving
Top