How to trace the path

I have written the code below for orbits of celetial bodies. However i am unable trace the path of the bodies. How would i do that? can anyone help please

function [p, v] = solarsystemM(p, v, mass, stop_time, hide_animation)
% SOLARSYSTEM Solve the gravitational dynamics of n bodies.
%
% SOLARSYSTEM(p, v, mass, stop_time) receives the initial position, initial
% velocity, and mass, and simulate for stop_time seconds. The inputs p and
% v are N-by-2 matrices where the Nth row refers to the Nth celestial
% body. The two columns are the X and Y values of position and velocity,
% respectively. The input mass is an N-by-1 vector where the Nth item is
% the mass of the Nth celestial body.
%
% [p, v] = SOLARSYSTEM(...) returns the final position and final velocity
% of each object in the same format as the input. This will be used during
% marking to test the accuracy of your program.
%
% SOLARSYSTEM(..., hide_animation) is an optional flag to indicate whether
% the graphical animation may be hidden. This is an advanced feature for
% students who are aiming for a high level of achievement. During marking,
% when the computation speed of your program is being tested, it will be
% run with the hide_animation flag set to true. This allows your program to
% skip the drawing steps and run more quickly.
%

if nargin < 5
    hide_animation = false;
end

% Below is example code showing how to read the initial position and
% velocity into the same notation as the assignment sheet, to help you get
% started on the two-body simulation:

% How to pull values out of p and v:
% p1 = p(1,:); % planet 1 position vector
% p2 = p(2,:); % planet 2 position vector
%
% v1 = v(1,:); % planet 1 velocity vector
% v2 = v(2,:); % planet 2 velocity vector
%
% m1 = mass(1); % planet 1 mass
% m2 = mass(2); % planet 2 mass


% Write your code here

% Constants
t = 0;
dT = 4800;
G = 6.673e-11;
N = size(p,1);
a = zeros(N,2);
hist = [0:dT:stop_time];
% L = zeros((N),numel(hist));
% v = v + (a*dT);
% p = p + (v*dT);
% r12 = p2 - p1
% F12 = (-G.*m1.*m2).*r12/norm(r12)^3

% Creating vectors for different color celestial bodies
C(1,:)=[1 0 0];
C(2,:)=[0 1 0];
C(3,:)=[0 0 1];

% Caculating speed velocity and and acceleration
% a1 = 1/m2.*F12
% V = v + a1*dT
% P = p2 + (V*dT)

% Plotting of the figures

if ~hide_animation
    figure (1);
    clf;
    hold on;
    xlim([-2e11 2e11]);
    ylim([-2e11 2e11]);
    xlabel('X axis')
    ylabel('Y axis')
    title ('Orbit simulations of clestial bodies')
    
    %Add for loop for the extra celestial bodies to be plotted
    
    for i = 1:N
        
        h(i) = plot(p(i,1),p(i,2), '*','markerfacecolor',C(i,:),'markersize',5);
        
    end
end


while t < stop_time
    for i = 1:N
        a(i,:) = 0;
        for j = 1:N
            if i~=j
                
                
                
                rij = p(j,:) - p(i,:);
                Fij = G.*((mass(i).*mass(j))./(norm(rij)^3)).*rij;
                Fji = -Fij;
                a(i,:) = a(i,:) +(1/mass(i)).*Fij;
                
            end
        end
    end
    
    %Updating my velocity and position
    
    v = v + (a*dT);
    p = p + (v*dT);
    
    
    if ~hide_animation
        
        for i = 1:N
            set(h(i),'Xdata', p(i,1),'Ydata', p(i,2));
            
        end
        
        
    end
    t = t + dT;
    drawnow limitrate
    
end

% How to pack the results back into the return values
p = p(1:N,:);
v = v(1:N,:);

end

[Adam Danz]

Not enough info.

Maybe if we had example input values we could run the function and see what you mean by "unable to trace the path".

[Juan Palacios]

If you run this test function with three celestial bodies it should run my function.

function test_moderate_level(unit_under_test)
% TEST_MODERATE_LEVEL Test the simulator against the moderate level of achievement.
%
% TEST_MODERATE_LEVEL(@unit) tests a function called "unit" instead of the
% default, "solarsystem".
%
% This is provided as a means for your to test your program's accuracy. We
% supply here high precision answers that can use as a benchmark against
% which to compare your code.
%
% A program similar to this one will be used during marking to test your
% program's accuracy and speed. We'll use different initial conditions, so
% don't try simply hard-coding these answers! :)
%

% Default to a function named "solarsystem"
if nargin < 1
    unit_under_test = @solarsystemM;
end

% Physical constants

% Sun
% Data from http://nssdc.gsfc.nasa.gov/planetary/factsheet/sunfact.html
data.sun.p = [0 0]; % put the sun at the origin
data.sun.v = [0 0]; % no velocity
data.sun.mass = 1988500e24; % kg

% Earth
% Data from http://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html
%
% Explanation: when the Earth is closest to the sun ("perihelion") its
% orbital velocity is maximised. Start the simulation at perihelion
% using the data from the above webpage.
data.earth.p = [147.09e9 0];
data.earth.v = [0 30.29*1000];
data.earth.mass = 5.9723e24;

