Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Prob with MATLAB code (Planetary motion)

  1. Mar 28, 2010 #1
    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 (Text):
    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;
     
  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted