Simulate a solar system

QUESTION/PROBLEM: I need to simulate multiple planets moving around my ellipse's drawn. Please use my code to solve this problem. My code only shows one planet in one ellipse, and I don't know what I'm doing wrong. I need my simulation to show multiple ellipses (each with their own planet)

clc
clear


x1 = input(" x value of sun (focus 1): ");
y1 = input("Y value of the sun(focus1): ");

np = input("How many Planets do you want to see? ");

xvals = zeros(np,1);
yvals = zeros(np,1);
mvals = zeros(np,1);
dvals = zeros(np,1);

for n=1:np
x2 = input(" x value of focus 2: ");
y2 = input(" Y value of focus 2: ");

m=input(" Enter the mass of the planet: ");
d=input(" Enter the density of the planet: ");


xvals(n) = x2;
yvals(n) = y2;
mvals(n) = m;
dvals(n) = d;

end



for f=1:np
    
    Sun(x1,y1,'yellow')
    
for i=1:np
hold on
drawellipse(x1,x2,y1,y2,m,d)
hold off
hold on
Planet(x(i),y(i),m,d,'green')
hold off
end
end


%% functions

function drawellipse(x1,x2,y1,y2,m,d)

%perihelion
Rp = sqrt ((x1)^2 + (y1)^2);
%aphelion
Ra = sqrt((x2)^2 + (y2)^2);
% eccentricity
eccentricity = (Ra-Rp)/(Ra+Rp);
numPoints = 500;
% semi-major axis
a = 0.5 * (Rp + Ra);
% c distance from the center to focus
c = eccentricity*a;
%semininor axis
b = sqrt(a^2 - c^2);



t = linspace(0, 2 * pi, numPoints); %from 0 radians to 2pi radians
X = a * cos(t);
Y = b * sin(t);

angles = atan2(y2 - y1, x2 - x1);
x = (x1 + x2) / 2 + X * cos(angles) - Y * sin(angles);
y = (y1 + y2) / 2 + X * sin(angles) + Y * cos(angles);





 for k=1:length(x)
     
     
    subplot(1, 1, 1);
    hPlot = plot(NaN,NaN);
    xlim([1.5*min(x) 1.5*max(x)]);
    ylim([1.5*min(y) 1.5*max(y)]);
    % update the data
    set(hPlot,'XData',x(1:k),'YData',y(1:k));
    Planet(x(k),y(k),m,d,'green')
    Sun(x1,y1,'yellow')
    drawnow
    % pause for 0.1 seconds
    pause(0.005);
    
    
 end
    grid on;
    axis equal
end

function Planet(s,v,m,d,color)
hold on
theta = linspace(0,2*pi);
r=((3*m)/(4*pi*d))^(1/3);
xunit = r * cos(theta) + s;
yunit = r * sin(theta) + v;
plot(polyshape(xunit,yunit),'FaceColor',color,'FaceAlpha',1.0);
hold off
end




function Sun(s,v,color)
hold on
color=color+".";
plot(s,v,color,'MarkerSize',75);
hold off
end

Answer 1

[Pavan Guntha]

Hello Jesus,

I understand that the issue faced here is that the simulation depicts only one planet at a time but the requirement is to simultaneously simulate multiple planets. One way to address this is by exploiting vectorization wherein the paths of all the planets are updated at every instant. Have a look at the following code snippet.

clc
clear

x1 = input(" x value of sun (focus 1): ");
y1 = input("Y value of the sun(focus1): ");
np = input("How many Planets do you want to see? ");

xvals = zeros(np,1);
yvals = zeros(np,1);
mvals = zeros(np,1);
dvals = zeros(np,1);

for n=1:np
x2 = input([' x value of focus ',  num2str(n), ': ']);
y2 = input([' y value of focus ',  num2str(n), ': ']);

m=input(" Enter the mass of the planet: ");
d=input(" Enter the density of the planet: ");
xvals(n) = x2;
yvals(n) = y2;
mvals(n) = m;
dvals(n) = d;
end

hold on
Sun(x1,y1,'yellow')
drawellipse(x1,xvals,y1,yvals,mvals,dvals) %passing entire vector data of all the planets at once (Vectorization).
hold off
%% functions

function drawellipse(x1,x2,y1,y2,m,d)

%perihelion
Rp = sqrt ((x1).^2 + (y1).^2);
%aphelion
Ra = sqrt((x2).^2 + (y2).^2);
% eccentricity
eccentricity = (Ra-Rp)./(Ra+Rp);
numPoints = 500;
% semi-major axis
a = 0.5 * (Rp + Ra);
% c distance from the center to focus
c = eccentricity.*a;
%semininor axis
b = sqrt(a.^2 - c.^2);

t = linspace(0, 2 * pi, numPoints); %from 0 radians to 2pi radians
X = a .* cos(t);
Y = b .* sin(t);

angles = atan2(y2 - y1, x2 - x1);
x = (x1 + x2) / 2 + X .* cos(angles) - Y .* sin(angles);
y = (y1 + y2) / 2 + X .* sin(angles) + Y .* cos(angles);

    for k=1:length(x)
        subplot(1, 1, 1);
        hPlot = plot(NaN,NaN);
        plot(x(:, 1:k),y(:, 1:k),'.b')
        xlim([1.5*min(x(:)) 1.5*max(x(:))]);
        ylim([1.5*min(y(:)) 1.5*max(y(:))]);
        % update the data
        Planet(x(:,k),y(:,k),m,d,'green')
        Sun(x1,y1,'yellow')
        drawnow
        % pause for 0.005 seconds
        pause(0.005);
    end
    grid on;
    axis equal
end

function Planet(s,v,m,d,color)
hold on
theta = linspace(0,2*pi);
r=((3*m)./(4*pi*d)).^(1/3);
xunit = r .* cos(theta) + s;
yunit = r .* sin(theta) + v;
sZ = size(xunit,1);
for i=1:sZ
    plot(polyshape(xunit(i,:),yunit(i,:),'simplify', false),'FaceColor',color,'FaceAlpha',1.0);
end
hold off
end

function Sun(s,v,color)
hold on
color=color+".";
plot(s,v,color,'MarkerSize',75);
hold off
end

You could even look at the documentation page of Vectorization for more details.

Hope this helps!