MATLAB模拟日地月系统。

Question

怎么在这部分代码加入判断月亮处于阴影和光照的情况，并且画出月亮的轨迹呢？
代码如下，注释部分是我的设想但是运行不出，求解

n=1000;
t=linspace(0,2*pi,n);
hold on;
[x,y,z]=sphere;
k=0;
view(3);
axis equal; 

axis([-17 17 -17 17 -10 10]);
title('Sunsystem') 

theta=[pi/2.3 pi/3];
r=[15 2];
w=[1 12]; 

trace=nan+[t;t;t];

for j=t 

k=k+1;
cla; %清除坐标区

surf(5*x,5*y,5*z,'AmbientStrength',1);%画半径5的球体：太阳
shading interp %实现对不平滑的图像进行平滑

i=1; 

T=[sin(theta(i)),0,cos(theta(i));0,1,0;-cos(theta(i)) 0 sin(theta(i))]; 

O=r(i)*T*[cos(t);sin(t);zeros(1,n)]; %放置地球轨迹线

o=r(i)*T*[cos(j*w(i));sin(j*w(i));0]; 

plot3(O(1,:),O(2,:),O(3,:),':'); %轨迹线

surf(x+o(1),y+o(2),z+o(3),'FaceLighting','gouraud','AmbientStrength',.5); %画大球体

i=2; 

T=[sin(theta(i)),0,cos(theta(i));0,1,0;-cos(theta(i)) 0 sin(theta(i))]; 

P=repmat(o,1,n)+r(i)*T*[cos(t);sin(t);zeros(1,n)]; %放置卫星轨迹线

p=o+r(i)*T*[cos(j*w(i));sin(j*w(i));0];

surf(x/i+p(1),y/i+p(2),z/i+p(3),'FaceLighting','gouraud','AmbientStrength',.5); %小球体
plot3(P(1,:),P(2,:),P(3,:),':'); %轨迹线

 % M=(O(2,:))/(O(1,:));
% m=(1+M^2)^(1/2);
% plot([O(1,:)+m,m],[O(2,:)+m,m])
% plot([O(1,:)-m,-m],[O(2,:)-m,-m])
% a=atan(M)+pi/2;
% b=a+pi;
% x1=O(1,:)+cos(a+(pi/2));
% y1=O(2,:)+sin(a+(pi/2));
% x2=2*O(1,:)-x1;
% y2=2*O(2,:)-y1;
% N=(P(2,:)-O(2,:))/(P(1,:)-O(1,:));
% c=atan(N);
trace(:,k)=p;
% if  (c>=a)&(c<=b)

plot3(trace(1,:),trace(2,:),trace(3,:),'r');

% end

shading interp;
light('position',[0 0 0],'style','local'); %设置光源

drawnow

end 

万分感谢

Answer 1

地球月球和太阳的公转可以用程序模拟一下，具体你可以看matlab的用法