天体.m

% 作者：IC设计者笔记
% 链接：https://www.zhihu.com/question/301995846/answer/3012429972
% 来源：知乎
% 著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

clc;clear;
h1=figure('name','太阳系模型','position',[0 25 2000 1200],'color',[0 0 0]);
xlabel('x');ylabel('y');zlabel('z');

[x,y,z]=sphere;
hold on;
axis equal;
set(gca,'xlim',[-13 13],'ylim',[-13 13],'zlim',[-8 8]);%设置坐标系范围

%设置球体大小
r0=2;r1=0.8;r2=1;r3=0.3;r4=0.1;%设置星球半径
x0=r0*x;y0=r0*y;z0=r0*z;%太阳
x01=r1*x;y01=r1*y;z01=r1*z;%金星
x02=r2*x;y02=r2*y;z02=r2*z;%地球
x03=r3*x;y03=r3*y;z03=r3*z;%月球
x04=r4*x;y04=r4*y;z04=r4*z;%金星一号

w=0:0.01*pi:2*pi;
pausetime=0.002;%暂停时间
t=0;uint64(t);%记录时间的参数
Q=pi/4;%月球旋转平面的倾斜角
step=0;%飞船太空旅行所经阶段
stop_time=0;%飞船停止时间
stop_n=100;%飞船到达金星后在进行表面停留stop_n个时间单位
rad_s=0.01*pi;%发射信号所占用角度
signal_rad=0;%信号传播方向角
signal_n=40;%信号时域压缩比
signal_q=0;%信号发射累积角度
signal_q1=0.02*pi;%信号采样间隔
signal_speed=0.1;%信号传播速度
signal_speed_n=int32(r2/signal_q1);%信号传至地球后消失
signal_long=0;%信号距地球距离
mark_4_m=0;%飞船围绕月球旋转标记
moon_s=1;%飞船飞向月球期间，月球转过的圈数
begin_moon=0;%飞船飞向月球标志
s_x=[];s_y=[];s_z=[];%信号采样点数组
tt=0;

long1=6;%太阳和金星间的距离
long2=10;%太阳和地球间的距离
long3=2;%月亮和地球间距
long4=1.5;%飞船和金星的距离
long4_m=0.6;%%飞船和月亮的距离
ship_long=long2-r2;%飞船起始点距太阳中心的距离
ship_long0=long4+long1;%飞船圆周运动时距太阳中心的距离；
ship_long_temp=long2-long1-r1;%飞船与地球的距离

speed1=0.004*pi;%金星运行速度
speed2=0.002*pi; %地球运行速度
speed3=0.03*pi;%月亮运行角速度
ship_speed1=0.02;%飞船直线飞行速度
ship_speed2=speed2;%飞船运行初始角速度
ship_speed_m=0.08*pi;%飞船围绕月亮运行速度
speed_add=0.0003*pi;%控制飞船的加速系数
ship_speed2_max=0.03*pi;%飞船运行最大角速度
ship_n=5;ship_m=5;%飞船围绕金星旋转 ship_n 周后着陆；在金星表面考察ship_m周后返回；
ship_drop=0.01;%飞船着陆速度
ship_Q=0;%飞船围绕月亮运行时的倾斜角
rad4_3=0;%飞船绕月时转过的角度
lamp=0;lamp_n=8;%灯光标记,时间t每隔n次灯光变化一次

ship_rad=pi/4;%地球转至此弧度处时发射飞船
rad1=0-(speed1-speed2)*(ship_rad/speed2+(ship_long-ship_long0)/ship_speed1);
%金星起始的弧度,使飞船刚好可以围绕金星运动
rad2=0;%地球起始弧度
rad3=0;%月球起始弧度
rad4=0;%飞船起始弧度
rad4_z=0;%飞船z向公转起始弧度
rad4_3=0;%飞船围绕月球旋转过的角度


x1=x01+long1*cos(rad1);
y1=y01+long1*sin(rad1);z1=z01;%金星所在处的坐标
%text(0,long1,1.5,'金星');
x2=x02+long2*cos(rad2);
y2=y02+long2*sin(rad2);z2=z02;%地球所在处的坐标
%text(0,long2,1.5,'地球');

x3=x03+long2*cos(rad2)+long3*cos(rad3)*sin(Q);
y3=y03+long2*sin(rad2)+long3*sin(rad3);
z3=z03+long3*cos(rad3)*sin(Q);%月球所在处的坐标

