MATLAB Answers

Simulate a solar system

9 views (last 30 days)
Jesus De Leon
Jesus De Leon on 4 Aug 2021
Answered: Pavan Guntha on 16 Aug 2021
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

Answers (1)

Pavan Guntha
Pavan Guntha on 16 Aug 2021
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!

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!