MATLAB三维绘图和相机图形术语玩到6起来会怎样？

一个摸鱼的上班日，一台电脑，加上MATLAB和一些奇思妙想，能得到什么？

答案是，我要给项目写个原创彩蛋，如下：

动起来的太阳系，还有运动的镜头。

代码来也！

---

function EasterEgg
% Draw several shapes rotating around the 
% parent star and the camera trajectory is an ellipse
planetsNum=4;
Sun_magnification=15;
Planets_magnification=10;
mainWin=figure(1);
set(mainWin, 'Toolbar', 'none','numbertitle','off');
set(mainWin,'menubar','none','Name', 'Easter Egg ')
axesmain=axes('Position', [0, 0, 1, 1]);
mainWin.WindowState = 'fullscreen';

%draw sun begin
[a,b,c]=sphere(99);
surf(axesmain,a*Sun_magnification,b*Sun_magnification,c*Sun_magnification);
newmap=[1 0 0];
colormap(newmap)
xx0=randn(3,999);       
xx0=(Sun_magnification+0.1)*xx0./vecnorm(xx0);
p=delaunay(xx0');

patch(axesmain,'faces',p,'vertices',xx0','FaceVertexCData',cool(size(p,1)),'FaceAlpha',.25);
set(mainWin,'color','k')
set(mainWin,'InvertHardCopy','off')
shading flat
%draw sun done


theta=[0 0 0 0 0 0 0 0];  %theta Planet_Sun_angle
phi=[7.01 3.39 0 1.85 1.31 2.49 0.77 1.77]; %phi Planet_Earth_angle
Ma=[0.4667,0.7282,1.0167,1.666,5.459,9.537,19.218,30.07]*100;
ecc_ratio=[0.2056,0.0068,0.0167,0.0934,0.0489,0.0565,0.0464,0.0097];
diametersOfPlanets = [4.88, 12.1, 12.742, 6.779, 139.82, 116.46, 50.72, 49.24]*Planets_magnification; % Relative size, unit: diameter of the earth
colorsOfPlanets = [0.5, 0.5, 0.5; 1, 0.8, 0; 0, 0, 1; 1, 0, 0; 1, 0.5, 0; 0.8, 0.8, 0.8; 0, 0.5, 1; 0, 0, 0.5]; % RGB of each planets
% period0=[0.24 0.62 1 1.88 11.86 29.46 84.01 164.79];
% period=[0.5 0.25 0.125 0.1 0.02 0.008 0.002 0.001];
period=[384 192 96 48 24 12 6 3];
if planetsNum<length(period)
    diametersOfPlanets(planetsNum+1:end) =[];
    colorsOfPlanets(planetsNum+1:end,:) =[];
    period(planetsNum+1:end)=[];
end


x0=[];y0=[];z0=[];
for i=1:planetsNum
    [x0(i,:),y0(i,:),z0(i,:)]=Orbit_coor_calcu(theta(i),phi(i),Ma(i),ecc_ratio(i));
    if i==1
        x=[x0(i,:)];
        y=[y0(i,:)];
        z=[z0(i,:)];
    elseif i>1
        x=[x,nan,x0(i,:)];
        y=[y,nan,y0(i,:)];
        z=[z,nan,z0(i,:)];
    end

end
% draw Planetary orbit
hold on;
plot3(axesmain,x, y, z, 'w', 'LineWidth', 0.1); 
% draw Planetary orbit done

% draw stars
tempBound=[x,y];
a0= max(abs(tempBound))*1.2;
farBound = a0*3;
r=@()rand(1,100);
arrayrand1 =(2*farBound) *r()-farBound;
arrayrand2 = (2*farBound) *r()-farBound;
arrayrand3 = (2*farBound/2) *r()-farBound/2;
strascoor=[arrayrand1;arrayrand2;arrayrand3];
index = strascoor(1, :) >= -a0 & strascoor(1, :) <= a0 & ...
        strascoor(2, :) >= -a0 & strascoor(2, :) <= a0 & ...
        strascoor(3, :) >= -a0/2 & strascoor(3, :) <= a0/2;
strascoor(:,index)=nan;

hold on
scatter3(axesmain,strascoor(1,:),strascoor(2,:),strascoor(3,:),'.w','SizeData', 50);
axis equal 
axis vis3d
axis off
% draw stars done

for j=1:planetsNum
    original_array = x0(j,:);
    interval=period(j)/period(end);
    expanded_new_array=Array_Thinning(original_array,interval);
    x1(j,:)=expanded_new_array;

    original_array = y0(j,:);
    expanded_new_array=Array_Thinning(original_array,interval);
    y1(j,:)=expanded_new_array;

    original_array = z0(j,:);
    expanded_new_array=Array_Thinning(original_array,interval);
    z1(j,:)=expanded_new_array;
    
end

[camtrack_x,camtrack_y,camtrack_z]=Orbit_coor_calcu(0,20,Ma(planetsNum),ecc_ratio(planetsNum));
camtrack_x=camtrack_x-max(camtrack_x)*0.8;

xlimset=size(x0,2);
m=xlimset;
step=1;
%draw Planets
hold on;
v=scatter3(axesmain,x1(:,m), y1(:,m), z1(:,m),diametersOfPlanets',colorsOfPlanets, 'filled'); 
%draw Planets done


%%%%%%%%%%%%% save gif
% filename = 'animation.gif';
% fps = 30; 
% gif_obj = [];
%%%%%%%%%%%%% save gif

while ishandle(1)
    m=m-step;
    if m==1
        m=xlimset;
    end
    % move plantes
    set(v,'xdata',x1(:,m),'ydata',y1(:,m),'zdata',z1(:,m))
    % move plantes

    camtarget(axesmain, [0, 0, 0]); 
    campos(axesmain,[camtrack_x(m)*8,camtrack_y(m)*8,camtrack_z(m)*5+250]); 
    pause(0.01)
    
    drawnow
%%%%%%%%%%%%%  save gif
%     frame = getframe(gcf);
%     img = frame2im(frame);
%     
%     if isempty(gif_obj)
%         [imind, cm] = rgb2ind(img, 256);
%         imwrite(imind, cm, filename, 'gif', 'Loopcount', inf, 'DelayTime', 1/fps);
%         gif_obj = false;
%     else
%         [imind, cm] = rgb2ind(img, 256);
%         imwrite(imind, cm, filename, 'gif', 'WriteMode', 'append', 'DelayTime', 1/fps);
%     end
%     
%     pause(0.01);
%%%%%%%%%%%%%  save gif
end


end

上面这些是主程序的内容，姑且叫“假”太阳系模拟吧，因为并不是按照真实的八大行星的运行速度模拟的。
大致就是①生成了一个太阳，②生成行星的轨道，③随机生成一些星星并且控制其在太阳系内不显示，④画出行星。这里用到了一点小心思，我并没有一个一个绘制所有行星，而是一次性绘制了所有行星，这样可以提高性能（尽管其实也没有必要提高），而且代码看着会精简一些。

---

另外动图嘛，当然要保存一哈，只需要将   %%%%%%%%%%%%% save gif   之内的代码取消注释就行。要想退出动画，那就按一下esc键就行了，因为窗口不关循环不停。

下面是两个自己写的函数，主程序会调用，建议将这俩函数分别存到不同的*.m文件中。

Orbit_coor_calcu函数是用法向量的方式生成三维椭圆用的；

Array_Thinning函数有点难以解释，大概是这个效果，例如：函数输入：original_array=[1 2 3 4 5 6 7 8 9]，interval=3，那么函数的输出就是：[1 4 7 1 4 7 1 4 7];

function [x_circle, y_circle, z_circle] = Orbit_coor_calcu(alpha, phi, major_axis, ecc_ratio)
num=768*4;
Ma = major_axis;
Sa = Ma * sqrt(1 - ecc_ratio^2);
x0 = sind(phi) * cosd(alpha);
y0 = sind(phi) * sind(alpha);
z0 = cosd(phi);
n = [x0, y0, z0]; % 法向量n

theta = linspace(0, 2*pi, num); % theta角从0到2*pi
a = cross(n, [1, 0, 0]); % n与i叉乘，求取a向量
if all(a == 0) % 如果a为零向量，将n与j叉乘
    a = cross(n, [0, 1, 0]);
end
b = cross(n, a); % 求取b向量
a = a / norm(a); % 单位化a向量
b = b / norm(b); % 单位化b向量

% 计算圆上各点的坐标，取实部确保结果是实数
x_circle = real(Ma * a(1) * cos(theta) + Sa * b(1) * sin(theta));
y_circle = real(Ma * a(2) * cos(theta) + Sa * b(2) * sin(theta));
z_circle = real(Ma * a(3) * cos(theta) + Sa * b(3) * sin(theta));

end



function expanded_new_array=Array_Thinning(original_array,interval)

    new_length = ceil(length(original_array) / interval);
    new_array = zeros(1, new_length);
    
    for i = 1:new_length
        index = (i - 1) * interval + 1;
        new_array(i) = original_array(index);
    end
    expanded_new_array = repmat(new_array, 1, interval);
end

哦对了，如果你看到这里了，那么我就说一嘴，主程序刚开始的位置

planetsNum=4;
Sun_magnification=15;
Planets_magnification=10;

这三个变量是可以更改的，planetsNum最大是8，要是你想让更多的行星出现，比如10，那么这些数组的维度都要改成10：

theta=[...];  %theta Planet_Sun_angle

phi=[...]; %phi Planet_Earth_angle

Ma=[...]*100;

ecc_ratio=[...];

diametersOfPlanets = [...]*Planets_magnification; % Relative size, unit: diameter of the earth

colorsOfPlanets = [...]; % RGB of each planets

period=[384 192 96 48 24 12 6 3];

要说一下数组period，里面的元素一定要是768*4的约数哟，至于为啥必须是768*4，是因为我闭着眼睛敲得三个数，写到最后发现：“哎！还可以哎”，但是代价就是，这个数不能改了（好像吧）

最后感谢 "Matlab矩阵研学社"，代码里的太阳就是模仿这个公众号里的内容画的。