col0=randn(size(z0));%随机取色
hand2=surf(x2,y2,z2,col0);%显示地球
hand1=surf(x1,y1,z1);%显示金星
set(hand1,'facecolor',[0.3 0.6 0.3],'edgecolor',[0.4 0.5 0.4]);
shading interp;%设置颜色为过渡
%get(hand1)
hand3=surf(x3,y3,z3);%显示月球
set(hand3,'facecolor',[0.2 0.2 0.2],'edgecolor','none','facelighting','none');
hand0=surf(2*x0,2*y0,2*z0);%显示太阳
%text(x0,y0,z0,'太阳');
set(hand0,'facecolor','w','edgecolor','w','facelighting','none','markersize',5,'MarkerFaceColor','w');%设置太阳颜色为红色
light('color','w','style','local','position',[0 0 0]);%在太阳中心添加光源

%绘制月球路径跟随线
line3_sx=long2*cos(rad2)+long3*cos(rad3)*sin(Q);
line3_sy=long2*sin(rad2)+long3*sin(rad3);
line3_sz=long3*cos(rad3)*sin(Q);
line3_hand_s=plot3(line3_sx,line3_sy,line3_sz);

line1_x=long1*cos(w);line1_y=long1*sin(w);line1_z=line1_x*0;
plot3(line1_x,line1_y,line1_z,'linewidth',3,'color',[0.4 0.4 0.4]);%绘制金星轨迹线
line2_x=long2*cos(w);line2_y=long2*sin(w);line2_z=line2_x*0;
plot3(line2_x,line2_y,line2_z,'linewidth',3,'color',[0.4 0.4 0.4]);%绘制地球轨迹线
line3_x=long3*cos(w)*cos(Q);line3_y=long3*sin(w);line3_z=long3*cos(w)*sin(Q);
l3_hand=plot3(line2_x,line2_y,line2_z,'linewidth',1.5,'color','b');%绘制月球轨迹线
axis off;