% A comet
% Data completely made up
phi = 45 * pi / 180;
R = [cos(phi) -sin(phi); sin(phi) cos(phi)]; % 2D rotation matrix
data.comet.p = (R * data.earth.p')';
data.comet.v = 0.4 * (R * data.earth.v')';
data.comet.mass = 5e21;

% Test 1
Test_2D_Sun_Earth_Comet(false);

% Test 2
Test_2D_Sun_Earth_Comet(true);

    function test_result(parameter, value, units, comparator, benchmark)
        if strcmp(units, '%')
            fprintf('  %28s :  %-15.6f', [parameter ' (' units ')'], value);
        else
            fprintf('  %28s :  %-15.6g', [parameter ' (' units ')'], value);
        end

        if nargin == 5
            if comparator(value, benchmark)
                fprintf('   ** PASS. Meets or exceeds the expectation of %g%s', benchmark, units);
            else
                fprintf('   ** FAIL. Does not meet the expectation of %g%s', benchmark, units);
            end
        end
        fprintf('\n');
    end


    function Test_2D_Sun_Earth_Comet(speed_test)
        if speed_test
            fprintf('<strong>*** [Moderate level] Three body simulation with Sun, Earth and a comet (execution speed test)</strong>\n');
        else
            fprintf('<strong>*** [Moderate level] Three body simulation with Sun, Earth and a comet</strong>\n');
        end
        
        p = [data.sun.p; data.earth.p; data.comet.p];
        v = [data.sun.v; data.earth.v; data.comet.v];
        mass = [data.sun.mass; data.earth.mass; data.comet.mass];

        % Run the program
        tic();
        [final_p, final_v] = unit_under_test(p, v, mass, 600*24*60*60, speed_test);
        t = toc();
        test_result('Execution time', t, 's');
        
        % Check the dimensions of the return values
        assert(all(size(final_p) == [3 2]), 'Expected size of return value "p" to be 3x2, received %ix%i instead. Did you implement 3 or more celestial bodies?', size(final_p, 1), size(final_p, 2));
        assert(all(size(final_v) == [3 2]), 'Expected size of return value "v" to be 3x2, received %ix%i instead. Did you implement 3 or more celestial bodies?', size(final_p, 1), size(final_p, 2));
        
        % Check the answers
        correct_p = [735284.879935314 5063046.74917656;-98070212939.0299 -115160878452.147;71150402679.2697 110734051521.905];
        correct_v = [-0.0680778317493048 0.147475054352529;22675.1355951941 -18802.003103355;-18556.0820979253 -3827.78340659555];
        
        expectations = [0.5 0.5 30]; % the comet is much harder to simulate due to its fast motion near the sun
        for i = 1:size(p,1)
            test_result(sprintf('Object %i position error', i), norm(final_p(i,:) - correct_p(i,:))/norm(correct_p(i,:))*100, '%', @le, expectations(i));
            test_result(sprintf('Object %i velocity error', i), norm(final_v(i,:) - correct_v(i,:))/norm(correct_v(i,:))*100, '%', @le, expectations(i));
        end  
    end

end

[Adam Danz]

When I run test_moderate_level(3)

*** [Moderate level] Three body simulation with Sun, Earth and a comet
Insufficient number of outputs from right hand
side of equal sign to satisfy assignment.
Error in jff/Test_2D_Sun_Earth_Comet (line 72)
        [final_p, final_v] = unit_under_test(p, v,
        mass, 600*24*60*60, speed_test);
Error in jff (line 42)
Test_2D_Sun_Earth_Comet(false); 

[Juan Palacios]

Not sure why, but that test simulation runs my code with 3 celestial bodies in orbit, im just trying to get a traced path for the orbits.

[Adam Danz]

We just need the inputs to solarsystemM(). So you can run your code to the point where solarsystemM() is called and tell us what the inputs are (or attach them as a mat file).

We've got to be able to see what you're talking about and figure out why it's not doing what you want it to do (none of that is clear yet).

[Juan Palacios]

the imputs are givem for sun earth and comet on that second code. But here is the actual file attached that you need to run against my function.

[Adam Danz]

That still throws an error.

test_moderate_level(3)

*** [Moderate level] Three body simulation with Sun, Earth and a comet
Insufficient number of outputs
from right hand side of equal
sign to satisfy assignment.
Error in
test_moderate_level/Test_2D_Sun_Earth_Comet
(line 84)
        [final_p, final_v] =
        unit_under_test(p, v,
        mass, 600*24*60*60,
        speed_test);
Error in test_moderate_level
(line 48)
Test_2D_Sun_Earth_Comet(false); 

I can't help without a clear description of the problem and seeing the relevant code that produces the problem. What does it mean " i am unable trace the path of the bodies" ?

[Juan Palacios]

do you have both the codes open when you run it?

My function simulates the orbit of two "planets" orbiting the sun. I just want to be able to add a trace to them, so it leaves a path of their orbit.

[Adam Danz]

Assuming you have the (x,y) coordinates of the orbits, why doesn't this work?

plot(x,y, '-')