Simulate a solar system

68 visualizzazioni (ultimi 30 giorni)
Jesus De Leon
Jesus De Leon il 4 Ago 2021
Risposto: Pavan Guntha il 16 Ago 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

Risposte (1)

Pavan Guntha
Pavan Guntha il 16 Ago 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!

Categorie

Scopri di più su Earth and Planetary Science in Help Center e File Exchange

Prodotti


Release

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by