while 1
if rad2>=ship_rad%判断是否发射飞船
    switch step
        case 0        %初始：发射飞船
            rad4=rad2;
            x4=x04+ship_long*cos(rad4);%定位飞船初始坐标
            y4=y04+ship_long*sin(rad4 );
            z4=z04;
            hand4=surf(x4,y4,z4);
            set(hand4,'facecolor','g','edgecolor','g','facelighting','none');%设置飞船颜色
            ship_rad=rad2;%记录飞船发射时地球转过的角度
            step=step+1;
        case 1         %第一步：飞向金星
            x4=x04+ship_long*cos(rad4);%定位飞船坐标
            y4=y04+ship_long*sin(rad4 );
            rad4=rad4+speed2;%定位飞船公转角度
            ship_long=ship_long-ship_speed1;%接近金星
            if ship_long<=ship_long0%到达距金星1.5长度的位置后，进入第三步：环绕金星运动
                step=step+1;
            end;
        case 2         %第二步：围绕金星公转
            x4=x04+long1*cos(rad1)+long4*cos(rad4);%定位飞船坐标
            y4=y04+long1*sin(rad1)+long4*sin(rad4);
            z4=z04;
            if (rad4<=(ship_n+ship_m)*2*pi)||(mod(rad4-rad1,2*pi)>=(ship_speed2_max-speed1)*pi)
            %当飞船围绕金星运行ship_n+ship_m圈，且背离太阳时停止公转
                if (ship_speed2<ship_speed2_max)&&(rad4<=ship_n*2*pi)
                    ship_speed2=ship_speed2+speed_add;%控制飞船加速，加速系数speed_add
                end;
                if rad4>=ship_n*2*pi%判断飞船是否开始着陆
                    if long4>r1+r4;%判断飞船是否到达金星地面，若否，则继续降落
                        long4=long4-ship_drop;
                    end;
                    if ship_speed2>speed_add;
                    ship_speed2=ship_speed2-speed_add;%控制飞船减速
                    elseif (long4<=r1+r4)&&(stop_time<=0)
                        ship_speed2=0;
                        stop_time=t+stop_n;%飞船在金星停留stop_n个时间单位
                    end;
                    if (stop_time>0)&&( t>stop_time)&&(ship_speed2<ship_speed2_max)
                        ship_speed2=ship_speed2+2*speed_add;%控制飞船加速，加速系数speed_add
                    end;
                end;
                rad4=rad4+ship_speed2;
            elseif (ship_speed2~=0)&&(rad4_z==0)
                ship_speed2=0;
            else
                if ship_speed2<ship_speed2_max
                    ship_speed2=ship_speed2+speed_add;
                end;
                if rad4_z<2*pi*3
                    rad4_z=rad4_z+ship_speed2;
                end;
            end;
            
            if rad1-rad2>=2*pi%当金星超越地球一周时，进行第三步
                step=step+1;
            end;
        case 3          %离开金星，并向地球传播信息
            rad4=rad4+speed2;%飞船与地球同步
            ship_long=long1+r1+r4; %找到此时飞船与太阳的距离
            if rad2<5*pi/2
                while signal_q<2*pi/signal_n%发射8个周期的正弦波
                    temp_z=sin(signal_q*signal_n)*0.2;%定位波形点
                    temp_y=ship_long*sin(rad4);
                    temp_x=ship_long*cos(rad4);
                    s_z=[s_z,temp_z];s_y=[s_y,temp_y];s_x=[s_x,temp_x];%扩展数组，显示正弦信号的产生过程
                    s_x=s_x+signal_q1*cos(rad4);s_y=s_y+signal_q1*sin(rad4);%正弦波向前运动
                    if (signal_q==0)&&(signal_long<=ship_long);
                        signal_hand=plot3(s_x,s_y,s_z);signal_rad=rad4;
                        set(signal_hand,'color',[1 0 0]);
                        signal_long=ship_long;
                    else
                        set(signal_hand,'XData',s_x,'YData',s_y,'ZData',s_z);signal_long=signal_long+signal_q1;
                    end;
                    signal_q=signal_q+signal_q1;
                end;
                signal_q=0;
            end;
                s_x=s_x+signal_speed*cos(rad4); 
                s_y=s_y+signal_speed*sin(rad4);
                set(signal_hand,'XData',s_x,'YData',s_y,'ZData',s_z);
                signal_long=signal_long+signal_speed;
            if signal_long>long2
                if length(s_x)>signal_speed_n
                s_x(1:signal_speed_n)=[];
                s_y(1:signal_speed_n)=[];
                s_z(1:signal_speed_n)=[];
                signal_long=signal_long-r2;
                else
                    set(signal_hand,'color',[0 0 0]);
                    step=step+1;
                end;
            end;
            x4=x04+ship_long*cos(rad4);y4=y04+ship_long*sin(rad4 );
        case 4      %飞向地球  
            if ship_long<long2
                ship_long=ship_long+ship_speed1;
            end;
            x4=x04+ship_long*cos(rad4);y4=y04+ship_long*sin(rad4 );
            rad4=rad4+speed2;
    end;
    set(hand4,'XData',x4,'YData',y4);
    if mod(t,lamp_n)==0%t每隔lamp_n次灯光变化一次
        lamp=~lamp;
    end;
    if lamp
         set(hand4,'facecolor','g','edgecolor','g','facelighting','none');
    else
         set(hand4,'facecolor','r','edgecolor','r','facelighting','none');
    end;
end;    
rad1=rad1+speed1;
rad2=rad2+speed2;
rad3=rad3+speed3;

x1=x01+long1*cos(rad1);y1=y01+long1*sin(rad1);
set(hand1,'XData',x1,'YData',y1);%移动金星位置

x2=x02+long2*cos(rad2);y2=y02+long2*sin(rad2);%地球所在处的坐标
set(hand2,'XData',x2,'YData',y2);%移动地球位置

x3=x03+long2*cos(rad2)+long3*cos(rad3)*cos(Q);
y3=y03+long2*sin(rad2)+long3*sin(rad3);
z3=z03+long3*cos(rad3)*sin(Q);%月球所在处的坐标
set(hand3,'XData',x3,'YData',y3,'ZData',z3);%移动月球位置
line3_x=long3*cos(w)*cos(Q)+long2*cos(rad2);line3_y=long3*sin(w)+long2*sin(rad2);
set(l3_hand,'XData',line3_x,'YData',line3_y,'ZData',line3_z);%绘制月球轨迹线

%扩充月亮跟随线，当线到达一定长度后停止增长
line3_sx=[line3_sx,long2*cos(rad2)+long3*cos(rad3)*cos(Q)];
line3_sy=[line3_sy,long2*sin(rad2)+long3*sin(rad3)];
line3_sz=[line3_sz,long3*cos(rad3)*sin(Q)];
if (length(line3_sx)>10*pi/speed3)
    line3_sx(1)=[];
    line3_sy(1)=[];
    line3_sz(1)=[];
end;
%set(line3_hand_s,'XData',line3_sx,'YData',line3_sy,'ZData',line3_sz,'color',[1 0 0]);

pause(pausetime);
t=t+1; 
drawnow;
end;