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
Hope this helps!
