在这个项目中,我们通过计算每个机构在太阳能系统上的系统的其他机构施加的引力。根据这些力量和身体的初始位置/速度,人们可以用牛顿第二定律预测他们的动作。我们假设以下信息给出:
- 参与了太阳系所有机构的质量。
- 在给定时间所有行星的位置和速度。
然后我们可以用牛顿万有引力定律来计算引力。每个身体都受到来自所有其他身体的重力的作用,因此身体上的总力量是所有其他身体所产生的力量的总和。一旦知道身体上的力量,就可以用牛顿的第二运动定律来确定每个身体的加速度。基于任何时刻t的加速度,我们可以计算出时间t +Δt处物体的新速度和位置。我们从时间t = 0处的已知位置和速度开始,然后在Δt处计算那些,然后在2Δt处等。
ASSIGNMENT:
编写一个MATLAB函数,接受作为输入一个结构阵列(其中每一个元素代表一个主体在太阳系),一个时间步长Δt和最后时间T.每个结构在该结构数组应具有以下网络的视场:
- 名称:持有星球的名称的字符串。
- x:行星初始位置的x坐标。
- y:行星初始位置的y坐标。
- z:行星初始位置的z坐标。
- vx:行星初始速度的x分量。
- vy:行星初始速度的y分量。
- vz:行星初始速度的z分量。
那么你的函数应该使用上述计算在任何时候的位置和所有行星的速度引力系统的模型。你的函数应该以适当的方式返回这些值,并产生一个三维图中所有行星位置的图表。情节应该至少包含一个标题和一个图例。图例应该使用现场名称来命名每个星球。此外,每个行星都应该以不同的颜色显示,其颜色是随机确定的。图1给出了一个示例图。
我的错误
有一个错误在我的最后一个循环,在这时候,有n=1:T/t;
关于矩阵的尺寸。我想利用每次迭代更新的位置来计算新的力量,但是我正在和错误,我的图显示了一堆的直线,而不是太阳系的轨道。**
function [x, y, z, vx, vy, vz] = solarsystemsim(F,t,T)
m = size(F,2);
G = 6.67384e-11;
j=1;
x = 1:10;
% r = zeros(m-1,1);
% gF = zeros(m-1,1);
for(i=1:m)
for(j=1:m)
r(i,j) = distform2(F(i).x,F(j).x,F(i).y,F(j).y,F(i).z,F(j).z);
gF(i,j) = (G*(F(i).mass)*(F(j).mass))/((r(i,j))^2);
gFx(i,j) = -((gF(i,j))*(F(i).x-F(j).x))/(r(i,j));
gFy(i,j) = -((gF(i,j))*(F(i).y-F(j).y))/(r(i,j));
gFz(i,j) = -((gF(i,j))*(F(i).z-F(j).z))/(r(i,j));
if(i==j)
gF(i,j) = 0;
gFx(i,j) = 0;
gFy(i,j) = 0;
gFz(i,j) = 0;
end
end
end
for(i=1:m)
gFxT(i) = sum(gFx(i,:));
gFyT(i) = sum(gFy(i,:));
gFzT(i) = sum(gFz(i,:));
end
tn = 0;
n = 1;
x = zeros(T/t,m);
y = zeros(T/t,m);
z = zeros(T/t,m);
vx = zeros(T/t,m);
vy = zeros(T/t,m);
vz = zeros(T/t,m);
for(i=1:m)
vx(1,i) = F(i).vx;
vy(1,i) = F(i).vy;
vz(1,i) = F(i).vz;
x(1,i) = F(i).x;
y(1,i) = F(i).y;
z(1,i) = F(i).z;
end
for(n=1:T/t)
for(i=1:m)
for(j=1:m)
r(i+1,j+1) = distform2(x(i,j),x(i,j),y(i,j),y(i,j),z(i,j),z(i,j));
gF(i+1,j+1) = (G*(F(i).mass)*(F(j).mass))/((r(i,j))^2);
gFx(i+1,j+1) = -((gF(i,j))*(x(i,j)-x(i,j+1)))/(r(i,j));
gFy(i+1,j+1) = -((gF(i,j))*(y(i,j)-y(i,j+1)))/(r(i,j));
gFz(i+1,j+1) = -((gF(i,j))*(z(i,j)-z(i,j+1)))/(r(i,j));
if(i==j)
gF(i,j) = 0;
gFx(i,j) = 0;
gFy(i,j) = 0;
gFz(i,j) = 0;
end
end
end
gFxT(i) = sum(gFx(i,:));
gFyT(i) = sum(gFy(i,:));
gFzT(i) = sum(gFz(i,:));
for(i=1:T/t)
for(j=1:m)
vx(i,j) = vx(i,j) + t*(gFxT(j))/(F(j).mass);
vy(i,j) = vy(i,j) + t*(gFyT(j))/(F(j).mass);
vz(i,j) = vz(i,j) + t*(gFzT(j))/(F(j).mass);
x(i,j) = x(i,j) + t*(vx(j));
y(i,j) = y(i,j) + t*(vy(j));
z(i,j) = z(i,j) + t*(vz(j));
end
end
plot3(x,y,z);
end
_Eppur SI muove_ – buzjwa